with(LinearAlgebra):
interface(rtablesize=50);#Maple returns degree(0,x) as -inf,
#we need it to return 0 for indexing
degree2 := proc(f,x)
local d;
d := degree(f,x);
if d = -infinity then return 0; fi;
return d;
end proc:expandSpace := proc(F,Space,dCur,dNew,k)
local i,A;
if(dNew+1 <= Space) then
return F,Space;
fi;
#create new array
A := Array(1..2^k);
#move elements over
for i from 1 to dCur+1 do
A[i] := F[i];
od;
return A,2^k;
end proc:#Calculates M,S polynomials for Diophantine Equations
#Verifies the initial n input polynomials are relatively prime
MultiEEA := proc(U::list,x::name,p::prime)
local n,M,sis,i,g,s,t,Mpolys,Spolys;
#local variables
n := nops(U);
Mpolys := Array(2..n);
Spolys := Array(1..n-1);
M := 1;
#Generate M
for i from n by -1 to 2 do
M := Expand(M*U[i]) mod p;
Mpolys[i] := M;
od;
Mpolys := convert(Mpolys,list);
sis := NULL;
#Generate S, verify polynomials relatively prime
for i from 1 to n-1 do
g := Gcdex(Mpolys[i],U[i],x,'Spolys[i]') mod p;
if g=1 then sis := sis,s; else return FAIL,FAIL,FAIL; fi;
od;
[sis],Mpolys,convert(Spolys,list);
end:#Solves the polynomial Diophantine Equation
DiophantineN := proc(U,c,M,S,p,x)
local n,q,g,ck,i,s,t,Sigmas;
n := nops(U);
ck := c;
Sigmas := Array(1..n);
for i from 1 to n-1 do
Sigmas[i] := Rem(ck*S[i],U[i],x) mod p;
ck := Quo(ck-Sigmas[i]*M[i],U[i],x) mod p;
end do;
Sigmas[n] := ck;
return(convert(Sigmas,list));
end proc:#calculates coeff(f_1*f_2*...*f_n,(y-alpha)^k) and stores the results in G
coeffExtract := proc(p,alpha,k,H,Space,Degrees,n,dy)
local i,j,m,s,t,L,D,MIN,MAX,Delta;
t := n; s := 1;
while(t > 1) do
for i from 0 to floor(t/2)-1 do
H[s+t+i],Space[s+t+i] := expandSpace(H[s+t+i],Space[s+t+i],Degrees[s+t+i],Degrees[s+2*i]+Degrees[s+2*i+1],ceil(log[2](Degrees[s+2*i]+Degrees[s+2*i+1]+1)));
Degrees[s+t+i] := Degrees[s+2*i] + Degrees[s+2*i+1];
if k <= Degrees[s+2*i] + Degrees[s+2*i+1] then
MIN := max(0,k-Degrees[s+2*i+1]);
MAX := min(k,Degrees[s+2*i]);
for L from MIN to MAX do
H[s+t+i][k+1] := H[s+t+i][k+1] + H[s+2*i][L+1]*H[s+2*i+1][k-L+1] mod p;
od;
fi;
#D[i] := D[2*i] + D[2*i+1];
od;
if t mod 2 = 1 then
H[s+t+i],Space[s+t+i] := expandSpace(H[s+t+i],Space[s+t+i],Degrees[s+t+i],Degrees[s+2*i],ceil(log[2](Degrees[s+2*i]+1)));
if Degrees[s+2*i] >= k then
Degrees[s+t+i] := Degrees[s+2*i];
H[s+t+i][k+1] := H[s+2*i][k+1] mod p;
fi;
fi;
s := s+t; t := ceil(t/2);
od;
if k+1 > Space[s] then
Delta := 0;
else
Delta := H[s][k+1];
fi;
return Delta,H,Space,Degrees;
end proc:#Fills in missing values in G
coeffUpdate := proc(p,alpha,k,H,Space,D,n)
local i,s,t,L,T;
s := 1; t := n;
while(t > 1) do
for i from 0 to floor(t/2)-1 do
if s=1 then
T[i] := 0;
if D[2*i+1] = k then T[i] := H[s+2*i][k+1]*H[s+2*i+1][0+1] mod p; fi;
if D[2*i+2] = k then T[i] := T[i] + H[s+2*i][0+1]*H[s+2*i+1][k+1] mod p; fi;
else
T[i] := H[s+2*i+1][0+1]*T[2*i] + H[s+2*i][0+1]*T[2*i+1] mod p;
fi;
H[s+t+i],Space[t+s+i] := expandSpace(H[s+t+i],Space[t+s+i],D[s+t+i],D[s+2*i] + D[s+2*i+1],ceil(log[2](D[s+2*i] + D[s+2*i+1]+1)));
D[s+t+i] := D[s+2*i] + D[s+2*i+1];
H[s+t+i][k+1] := H[s+t+i][k+1] + T[i] mod p;
od;
if(t mod 2 = 1) then
T[i] := H[s+2*i][k+1];
H[s+t+i],Space[t+s+i] := expandSpace(H[s+t+i],Space[t+s+i],D[s+t+i],D[s+2*i],ceil(log[2](D[s+2*i]+1)));
D[s+t+i] := D[s+2*i];
H[t+s+i][k+1] := T[i];
fi;
s := s+t; t := ceil(t/2);
od;
return H,Space,D;
end proc:CubicBivariateHensel := proc(A::polynom,F0::list,x::name,y::name,alpha::integer,p::prime)
#A(x,y) - polynomial to factor
#F0 - list of monic, relatively prime polynomials in x s.t. the product is equal to...
#x - variable 1 (usually x)
#y - variable 2 (usually y)
#alpha - integer to use for calculating the taylor coeff.
#p - prime
#local variables
local n,m,B,F,W,G,i,j,k,t,ck,sigmas,deltas,Delta,Coeffs,evalPoints,dx,dy,T,M,S,Space,Degrees,numPolys;
n := nops(F0);
F := F0;
dx := degree(A,x);
dy := degree(A,y);
#find the number of polynomials we will calculate
t := n;
numPolys := 1;
while t > 1 do
numPolys := numPolys + t;
t := ceil(t/2);
od;
#initial arrays/lists
deltas := [seq(0,i=1..dx)];
evalPoints := [seq(i,i=0..dx-1)];
W := Array(1..dx);
G := Array(1..dx);
Degrees := Array(1..dx);
Space := Array(1..dx);
for i from 1 to dx do
G[i] := Array(1..numPolys);
for j from 1 to numPolys do
G[i][j] := Array(1..1);
od;
Space[i] := Array([seq(1,i=1..numPolys)]);
Degrees[i] := Array(1..numPolys);
od;
#get a taylor series around (y-alpha) (does NOT use Shaw and Traub's method)
B := taylor(A,y=alpha,dy+1);
#Solve for M polynomials to use for the Diophantine Equation (Optimization)
T,M,S := MultiEEA(F0,x,p); # Solve this once for re-use
#do initial evaluations
for i from 1 to dx do
for j from 1 to n do
G[i][j][1] := Eval(F0[j],x=evalPoints[i]) mod p;
od;
deltas[i],G[i],Space[i],Degrees[i] := coeffExtract(p,alpha,0,G[i],Space[i],Degrees[i],n);
od;
#print(G[1]);
#if the EEA failed
if (T,M) = (FAIL,FAIL) then
return "Initial Factors are not relatively prime";
fi;
#main loop
for k from 1 to dy do
#printf("k=%d\134n",k);
#CoefficientExtraction
for i from 1 to dx do
deltas[i],G[i],Space[i],Degrees[i] := coeffExtract(p,alpha,k,G[i],Space[i],Degrees[i],n);
od;
#Interpolation
Delta := Interp(evalPoints,deltas,x) mod p;
#print("got here");
ck := Expand(coeff(B,(y-alpha),k) - Delta) mod p;
if add(degree(F[i],y),i=1..n) = dy and ck <> 0 then
return (FAIL);
fi;
if ck <> 0 then
#Solve Diophantine Equation for coefficients
sigmas := DiophantineN(F0,ck,M,S,p,x);
#Update the values of F
for i from 1 to n do
F[i] := F[i] + sigmas[i]*(y-alpha)^k;
od;
#Update Evaluations
for i from 1 to dx do
for j from 1 to n do
t := Eval(sigmas[j],x=evalPoints[i]) mod p;
if t <> 0 then
G[i][j],Space[i][j] := expandSpace(G[i][j],Space[i][j],Degrees[i][j],k,ceil(log[2](k+1)));
Degrees[i][j] := k;
G[i][j][k+1] := t;
fi;
od;
od;
#Perform CoefficientUpdate
for i from 1 to dx do
G[i],Space[i],Degrees[i] := coeffUpdate(p,alpha,k,G[i],Space[i],Degrees[i],n);
od;
fi;
od;
#return bivar polynomials or fail
if add(degree(F[i],y),i=1..n) = dy then
return(F);
else
return(FAIL);
fi;
end proc:#`mod` := mods;p := 1009;alpha := 3;f1 := x^4 + randpoly([x,y],dense,degree=3);
f2 := x^4 + randpoly([x,y],dense,degree=3);
f3 := x^4 + randpoly([x,y],dense,degree=3);
f4 := x^4 + randpoly([x,y],dense,degree=3);A := Expand(f1*f2*f3*f4) mod p;f10 := Eval(f1,y=alpha) mod p;
f20 := Eval(f2,y=alpha) mod p;
f30 := Eval(f3,y=alpha) mod p;
f40 := Eval(f4,y=alpha) mod p;F0 := [f10,f20,f30,f40];Rem(Expand(A - mul(F0)) mod p,(y-3),y) mod p;C := CubicBivariateHensel(A,F0,x,y,alpha,p):
Expand(C) mod p;Expand(A-mul(C)) mod p;