#### Modular Arithmetic Without Integer Overflow

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,
multmod_1(B1, A1, N, C1, C).
multmod_1(B, A, N, C0, C):-
B1 is B // 2,
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.
A < N - B, !,
C is A + B.
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).
```