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. p_addmod(Xs, Ys, P, Zs):- reverse(Xs, ReversedXs), reverse(Ys, ReversedYs), p_addmod_1(ReversedXs, ReversedYs, P, ReversedZs), reverse(ReversedZs, Zs0), trim_left(Zs0, Zs). p_addmod_1([], [], _, []):-!. p_addmod_1([], [X|Xs], _, [X|Xs]):-!. p_addmod_1([X|Xs], [], _, [X|Xs]):-!. p_addmod_1([X|Xs], [Y|Ys], P, [Z|Zs]):- addmod(X, Y, P, Z), p_addmod_1(Xs, Ys, P, 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_addmod_1([0|Zs0], Zs1, P, Zs2), 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).