#### Square Roots Modulo the Product of Two Primes

This code requires the Square Roots Modulo a Prime and Modular Arithmetic packages.

```% sqrt_pq(A, P, Q, R1, R2, R3, R4) is true if P and Q are distinct prime
%   numbers and R1, R2, R3 and R4 are the square roots of A modulo P*Q.
%   (Calculate x0 from x0^2=a mod p and calculate x1 from x1^2=a mod q.
%   Then solve x=x0 mod p and x=x1 mod q (and also x'=(n-x0) mod p and
%   x'=x1 mod q) using the Chinese Remainder Theorem.)
% e.g. sqrt_pq(16, 13, 29, R1, R2, R3, R4)
%   gives R1=373, R2=4, R3=199, R4=178
sqrt_pq(A, P, Q, R1, R2, R3, R4):-
N is P * Q,
X0_2 is A mod P,
X1_2 is A mod Q,
sqrt_prime(X0_2, P, X0),
sqrt_prime(X1_2, Q, X1),
N_X0 is (N - X0) mod P,
crt([X0,  X1], [P,Q], R1),
crt([N_X0,X1], [P,Q], R3),
R2 is N - R1,
R4 is N - R3.

% crt(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 is the "Chinese Remainder Theorem".
% e.g. crt([3,4,5], [125,17,12], X) gives X=23753
crt(As, Ms, X):-
are_relatively_prime(Ms),
mod_mixed(As, Ms, Vs),
mixed_int(Vs, Ms, X).

% are_relatively_prime(Xs) is true if the elements of Xs are relatively
%   prime to each other in pairs.
% e.g. are_relatively_prime([125,17,12]).
are_relatively_prime([]).
are_relatively_prime([X|Xs]):-
are_relatively_prime_1(Xs, X),
are_relatively_prime(Xs).

are_relatively_prime_1([], _).
are_relatively_prime_1([X2|Xs], X1):-
gcd(X1, X2, 1),
are_relatively_prime_1(Xs, 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).

% mod_mixed(Us, Ms, Vs) is true if the list Vs is a mixed-radix
%   representation of the modular representation Us with respect to the
%   list of moduli Ms.
%   [See Knuth 4.3.2(24)]
% e.g. mod_mixed([9,25], [13,29], Vs) gives Vs=[9,28]
mod_mixed([], [], []).
mod_mixed([Uj|Us], [Mj|Ms], [Uj|Vs]):-
mod_mixed_1(Us, Ms, Uj, Mj, Us1),
mod_mixed(Us1, Ms, Vs).

mod_mixed_1([], [], _, _, []).
mod_mixed_1([Uk|Us], [Mk|Ms], Uj0, Mj, [Vk|Vs]):-
inversemod(Mj, Mk, Cjk),
Uj is Uj0 mod Mk,
subtmod(Uk, Uj, Mk, Temp2),
multmod(Temp2, Cjk, Mk, Vk),
mod_mixed_1(Us, Ms, Uj0, Mj, Vs).

% mixed_int(Vs, Ms, U) is true if U is the integer value corresponding to
%   the mixed-radix representation Vs with respect to the list of
%   moduli Ms.
%   [See Knuth 4.3.2(25)]
% e.g. mixed_int([4,15], [13,29], U) gives U=199
mixed_int(Vs, Ms, U):-mixed_int_1(Vs, Ms, 1, 0, U).

mixed_int_1([], [], _, U, U).
mixed_int_1([V|Vs], [M|Ms], Product0, U0, U):-
Product1 is Product0 * M,
U1 is U0 + Product0 * V,
mixed_int_1(Vs, Ms, Product1, U1, U).
```