// Factorization of polynomials over finite fields function modexp(a,n,f) // compute a^n mod f, n >= 1 R := quo; return Parent(f)!((R!a)^n); // if n eq 1 then // return a mod f; // elif IsOdd(n) then // return (a*(modexp(a, n div 2, f)^2 mod f)) mod f; // else // return modexp(a, n div 2, f)^2 mod f; // end if; end function; // Distinct degree factorization function ddf(f) // f in F_q[X] monic and squarefree. // returns [, ..., ] with u_j product of >= 1 irreducible // polynomials of degree d_j, f = u_1*...*u_k. P := Parent(f); // F_q[X] F := CoefficientRing(P); // F_q X := P.1; q := #F; u := []; // list for the result d := 0; // d = degree of irreducible factors currently considered a := X mod f; // a = X^(q^d) mod f while Degree(f) gt 0 do d +:= 1; if Degree(f) lt 2*d then // current f must be irreducible return Append(u, ); end if; a := modexp(a, q, f); // compute a = X^(q^d) mod f h := GCD(f, a-X); if Degree(h) gt 0 then // found a contribution from irreducibles of degree d Append(~u, ); f := ExactQuotient(f, h); // remove the factors found end if; end while; return u; end function; // Produce a random monic polynomial of given degree in the polynomial ring "ring". function randpoly(ring, degree) F := CoefficientRing(ring); return ring!([Random(F) : i in [0..degree-1]] cat [F!1]); end function; // Find roots of f in F_q, q odd. function roots(f : level := -1, reduce := true) // f in F_q[X], q odd. // returns {r_1, ..., r_k}, where r_j are the roots of F in F_q. // set reduce := false if it is known that all roots are in F_q and simple. P := Parent(f); // F_q[X] F := CoefficientRing(P); // F_q q := #F; X := P.1; if reduce then // replace f by product of linear factors of f, each taken once f := GCD(f, modexp(X, q, f) - X); end if; assert LeadingCoefficient(f) eq 1; if level ge 0 then printf " "^level*"degree %o\n", Degree(f); end if; if Degree(f) eq 0 then // constant polynomial return {F|}; elif Degree(f) eq 1 then // f = X - r: return r return {-Coefficient(f,0)}; end if; a := Random(F); h := GCD(f, modexp(X-a, ExactQuotient(q-1, 2), f) - 1); // GCD(f, (X-a)^((q-1)/2) - 1) newlevel := level ge 0 select level+1 else level; return roots(h : level := newlevel, reduce := false) join roots(ExactQuotient(f, h) : level := newlevel, reduce := false); end function; p := NextPrime(314159); F := FiniteField(p); P := PolynomialRing(F); // Produce a random monic polynomial of degree n with only simple roots, // all in F_q. function randpoly1(ring, n) F := CoefficientRing(ring); X := ring.1; zeros := {F|}; while #zeros lt n do Include(~zeros, Random(F)); end while; return &*[ring | X-a : a in zeros]; end function; // Test roots function procedure testroots(P, n) t0 := Cputime(); f := randpoly1(P, n); t1 := Cputime(t0); t0 := Cputime(); roots1 := roots(f); t2 := Cputime(t0); t0 := Cputime(); roots2 := Roots(f); t3 := Cputime(t0); assert roots1 eq {r[1] : r in roots2}; printf "f in F_%o[X], degree = %o:\n", #CoefficientRing(P), n; printf " generate f: %o sec\n", t1; printf " roots: %o sec\n", t2; printf " Roots: %o sec\n\n", t3; end procedure; // > testroots(P, 5000); // f in F_314161[X], degree = 5000: // generate f: 0.150 sec // roots: 1.520 sec // Roots: 2.240 sec // Equal degree factorization function edf(f, d) // f in F_q[X], q odd, f product of distinct irreducible polynomials of degree d assert IsDivisibleBy(Degree(f), d); P := Parent(f); // F_q[X] if Degree(f) eq 0 then return [P|]; elif Degree(f) eq d then return [f]; end if; F := CoefficientRing(P); // F_q repeat a := P![Random(F) : i in [0..Degree(f)-1]]; h := GCD(f, modexp(a, ExactQuotient(#F^d-1,2), f)-1); until Degree(h) gt 0 and Degree(h) lt Degree(f); return edf(h, d) cat edf(ExactQuotient(f, h), d); end function; // Complete factorization for squarefree f function sqf(f) // f in F_q[X], q odd dd := ddf(f); // distinct degree factorization return &cat[edf(pair[2], pair[1]) : pair in dd]; end function; // Complete factorization function fact(f) // f in F_q[X] monic, q odd. // returns a list [, ...] with h_j irreducible, e_j >= 1, such that // f = h_1^(e_1)*... assert LeadingCoefficient(f) eq 1; P := Parent(f); // F_q[X] F := CoefficientRing(P); // F_q X := P.1; q := #F; u := []; // list for the result d := 0; // d = degree of irreducible factors currently considered a := X mod f; // a = X^(q^d) mod f while Degree(f) gt 0 do // do distinct degree factorization d +:= 1; // printf "[%o]", d; if Degree(f) lt 2*d then // current f must be irreducible // printf "[%o]*\n", Degree(f); return Append(u, ); end if; a := modexp(a, q, f); // compute a = X^(q^d) mod f g := GCD(f, a-X); if Degree(g) gt 0 then // found a contribution from irreducibles of degree d // now do equal degree factorization to find all the irreducibles of degree d ed := edf(g, d); // printf "*"^#ed; f := ExactQuotient(f, g); // divide out each factor once // for each irreducible factor, find its multiplicity in f for h in ed do e := 1; // exponent is at least 1 qu, r := Quotrem(f, h); while r eq 0 do e +:= 1; f := qu; qu, r := Quotrem(f, h); end while; Append(~u, ); // append information to u end for; // here all the irreducible factors of degree <= d are removed from f end if; end while; return u; end function; // Test fact function procedure testfact(P, n) t0 := Cputime(); f := randpoly(P, n); t1 := Cputime(t0); t0 := Cputime(); fact1 := fact(f); t2 := Cputime(t0); t0 := Cputime(); fact2 := Factorization(f); t3 := Cputime(t0); assert {a : a in fact1} eq {a : a in fact2}; printf "f in F_%o[X], degree = %o:\n", #CoefficientRing(P), n; printf " generate f: %o sec\n", t1; printf " fact: %o sec\n", t2; printf " Factorization: %o sec\n", t3; printf "Type: %o\n\n", [ : a in fact2]; end procedure;