Clear["Factor`*"]
BeginPackage["Factor`", "Numbers`", "Polynom`", "Euclid`"]
(* Copyright: Erich Kaltofen, 4/27/1990.
   Permission to use provided the copyright notice is not removed. *)

Factor`FTrace::usage =
"Type \"Factor`FTrace[s__]:=Print[s]\" to trace this package."

FactorPolyMod::usage =
"FactorPolyMod[f_, x_, p_] returns a list of {g, e} such that g
 is an irreducible factor of f and e is its multiplicity;
 Mathematica's Factor[f, Modulus->p]."

SquarefreeMod::usage =
"SquarefreeMod[f_, x_, p_] returns a list of {g, e} such that
 g is the factor of multiplicity e in f modulo p."
 
CantorZassenhaus::usage =
"CantorZassenhaus[f_, x_, d_, p_] returns a list of factors
 of the integer polynomial f in x modulo p under the assumption
 that each factor has degree d, p an odd prime."

DistinctDegree::usage =
"DistinctDegree[f_, x_, p_] returns a list of pairs
 {g, d} where g is the product of all irreducible factors
 of the integer polynomial f in x modulo p of degree d."

PolynomialPowerMod::usage =
"PolynomialPowerMod[g_, x_, e_, f_, p_] returns
 g^e (modulo f, p)."

PolynomialGCDMod::usage =
"PolynomialGCDMod[f_, g_, x_, p_] has the same functionality
 as Mathematica's PolynomialGCD[f, g, Modulus->p] but does not
 restrict the size of p."

FactorSqfPrimPoly::usage =
"FactorSqfPrimPoly[f_, x_] computes all irreducible factors of the
 squarefree primitive integer polynomial f in x."

MultiComplement::usage =
"MultiComplement[list1_, list2_, equal_:Equal] computes the multi-set
 complement of list1 minus list2."

Begin["`private`"]

FTrace[s__]:=Print[s]

FactorPolyMod[f_, x_, p_]:=
Block[{F, (* monic modular version of f *)
       sqflist, (* list of squarefree factors *)
       sqflength, (* length of sqflist *)
       sqfpair, (* individual entry in sqflist *)
       i, j, (* loop indices *)
       g, (* squarefree factor *)
       e, (* multiplicity *)
       ddlist, (* distinct degree factor list *)
       ddlength, (* lenght of ddlist *)
       ddpair, (* individual entry in ddlist *)
       h, (* distinct degree factor *)
       d, (* distinct degree *)
       czlist, (* irreducible factors list *)
       czlength, (* lenght of czlist *)
       factorlist = {}
      },
      (* Make f monic *)
      F = PolynomialMod[f, p]; If[Exponent[F, x] < 1, Return[ {{F, 1}} ] ];
      F = Polynom`MakeMonicMod[F, x, p];

      (* Squarefree decomposition *)
      sqflist = SquarefreeMod[F, x, p];
      FTrace[StringForm["sqflist=``", sqflist] ];

      (* Factor each squarefree factor *)
      sqflength = Length[sqflist];
      For[i=1, i<=sqflength, i++,
          sqfpair = First[sqflist]; sqflist = Rest[sqflist];
          g = First[sqfpair]; e = sqfpair[[2]];
          ddlist = DistinctDegree[g, x, p];
          FTrace[StringForm["ddlist=``", ddlist] ];

          (* Factor each distinct degree product *)
          ddlength = Length[ddlist];
          For[j=1, j<=ddlength, j++,
              ddpair = First[ddlist]; ddlist = Rest[ddlist];
              h = First[ddpair]; d = ddpair[[2]];
              czlist = CantorZassenhaus[ h, x, d, p];  
              FTrace[StringForm["czlist=``", czlist] ];

              (* Append factors and corresponding multiplicities to output *)
              czlength = Length[czlist];
              For[k=1, k<=czlength, k++,
                  AppendTo[factorlist, List[First[czlist], e] ];
                  czlist = Rest[czlist]
                 ]
             ] (* end process distinct degree factors *)
         ]; (* end process squarefree factors *)
     Return[factorlist]
] (* end FactorPolyMod *)

SquarefreeMod[f_, x_, p_]:=
(* This function is not as efficient as can be *)
Block[{fp, (* derivative of f in x *)
       g, (* GCD[f, fp] modulo p *)
       glist, (* squarefree decomposition of g *)
       h, (* product of factors of higher multiplicities in f *)
       foverh, (* f/h *)
       i, (* loop index *) j, (* Product index *)
       sqflist
      },
      fp = PolynomialMod[ D[f, x], p ];
      If[fp===0,
         If[ Exponent[f, x] < 1, Return[ {{f, 1}} ] ];
         (* select coefficients of x^(p*i) *)
         g = 0;
         For[i=0, i<=Exponent[f, x], i++,
             If[Mod[i, p]==0, g += Coefficient[f, x, i] * x^(i/p)]
            ];
         FTrace[StringForm["Higher multiplicities, g=``", g] ];
         glist = SquarefreeMod[g, x, p];
         (* multiply p into each multiplicity of the squarefree factors of g *)
         sqflist = Map[ Function[ {#[[1]], #[[2]] * p} ],
                        glist];
         Return[sqflist]
        ];
      g = PolynomialGCDMod[f, fp, x, p];
      If[Exponent[g, x]>0,
         FTrace[StringForm["Higher multiplicities, g=``", g] ];
         glist = SquarefreeMod[g, x, p];
         (* add 1 to each multiplicity of the squarefree factors of g *)
         sqflist = Map[ Function[ {#[[1]], #[[2]]+1} ],
                        glist];
         h = Product[sqflist[[j,1]]^sqflist[[j,2]], {j, 1, Length[sqflist]} ];
         foverh = PolynomialMod[ PolynomialQuotient[f, PolynomialMod[Expand[h], p], x], p];
         If[Exponent[foverh, x] > 0, PrependTo[sqflist, List[foverh, 1] ] ];
         Return[sqflist]
        ];
      Return[ List[List[f, 1] ] ]      
] (* end SquarefreeMod *)

DistinctDegree[f_, x_, p_]:=
(* Compute the distinct degree factorization of f modulo p *)
Block[{d, (* degree tried *)
       h, (* polynomial remaining *)
       list, (* accumulated answer *)
       xpower, (* x^(p^d) (modulo (f, p)) *)
       g (* GCD[ x^(p^d) - x, h ] *)
      },
      h = f; list = {}; 
      For[d=1, Exponent[h, x] > 0, d++,
          xpower = PolynomialPowerMod[x, x, p^d, h, p];
          g = PolynomialGCDMod[xpower-x, h, x, p];
          If[Exponent[g, x] > 0,
             FTrace[StringForm["Factor-product of degree `` is ``", d, g] ];
             AppendTo[list, List[g, d] ];
             h = PolynomialMod[ PolynomialQuotient[h, g, x], p ];
            ]
         ];
      Return[list]
] (* end DistinctDegree *)

CantorZassenhaus[f_, x_, d_, p_]:=
Block[{deg, (* degree of f in x *)
       r, (* number of irreducible factors of f *)
       t, (* random polynomial *)
       tpower, (* t^(p^d-1) (modulo f) *)
       h, (* possible split *)
       foverh (* f/h *)
      },
      deg = Exponent[f, x];
      If[deg == d, Return[{f}] ];
      r = deg/d; FTrace[StringForm["r=``", r] ];
      h = 1; (* to enter loop once *)
      While[Exponent[h, x]==0 || Exponent[h, x]==deg,
            t = Polynom`RandomPoly[2*d-1, (p-1)/2, x];
            FTrace[StringForm["Random choice=``", t] ];
            tpower = PolynomialPowerMod[t, x, (p^d-1)/2, f, p];
            h = PolynomialGCDMod[tpower-1, f, x, p];
            FTrace[StringForm["GCD=``", h] ];
           ];
      foverh = PolynomialMod[ PolynomialQuotient[f, h, x], p];
      FTrace[StringForm["f/GCD=``", foverh] ];
      Return[Join[CantorZassenhaus[h, x, d, p],
                  CantorZassenhaus[foverh, x, d, p ]
            ]    ]
] (* end CantorZassenhaus *)

FactorSqfPrimPoly[f_, x_]:=
Block[{prime, (* modulus used *)
       modfactorlist, (* factors modulo prime *)
       combinlist, (* factors for combinations without multiplicities *)
       g, (* individual modular factor *) e, (* its multiplicity *)
       i, k, (* loop indices *)
       TryCombinations, (* local recursive procedure *)
       LC, (* leading coefficient of f and remf *)
       fLC, (* f * LC and remf * LC *)
       TC, (* trailing coefficient of f*LC and remf * LC *)
       trial, (* trial combination result *)
       combinations, (* list of modular factors leading to integer factor *)
       factorlist, (* accumulated list of integer factors *)
       remf (* remaining co-factor of f *)
      }, 

   prime = Numbers`FindPrime[ 2 * Polynom`FactorCoefficientBound[f, x] ];
   FTrace[ StringForm["Modulus chosen = ``", prime] ];

   (* Factor polynomial modulo prime *)
   modfactorlist = FactorPolyMod[f, x, prime];
   FTrace[ StringForm["Modular factors = ``", modfactorlist] ];

   (* List multiple factors in duplicates *)
   combinlist = {};
   While[Length[modfactorlist] > 0,
         g = modfactorlist[[1,1]]; e = modfactorlist[[1,2]];
         modfactorlist = Rest[modfactorlist]; 
         For[i=1, i<=e, i++,
             AppendTo[combinlist, g]
            ]
        ];
    FTrace[ StringForm["cominlist = ``", combinlist] ];

    (* Try all combinations (see Polynom`FactorCoefficientBound) *)

    TryCombinations[g_, list_, k_, h_, hlist_]:=
    (* g......remaining unfactored polynomial (with leading coeff. squared),
       list...remaining modular factor list,
       k......cardinality of combinations to be tried,
       h......pre-factor for combinations,
       hlist..list of modular factors in h (for debugging only)
          Returns a pair {fac, sublist} such that fac is an integer factor of g,
       and sublist is the list of k modular factors needed to build fac;
       if no such list exists, {} is returned. *)
       Block[{f1, (* first element on list *)
              rest, (* remaining elements on list *)
              trial1,
              trial2,
              testfac
             },
             If[Length[list] < k,
                Return[{}]
               ];
             If[k==0,
                (* h holds accumulated modular factor product;
                   perform a trial division *)
                testfac = Polynom`ModBalanced[LC * h, x, prime];
                testTC = Coefficient[testfac, x, 0];
                If[testTC == 0, testTC = 1]; (* to prevent 0-division *)
                FTrace[StringForm["Test-division by ``", testfac] ];
                FTrace[StringForm["Combined factors list ``", hlist] ];
                If[Mod[TC, Abs[testTC] ] == 0,
                   If[PolynomialRemainder[g, testfac, x] === 0,
                      Return[ List[testfac, {}] ],
                      FTrace[StringForm["Test polynomial division failed"] ];
                      Return[ {} ]
                     ],
                   FTrace[StringForm["Test trailing coefficient division failed"] ];
                   Return[ {} ]
                  ]
               ];
             f1 = list[[1]]; rest = Rest[list];

             (* try combinations of cardinality k without f1 *)
             trial1 = TryCombinations[g, rest, k, h, hlist];
             If[trial1 == {}, (* no success *)
                (* try combinations of cardinality k-1 with f1 *)
                trial2 = TryCombinations[g, rest, k-1, PolynomialMod[h * f1, prime], Append[hlist, f1] ];
                If[trial2 != {},
                   Return[{trial2[[1]], Prepend[ trial2[[2]], f1 ]} ]
                  ];
                Return[trial2], (* {} *)
                (* else *) Return[trial1] (* {} *)
               ]
            ]; (* end local procedure TryCombinations *)
    
    (* prepare initial polynomial for trial division *)
    LC = Part[Polynom`MonicLC[f, x], 2];
    fLC = Expand[f * LC];
    TC = Coefficient[fLC, x, 0];

    remf = f; factorlist = {};
    For[k=1, k<=Length[combinlist]-1, k++,
        FTrace[ StringForm["Trying combinations of `` factors", k] ];
        trial = TryCombinations[fLC, combinlist, k, 1, {} ];
        If[trial != {},
           (* factor found with k combinations *)
           factor = First[trial];
           factor = Expand[2 * factor]; (* to have a content;
                                           otherwise, Part selects wrong thing *)
           factor = Part[FactorTerms[factor], 2]; (* take content *)
           AppendTo[factorlist, factor];

           (* remove combinations leading to factor from combinlist *)
           combinations = trial[[2]];
           FTrace[StringForm["Factor `` found, combinations = ``",
                            factor, combinations ] ];
           combinlist = MultiComplement[ combinlist, combinations ];
           FTrace[StringForm["Remaining combinations = ``", combinlist] ];

           (* prepare co-factor for trial division *)
           remf = PolynomialQuotient[remf, factor, x];
           LC = Part[Polynom`MonicLC[remf, x], 2];
           fLC = Expand[remf * LC];
           TC = Coefficient[fLC, x, 0];

           (* dirty trick: let loop try k combinations
              again for remaining combinlist *)
           k -= 1
          ]
       ];
          
   AppendTo[factorlist, remf];
   Return[factorlist]
] (* end FactorSqfPrimPoly *)

PolynomialPowerMod[g_, x_, e_, f_, p_]:=
Block[{localTimes},
      localTimes[a_, b_]:=PolynomialMod[ PolynomialRemainder[a*b, f, x], p ];
      Return[ Numbers`GenericExponentiation[g, e, localTimes] ]
] (* end PolynomialPowerMod *)

PolynomialGCDMod[f_, g_, x_, p_]:=
Block[{localremainder, monicg, h1, h2},
   localremainder[a_, b_]:=MakeMonicMod[ PolynomialRemainder[a, b, x], x, p];
   monicg = MakeMonicMod[g, x, p]; (* to guarantee first remainder to be
                                      fraction free *)
   Return[ Euclid`PlainGeneric[f, monicg, localremainder] ];
] (* end PolynomialGCDMod *)                             

MultiComplement[list1_, list2_, equal_:Equal]:=
Block[{newlist1, newlist2},
   newlist1 = list1; newlist2  = list2;
   While[newlist2 != {},
         newlist1 = RemoveSingleOccurrence[newlist1, First[newlist2], equal ];
         newlist2 = Rest[newlist2]
        ];
   Return[newlist1] 
]

RemoveSingleOccurrence[list_, elem_, equal_:Equal]:=
Block[{first, rest}, 
  If[list == {}, Return[{}] ];
  first = First[list]; rest = Rest[list];
  If[equal[first, elem], Return[rest] ];
  Return[ Prepend[ RemoveSingleOccurrence[rest, elem, equal],
                   first
        ]        ]
]

End[] (* `private` *)
EndPackage[] (* Factor` *)
Print["Package Factor` loaded"]