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