(*:Version: Mathematica 2.0 *)

(*:Name: NumberTheory`PolynomialMod` *)

(* :Title:  PolynomialMod *)

(* :Author: Ilan Vardi *)

(*:Keywords: Polynomial operations modulo a prime > 2^16.

*)


(* :Source: D.E. Knuth, ``Seminumerical Algorithms,'' Second Edition,
            Addison Wesley 1981.

*)

(* :Warning: This top level code is much slower than the built in 
             functions for primes < 2^16. For p = 2, this code is
             expecially inefficient.
             
*)


(* :Limitation: FactorMod can factor degree 20 polynomials modulo hundred
                digit primes in reasonable time. 
                These algorithms are for polynomials in one varable only.
*)  

(* :Discussion:  This package supplants the built-in polynomial functions
                 PolynomialGCD, FactorList, and Factor to primes larger 
                 than 2^16. As well it generalizes PolynomialQuotient and
                 PolynomialRemainder to arithmetic modulo a prime, and 
                 introduces the function PolynomialPowerMod, which computes
                 powers of a polynomial modulo a prime and another polynomial
                 efficiently.

    The factoring algorithm is described in the Knuth reference,
    Section 4.6.2: One first uses the derivative to factor the polynomial
    f(x) into squarefree factors. One computes the degree d factorization
    of f by taking GCD[f(x), x^(p^d) - x], for d = 1,2,.... This is 
    computed quickly by using the power mod algorithm, i.e., repeated
    squaring. Finally, one uses the probabilistic factoring algorithm of
    Cantor-Zassenhaus to factor these into irreducible degree d factors.

*)

BeginPackage["NumberTheory`PolynomialMod`",
	"NumberTheory`NumberTheoryFunctions`"]  

PolynomialQuotientMod::usage := "PolynomialQuotientMod[f, g, x, p]
computes PolynomialMod[PolynomialQuotient[f, g, x], p]."

PolynomialRemainderMod::usage := "PolynomialRemainderMod[f, g, x, p]
computes PolynomialMod[PolynomialRemainder[f, g, x], p], i.e., 
f modulo g and p, where g is a polynomial and p is a prime."

PolynomialGCDMod::usage := "PolynomialGCDMod[f, g, h,..., p] computes
PolynomialGCD[f, g, h,..., Modulus -> p], when p > 2^16."

PolynomialPowerMod::usage := "PolynomialPowerMod[f, n, {g, p}]
computes f^n modulo g and p, where g is a polynomial and p is a prime."

FactorMod::usage := "FactorMod[f, p] computes Factor[f, Modulus -> p], 
where p > 2^16."

FactorListMod::usage := "FactorListMod[f, p] computes 
FactorList[f, Modulus -> p] where p > 2^16."

Begin["`Private`"]
                        
PolynomialQuotientMod[f_, g_, x_, p_]:= 
Block[{fp = PolynomialMod[f, p], gp = PolynomialMod[g, p]}, 
       {fp, gp} = 
       PolynomialMod[
        PowerMod[Coefficient[gp, x, Exponent[gp, x]], -1, p] {fp, gp}, 
        p];
       PolynomialMod[PolynomialQuotient[fp, gp, x], p]
      ]

PolynomialRemainderMod[f_, g_, x_, p_]:= 
    PolynomialMod[PolynomialRemainder[PolynomialMod[f, p], 
                  PolynomialMod[g, p], x], p]

PolynomialGCDMod[f_, p_]:= PolynomialMod[f, p]

PolynomialGCDMod[f_, g_, p_]:= 
    polynomialgcdmod[PolynomialMod[f, p], PolynomialMod[g, p],
                     If[Variables[f] != {}, First[Variables[f]],
                                            First[Variables[g]]],
                     p]

PolynomialGCDMod[f_, g__, p_]:= 
     PolynomialGCDMod[f, PolynomialGCDMod[g, p], p]

polynomialgcdmod[f_, g_, x_, p_Integer]:= 
   polynomialgcdmod[g, f, x, p] /; Exponent[f, x] < Exponent[g, x] 

polynomialgcdmod[f_, g_, x_, p_Integer]:= f  /;  Exponent[g,x] == -Infinity

polynomialgcdmod[f_, g_, x_, p_Integer] := 
          Block[{fp = PolynomialMod[f, p], gp = PolynomialMod[g, p], 
                 q, r, monic},
                     monic = Coefficient[gp, x, Exponent[gp, x]];
                     If[monic != 1, 
                        gp = PolynomialMod[PowerMod[monic, -1, p] gp, p]];
                     q = PolynomialMod[PolynomialQuotient[fp, gp, x], p];
                     r = PolynomialMod[PolynomialRemainder[fp, gp, x], p];
                     If[Exponent[r, x] <= 0,
                        If[r == 0, Return[gp], Return[1]],
                     Return[polynomialgcdmod[gp, r, x, p]]
                     ]
             ]   /; Exponent[g, x] >= 0  


PolynomialPowerMod[f_, n_Integer, {g_, p_Integer}]:=
PolynomialPowerMod[f, n, First[Variables[g]], {g, p}]

PolynomialPowerMod[f_, n_Integer, x_, {g_, p_Integer}]:=
Fold[PolynomialRemainderMod[#1^2 #2, g, x, p] &,
     1, 
     Block[{prm = PolynomialRemainderMod[f, g, x, p]},
            If[# == 0, 1, prm]& /@ IntegerDigits[n, 2]]] /; n > 0


FactorMod[f_, p_]:= Times @@ (#[[1]]^#[[2]]& /@ FactorListMod[f, p])

FactorListMod[f_, p_]:= 
Block[{flm = factorlistmod[f, p]},
      {#, Count[flm, #]}& /@ Union[flm]]

factorlistmod[f_, p_]:= factorlistmod[f, First[Variables[f]], p]

factorlistmod[f_, x_, p_]:= 
Block[{fp = PolynomialMod[f, p], monic, factor, df},
       monic = Coefficient[fp, x, Exponent[fp, x]];
       If[monic != 1, 
          Return[
          Join[{monic}, 
               factorlistmod[
                 PolynomialMod[PowerMod[monic, -1, p] fp, p], x, p]]
                ]];
       df = PolynomialMod[D[fp, x], p]; 
       If[df === 0, 
          Return[Flatten[
                  Table[factorlistmod[f /. x -> x^(1/p), x, p], {p}]]]];
       factor = PolynomialGCDMod[fp, df, p];
       If[Exponent[factor, x] == 0, 
          factorsquarefree[fp, x, 1, x, p], 
          Join[factorlistmod[PolynomialQuotientMod[fp, factor, x, p], x,p],
               factorlistmod[factor, x, p]]]]

factorsquarefree[f_, x_, d_, power_, p_]:= {f} /; Exponent[f, x] == d

factorsquarefree[f_, x_, d_, xpd_, p_]:= 
Block[{power = PolynomialPowerMod[xpd, p, x, {f, p}], factor, q}, 
       factor = PolynomialGCDMod[f, power - x, p];
       If[Exponent[factor, x] == 0, 
          Return[factorsquarefree[f, x, d + 1, power, p]]];
       If[Degree[factor, x] == Degree[f, x], 
          Return[factordegree[f, x, d, p]]];
       q = PolynomialQuotientMod[f, factor, x, p];
       Join[factorsquarefree[q, x, d + 1,
                             PolynomialRemainderMod[power, q, x, p],
                             p], 
            factordegree[factor, x, d, p]]]

     
factordegree[f_, x_, d_, p_]:= {f} /; Exponent[f, x] == d

factordegree[f_, x_, d_, p_]:= 
  Block[{exp = Exponent[f, x], fp = PolynomialMod[f, p], 
         s, t, i, q, pd = p^d},
  While[True, 
        s = Plus @@ Table[Random[Integer,{0, p}] x^i, {i,0,exp-2}] + 
            x^(exp - 1); 
        t = PolynomialGCDMod[s, fp, p]; 
       If[0 < Exponent[t,x] < exp, Break[],
         t =  PolynomialGCDMod[fp, 
                PolynomialPowerMod[s, (pd-1)/2, x, {fp, p}] - 1, p];
            If[ 0 < Exponent[t, x] < exp, Break[]]
        ]]; 
    q = PolynomialQuotientMod[fp, t, x, p];
    t = If[Exponent[q, x] < Exponent[t, x], q, t];
    Join[factordegree[t, x, d, p], 
         factordegree[PolynomialQuotientMod[f, t, x, p], x, d, p]]
   ]

Protect[PolynomialQuotientMod, PolynomialRemainderMod, PolynomialGCDMod,
        PolynomialPowerMod, FactorMod, FactorListMod]

End[]


EndPackage[]


(*:Tests:    (* Timings on DEC3100 *)

Test[PolynomialQuotientMod[1 + 3 x^2, 1 + x, x, 3], 
     0
]

Test[PolynomialQuotientMod[1 + 2 x + 4 x^2 + 2 x^5, 1 + x, x, 5]
     2 x + 2 x^2 + 3 x^3 + 2 x^4
]

Test[PolynomialQuotientMod[(1 + x)^3 (1+ x^3), x^2 + 2, x, 11] == 
     PolynomialMod[PolynomialQuotient[(1 + x)^3 (1+ x^3), x^2 + 2, x], 11],
     True
]

Test[PolynomialRemainderMod[1 + x + x^2 + x^3, 2 + x, x, 13],
     8
]

Test[PolynomialRemainderMod[(1 + x)^3 (1+ x^3), x^2 + 2, x, 11] == 
     PolynomialMod[PolynomialRemainder[(1 + x)^3 (1+ x^3), x^2 + 2, x], 11],
     True
]

Test[PolynomialGCDMod[(1 + 14 x) (1 + x + x^3), 1 + x, 13],
     1 + x
]

Test[PolynomialGCDMod[x^3 (1 + x) (1 + x^2)^3, x^2 (1 + x)^2, x^5 (1 + x), 13],
     x^2 + x^3
]

Test[PolynomialGCDMod[x^2 -1, x - 1, Prime[10^8]],
     2038074742 + x
]

Test[PolynomialPowerMod[1 + x, 5, {x^2 + 1, 7}], 
     3 + 3 x
]

Test[PolynomialPowerMod[1 + x, 5, {x^2 + 1, 7}] == 
     PolynomialRemainderMod[(1 + x)^5, x^2 + 1, x, 7], 
     True
]

Test[PolynomialPowerMod[1 + x, 10^4, {x^3 + x^2 + 1, Prime[10^9]}],
     8076973906 + 12198694781 x + 13390038726 x^2  
]                              (* 1 second *)


Test[FactorMod[1 + 2 x^3, 3], 
     2 (2 + x)^3
]

Test[FactorMod[Expand[(x + 1)^9 (x^4 + x + 1)], 3],
     (1 + x)^9 (2 + x) (2 + x + x^2 + x^3)
]

Test[FactorMod[Expand[(x + 1)^9 (x^4 + x + 1)], 3] == 
     Factor[Expand[(x + 1)^9 (x^4 + x + 1)], Modulus -> 3],
     True
]

Test[FactorMod[Expand[(x+1)^2 (x^2 + 5)^2 (x^2 + 13)], 10^64 + 57],
     (1 + x)^2 (5 + x^2)^2 (13 + x^3)
]                                      (* 100 seconds *)

Test[FactorListMod[1 + 2 x^3, 3],
     {{2, 1}, {2 + x, 3}}
]

Test[FactorListMod[Expand[(x + 1)^9 (x^4 + x + 1)], 3] == 
     FactorList[Expand[(x + 1)^9 (x^4 + x + 1)], Modulus -> 3],
     True

*)




