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

LPA Index     Home Page

Valid HTML 4.0!