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