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),
adjust_zeros(Decimal0, Decimal),
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.
add(Xs, Ys, Zs):-
reverse(Xs, ReversedXs),
reverse(Ys, ReversedYs),
add_1(ReversedXs, ReversedYs, 0, ReversedZs),
reverse(ReversedZs, Zs).
add_1([], Ys, 0, Ys):-!.
add_1([], Ys, 1, Zs):-
!, add_1([1], Ys, 0, Zs).
add_1(Xs, [], 0, Xs):-!.
add_1(Xs, [], 1, Zs):-
!, add_1(Xs, [1], 0, Zs).
add_1([X|Xs], [Y|Ys], Carry, [Z|Zs]):-
Z1 is X + Y + Carry,
Z1 >= 32767/*BASE*/, !,
Z is Z1 - 32767/*BASE*/,
add_1(Xs, Ys, 1, Zs).
add_1([X|Xs], [Y|Ys], Carry, [Z|Zs]):-
Z is X + Y + Carry,
add_1(Xs, Ys, 0, Zs).
% 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),
adjust_zeros(Zs0, Zs).
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),
adjust_zeros(Rs2, Rs),
adjust_zeros(Qs1, Qs).
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):-
pad_if_odd_length(Xs, [X1,X2|Xs1]),
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):-
add(Ys0, Ys0, Ts1),
append(Rs, [X1], Ts2),
divide(Ts2, Ts1, Ds, _),
Base1 = 32767/*BASE*/ - 1,
minimum(Ds, [Base1], [D0]),
append(Ts2, [X2], Ts3),
adjust_zeros(Ts3, Ts4),
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),
adjust_zeros([B|Bs], Ys).
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),
adjust_zeros(Ys1, Ys).
bigint_decimal_1([], _, Ys, Ys).
bigint_decimal_1([X|Xs], Base, Ys0, Ys):-
int_decimal(X, Ys1),
multiply10(Ys0, Base, Ys2),
add10(Ys1, Ys2, Ys3),
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.
add10(Xs, Ys, Zs):-
reverse(Xs, ReversedXs),
reverse(Ys, ReversedYs),
add10_1(ReversedXs, ReversedYs, 0, ReversedZs),
reverse(ReversedZs, Zs).
add10_1([], Ys, 0, Ys):-!.
add10_1([], Ys, 1, Zs):-
!, add10_1([1], Ys, 0, Zs).
add10_1(Xs, [], 0, Xs):-!.
add10_1(Xs, [], 1, Zs):-
!, add10_1(Xs, [1], 0, Zs).
add10_1([X|Xs], [Y|Ys], Carry, [Z|Zs]):-
Z1 is X + Y + Carry,
Z1 >= 10, !,
Z is Z1 - 10,
add10_1(Xs, Ys, 1, Zs).
add10_1([X|Xs], [Y|Ys], Carry, [Z|Zs]):-
Z is X + Y + Carry,
add10_1(Xs, Ys, 0, Zs).
% 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),
adjust_zeros(Qs1, Qs).
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.
adjust_zeros([], [0]):-!.
adjust_zeros([0], [0]):-!.
adjust_zeros([0|Xs], Ys):-!, adjust_zeros(Xs, Ys).
adjust_zeros(Xs, Xs).
% 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),
adjust_zeros(Qs1, Qs).
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.
pad_if_odd_length(Xs, Ys):-pad_if_odd_length_1(Xs, Xs, Ys).
pad_if_odd_length_1([], Xs, Xs).
pad_if_odd_length_1([_], Xs, [0|Xs]):-!.
pad_if_odd_length_1([_,_|Xs], Ys0, Ys):-pad_if_odd_length_1(Xs, Ys0, Ys).
% 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
%
addmod(A, B, N, C):-
add(A, B, T),
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),
add(T1, U2, W2), % W2=U2+V2*Q
multiply(V1, Q, T2),
add(T2, U1, W1), % W1=U1+V1*Q
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):-
add(P, [1], P_1),
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
% r is the answer:
% 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
addmod(H1, H1, P, TwoH1),
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).