Discrete Logarithm Problem

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

LPA Index     Home Page

Valid HTML 4.01 Strict