#### Polynomial Operations Modulo a Prime

These procedures require the Modular arithmetic package.

```/*
A polynomial is represented by a list of its coefficients, high-order first
e.g. x^3 + 3*x^2 - 2*x + 4 is represented by [1,3,-2,4].
*/

% p_addmod(Xs, Ys, P, Zs) is true if adding the polynomial Xs to the
%   polynomial Ys modulo P gives the polynomial Zs.
reverse(Xs, ReversedXs),
reverse(Ys, ReversedYs),
reverse(ReversedZs, Zs0),
trim_left(Zs0, Zs).

% p_subtmod(Xs, Ys, P, Zs) is true if subtracting the polynomial Ys from
%   the polynomial Xs modulo P gives the polynomial Zs.
p_subtmod(Xs, Ys, P, Zs):-
reverse(Xs, ReversedXs),
reverse(Ys, ReversedYs),
p_subtmod_1(ReversedXs, ReversedYs, P, ReversedZs),
reverse(ReversedZs, Zs0),
trim_left(Zs0, Zs).

p_subtmod_1([], [], _, []):-!.
p_subtmod_1([], [Y|Ys], P, [Z|Zs]):-!,
Z is P - Y,
p_subtmod_1([], Ys, P, Zs).
p_subtmod_1([X|Xs], [], _, [X|Xs]):-!.
p_subtmod_1([X|Xs], [Y|Ys], P, [Z|Zs]):-
subtmod(X, Y, P, Z),
p_subtmod_1(Xs, Ys, P, Zs).

% p_multmod(Xs, Ys, P, Zs) is true if multiplying the polynomial Xs by the
%   polynomial Ys modulo P gives the polynomial Zs.
p_multmod(Xs, Ys, P, Zs):-
reverse(Ys, ReversedYs),
p_multmod_1(Xs, ReversedYs, P, [], ReversedZs),
reverse(ReversedZs, Zs0),
trim_left(Zs0, Zs).

p_multmod_1([], _, _, Zs, Zs).
p_multmod_1([X|Xs], Ys, P, Zs0, Zs):-
p_multmod_2(Ys, X, P, Zs1),
p_multmod_1(Xs, Ys, P, Zs2, Zs).

p_multmod_2([], _, _, []):-!.
p_multmod_2(_, 0, _, [0]):-!.
p_multmod_2([Y|Ys], X, P, [Z|Zs]):-
multmod(X, Y, P, Z),
p_multmod_2(Ys, X, P, Zs).

% p_divmod(Xs, Ys, P, Qs, Rs) is true if Qs and Rs are the quotient and
%   remainder when dividing the polynomial Xs by the polynomial Ys
%   modulo P.
% [Knuth 4.6.1 Algorithm D]
% e.g. p_divmod([1,0,1,0,10,10,8,2,8], [3,0,5,0,9,4,8], 13, Qs, Rs) gives
%   Qs=[9,0,7], Rs=[11,0,3,0,4]
p_divmod(_, [0|_], _, [0], [0]):-!,
write(`Error - attempted division by zero`), nl, fail.
p_divmod(Xs, [Y|Ys], P, Qs, Rs):-
inversemod(Y, P, InvY),
p_divmod_1(Xs, InvY, Ys, P, Qs0, Rs0),
trim_left(Qs0, Qs),
trim_left(Rs0, Rs).

p_divmod_1([X|Xs], InvY, Ys, P, [Q|Qs], Rs):-
multmod(X, InvY, P, Q),
p_mul_1(Ys, Q, P, Ys1),
p_subt_1(Ys1, Xs, P, Xs1), !,
p_divmod_1(Xs1, InvY, Ys, P, Qs, Rs).
p_divmod_1(Xs, _, _, _, [], Xs).

p_mul_1([], _, _, []).
p_mul_1([Y|Ys], X, P, [Z|Zs]):-
multmod(X, Y, P, Z),
p_mul_1(Ys, X, P, Zs).

% Fails if Xs is longer than Ys
p_subt_1([], Ys, _, Ys).
p_subt_1([X|Xs], [Y|Ys], P, [Z|Zs]):-
subtmod(Y, X, P, Z),
p_subt_1(Xs, Ys, P, Zs).

% p_powmod(Xs, I, P, Ys) is true if raising the polynomial Xs to the
%   integral power I modulo P gives the polynomial Ys.
% e.g. p_powmod([3,1,4], 5, 7, Ys) gives Ys=[5,6,0,3,4,4,3,3,0,6,2]
p_powmod(Xs, I, P, Ys):-p_powmod_1(I, Xs, P, [1], Ys).

p_powmod_1(0, _, _, Ys, Ys):-!.
p_powmod_1(I, Xs, P, Ys0, Ys):-
I1 is I - 1,
p_multmod(Ys0, Xs, P, Ys1),
p_powmod_1(I1, Xs, P, Ys1, Ys).

% p_powermod(Xs, N, Ys, P, Zs) is true if Zs is the result of raising the
%   polynomial Xs to the power N modulo the polynomial Ys, arithmetic being
%   performed modulo P.
% e.g. p_powermod([1,6], 4, [1,9,4,6,5], 7, Zs) gives Zs=[1,2,4,3]
p_powermod(Xs, N, Ys, P, Zs):-p_powermod_1(N, Xs, Ys, P, [1], Zs).

p_powermod_1(0, _, _, _, Zs, Zs):-!.
p_powermod_1(N, Xs, Ys, P, Zs0, Zs):-
1 =:= N mod 2, !,
N1 is N // 2,
p_multmod(Xs, Xs, P, X2s),
p_divmod(X2s, Ys, P, _, Xs1),
p_multmod(Zs0, Xs, P, ZXs),
p_divmod(ZXs, Ys, P, _, Zs1),
p_powermod_1(N1, Xs1, Ys, P, Zs1, Zs).
p_powermod_1(N, Xs, Ys, P, Zs0, Zs):-
N1 is N // 2,
p_multmod(Xs, Xs, P, X2s),
p_divmod(X2s, Ys, P, _, Xs1),
p_powermod_1(N1, Xs1, Ys, P, Zs0, Zs).

% p_gcdmod(Xs, Ys, P, Zs) is true if the polynomial Zs is the greatest
%   common divisor of the polynomials Xs and Ys modulo P.
% e.g. p_gcdmod([1,0,1,0,10,10,8,2,8], [3,0,5,0,9,4,8], 13, Zs) gives Zs=[1]
%      p_gcdmod([1,0,5], [1,0,0,6], 7, Zs) gives Zs=[1,3]
p_gcdmod(Xs, [0], P, Ys):-!,
monic(Xs, P, Ys).
p_gcdmod(_, [_], _, [1]):-!.
p_gcdmod(Xs, Ys, P, Zs):-
p_divmod(Xs, Ys, P, _, Rs),
p_gcdmod(Ys, Rs, P, Zs).

% p_gcdmodext(X, Y, P, A, B, Z) is true if X*A + Y*B = Z = gcd(X, Y),
%   where X, Y, A, B and Z are polynomials, all operations being modulo P.
% e.g. p_gcdmodext([1,0,1,0,10,10,8,2,8], [3,0,5,0,9,4,8], 13, A, B, Z).
%        gives A=[9,1,5,2,9,0], B=[10,4,9,1,7,4,8,5], Z=[1]
%      p_gcdmodext([1,0,5], [1,0,0,6], 7, A, B, Z).
%        gives A=[3,0], B=[4], Z=[1,3]
p_gcdmodext(X, Y, P, A, B, Z):-
p_gcdmodext_1(Y, [1], [0], P,  X, [0], [1],  Z0, B0, A0),
Z0=[W|_],
dividemod(A0, W, P, A),
dividemod(B0, W, P, B),
monic(Z0, P, Z).

p_gcdmodext_1([0], _, _, _,  U3, U2, U1,  V3, V2, V1):-!,
V3=U3, V2=U2, V1=U1.
p_gcdmodext_1(V3, V2, V1, P,  U3, U2, U1, U3P, U2P, U1P):-
p_divmod(U3, V3, P, Q, W3),
p_multmod(V2, Q, P, T1),
p_subtmod(U2, T1, P, W2),
p_multmod(V1, Q, P, T2),
p_subtmod(U1, T2, P, W1),
p_gcdmodext_1(W3, W2, W1, P,  V3, V2, V1,  U3P, U2P, U1P).

% Divides the elements of Xs by Y mod P
dividemod([], _, _, []).
dividemod([X|Xs], Y, P, [Z|Zs]):-
divmod(X, Y, P, Z),
dividemod(Xs, Y, P, Zs).

% normalise(Ux, P, Wx) is true if the polynomial Wx contains the elements
%   of the polynomial Ux modulo P.
normalise(Ux, P, Wx):-
normalise_1(Ux, P, Wx0),
trim_left(Wx0, Wx).

normalise_1([], _, []).
normalise_1([U|Ux], P, [W|Wx]):-
W is U mod P,
normalise_1(Ux, P, Wx).

% monic(Wx, P, Vx) is true if Vx is the polynomial Wx divided (modulo P) by
%   the leading coefficient of Wx.
monic([W|Wx], P, [1|Vx]):-
P > 0,
monic_1(Wx, W, P, Vx).

monic_1([], _, _, []).
monic_1([W|Wx], W0, P, [V|Vx]):-
divmod(W, W0, P, V),
monic_1(Wx, W0, P, Vx).

% trim_left(Xs, Zs) is true if Zs is the list Xs with all leading zeroes
%   removed, except that [] is replaced by [0].
trim_left([], [0]):-!.
trim_left(Xs, Zs):-
trim_left_1(Xs, Zs).

trim_left_1([0|Xs], Zs):-
Xs=[_|_], !,
trim_left_1(Xs, Zs).
trim_left_1(Xs, Xs).

% reverse(Xs, Ys) is true if Ys is the result of reversing the order of
%   the elements in the list Xs.
%reverse(Xs, Ys):-reverse_1(Xs, [], Ys).

%reverse_1([], As, As).
%reverse_1([X|Xs], As, Ys):-reverse_1(Xs, [X|As], Ys).
```