Factorization of Polynomials Modulo a Prime

These procedures require the Modular arithmetic and Polynomial Operations Modulo a Prime packages.

/* fact_m(Ux, P, Fxs) is true if the elements of Fxs are the factors modulo  */
/*   the prime P of the polynomial Ux.  See Knuth 4.6.2 Algorithm D.         */
/* e.g. fact_m([1,0,1,0,10,10,8,2,8], 13, X) gives                           */
/*   X=[[1,2,3,4,6],[1,8,4,12],[1,3]], that is,                              */
/*   x^8+x^6+10x^4+10x^3+8x^2+2x+8=(x^4+2x^3+3x^2+4x+6)(x^3+8x^2+4x+12)(x+3) */
/*     modulo 13.                                                            */
fact_m([U], P, [[Fx]]):-!,
  is_prime(P),
  Fx is U mod P,
  Fx > 0.
fact_m(Ux, P, Fxs):-
  normalise(Ux, P, Wx),
  Wx \= [0],
  monic(Wx, P, Vx),
  is_prime(P),
  squarefree(Vx, P, [], Fxs0),
  Wx=[Fx|_],
  avoid_1(Fx, Fxs0, Fxs),
  product(Fxs, P, [1], Wx). % Sanity check

% Avoids returning [1] as a factor as a result of making the polynomial monic
avoid_1(1, Fxs, Fxs):-!.
avoid_1(Fx, Fxs, [[Fx]|Fxs]).

% Let g(x)=gcd(u(x),u'(x)). Then
%   g(x)=1 => u(x) squarefree => factorize u(x)
%   g(x)=u(x) => u'(x)=0 => the coefficient of x^k is nonzero only when k is a
%     multiple of p, so that u(x) can be rewritten as v(x^p). Then 
%     u(x)=(v(x))^p. Factorize v(x) and replicate its factors p times.
%   g(x)!=u(x) => g(x) is a factor of u(x) => factorize g(x) and v(x)=u(x)/g(x)
squarefree(Ux, P, Fxs0, Fxs):-
  derivative(Ux, P, Dx),
  p_gcdmod(Ux, Dx, P, Gx),                % gcd(u(x), u'(x))
  squarefree_1(Gx, Ux, P, Fxs0, Fxs).

squarefree_1([1], Ux, P, Fxs0, Fxs):-!,   % u(x) is squarefree - factorize u(x)
  fact_m_1(0, Ux, [1,0], P, Fxs0, Fxs).
squarefree_1(Ux, Ux, P, Fxs0, Fxs):-!,    % g(x)=u(x) - rewrite u(x) as v(x^p)
  rewrite(Ux, P, P, Vx),
  squarefree(Vx, P, Fxs0, Fxs1),          % Factorize v(x)
  replicate(Fxs1, P, [], Fxs).
squarefree_1(Gx, Ux, P, Fxs0, Fxs):-      % g(x)!=u(x) - g(x) divides u(x)
  squarefree(Gx, P, Fxs0, Fxs1),          % Factorize g(x)
  p_divmod(Ux, Gx, P, Vx, [0]),
  squarefree(Vx, P, Fxs1, Fxs).           % Factorize v(x)=u(x)/g(x)

% e.g. rewrite([1,0,0,2,0,0,2], 3, 3, [1,2,2]).
rewrite([], _, _, []):-!.
rewrite([U|Us], P, P, [U|Vs]):-!,
  rewrite(Us, 1, P, Vs).
rewrite([_|Us], I, P, Vs):-
  I1 is I + 1,
  rewrite(Us, I1, P, Vs).
  
% e.g. replicate([[1,2],[3,4]], 3, [], [[3,4],[3,4],[3,4],[1,2],[1,2],[1,2]]).
replicate([], _, Yss, Yss).
replicate([Xs|Xss], P, Yss0, Yss):-
  replicate_1(P, Xs, Yss0, Yss1),
  replicate(Xss, P, Yss1, Yss).

replicate_1(0, _, Yss, Yss):-!.
replicate_1(I, Xs, Yss0, Yss):-
  I1 is I - 1,
  replicate_1(I1, Xs, [Xs|Yss0], Yss).

% Continue the factorization process  
fact_m_1(_, [1], _, _, Fxs, Fxs):-!.      % v(x)=1
fact_m_1(D, Vx, _, _, Fxs, [Vx|Fxs]):-
  length_1(Vx, -1, Deg),
  D + 1 > Deg // 2, !.                    % d+1>deg(v(x))/2
fact_m_1(D, Vx, Ws, P, Fxs0, Fxs):-
  p_powermod(Ws, P, Vx, P, Ws1),          % w(x)=w(x)^p modulo v(x)
  p_subtmod(Ws1, [1,0], P, Ws2),          % w(x)-x
  p_gcdmod(Ws2, Vx, P, Gx),               % g(x)=gcd(w(x)-x,v(x))
  D1 is D + 1,
  fact_m_2(Gx, D1, Vx, Ws1, P, Fxs0, Fxs).
  
% Continue the factorization process  
fact_m_2([1], D, Vx, Ws, P, Fxs0, Fxs):-!,% g(x)=1
  fact_m_1(D, Vx, Ws, P, Fxs0, Fxs).
fact_m_2(Gx, D, Vx, Ws, P, Fxs0, Fxs):-   % g(x)!=1
  length_1(Gx, -1, Deg),
  Deg =< D, !,                            % deg(g(x))<=d
  p_divmod(Vx, Gx, P, Vx1, [0]),          % v(x)=v(x)/g(x)
  p_divmod(Ws, Vx1, P, _, Ws1),           % w(x)=w(x) modulo v(x)
  fact_m_1(D, Vx1, Ws1, P, [Gx|Fxs0], Fxs).
fact_m_2(Gx, D, Vx, Ws, 2, Fxs0, Fxs):-!, % g(x)!=1, deg(g(x))>d
  irred_2(Gx, D, Fxs0, Fxs1),
  p_divmod(Vx, Gx, 2, Vx1, [0]),          % v(x)=v(x)/g(x)
  p_divmod(Ws, Vx1, 2, _, Ws1),           % w(x)=w(x) modulo v(x)
  fact_m_1(D, Vx1, Ws1, 2, Fxs1, Fxs).
fact_m_2(Gx, D, Vx, Ws, P, Fxs0, Fxs):-   % g(x)!=1, deg(g(x))>d, p>2
  irred_p(Gx, D, P, Fxs0, Fxs1),
  p_divmod(Vx, Gx, P, Vx1, [0]),          % v(x)=v(x)/g(x)
  p_divmod(Ws, Vx1, P, _, Ws1),           % w(x)=w(x) modulo v(x)
  fact_m_1(D, Vx1, Ws1, P, Fxs1, Fxs).
  
% Gx is the product of all the irreducible factors of degree D of u(x)
% Returns all irreducible factors of degree D of Gx for P=2 
irred_2(Gx, D, Fxs, [Gx|Fxs]):-    
  length_1(Gx, -1, D), !.                 % g(x) is an irreducible factor
irred_2(Gx, D, Fxs0, Fxs):-
  D1 is D - 1,
  D2 is 2 * D - 1,
  repeat,
    rand_poly(D2, 2, Tx),                 % t(x)
    irred_2a(D1, Gx, Tx, Tx, Tx1),        % T(t(x)) modulo g(x)
    p_gcdmod(Gx, Tx1, 2, Jx),             % gcd(g(x), T(t(x)))
  Jx \= [1], Jx \= Gx, !,                 % j(x) is a factor of g(x)
  p_divmod(Gx, Jx, 2, Kx, [0]),           % k(x) is a factor of g(x)
  irred_2(Jx, D, Fxs0, Fxs1),
  irred_2(Kx, D, Fxs1, Fxs).

% Calculates T(t(x)) modulo g(x) where T(x)=x+x^p+...+x^(2^(p-1)). Here p=2.
irred_2a(0, _, _, Us, Us):-!.
irred_2a(I, Gx, Tx, Us0, Us):-
  p_multmod(Us0, Us0, 2, Us1),            % u(x)^2 modulo 2
  p_addmod(Tx, Us1, 2, Us2),              % t(x)+u(x)^2 modulo 2
  p_divmod(Us2, Gx, 2, _, Us3),           % (t(x)+u(x)^2) modulo g(x)
  I1 is I - 1,
  irred_2a(I1, Gx, Tx, Us3, Us).  

% Gx is the product of all the irreducible factors of degree D of u(x)
% Returns all irreducible factors of degree D of Gx for P>2 
irred_p(Gx, D, _, Fxs, [Gx|Fxs]):-    
  length_1(Gx, -1, D), !.                 % g(x) is an irreducible factor
irred_p(Gx, D, P, Fxs0, Fxs):-
  power(P, D, PD),
  PD12 is (PD - 1) // 2,                  % (p^d-1)/2
  D1 is 2 * D - 1,
  repeat,
    rand_poly(D1, P, Tx),                 % t(x)
    p_powermod(Tx, PD12, Gx, P, Tx1),     % t(x)^((p^d-1)/2) modulo g(x)
    p_subtmod(Tx1, [1], P, Tx2),          % t(x)^((p^d-1)/2)-1 modulo g(x)
    p_gcdmod(Gx, Tx2, P, Jx),             % gcd(g(x), t(x)^((p^d-1)/2)-1)
  Jx \= [1], Jx \= Gx, !,                 % j(x) is a factor of g(x)
  p_divmod(Gx, Jx, P, Kx, [0]),           % k(x) is a factor of g(x)
  irred_p(Jx, D, P, Fxs0, Fxs1),
  irred_p(Kx, D, P, Fxs1, Fxs).

% Calculates the product modulo P of the given polynomials 
product([], _, Wx, Wx):-!.
product([Fx|Fxs], P, Wx0, Wx):-!,
  normalise(Fx, P, Fx1),
  p_multmod(Wx0, Fx1, P, Wx1),
  product(Fxs, P, Wx1, Wx).
product(_, _, _, _):-  
  write(`Sanity check failed`), nl,
  fail.

% Generates a pseudo-random monic, non-constant polynomial of degree at most D
%   with coefficients less than P.
rand_poly(0, _, [1]):-!.
rand_poly(D, P, [1|Xs]):-
  D > 0,
  rand_int(D, D1),
  D2 is D1 + 1,
  rand_poly_1(D2, P, Xs).

rand_poly_1(0, _, []):-!.
rand_poly_1(D, P, [X|Xs]):-
  rand_int(P, X),
  D1 is D - 1,
  rand_poly_1(D1, P, Xs).

/* rand_int(I, J) is true if J is a pseudo-random integer in the range       */
/*   0 to I-1.                                                               */
rand_int(I, J):-J is int(rand(I)).

/* derivative(Vx, P, Dx) is true if Dx is the result of differentiating      */
/*   the polynomial Vx modulo P.                                             */
/* e.g. derivative([8,7,6,5,4,3,2,1,0], 10, [4,9,6,5,6,9,4,1]).              */
derivative([V|Vx], P, Dx):-
  length(Vx, Length),
  derivative_1(Vx, V, Length, P, Dx0),
  trim_left(Dx0, Dx).
  
derivative_1([], _, _, _, []).
derivative_1([V|Vx], V0, I, P, [D|Dx]):-
  ImodP is I mod P,
  multmod(ImodP, V0, P, D),
  I1 is I - 1,
  derivative_1(Vx, V, I1, P, Dx).

/* is_prime(P) is true if P is a prime number.                               */
is_prime(2):-!.
is_prime(P):-
  P > 2,
  P mod 2 =:= 1,
  SqrtP is sqrt(P),
  is_prime_1(3, P, SqrtP).

is_prime_1(I, _, SqrtP):-
  I > SqrtP, !.
is_prime_1(I, P, SqrtP):-
  P mod I > 0,
  I1 is I + 2,
  is_prime_1(I1, P, SqrtP).

/* length(Xs, L) is true if L is the number of elements in the list Xs.      */
%length(Xs, L):-length_1(Xs, 0, L).

length_1([], L, L).
length_1([_|Xs], L0, L):-L1 is L0 + 1, length_1(Xs, L1, L).

/* power(X, N, Y) is true if Y is the result of raising X to the power N.    */
power(X, N, Y):-
  N >= 0, 
  power_1(N, X, 1, Y).

power_1(0, _, Y, Y):-!.
power_1(1, X, Y0, Y):-!,
  Y is Y0 * X.
power_1(N, X, Y0, Y):-
  1 =:= N mod 2, !,
  N1 is N // 2,
  X1 is X * X,
  Y1 is Y0 * X,
  power_1(N1, X1, Y1, Y).
power_1(N, X, Y0, Y):-
  N1 is N // 2,
  X1 is X * X,
  power_1(N1, X1, Y0, Y).

/* repeat/0 always succeeds.                                                 */
%repeat.
%repeat:-repeat.

/*
 * Procedures to display polynomials prettily
 */  

% Display all the polynomials in the list
% e.g. show_polys([[1,-1],[1,1,1],[1,1,1,1,1],[1,-1,0,1,-1,1,0,-1,1]]) displays
%   x-1
%   x^2+x+1
%   x^4+x^3+x^2+x+1
%   x^8-x^7+x^5-x^4+x^3-x+1
show_polys([]).
show_polys([P|Ps]):-
  show_poly(P),
  show_polys(Ps).

% Display the polynomial
show_poly([V|Vs]):-
  length(Vs, L), 
  show_sign_1(V), 
  show_value(V), 
  show_exponent(L), 
  L1 is L - 1, 
  show_poly(Vs, L1),
  nl.
  
show_poly([], _):-!.
show_poly([0], _):-!.
show_poly([V], _):-!, 
  show_sign(V), AbsV is abs(V), write(AbsV).
show_poly([0|Vs], L):-!,
  L1 is L - 1, show_poly(Vs, L1).
show_poly([V|Vs], L):-
  show_sign(V), show_value(V), show_exponent(L), L1 is L-1, show_poly(Vs, L1).

show_sign_1(V):-V >= 0, !.  
show_sign_1(_):-write(`-`).

show_sign(0):-!.  
show_sign(V):-V > 0, write(`+`), !.  
show_sign(_):-write(`-`).
  
show_value(1):-!.
show_value(-1):-!.
show_value(V):-AbsV is abs(V), write(AbsV).
  
show_exponent(0):-!.
show_exponent(1):-!, write(`x`).
show_exponent(L):-write(`x^`), write(L).  

LPA Index     Home Page

Valid HTML 4.01 Strict