This code requires the Modular Arithmetic package.
% ph(A, B, P, X) is true if A is a primitive root of the prime P, P is % less than 10^15, B is a value between 1 and P-1, and X is the % solution to A^X = B mod P. % This is the Pohlig-Hellman Discrete Logarithm algorithm. % e.g. ph(31415, 926535897, 93238462643, X) % gives X=56545841881 % ph(12747889823615, 13264032151410, 29176179583009, X) % gives X=19989561017406 % Solutions may be obtained for some values of A which are not primitive % roots of P. For example, ph(3332, 22333, 23333, X) gives X=1799 ph(A, B, P, X):- 0 < A, A < P, 0 < B, B < P, is_probably_prime(P), N is P - 1, prime_factorization(N, QEs), ph(A, B, P, QEs, X). % ph(A, B, P, QEs, X) is true if A is a primitive root of the prime P, % B is a value between 1 and P-1, the elements [Q,E] of QEs give the % prime factorization of P-1, and X is the solution to A^X = B mod P. % e.g. ph(31415, 92653, 589793, [[2,5],[7,1],[2633,1]], X) gives X=180541 ph(A, B, P, QEs, X):- N is P - 1, ph_1(QEs, A, B, P, N, Xs, QE1s), art(Xs, QE1s, X), powmod(A, X, P, B). % Sanity check % For each [Q,E] do... ph_1([], _, _, _, _, [], []). ph_1([[Q,E]|QEs], A, B, P, N, [X|Xs], [QE|QE1s]):- NQ is N // Q, powmod(A, NQ, P, A1), % a'=a^(n/q) mod p ph_2(0, E, Q, 0, 1, A, A1, B, P, N, 1, 0, 0, X0), power(Q, E, QE), X is X0 mod QE, ph_1(QEs, A, B, P, N, Xs, QE1s). % For J from 0 to E-1 do... ph_2(E, E, _, _, _, _, _, _, _, _, _, _, X, X):-!. ph_2(J, E, Q, QJm, QJ, A, A1, B, P, N, G, Prev, X0, X):- multmod(QJ, Q, P, QJp), multmod(Prev, QJm, P, T1), powmod(A, T1, P, T2), multmod(G, T2, P, G1), % g=g*a^(prev*q^(j-1)) inversemod(G1, P, T3), multmod(B, T3, P, T4), NQ is N // QJp, powmod(T4, NQ, P, B1), % b'=(b*g^(-1))^(n/q^(j+1)) rho(A1, B1, P, Q, Prev1), % prev=log_a'(b') in the subgroup of order q multmod(Prev1, QJ, P, T6), X1 is X0 + T6, % x=x+prev*q^j J1 is J + 1, ph_2(J1, E, Q, QJ, QJp, A, A1, B, P, N, G1, Prev1, X1, X). % power(X, N, Y) is true if Y is the result of raising X to the power N. power(X, N, Y):-N >= 0, power_1(N, X, 1, Y). power_1(0, _, Y, Y):-!. power_1(1, X, Y0, Y):-!, Y is Y0 * X. power_1(N, X, Y0, Y):- 1 =:= N mod 2, !, N1 is N // 2, X1 is X * X, Y1 is Y0 * X, power_1(N1, X1, Y1, Y). power_1(N, X, Y0, Y):- N1 is N // 2, X1 is X * X, power_1(N1, X1, Y0, Y). % art(As, Ms, X) is true if each two elements of Ms are coprime, and X is % the smallest solution of the congruences X = A mod M for each element % M of Ms and the corresponding element A of As. This algorithm has the % advantage of performing all operations modulo the values of M, whereas % other algorithms perform some operations modulo the product of the values % of M. See "Aryabhata Remainder Theorem: Relevance to public-key crypto % algorithms" by T.R.N. Rao and Chung-Huang Yang. % e.g. art([3,1,4,1,59], [265,3,589,79,3238462], X) gives X=96506481730873 art([A|As], [M|Ms], X):- A >= 0, M > 1, art_1(As, Ms, M, 1, A, X). % n_1=1, x_1=a_1 art_1([], [], _, _, X, X). art_1([A|As], [M|Ms], M0, N0, X0, X):- A >= 0, M > 1, N1 is N0 * M0, % n_i=n_(i-1)*m_(i-1) inversemod(N1, M, C), % c=n_i^-1 mod m_i X0M is X0 mod M, subtmod(A, X0M, M, AX0), multmod(AX0, C, M, U), % u=((a_i-x_(i-1))*c) mod m_i X1 is X0 + U * N1, % x_i=x_(i-1)+u*n_i art_1(As, Ms, M, N1, X1, X). % e.g. prime_factorization(314159265, PEs) % gives PEs=[[3,2],[5,1],[7,1],[127,1],[7853,1]] prime_factorization(N, PEs):- prime_factors(N, [P|Ps]), prime_factorization_1(Ps, P, 1, PEs). prime_factorization_1([], P, E, [[P,E]]). prime_factorization_1([P|Ps], P, E, PEs):-!, E1 is E + 1, prime_factorization_1(Ps, P, E1, PEs). prime_factorization_1([P|Ps], P0, E, [[P0,E]|PEs]):- prime_factorization_1(Ps, P, 1, PEs). % distinct_prime_factors(N, Qs) is true if Qs contains the distinct prime % factors of N in ascending order of sequence. % e.g. distinct_prime_factors(314159265, Qs) gives Qs=[3,5,7,127,7853] distinct_prime_factors(N, Qs):- prime_factors(N, [P|Ps]), distinct_prime_factors_1(Ps, P, Qs). distinct_prime_factors_1([], P, [P]). distinct_prime_factors_1([P|Ps], P, Qs):-!, distinct_prime_factors_1(Ps, P, Qs). distinct_prime_factors_1([P|Ps], Q, [Q|Qs]):- distinct_prime_factors_1(Ps, P, Qs). % prime_factors(N, Ps) is true if the elements of the list Ps are the prime % factors of N in ascending order of sequence, and N is less than 10^15. % If N cannot be factored, [0] is returned. % e.g. prime_factors(31415926535897, Ps) gives Ps=[163,44777,4304347] prime_factors(N, Ps):- N < 1.0E15, prime_factors_1(N, [], Ps). prime_factors_1(N, Ps0, Ps):- prime_factors_3(1, N, P), % P is a factor of N prime_factors_2(P, N, Ps0, Ps). prime_factors_2(0, _, _, [0]):-!. % Unable to factor - return 0 prime_factors_2(N, N, Ps0, Ps):-!, % N is prime insert(Ps0, N, Ps). % Insert N prime_factors_2(P, N, Ps0, Ps):- % N is composite Q is N // P, % N = P * Q prime_factors_1(P, Ps0, Ps1), % Insert the factors of P prime_factors_1(Q, Ps1, Ps). % Insert the factors of Q prime_factors_3(K, N, P):- squfof_1(N, K, P), !. % P is a factor of N prime_factors_3(K, N, P):- K1 is K + 1, % Try the next value of K prime_factors_3(K1, N, P). % insert(Ys, X, Zs) is true if Zs is the result of inserting X into the % ordered list Ys. insert([Y|Ys], X, [Y|Zs]):-Y < X, !, insert(Ys, X, Zs). insert(Ys, X, [X|Ys]). % squfof(N, F) is true if F is a factor, not necessarily prime, of N, and % is less than 10^15. % This is Shanks' square forms factorization algorithm. % e.g. squfof(31415926535897, F) gives F=44777 squfof(N, F):- for(1, 100, K), % Find a suitable multiplier squfof_1(N, K, F), !. squfof_1(N, _, N):- N < 1.0E15, is_probably_prime(N), !. % Correct for all N < 341550071728321 squfof_1(N, _, F):- F is sqrt(N), F =:= ip(F), !. % N is a perfect square squfof_1(N, K, F):- KN is K * N, KN < 1.0E15, P0 is ip(sqrt(KN)), Q1 is KN - P0 * P0, reduce_1(P0, P0, 1, Q1, Pi, S0), % S0 = sqrt(Qi) B is (P0 + Pi) // S0, % Smallest B such that P0+(Pi-B*S0) < S0 R0 is B * S0 - Pi, S1 is (KN - R0 * R0) // S0, % S1 = (KN-R0*R0) div S0 reduce_2(P0, R0, S0, S1, Pj), gcd(N, Pj, F), 1 < F, F < N, % F is a non-trivial factor of N N mod F =:= 0. % Sanity check % e.g. N=1353,X=sqrt(N),reduce_1(X,36,1,57,P,Q) gives P=27, Q=4 reduce_1(X, P0, Q0, Q1, P, Q):- Q1 > 0, % Avoids division by 0 if KN is a square reduce(X, P0, Q0, Q1, P1, Q2), B2 is (X + P1) // Q2, P2 is B2 * Q2 - P1, Q3 is Q1 + B2 * (P1 - P2), reduce_1_1(X, P1, P2, Q2, Q3, P, Q). reduce_1_1(_, _, P, Q2, _, P, Q):- Q is sqrt(Q2), Q =:= ip(Q), !. % Q2 is a perfect square reduce_1_1(X, P1, P2, Q2, Q3, P, Q):- P2 \= P1, reduce_1(X, P2, Q2, Q3, P, Q). reduce_2(X, P0, Q0, Q1, P):- reduce(X, P0, Q0, Q1, P1, Q2), !, reduce_2(X, P1, Q1, Q2, P). reduce_2(_, P, _, _, P). reduce(X, P0, Q0, Q1, P1, Q2):- B1 is (X + P0) // Q1, P1 is B1 * Q1 - P0, P1 \= P0, Q2 is Q0 + B1 * (P0 - P1). % is_probably_prime(N) is true if N is probably prime. For all N less than % 341550071728321, if N is prime the procedure always succeeds, and if N % is composite the procedure always fails. % e.g. is_probably_prime(314159) succeeds % is_probably_prime(31415926535897) fails is_probably_prime(2):-!. is_probably_prime(N):- 1 =:= N mod 2, is_probably_prime_1([2,3,5,7,11,13,17], N). % For each prime in As... is_probably_prime_1([], _). is_probably_prime_1([N|As], N):-!, % Otherwise the primes in As will not be found as primes is_probably_prime_1(As, N). is_probably_prime_1([A|As], N):- is_sprp(N, A), is_probably_prime_1(As, N). % is_sprp(N, B) is true if N is odd, and N is a strong probable-prime base B. % [If we write n-1 = 2^s.d where d is odd and s is nonnegative, then % n is a strong probable-prime base b (a b-SPRP) if either b^d=1(mod n) or % (b^d)^(2r)=-1(mod n) for some nonnegative r less than s.] % e.g. is_sprp(2047, 2). is_sprp(N, B):- N > 1, 1 =:= N mod 2, N_1 is N - 1, is_sprp_1(N_1, D, 0, S), % N = 2^S*D+1 B0 is B mod N, powmod(B0, D, N, X), % X = B^D mod N is_sprp_2(0, S, N, X). % is_sprp_1(N, D, 0, S) is true if N = 2^S*D, and D is odd. is_sprp_1(N, D, S0, S):- 0 =:= N mod 2, !, N1 is N // 2, S1 is S0 + 1, is_sprp_1(N1, D, S1, S). is_sprp_1(N, N, S, S). % Succeeds if N is probably prime. % For each R in the range [0, S-1]... is_sprp_2(0, _, _, 1):-!. % B^D = 1(mod N) is_sprp_2(_, _, N, X):- X is N - 1, !. % (B^D)^(2*R) = -1(mod N) is_sprp_2(R, S, N, X):- R1 is R + 1, R1 < S, X0 is X mod N, multmod(X0, X0, N, X1), is_sprp_2(R1, S, N, X1). % gcd(X, Y, Z) is true if Z is the greatest common divisor of X and Y. gcd(X, 0, X):-!. gcd(X, Y, Z):-X1 is abs(X), Y1 is abs(Y), gcd_1(X1, Y1, Z). gcd_1(X, 0, Z):-!, Z=X. gcd_1(X, Y, Z):-R is X mod Y, gcd_1(Y, R, Z). % for(I, J, K) is true if K is an integer between I and J inclusive. for(I, J, I):-I =< J. for(I, J, K):-I < J, I1 is I + 1, for(I1, J, K). % rho(A, B, P, X) is true if A is a primitive root of the prime P, B is a % value between 1 and P-1, and X is the solution to A^X = B mod P. % This is Pollard's Rho Discrete Logarithm algorithm. % e.g. rho(31415, 92653, 589793, X) gives X=180541 rho(A, B, P, X):- 0 < A, A < P, 0 < B, B < P, Q is P - 1, rho_0(A, B, P, Q, X). % rho(A, B, P, Q, X) is true if X is the solution to A^X = B mod P, where % A is a generator of a subgroup of prime order Q of z_star(P), and B is an % element of the subgroup. % e.g. rho(2, 228, 383, 191, X) gives X=110 rho(A, A, _, _, X):-!, % Required when called from ph_2/14 X=1. rho(P_1, 1, P, _, X):- % Required when called from ph_2/14 P_1 is P - 1, !, X=2. rho(A, B, P, Q, X):- 0 < A, A < P, 0 < B, B < P, rho_0(A, B, P, Q, X). % Initialisation... rho_0(A, B, P, Q, X):- % Solutions are not always obtained if A0=B0=0, so... for(0, 10, A0), for(0, 10, B0), powmod(A, A0, P, C), powmod(B, B0, P, D), multmod(C, D, P, X0), rho_1(0, A, B, P, Q, X0, A0, B0, X0, A0, B0, X), powmod(A, X, P, B), !. % Sanity check % Loops until Xi=X2i, but does not succeed at the first call. % e.g. rho_1(0, 2, 228, 383, 191, 1, 0, 0, 1, 0, 0, X) gives X=110 rho_1(1, A, B, P, Q, Xi, Ai, Bi, Xi, A2i, B2i, X):-!, subtmod(B2i, Bi, Q, F), % f=b2i-bi mod p-1 F > 0, subtmod(Ai, A2i, Q, E), % e=ai-a2i mod p-1 gcd_extended(F, Q, U0, _, D), % gcd(f,p-1)=u*f+v*(p-1)=d make_positive(U0, Q, U), multmod(U, E, Q, UE), UED is UE // D, powmod(A, UED, P, S), % s=a^(u*e/d) mod p QD is Q // D, powmod(A, QD, P, T), % t=a^((p-1)/d) mod p find_j(S, B, 0, D, T, P, J), % Finds smallest j such that s*t^j=b mod p X is UED + J * QD. % x=(u*e/d)+j*((p-1)/d) rho_1(_, A, B, P, Q, Xi, Ai, Bi, X2i, A2i, B2i, X):- rho_2(A, B, P, Q, Xi, Ai, Bi, Xi1, Ai1, Bi1), rho_2(A, B, P, Q, X2i, A2i, B2i, X2i0, A2i0, B2i0), rho_2(A, B, P, Q, X2i0, A2i0, B2i0, X2i1, A2i1, B2i1), rho_1(1, A, B, P, Q, Xi1, Ai1, Bi1, X2i1, A2i1, B2i1, X). rho_2(_, B, P, Q, Xi, Ai, Bi, Xi1, Ai, Bi1):- % ai1=ai Xi mod 3 =:= 1, !, multmod(B, Xi, P, Xi1), % xi1=xi*b mod p addmod(Bi, 1, Q, Bi1). % bi1=bi+1 mod p-1 rho_2(_, _, P, Q, Xi, Ai, Bi, Xi1, Ai1, Bi1):- Xi mod 3 =:= 0, !, multmod(Xi, Xi, P, Xi1), % xi1=xi^2 mod p addmod(Ai, Ai, Q, Ai1), % ai1=2*ai mod p-1 addmod(Bi, Bi, Q, Bi1). % bi1=2*bi mod p-1 rho_2(A, _, P, Q, Xi, Ai, Bi, Xi1, Ai1, Bi):- % bi1=bi multmod(A, Xi, P, Xi1), % xi1=xi*a mod p addmod(Ai, 1, Q, Ai1). % ai1=ai+1 mod p-1 % Finds the smallest j, where 0<=j<=d-1, such that s*t^j = b mod p, that is, % a^(u*e/d)*(a^((p-1)/d))^j = b mod p find_j(B, B, J, _, _, _, J):-!. find_j(S, B, J0, D, T, P, J):- multmod(S, T, P, S1), J1 is J0 + 1, J1 < D, find_j(S1, B, J1, D, T, P, J). make_positive(L0, Q, L):-L0 < 0, !, L is Q + L0. make_positive(L, _, L).