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