#### Multi-Precision Integer Arithmetic

Procedures for performing arithmetic on non-negative integers of arbitrary length, which we shall call bigints. Internally, a bigint is represented in base 32767 by a list of long integers. A bigint will not have any leading zeros, except in the case of the bigint 0.

Examples using this package can be found in Factoring (Pollard's p-1 Method) and Factoring (Pollard's rho Method).

```%
% User procedures
%

% int_bigint(Int, Bigint) is true if Bigint contains the value of Int.
int_bigint(Int, Bigint):-
int_decimal(Int, Decimal),
decimal_bigint(Decimal, Bigint).

% bigint_int(Bigint, Int) is true if Int contains the value of Bigint.
bigint_int(Bigint, Int):-
bigint_decimal(Bigint, Decimal),
decimal_int(Decimal, Int).

% string_bigint(String, Bigint) is true if String contains a sequence
%   of decimal digits, and Bigint contains the equivalent bigint
%   with leading zeros removed if present.
string_bigint(String, Bigint):-
name(String, Chars),
chars_decimal(Chars, Decimal0),
decimal_bigint(Decimal, Bigint).

% bigint_string(Bigint, String) is true if String contains the sequence
%   of decimal digits corresponding to the Bigint value.
bigint_string(Bigint, String):-
bigint_decimal(Bigint, Decimal),
decimal_chars(Decimal, Chars),
string_chars(String, Chars).

% Converts a list of bigints to strings.
bigints_strings([], []).
bigints_strings([Bigint|Bigints], [String|Strings]):-
bigint_string(Bigint, String),
bigints_strings(Bigints, Strings).

% equal(Xs, Xs) is true if the bigint Xs is equal to the bigint Ys.
equal(Xs, Xs).

% less_than(Xs, Xs) is true if the bigint Xs is less than the bigint Ys.
less_than(Xs, Ys):-
shorter_than(Xs, Ys), !.
less_than(Xs, Ys):-
same_length(Xs, Ys),
less_than_1(Xs, Ys).

less_than_1([X|_], [Y|_]):-
X < Y, !.
less_than_1([X|Xs], [X|Ys]):-
less_than_1(Xs, Ys).

% less_than_or_equal(Xs, Xs) is true if the bigint Xs is less than or
%   equal to the bigint Ys.
less_than_or_equal(Xs, Ys):-
\+ less_than(Ys, Xs).

% add(Xs, Ys, Zs) is true if adding the bigint Xs to the bigint Ys gives
%   the bigint Zs.
reverse(Xs, ReversedXs),
reverse(Ys, ReversedYs),
reverse(ReversedZs, Zs).

Z1 is X + Y + Carry,
Z1 >= 32767/*BASE*/, !,
Z is Z1 - 32767/*BASE*/,
Z is X + Y + Carry,

% subtract(Xs, Ys, Zs) is true if subtracting the bigint Ys from the
%   bigint Xs gives the bigint Zs.
subtract(Xs, Ys, Zs):-
less_than_or_equal(Ys, Xs),
reverse(Xs, ReversedXs),
reverse(Ys, ReversedYs),
subtract_1(ReversedXs, ReversedYs, 0, ReversedZs),
reverse(ReversedZs, Zs0),

subtract_1(Xs, [], 0, Xs):-!.
subtract_1(Xs, [], -1, Zs):-
!, subtract_1(Xs, [1], 0, Zs).
subtract_1([X|Xs], [Y|Ys], Carry, [Z|Zs]):-
Z1 is X - Y + Carry,
Z1 < 0, !,
Z is Z1 + 32767/*BASE*/,
subtract_1(Xs, Ys, -1, Zs).
subtract_1([X|Xs], [Y|Ys], Carry, [Z|Zs]):-
Z is X - Y + Carry,
subtract_1(Xs, Ys, 0, Zs).

% absolute_difference(Xs, Ys, Zs) is true if Zs is the absolute
%   difference between the bigints Xs and Ys.
absolute_difference(Xs, Ys, Zs):-
less_than(Xs, Ys),
!,
subtract(Ys, Xs, Zs).
absolute_difference(Xs, Ys, Zs):-
subtract(Xs, Ys, Zs).

% multiply(Xs, Ys, Zs) is true if multiplying the bigint Xs by the bigint
%   Ys gives the bigint Zs.
multiply([0], _, [0]):-!.
multiply(_, [0], [0]):-!.
multiply(Xs, Ys, Zs):-
reverse(Ys, ReversedYs),
multiply_1(Xs, ReversedYs, [], ReversedZs),
reverse(ReversedZs, Zs).

multiply_1([], _, Zs, Zs).
multiply_1([X|Xs], Ys, Zs0, Zs):-
multiply_2(Ys, X, 0, Zs1),
add_1([0|Zs0], Zs1, 0, Zs2), % The prefixed zero is the low-order digit
multiply_1(Xs, Ys, Zs2, Zs).

multiply_2([], _, 0, []):-!.
multiply_2([], _, Carry, [Carry]):-
Carry > 0.
multiply_2([Y|Ys], X, Carry, [Z|Zs]):-
Z1 is X * Y + Carry,
Z  is Z1 mod 32767/*BASE*/,
Carry1 is Z1 // 32767/*BASE*/,
multiply_2(Ys, X, Carry1, Zs).

% divide(Xs, Ys, Qs, Rs) is true if the result of dividing the bigint Xs
%   by the bigint Ys is the bigint quotient Qs and the bigint remainder Rs.
divide(_, [0], [0], [0]):-
!, write(`Error - attempted division by zero`), nl,
fail.
divide(Xs, Ys, [0], Xs):-
less_than(Xs, Ys), !.
divide(Xs, [Y], Qs, [R]):-
!, divide_short(Xs, Y, Qs, R).
divide(Xs, Ys, Qs, Rs):-
Ys = [Y|_],
D is 32767/*BASE*/ // (Y + 1),
mult(Xs, D, Xs0),
mult(Ys, D, [0|Ys0]),
split(Ys, Xs0, Xs1, Xs2), % Xs1 has same length as Ys
divide_1(Xs2, Ys0, Xs1, Qs1, Rs1),
% Qs1 may have a leading zero e.g. 20/4
div(Rs1, D, Rs2, 0),

divide_1([], _, Rs, [], Rs).
divide_1([X|Xs], Ys, Zs, [Q|Qs], Rs):-
Zs = [Z1,Z2|_],
Ys = [Y1|_],
Q1 is (Z1 * 32767/*BASE*/ + Z2) // Y1,
Base1 is 32767/*BASE*/ - 1,
min(Base1, Q1, Qmax),
append(Zs, [X], Xs1),
divide_2(Qmax, Xs1, Ys, Q, [0|Zs1]),
% Here Xs1 has n+1 digits, and Ys and Zs1 have n digits
divide_1(Xs, Ys, Zs1, Qs, Rs).

% divide_2(Qmax, Xs, Ys, Q, Rs) is true if the result of dividing the
%   (n+1)-place bigint Xs by the n-place bigint Ys is the integer
%   quotient Q and the n-place bigint remainder Rs.  It is required, but
%   not checked, that 0 <= Xs/Ys < Base, and that Qmax is an upper limit
%   of the value of Q.
divide_2(Q, Xs, Ys, Q, Rs):-
mult(Ys, Q, QYs),
le(QYs, Xs), !,
subt(Xs, QYs, Rs).
divide_2(Q0, Xs, Ys, Q, Rs):-
Q1 is Q0 - 1,
divide_2(Q1, Xs, Ys, Q, Rs).

% square_root(Xs, Ys) is true if the bigint Ys is the floor of the square
%   root of the bigint Xs.
square_root(Xs, Ys):-
T is X1 * 32767/*BASE*/ + X2,
Y is int(sqrt(T)),
R is T - Y * Y,
int_decimal(R, Rs0),
decimal_bigint(Rs0, Rs),
square_root_1(Xs1, Rs, [Y], Ys).

square_root_1([], _, Ys, Ys).
square_root_1([X1,X2|Xs], Rs, Ys0, Ys):-
append(Rs, [X1], Ts2),
divide(Ts2, Ts1, Ds, _),
Base1 = 32767/*BASE*/ - 1,
minimum(Ds, [Base1], [D0]),
append(Ts2, [X2], Ts3),
square_root_2(D0, Ts1, Ts4, D, Rs1),
append(Ys0, [D], Ys1),
square_root_1(Xs, Rs1, Ys1, Ys).

square_root_2(0, _, Bs, 0, Bs):-!.
square_root_2(D, As, Bs, D, Rs):-
append(As, [D], Ts1),
multiply_short(Ts1, D, Ts2),
less_than_or_equal(Ts2, Bs), !,
subtract(Bs, Ts2, Rs).
square_root_2(D0, As, Bs, D, Rs):-
D1 is D0 - 1,
square_root_2(D1, As, Bs, D, Rs).

% power(Xs, Ns, Ys) is true if Ys is the result of raising Xs to the
%   power Ns, where all the arguments are bigints.
power(Xs, Ns, Ys):-power_1(Ns, Xs, [1], Ys).

power_1([0], _, Y, Y):-!.
power_1(N, X, Y0, Y):-
divide_short(N, 2, N1, R),
R = 1, !,
multiply(Y0, X, Y1),
multiply(X, X, X1),
power_1(N1, X1, Y1, Y).
power_1(N, X, Y0, Y):-
divide_short(N, 2, N1, _),
multiply(X, X, X1),
power_1(N1, X1, Y0, Y).

% gcd(Xs, Ys, Zs) is true if the bigint Zs is the Greatest Common Divisor
%   of the bigints Xs and Ys.
gcd(Xs, [0], Zs):-!,
Zs=Xs.
gcd(Xs, Ys, Zs):-
divide(Xs, Ys, _, Rs),
gcd(Ys, Rs, Zs).

% random_bigint(Xs, Ys) is true if Ys is a pseudo-random bigint greater
%   than or equal to 0 and less than the bigint Xs.
random_bigint(Xs, Ys):-
\+ equal(Xs, [0]),
Xs=[A|As],
B is ip(rand(A)),
random_bigint_1(As, Bs),

random_bigint_1([], []).
random_bigint_1([_|As], [B|Bs]):-
B is ip(rand(32767/*BASE*/)),
random_bigint_1(As, Bs).

%
% Utility procedures
%

% e.g. chars_decimal("1234567", X) gives X=[1,2,3,4,5,6,7]
chars_decimal([], []).
chars_decimal([C|Cs], [D|Ds]):-
C >= "0", C =< "9",
D is C - "0",
chars_decimal(Cs, Ds).

% e.g. decimal_chars([1,2,3,4,5,6,7], X) gives X=[49,50,51,52,53,54,55]
decimal_chars([], []).
decimal_chars([D|Ds], [C|Cs]):-
C is D + "0",
decimal_chars(Ds, Cs).

% e.g. decimal_bigint([1,1,1,1,1,1], X) gives X=[3,12810]
decimal_bigint([0], [0]):-!.
decimal_bigint([X|Xs], Ys):-
decimal_bigint_1([X|Xs], [], Ys).

decimal_bigint_1([0], Ys, Ys):-!.
decimal_bigint_1(Xs, Ys0, Ys):-
divide10(Xs, 32767/*BASE*/, Qs, R),
decimal_bigint_1(Qs, [R|Ys0], Ys).

% e.g. bigint_decimal([3,12810], X) gives X=[1,1,1,1,1,1]
bigint_decimal([X|Xs], Ys):-
int_decimal(X, Ys0),
int_decimal(32767/*BASE*/, Base1),
bigint_decimal_1(Xs, Base1, Ys0, Ys1),

bigint_decimal_1([], _, Ys, Ys).
bigint_decimal_1([X|Xs], Base, Ys0, Ys):-
int_decimal(X, Ys1),
multiply10(Ys0, Base, Ys2),
bigint_decimal_1(Xs, Base, Ys3, Ys).

% add10(Xs, Ys, Zs) is true if adding the decimal Xs to the decimal Ys
%   gives the decimal Zs.
reverse(Xs, ReversedXs),
reverse(Ys, ReversedYs),
reverse(ReversedZs, Zs).

Z1 is X + Y + Carry,
Z1 >= 10, !,
Z is Z1 - 10,
Z is X + Y + Carry,

% multiply10(Xs, Ys, Zs) is true if multiplying the decimal Xs by the
%   decimal Ys gives the decimal Zs.
multiply10(Xs, Ys, Zs):-
reverse(Ys, ReversedYs),
multiply10_1(Xs, ReversedYs, [], ReversedZs),
reverse(ReversedZs, Zs).

multiply10_1([], _, Zs, Zs).
multiply10_1([X|Xs], Ys, Zs0, Zs):-
multiply10_2(Ys, X, 0, Zs1),
add10_1([0|Zs0], Zs1, 0, Zs2), % The prefixed zero is the low-order digit
multiply10_1(Xs, Ys, Zs2, Zs).

multiply10_2([], _, 0, []):-!.
multiply10_2([], _, Carry, [Carry]):-Carry > 0, !.
multiply10_2([Y|Ys], X, Carry, [Z|Zs]):-
Z1 is X * Y + Carry,
Z  is Z1 mod 10,
Carry1 is Z1 // 10,
multiply10_2(Ys, X, Carry1, Zs).

% divide10(Xs, Y, Qs, R) is true if the result of dividing the decimal
%   Xs by the integer Y is the decimal quotient Qs and the integer
%   remainder R.
divide10(_, 0, [0], 0):-
!, write(`Error - attempted division by zero`), nl,
fail.
divide10([X], Y, [0], X):-
X < Y, !.
divide10(Xs, Y, Qs, R):-
Xs = [_|_],
divide10_1(Xs, Y, 0, Qs1, R),

divide10_1([], _, R, [], R).
divide10_1([X|Xs], Y, R0, [Q|Qs], R):-
T  is R0 * 10 + X,
Q  is T // Y,
R1 is T mod Y,
divide10_1(Xs, Y, R1, Qs, R).

% adjust_zeros(Xs, Ys) is true if the list Ys is the result of removing
%   leading zeros from the list Xs, except that the list consisting of a
%   single zero is not changed, and an empty list is converted to a list
%   consisting of a single zero.

% same_length(Xs, Ys) is true if the lists Xs and Ys contain the same
%   number of elements.
same_length([], []).
same_length([_|Xs], [_|Ys]):-same_length(Xs, Ys).

% shorter_than(Xs, Ys) is true if the list Xs contains fewer elements
%   than the list Ys.
shorter_than([], [_|_]).
shorter_than([_|Xs], [_|Ys]):-shorter_than(Xs, Ys).

% divide_short(Xs, Y, Qs, R) is true if the result of dividing the bigint
%   Xs by the integer Y (in the same base) is the bigint quotient Qs and
%   the integer remainder R.
divide_short(_, 0, [0], 0):-
!, write(`Error - attempted division by zero`), nl,
fail.
divide_short([X], Y, [0], X):-
X < Y, !.
divide_short(Xs, Y, Qs, R):-
Xs = [_|_],
divide_short_1(Xs, Y, 0, Qs1, R),

divide_short_1([], _, R, [], R).
divide_short_1([X|Xs], Y, R0, [Q|Qs], R):-
T  is R0 * 32767/*BASE*/ + X,
Q  is T // Y,
R1 is T mod Y,
divide_short_1(Xs, Y, R1, Qs, R).

% subt(Xs, Ys, Zs) is true if subtracting the n-place bigint Ys from the
%   n-place bigint Xs gives the n-place bigint Zs.
subt(Xs, Ys, Zs):-
reverse(Xs, ReversedXs),
reverse(Ys, ReversedYs),
subtract_1(ReversedXs, ReversedYs, 0, ReversedZs),
reverse(ReversedZs, Zs).

% mult(Xs, Y, Zs) is true if multiplying the n-place bigint Xs by the
%   integer Y gives the (n+1)-place bigint Zs.
mult(Xs, Y, Zs):-
reverse(Xs, Xs1),
mult_1(Xs1, Y, 0, Zs1),
reverse(Zs1, Zs).

mult_1([], _, Carry, [Carry]).
mult_1([Y|Ys], X, Carry, [Z|Zs]):-
Z1 is X * Y + Carry,
Z  is Z1 mod 32767/*BASE*/,
Carry1 is Z1 // 32767/*BASE*/,
mult_1(Ys, X, Carry1, Zs).

% split(As, Bs, Cs, Ds) is true if the list Cs has the same number of
%   elements as the list As and it contains the leading elements of the
%   list Bs, and the list Ds contains the remaining elements of the
%   list Bs.
split([], Bs, [], Bs).
split([_|As], [B|Bs], [B|Cs], Ds):-split(As, Bs, Cs, Ds).

% div(Xs, Y, Qs, R) is true if the result of dividing the n-place bigint
%   Xs by the integer Y (in the same base) is the n-place bigint quotient
%   Qs and the integer remainder R.
div([X], Y, [0], X):-
X < Y, !.
div(Xs, Y, Qs, R):-
Xs = [_|_],
divide_short_1(Xs, Y, 0, Qs, R).

% lt(Xs, Xs) is true if the n-place bigint Xs is less than the n-place
%   bigint Ys.
lt([X|_], [Y|_]):-
X < Y, !.
lt([X|Xs], [X|Ys]):-
lt(Xs, Ys).

% le(Xs, Xs) is true if the n-place bigint Xs is less than or equal to
%   the n-place bigint Ys.
le(Xs, Ys):-
\+ lt(Ys, Xs).

% min(X, Y, Z) is true if Z is the minimum of the numbers X and Y.
min(X, Y, Z):-X =< Y, Z = X, !.
min(_, Y, Y).

% pad_if_odd_length(Xs, Ys) is true if the list Ys is identical to the
%   list Xs except that it is prefixed with a 0 if Xs has an odd number
%   of elements.

% int_decimal(X, Ys) is true if Ys is the base-10 bigint representation
%   of the integer X.
% e.g. int_decimal(12345, X) gives X=[1,2,3,4,5]
int_decimal(0, [0]):-!.
int_decimal(X, Ys):-
X > 0,
int_decimal_1(X, [], Ys).

int_decimal_1(0, Ys, Ys):-!.
int_decimal_1(X, Ys0, Ys):-
Y  is X mod 10,
X1 is X // 10,
int_decimal_1(X1, [Y|Ys0], Ys).

% decimal_int(Ints, X) is true if X is the int represented by the base-10
%   bigint Ints.
% e.g. decimal_int([1,2,3,4,5], X) gives X=12345
decimal_int([Int|Ints], X):-decimal_int_1(Ints, Int, X).

decimal_int_1([], X, X):-!.
decimal_int_1([Int|Ints], X0, X):-
Int >= 0,
X1 is X0 * 10 + Int,
decimal_int_1(Ints, X1, X).

% minimum(Xs, Ys, Zs) is true if the bigint Zs is the smaller of the
%   bigints Xs and Ys.
minimum(Xs, Ys, Zs):-less_than_or_equal(Xs, Ys), !, Zs = Xs.
minimum(_, Ys, Ys).

% multiply_short(Xs, Y, Zs) is true if multiplying the bigint Xs by the
%   integer Y gives the bigint Zs.
multiply_short(_, 0, [0]):-!.
multiply_short(Xs, Y, Zs):-
reverse(Xs, Xs1),
multiply_2(Xs1, Y, 0, Zs1),
reverse(Zs1, Zs).

%
% Modular operations
%
divide(T, N, _, C).

subtmod(A, B, N, C):-
less_than_or_equal(B, A), !,
subtract(A, B, T),
divide(T, N, _, C).
subtmod(A, B, N, C):-
subtract(B, A, T),
divide(T, N, _, C0),
subtract(N, C0, C).

multmod(A, B, N, C):-
multiply(A, B, T),
divide(T, N, _, C).

divmod(A, B, N, C):-
inversemod(B, N, T),
multmod(A, T, N, C).

% inversemod(B, N, X) is true if X is the inverse of B modulo N, that is,
%   B * X = (1 modulo N). This requires that B and N be mutually prime.
%   The variable S is used in order to avoid negative values.
% e.g. inversemod([60], [103], X) gives X=[91]
%      inversemod([59], [103], X) gives X=[7]
inversemod(B, N, X):-
less_than([1], N),
gcd_extended(B, N, S, X0, _, [1]),
inversemod_1(S, N, X0, X).

inversemod_1(-1, N, X0, X):-!,
subtract(N, X0, X).
inversemod_1(_, _, X, X).

% gcd_extended(X, Y, S, A, B, G) is true
%   if S=1  and X*A - Y*B = G = gcd(X, Y),
%   or S=-1 and Y*B - X*A = G = gcd(X, Y).
% e.g. gcd_extended([23],  [69],  S, A, B, G) gives S=1, A=[1],B=[0],G=[23]
%      gcd_extended([112], [126], S, A, B, G) gives S=-1,A=[1],B=[1],G=[14]
gcd_extended(X, Y, S, A, B, G):-
gcd_extended_1(1,  Y, [1], [0],  X, [0], [1],  S, G, B, A).

gcd_extended_1(S0,  [0], _, _,  U3, U2, U1,  S, G, B, A):-!,
S = S0,
G = U3,
B = U2,
A = U1.
gcd_extended_1(S,  V3, V2, V1,  U3, U2, U1,  N, U3P, U2P, U1P):-
divide(U3, V3, Q, W3), % Q=U3/V3 W3=U3%V3
multiply(V2, Q, T1),
multiply(V1, Q, T2),
S1 is -S,
gcd_extended_1(S1,  W3, W2, W1,  V3, V2, V1,  N, U3P, U2P, U1P).

% powermod(Xs, Ns, Ms, Ys) is true if Ys is the result of raising Xs to
%   the power Ns, mod Ms, where all the arguments are bigints.
powermod(Xs, Ns, Ms, Ys):-powermod_1(Ns, Xs, Ms, [1], Ys).

powermod_1([0], _, _, Y, Z):-!,
Z=Y.
powermod_1(N, X, M, Y0, Y):-
divide_short(N, 2, N1, R),
R = 1, !,
multiply(Y0, X, Y1),
divide(Y1, M, _, Y2),
multiply(X, X, X1),
divide(X1, M, _, X2),
powermod_1(N1, X2, M, Y2, Y).
powermod_1(N, X, M, Y0, Y):-
divide_short(N, 2, N1, _),
multiply(X, X, X1),
divide(X1, M, _, X2),
powermod_1(N1, X2, M, Y0, Y).

% sqrtmod(N, P, R) is true if P is a prime number (not verified), N is a
%   quadratic residue modulo P, and +R and -R are the square roots of N
%   modulo P. All the arguments are bigints.
sqrtmod(N, [2], [R]):-!,
divide_short(N, 2, _, R).
sqrtmod(N, P, R):-
divide_short(P, 4, _, 3), !,
sqrt_3_4(N, P, R).
sqrtmod(N, P, R):-
divide_short(P, 8, _, 5), !,
sqrt_5_8(N, P, R).
sqrtmod(N, P, R):-
divide_short(P, 8, _, 1), !,
sqrt_1_8(N, P, R).

% sqrt_3_4(N, P, R) is true if P is a prime number congruent to 3 modulo 4,
%   N is a quadratic residue modulo P, and +R and -R are the square roots
%   of N modulo P. All the arguments are bigints.
% e.g. sqrt_3_4([105], [107], X) gives X=[76]
sqrt_3_4(N, P, R):-
divide_short(P_1, 4, P_1_4, _),
powermod(N, P_1_4, P, R), % r=n^((p+1)/4) mod p
multmod(R, R, P, N).      % Sanity check

% sqrt_5_8(N, P, R) is true if P is a prime number congruent to 5 modulo 8,
%   N is a quadratic residue modulo P, and +R and -R are the square roots
%   of N modulo P. All the arguments are bigints.
% e.g. sqrt_5_8([102], [109], X) gives X=[50]
sqrt_5_8(N, P, R):-
subtract(P, [5], P_5),
divide_short(P_5, 8, P_5_8, _),
addmod(N, N, P, N2),       % Atkin's algorithm:
powermod(N2, P_5_8, P, V), % v=(2n)^((p-5)/8) mod p
multmod(N2, V, P, I0),
multmod(I0, V, P, I),      % i=(2nv^2) mod p
subtmod(I, [1], P, I_1),
multmod(N, V, P, NV),
multmod(NV, I_1, P, R),    % r=(nv(i-1)) mod p
multmod(R, R, P, N).       % Sanity check

% sqrt_1_8(N, P, R) is true if P is a prime number congruent to 1 modulo 8,
%   N is a quadratic residue modulo P, and +R and -R are the square roots
%   of N modulo P. All the arguments are bigints.
% Shank's algorithm [Knuth Exercise 4.6.2-15]: Find e and q such that
%   p=(2^e)q+1 and e>2. Choose x at random in the range 1 < x < p, and set
%   z=x^q mod p. If z^(2^(e-1)) mod p=1, repeat this step with a different
%   value of x. Otherwise initialize y=z, s=e, x=a^((q-1)/2) mod p,
%   r=ax mod p and w=ax^2 mod p, and repeat the following until w=1, when
%     Find the smallest k such that w^(2^k) mod p=1. If k=s there is no
%     solution. Otherwise set y=y^(2^(s-k)), s=k, r=ry^(2^(s-k-1)) and
%     w=wy^(2^(s-k)).
% e.g. sqrt_1_8([102], [113], X) gives X=[21]
sqrt_1_8([0], _, [0]):-!.
sqrt_1_8(N, P, R):-
divide(P, [8], _, [1]),      % p==1(mod 8)
subtract(P, [1], P1),
sqrt_1_8_e_q(P1, Q, 0, E),   % p=(2^e)q+1
sqrt_1_8_z(1, P, E, Q, Z),   % z=x^q mod p such that z^(2^(e-1)) mod p=1
subtract(Q, [1], QQ),
divide_short(QQ, 2, Q1, _),
powermod(N, Q1, P, X),       % x=a^((q-1)/2) mod p
multmod(N, X, P, R0),        % r=ax mod p
multmod(R0, X, P, W),        % w=ax^2 mod p
sqrt_1_8_iterate(W, P, X, Z, E, R0, R),
multmod(R, R, P, N).         % Sanity check

sqrt_1_8_iterate([1], _, _, _, _, R, R):-!.
sqrt_1_8_iterate(W, P, X, Y, S, R0, R):-
sqrt_1_8_k(1, S, W, P, K),   % Smallest k < s such that w^(2^k) mod p=1
S_K is S - K,
S_K_1 is S - K - 1,
int_bigint(S_K, S_Ks),
int_bigint(S_K_1, S_K_1s),
powermod([2], S_Ks, P, S_K_2),
powermod([2], S_K_1s, P, S_K_1_2),
powermod(Y, S_K_2, P, Y1),   % y1=y^(2^(s-k))
S1 is K,                     % s1=k
powermod(Y, S_K_1_2, P, T),
multmod(R0, T, P, R1),       % r1=ry^(2^(s-k-1))
multmod(W, Y1, P, W1),       % w1=wy^(2^(s-k))
sqrt_1_8_iterate(W1, P, X, Y1, S1, R1, R).

% sqrt_1_8_k(1, S, W, P, K) is true if k is the smallest integer less than
%   s such that w^(2^k) mod p=1.
sqrt_1_8_k(K, S, W, P, K):-
K < S,
int_bigint(K, Kb),
powermod([2], Kb, P, T),
powermod(W, T, P, [1]), !.
sqrt_1_8_k(I, S, W, P, K):-
I < S,
I1 is I + 1,
sqrt_1_8_k(I1, S, W, P, K).

% sqrt_1_8_e_q(Q0, Q, 0, E) is true if q0 = (2^e)q.
% e.g. sqrt_1_8_e_q([88], [11], 0, X) gives X=3
sqrt_1_8_e_q(Q0, Q, E0, E):-
divide_short(Q0, 2, Q1, 0), !,
E1 is E0 + 1,
sqrt_1_8_e_q(Q1, Q, E1, E).
sqrt_1_8_e_q(Q, Q, E, E).

% sqrt_1_8_z(2, P, E, Q, Z) is true if x is the smallest integer [it should
%   be a random integer] such that z=x^q mod p and z^(2^(e-1)) mod p != 1.
sqrt_1_8_z(X, P, E, Q, Z):-
int_bigint(X, Xb),
powermod(Xb, Q, P, Z),
E_1 is E - 1,
int_bigint(E_1, E_1b),
powermod([2], E_1b, P, T),
powermod(Z, T, P, S),
S \= [1], !.
sqrt_1_8_z(X, P, E, Q, Z):-
X1 is X + 1,
sqrt_1_8_z(X1, P, E, Q, Z).

% sqrt_prime_squared(N, P, R1, R2) is true if P is a prime (not verified),
%   N is a quadratic residue modulo P, and R1 and R2 are the square roots
%   of N modulo P*P. All the arguments are bigints.
% e.g. sqrt_prime_squared([39], [7], X, Y) gives X=[23], Y=[26]
sqrt_prime_squared(N, P, R1, R2):-
divide(N, P, _, NP),
sqrtmod(NP, P, H1),          % h1=n^0.5 mod p
multiply(P, P, P2),
divide(N, P2, _, NP2),
multiply(H1, H1, H12),
sqrt_prime_squared_1(NP2, H12, P, Temp),
divmod(Temp, TwoH1, P, H2),  % h2=((n-h1^2)/p)/(2h1) mod p
multiply(H2, P, H2P),
addmod(H1, H2P, P2, R1),     % r1=h1+h2p mod p^2
subtract(P2, R1, R2),        % r2=p^2-r1
multmod(R1, R1, P2, NP2),    % Sanity check
multmod(R2, R2, P2, NP2).    % Sanity check

sqrt_prime_squared_1(NP2, H12, P, R):-
less_than_or_equal(H12, NP2), !,
subtract(NP2, H12, Temp),
divide(Temp, P, R, _).
sqrt_prime_squared_1(NP2, H12, P, R):-
subtract(H12, NP2, Temp1),
divide(Temp1, P, Temp2, _),
subtract(P, Temp2, R).

% append(Xs, Ys, Zs) is true if Zs is the result of appending the list Xs to
%   the list Ys.
%append([], Ys, Ys).
%append([X|Xs], Ys, [X|Zs]):-append(Xs, Ys, Zs).
```