/* This module implements factorizing integers and polynomials */ /* This is Pollard's Rho method of factorizing, as described in * "Modern Computer Algebra". It is a rather fast algorithm for * factoring, but doesn't scale to polynomials regrettably. * * It acts 'by chance'. This is the Floyd cycle detection trick, where * you move x(i+1) = f(x(i)) and y(i+1) = f(f(y(i))), so the y goes twice * as fast as x, and for a certain i x(i) will be equal to y(i). * * "Modern Computer Algebra" reasons that if f(x) = (x^2+1) mod n for * the value n to be factored, then chances are good that gcd(x-y,n) * is a factor of n. The function x^2+1 is arbitrary, a higher order * polynomial could have been chosen also. * */ PollardRhoRandom(n_IsInteger) <-- MathFloor(Random()*n); 10 # PollardRhoFactorize(n_IsPrimePower) <-- {GetPrimeFactors(n)}; 15 # PollardRhoFactorize(0) <-- {}; 20 # PollardRhoFactorize(_n) <-- [ Local(x,y,i,gcd); gcd:=n; While(gcd = n) [ /* Pick a random value between 0 and n-1 */ x:=PollardRhoRandom(n); /* Initialize loop */ gcd:=1; y:=x; i:=0; /* loop until failure or success found */ While(gcd = 1) [ x:= Mod( ((x^2)+1), n); y:= Mod( ((y^2)+1)^2+1, n); i++; If(x-y = 0, gcd := 1, gcd:=Gcd(x-y,n) ); ]; ]; /* Return result found */ PollardCombineLists(PollardRhoFactorize(gcd),PollardRhoFactorize(Div(n,gcd))); ]; /* PollardCombineLists combines two assoc lists used for factoring. the first element in each item list is the factor, and the second the exponent. Thus, an assoc list of {{2,3},{3,5}} means 2^3*3^5. */ 5 # PollardMerge(_list,{1,_n}) <-- True; 10 # PollardMerge(_list,_item)_(Assoc(item[1],list) = Empty) <-- DestructiveInsert(list,1,item); 20 # PollardMerge(_list,_item) <-- [ Local(assoc); assoc := Assoc(item[1],list); assoc[2]:=assoc[2]+item[2]; ]; PollardCombineLists(_left,_right) <-- [ ForEach(item,right) [ PollardMerge(left,item); ]; left; ]; /* TODO for Bk: - ModDivide, ModGcd, make an interface for in the scripts. - The cases non-monic polynomials and polys with rational coefficients are not handled very well yet. These two are related: if the non-monic goes well then the rational coeffs goes well. Problem is that Knuth describes an order n! algorithm, for the n possible options! This cannot be right... When you have many factors this might go wrong. The algorithm will terminate (a great improvement over the previous) but might not fully factor the polynomial. */ 10 # FindMod(_poly,_item,_p,_n)_(Degree(Gcd(poly,item)) > 0) <-- Expand(item); 11 # FindMod(_poly,_item,_p,_n) <-- FindModAux(poly,item,p,n); 11 # FindModAux(_poly,_item,_p,_n)_(Mod(poly,item+n*p) = 0) <-- Expand(item+n*p); 12 # FindModAux(_poly,_item,_p,_n)_(Mod(poly,item-n*p) = 0) <-- Expand(item-n*p); 13 # FindModAux(_poly,_item,_p,_n)_(Mod(poly,(LeadingCoef(poly))*item+n*p) = 0) <-- Expand((LeadingCoef(poly))*item+n*p); 14 # FindModAux(_poly,_item,_p,_n)_(Mod(poly,(LeadingCoef(poly))*item-n*p) = 0) <-- Expand((LeadingCoef(poly))*item-n*p); 20 # FindModAux(_poly,_item,_p,_n)_(n<10) <-- [ /* Echo({item,p,"n = ",n}); Echo({(Mod(LeadingCoef(poly),p))*item+n*p}); */ FindModAux(poly,item,p,n+1); ]; 30 # FindModAux(_poly,_item,_p,_n) <-- 1; MergeTerm(_f,_i) <-- MergeTerm(f,i,Assoc(i[1],f)); 10 # MergeTerm(_f,_i,Empty) <-- DestructiveAppend(f,i); 20 # MergeTerm(_f,_i,list_IsList) <-- [ list[2] := list[2] + i[2]; ]; 10 # MergeFactors(_f,_s) <-- [ ForEach(i,s) MergeTerm(f,{Monic(i[1]),i[2]}); f; ]; FactorUniVarSec(_poly) <-- [ Local(checked,i); checked := Bk(poly); Local(pre); pre:=Div(poly,FW(checked)); // Rescan the solutions, because there may be some that still // divide the remaining polynomial. For(i:=1,i<=Length(checked),i++) [ While (Mod(pre,checked[i][1]) = 0) [ checked[i][2] := checked[i][2] + 1; pre := Div(pre, checked[i][1]); ]; ]; if (pre != 1) [ Local(lc,monic); lc := LeadingCoef(pre); monic := Monic(pre); if (monic != 1) [checked := {monic,1}:checked;]; if (lc != 1) [checked := {lc,1}:checked;]; ]; checked; ]; MakeInteger(_poly) <-- [ Local(cf,t1,t2); cf := MapSingle("Denom",Abs(Coef(poly,(0 .. Degree(poly))))); cf := Select({{n},Not IsZero(n)},cf); While (Length(cf) > 1) [ t1:=Head(cf); cf:=Tail(cf); t2:=Head(cf); cf:=Tail(cf); cf := Lcm(t1,t2):cf; ]; poly := poly*(cf[1]); poly; ]; Bk(_poly) <-- [ Local(p,checked,factored,result); // result will contain the final result returned to the user. checked := {}; result:= {}; // Create a internal representation, to speed up processing poly := MakeUni(poly,VarList(poly)); // Convert from polynomial containing rational coefficients // to one containing integer coefficients poly := MakeInteger(poly); // Trivial factors if (poly[2] > 0) [ poly := UniVariate(poly[1],0,poly[3]); ]; // poly := poly*LeadingCoef(poly); /*Echo({poly," should have integer coefficients"});*/ // square-free factorization [ Local(dv,dd); dd := Apply("D",{x,poly}); dv:=0; if (Degree(dd) > 0) [ dv:=Gcd(poly,dd); ]; if (Degree(dv)>0) [ Local(new); dv := Monic(dv); poly := Div(poly,dv); poly := MakeInteger(poly); /*Echo({poly," should have integer coefficients"});*/ new := Bk(dv); ForEach(item,new) [ While(Mod(poly,item[1]) = 0) [ item[2] := item[2] + 1; poly := Div(poly,item[1]); ]; ]; poly := MakeInteger(poly); MergeFactors(checked, new); /*Echo({"squarefree: ",checked}); */ ]; ]; p:=Max(1+2*Abs(Coef(poly,(0 .. Degree(poly))))); if (p<5) [p := 5;]; if (p > 100) [ p := 100; ]; While (Not IsPrime(p)) p++; factored := Bk(poly,p); if (Length(factored) < 2) [ p++; While (Not IsPrime(p)) p++; factored := Bk(poly,p); ]; ForEach(item,factored) [ if (Degree(item[1]) = 1) MergeFactors(result,{{Monic(FindMod(poly,item[1],p,0)),item[2]}}); ]; MergeFactors(checked,result); checked; ]; Bk(_poly,_prime) <-- [ Local(coefs,result); if (Degree(poly)<2) [ result := {{poly,1}}; ] else [ coefs := MakeUni(poly)[3]; result := Berlekamp(coefs,prime); Local(i); For(i:=1,i<=Length(result),i++) [ result[i] := {Monic(NormalForm(UniVariate(x,0,result[i]))),1}; ]; ]; result; ]; /* New factorization : split between integers and polynomials. */ 10 # Factors(p_IsInteger) <-- PollardRhoFactorize(p); 20 # Factors(p_CanBeUni)_(Length(VarList(p)) = 1) <-- FactorUniVarSec(p); Factor(_p) <-- FW(Factors(p)); /* FW: pass FW the result of Factors, and it will show it in the * form of p0^n0*p1^n1*... */ 10 # FWatom({_a,1}) <-- a; 20 # FWatom({_a,_n}) <-- UnList({Atom("^"),a, n}); 5 # FW(_list)_(Length(list) = 0) <-- 1; 10 # FW(_list)_(Length(list) = 1) <-- FWatom(list[1]); 20 # FW(_list) <-- [ Local(result); result:=FWatom(Head(list)); ForEach(item,Tail(list)) [ result := UnList({ Atom("*"),result,FWatom(item)}); ]; result; ]; 10 # Roots(poly_CanBeUni) <-- [ Local(factors,result,uni,root,i); factors:=Factors(poly); result:={}; ForEach(item,factors) [ uni:=MakeUni(item[1]); If(Degree(uni) = 1, [ root:= -Coef(uni,0)/Coef(uni,1); For(i:=0,i