An example using this package can be found in Square Roots Modulo a Prime.
% powmod(A, B, N, C) is true if 0 <= A < N and C = A^B mod N. powmod(A, B, N, C):-B >= 0, N > 0, powmod_1(B, A, N, 1, C). powmod_1(0, _, _, C0, C):-!, C = C0. powmod_1(B, A, N, C0, C):- B mod 2 =:= 1, !, B1 is B // 2, multmod(A, A, N, A1), multmod(A, C0, N, C1), powmod_1(B1, A1, N, C1, C). powmod_1(B, A, N, C0, C):- B1 is B // 2, multmod(A, A, N, A1), powmod_1(B1, A1, N, C0, C). % multmod(A, B, N, C) is true if 0 <= A < N and C = A*B mod N. multmod(A, B, N, C):-B >= 0, N > 0, multmod_1(B, A, N, 0, C). multmod_1(0, _, _, C0, C):-!, C = C0. multmod_1(B, A, N, C0, C):- B mod 2 =:= 1, !, B1 is B // 2, addmod(A, A, N, A1), addmod(A, C0, N, C1), multmod_1(B1, A1, N, C1, C). multmod_1(B, A, N, C0, C):- B1 is B // 2, addmod(A, A, N, A1), multmod_1(B1, A1, N, C0, C). % addmod(A, B, N, C) is true if 0 <= A,B < N and C = A+B mod N. addmod(A, B, N, C):- A < N - B, !, C is A + B. addmod(A, B, N, C):- C is A - (N - B). % subtmod(A, B, N, C) is true if 0 <= A,B < N and C = A-B mod N. subtmod(A, B, _, C):- A >= B, !, C is A - B. subtmod(A, B, N, C):- C is A + (N - B). % divmod(A, B, N, C) is true if 0 <= A,B < N and C = A/B mod N, where % B and N are mutually prime. % e.g. divmod(3, 11, 17, C) gives C=8 divmod(A, B, N, C):- inversemod(B, N, X), multmod(A, X, N, C). % inversemod(B, N, X) is true if X is the inverse of B modulo N, that is, % B * X = (1 modulo N), where B and N are mutually prime. % e.g. inversemod(1234, 567, X) gives X=550 inversemod(B, N, X):- B > 0, N > 0, gcd_extended(B, N, X0, _, 1), inversemod_1(X0, N, X). inversemod_1(X0, N, X):- X0 < 0, !, X is X0 + N. inversemod_1(X, _, X). % gcd_extended(X, Y, U1, U2, Z) is true if X*U1 + Y*U2 = Z = gcd(X, Y). % e.g. gcd_extended(24140, 16762, U1, U2, Z) gives U1=-234, U2=337, Z=34 gcd_extended(X, Y, U1, U2, Z):- gcd_extended_1(Y, 1, 0, X, 0, 1, Z, U2, U1). gcd_extended_1(0, _, _, U3, U2, U1, V3, V2, V1):-!, V3=U3, V2=U2, V1=U1. gcd_extended_1(V3, V2, V1, U3, U2, U1, U3P, U2P, U1P):- Q is U3 // V3, W3 is U3 mod V3, W2 is U2 - V2 * Q, W1 is U1 - V1 * Q, gcd_extended_1(W3, W2, W1, V3, V2, V1, U3P, U2P, U1P).