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