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