(*:Mathematica Version: 2.0 *)

(*:Package Version: 1.1 *)

(*:Name: Algebra`PolynomialPowerMod` *)

(* :Title:  PolynomialPowerMod *)

(* :Author: Ilan Vardi *)

(*:Copyright: Copyright 1988-1999, Wolfram Research, Inc. *)

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

(* :History: Version 1.0 by Ilan Vardi, 1988.
	     Functionality shifted to kernel functions Factor, FactorList,
		PolynomialGCD, PolynomialQuotient, and PolynomialRemainder,
		ECM, 1996.
*)

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

(*:Summary:
This package provides PolynomialPowerMod (similar to the kernel
function PowerMod) for efficiently computing powers of a polynomial modulo
a prime and another polynomial.  It also extends Factor, FactorList, and
PolynomialGCD to handle prime modulii larger than 2^15.  It adds support
for the Modulus option to PolynomialQuotient and PolynomialRemainder.
*)



(* :Limitation:
	The Modulus extension to Factor can factor degree 20 polynomials
		modulo hundred digit primes in reasonable time.
        These algorithms are for polynomials in one variable only.
*)  


(* :Discussion:

This package extends the built-in polynomial functions
PolynomialGCD, FactorList, and Factor to primes larger 
than 46337. It also adds support for the Modulus option to PolynomialQuotient 
and PolynomialRemainder.  Most importantly, it
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.

*)

(* Note: This top level code is much slower than the built-in 
             functions for primes < 2^16. For p = 2, this code is
             especially inefficient.  Thus the top level code for
	     Factor, FactorList, PolynomialGCD, PolynomialQuotient,
	     and PolynomialRemainder is only called for large prime
	     modulii (above the max modulus 46337).  Note also
	     that the internal functions iPolynomialQuotientMod
	     iPolynomialRemainderMod, and iPolynomialGCDMod are
	     used to minimize the timing hit caused by attaching rules
	     to built-in symbols.
*)

BeginPackage["Algebra`PolynomialPowerMod`",
	"NumberTheory`NumberTheoryFunctions`"]  

PolynomialPowerMod::usage :=
"PolynomialPowerMod[poly1, n, {poly2, p}] gives the polynomial poly1^n reduced
modulo poly2 and p, where in n is an integer and p is a prime."

Begin["`Private`"]

$MaxModulus = 46337 (* Prime[4792] *) (* between 2^15 and 2^16 *)
                        
Unprotect[PolynomialPowerMod] (* new function *)
Unprotect[PolynomialQuotient, PolynomialRemainder, PolynomialGCD,
		Factor, FactorList]

Options[PolynomialQuotient] = {Modulus -> 0}
Options[PolynomialRemainder] = {Modulus -> 0}


(* ======================================================================== *)
PolynomialQuotient[f_, g_, x_, Modulus -> p_?PrimeQ]:= 
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]
      ] (* p any value *)

(* ================== (internal) iPolynomialQuotientMod ==================== *)
iPolynomialQuotientMod[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]
 ] 

(* ======================================================================== *)
PolynomialRemainder[f_, g_, x_, Modulus -> p_?PrimeQ]:= 
	(
    PolynomialMod[PolynomialRemainder[PolynomialMod[f, p], 
                  PolynomialMod[g, p], x], p]
	)	(* p any value *)

(* ==================== (internal) iPolynomialRemainderMod ================= *)
iPolynomialRemainderMod[f_, g_, x_, p_]:= 
    PolynomialMod[PolynomialRemainder[PolynomialMod[f, p], 
                  PolynomialMod[g, p], x], p]

(* ======================================================================== *)
PolynomialGCD[f_, Modulus -> p_?PrimeQ]:= PolynomialMod[f, p]

PolynomialGCD[f_, g_, Modulus -> p_?PrimeQ]:= 
  (
     polynomialgcdmod[PolynomialMod[f, p], PolynomialMod[g, p],
                     If[Variables[f] != {}, First[Variables[f]],
                        If[Variables[g] != {}, First[Variables[g]],
				Null]],
                     p]
  ) ; p > $MaxModulus


PolynomialGCD[f_, g__, Modulus -> p_?PrimeQ]:= 
  (
     iPolynomialGCDMod[f, iPolynomialGCDMod[g, p], p]
  ) /; p > $MaxModulus

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  

(* ==================== (internal) iPolynomialGCDMod ====================== *)
iPolynomialGCDMod[f_, p_]:= PolynomialMod[f, p]

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

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


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

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


PolynomialPowerMod[f_, 0, x_, {g_, p_Integer}]:= 1

PolynomialPowerMod::notaunit = "Polynomial `1` in the variable `2` cannot
be inverted modulo `3` and `4`."

PolynomialPowerMod[f_, n_Integer, x_, {g_, p_Integer}] :=
    Block[{finv, gb, soln, gcd, vars},
        gcd = iPolynomialGCDMod[f, g, p];
	vars = Variables[gcd];
	If [Select[vars, #==x&] != {}, 

	    Message[PolynomialPowerMod::notaunit, f, x, g, p];
	    Return[Null],

	    gb = GroebnerBasis[{finv*f-1, g}, {finv, x},
	      Modulus->p, CoefficientDomain->RationalFunctions];
	    soln = Solve[gb[[2]]==0, finv][[1]];
	    finv = finv /. soln;
	    PolynomialPowerMod[finv, -n, x, {g, p}]
	    ]
    ] /; n < 0

(* ======================================================================== *)
Factor[f_, Modulus -> p_?PrimeQ] :=
   (
	Times @@ (#[[1]]^#[[2]]& /@ FactorList[f, Modulus -> p])
   ) /; p > $MaxModulus

(* ======================================================================== *)
FactorList[f_, Modulus -> p_?PrimeQ]:= 
  Block[{flm = factorlistmod[f, p]},
      {#, Count[flm, #]}& /@ Union[flm]
  ] /; p > $MaxModulus

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 = iPolynomialGCDMod[fp, df, p];
       If[Exponent[factor, x] == 0, 
          factorsquarefree[fp, x, 1, x, p], 
          Join[factorlistmod[iPolynomialQuotientMod[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 = iPolynomialGCDMod[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 = iPolynomialQuotientMod[f, factor, x, p];
       Join[factorsquarefree[q, x, d + 1,
                             iPolynomialRemainderMod[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 = iPolynomialGCDMod[s, fp, p]; 
       If[0 < Exponent[t,x] < exp, Break[],
         t =  iPolynomialGCDMod[fp, 
                PolynomialPowerMod[s, (pd-1)/2, x, {fp, p}] - 1, p];
            If[ 0 < Exponent[t, x] < exp, Break[]]
    ]]; 
    q = iPolynomialQuotientMod[fp, t, x, p];
    t = If[Exponent[q, x] < Exponent[t, x], q, t];
    Join[factordegree[t, x, d, p], 
         factordegree[iPolynomialQuotientMod[f, t, x, p], x, d, p]]
   ]

Protect[PolynomialQuotient, PolynomialRemainder, PolynomialGCD,
	Factor, FactorList]

Protect[PolynomialPowerMod]

End[]


EndPackage[]


(*:Tests:

Test[PolynomialQuotient[1 + 3 x^2, 1 + x, x, Modulus -> 3], 
     0
] (* package function *)

Test[PolynomialQuotient[1 + 2 x + 4 x^2 + 2 x^5, 1 + x, x, Modulus -> 5],
     2 x + 2 x^2 + 3 x^3 + 2 x^4
] (* package function *)

Test[PolynomialQuotient[(1 + x)^3 (1+ x^3), x^2 + 2, x, Modulus -> 11],
	1 + 7*x + x^2 + 3*x^3 + x^4
] (* package function *)

Test[PolynomialRemainder[1 + x + x^2 + x^3, 2 + x, x, Modulus -> 13],
     8
] (* package function *)

Test[PolynomialRemainder[(1 + x)^3 (1+ x^3), x^2 + 2, x, Modulus -> 11],
     10
] (* kernel function *)

Test[PolynomialGCD[(1 + 14 x) (1 + x + x^3), 1 + x, Modulus -> 13],
     1 + x
] (* kernel function *)

Test[PolynomialGCD[x^3 (1 + x) (1 + x^2)^3, x^2 (1 + x)^2, x^5 (1 + x),
	Modulus -> 13],
     x^2 + x^3
] (* kernel function *)

Test[PolynomialGCD[x^2 -1, x - 1, Modulus -> Prime[10^8]],
     2038074742 + x
] (* package function *)

Test[PolynomialPowerMod[1 + x, 5, {x^2 + 1, 7}], 
     3 + 3 x
] (* package function *)

Test[PolynomialPowerMod[1 + x, 5, {x^2 + 1, 7}] == 
     PolynomialRemainder[(1 + x)^5, x^2 + 1, x, Modulus -> 7], 
     True
] (* package function *)

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

Test[Factor[1 + 2 x^3, Modulus -> 3], 
     2 (2 + x)^3
] (* kernel function *)

Test[Factor[Expand[(x + 1)^9 (x^4 + x + 1)], Modulus -> 3],
     (1 + x)^9 (2 + x) (2 + x + x^2 + x^3)
] (* kernel function *)

Test[Factor[Expand[(x+1)^2 (x^2 + 5)^2 (x^2 + 13)], Modulus -> 10^64 + 57],
     (1 + x)^2 (5 + x^2)^2 (13 + x^3)
] (* package function *)                  (* 100 seconds on DEC3100 *)
	(* 282 seconds on SPARC, but only 188 seconds if algorithm
		implemented w/out attaching rules to system symbols *)

Test[FactorList[1 + 2 x^3, Modulus -> 3],
     {{2, 1}, {2 + x, 3}}
] (* kernel function *)


*)




