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),
  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).

LPA Index     Home Page

Valid HTML 4.01 Strict