This Maple worksheet generates the data appearing in Tables 1, 2, 3, 4 and 5 in the paper "The complexity of sparse Hensel lifting and sparse polynomial factorization" by Michael Monagan and Baris Tuncer and examples for estimates in the paper.Michael Monagan and Baris Tuncer, March 2018.restart;###Create a random polynomial for examples:
n:=7; ## the number of variables
deg:=15; ## the total degree
Tf:=17161; ## the number of terms
X:=[seq(x[i],i=1..n)]: ## variables
p:=2^31-1; ### prime number
F:=randpoly(X,degree=deg,terms=Tf,coeffs = rand(1..1000)): ## generate a random polynomial
Roll:=rand(1..p): ## random number generator
sd:=binomial(n+deg,deg): ## possible number of monomials
tF:=evalf(Tf/sd); ## the density ratio of F ###to compute the number of terms of a polynomial
numterms := proc(a) if type(a,`+`) then nops(a) else 1 fi; end;### Eqn (3) in the paperTEstimate:=proc(Tf,n,d)
local sd,s1,k;
sd:=binomial(n+d,d);
s1:=sum(binomial(n+k-2,k)*(1-binomial(sd-(d-k+1),Tf)/binomial(sd,Tf)),k=0..d);
return evalf(s1);
end proc:##### An example
G:=eval(F,x[n]=Roll()) mod p: ##evaluate F mod p at a random point
print("the actual numterms of G= ",numterms(G));
print("the estimated numterms of G by (6.2) = ",TEstimate(Tf,n,deg));#### Theoretical estimate after m random evaluations Eqn(8)TEAftermEvaluationsEstimate:=proc(Tf,n,m,d)
local sd,s1,k;
sd:=binomial(n+d,d);
s1:=sum(binomial(n-m-1+k,k)*(1-binomial(sd-binomial(d-k+m,m),Tf)/binomial(sd,Tf)),k=0..d);
return evalf(s1);
end proc:###An exampleG:=eval(F,{x[n-2]=Roll(), x[n-1]=Roll(), x[n]=Roll()}) mod p:
print("the actual numterms after 3 evals = ",numterms(G));
print("the estimated numterms after 3 evals by (6.7)= ",TEAftermEvaluationsEstimate(Tf,n,3,deg));##### Approximation to the theoretical estimate Eqn (9)FirstApproxEstimate:=proc(tf,n,d)
local temp1,temp2,k;
temp1:=binomial(n+d-1,d);
temp2:=sum(binomial(n+k-2,k)*(1-tf)^(d-k+1),k=0..d);
return 1.0*(temp1-temp2);
end proc:###An exampleG:=eval(F,x[n]=Roll()) mod p:
print("the actual numterms of G = ",numterms(G));
print("the estimated numterms of G by (6.2) = ",TEstimate(Tf,n,deg));
print("the approx to estimate to (6.2) by (6.4) = ",FirstApproxEstimate(tF,n,deg));
############ Example 8, Table (1) in the paper:
n:=7; ## the number of variables
deg:=15; ## the total degree
Tf:=17161; ## the number of terms
X:=[seq(x[i],i=1..n)]: ## variables
p:=2^31-1; ### prime number
F:=randpoly(X,degree=deg,terms=Tf,coeffs = rand(1..1000)):
Roll:=rand(1..p): ## random number generator
sd:=binomial(n+deg,deg): ## possible number of monomials
tF:=evalf(Tf/sd); ## the density ratio of F
G:=eval(F,x[n]=Roll()) mod p:
nt:=numterms(F):
tf:=evalf(nt/sd):
Tg:=numterms(G):
ETg:= TEstimate(Tf,n,deg):
eTg := FirstApproxEstimate(tF,n,deg):
tg:=evalf(Tg/binomial(n-1+degree(G),n-1)):
Etg:=ETg/binomial(n-1+degree(G),n-1):
etg:=eTg/binomial(n-1+degree(G),n-1):
print("the numterms of f, Tf = ", nt);
print("the density ratio of f, tf = ",tf);
print("the actual numterms of g, Tg= ",Tg);
print("the estimated numterms by (6.2) E[Tg]= ",ETg);
print("the approx to est numterms by (6.4) eTg = ", eTg);
print("the actual density ratio of g, tg = ", tg);
print("the estimated density ratio of g by (6.5), E[tg] = ",etg);
print("the approx to estimated density ratio of g by (6.5), etg =", etg);###################Example 9, Table (2) in the Paper##############n:=7;
n:=10;
#deg:=15;
deg:=20;
#Tf:=17958; ##The number of terms to be guessed
Tf:=20000;
X:=[seq(x[i],i=1..n)]:
Y:=X[1..n-1]:
p:=2^31-1;
f:=randpoly(X,degree=deg,terms=Tf,coeffs = rand(1..1000)):
Roll:=rand(1..p):
sd:=binomial(n+deg,deg); ### upper bound
tf:=evalf(Tf/sd);### m random evaluations....
m:=4:
EvalPoints:={seq(x[n-i]=Roll(), i=1..m)}:
g:=eval(f,EvalPoints) mod p: ### What we are given if the numterms Tg of g
Tg:=numterms(g):
tg:=evalf(Tg/binomial(n-m+degree(g),degree(g))):
gammam:=binomial(n-m+deg,deg): ### Here we assume that there is no degree loss
h:=1-(1/gammam)*sum(binomial(n-m-1+k,n-m-1)*(1-y)^binomial(deg-k+m,m),k=0..deg)-tg:
h2:=subs(y=1+z,h):
solution:=fsolve(h2,z,-1..0):
solution := solution +1:
E[t_f]:=solution:
E[T_f]:=E[t_f]*binomial(n+deg,deg):
print("The numterms of g, Tg = ", Tg);
print("The actual density ratio of g, tg = ", tg);
print("The guessed density ratio of f ,E[t_f] =", E[t_f]);
print("The actual density ratio of f ,t_f =", tf);
print("The guesses numterms of f, E[T_f] =", E[T_f]);
print("The actual numterms of f, T_f =", Tf);
############# Example 11, Table (3) in the paper (Sparse case, on Zippel's assumption)n:=12;
deg:=20;
T:=10^4:
X:=[seq(x[i],i=1..n)]:
p:=2^31-1;
F:=randpoly(X,degree=deg,terms=T,coeffs = rand(1..1000)):
T:=numterms(F);
Roll:=rand(1..p):
sd:=binomial(n+deg,deg):
tF:=evalf(T/sd);
f[0]:=F: tf[0]:=tF: Tf[0]:=T:print("i, Tf_i, E[Tf_i], tf_i, tf_i(1+d/n-i), Tf_i/Tf_(i+1), 1+d/(n-i+1)");
for i from 1 to n-1 do
alpha:=Roll();
f[i]:=eval(f[i-1],x[n-i+1]=alpha) mod p;
Tf[i]:=numterms(f[i]);
tf[i]:=1.0*numterms(f[i])/binomial(n-i+degree(f[i]),degree(f[i]));
print(i-1,Tf[i-1],TEAftermEvaluationsEstimate(Tf[0],n,i-1,deg),tf[i-1],tf[i-1]*(1+degree(f[0])/(n-i+1)),1.0*Tf[i-1]/Tf[i],(1+1.0*degree(f[0])/(n-i)));
od:############# Example 12, Table (4) in the paper (Dense case on Zippel's assumption)
n:=7;
deg:=13;
T:=7752:
X:=[seq(x[i],i=1..n)]:
p:=2^31-1;
F:=randpoly(X,degree=deg,terms=T,coeffs = rand(1..1000)):
T:=numterms(F);
Roll:=rand(1..p):
sd:=binomial(n+deg,deg):
tF:=evalf(T/sd);
f[0]:=F: tf[0]:=tF: Tf[0]:=T:for i from 1 to n-1 do
alpha:=Roll();
f[i]:=eval(f[i-1],x[n-i+1]=alpha) mod p;
Tf[i]:=numterms(f[i]);
tf[i]:=1.0*numterms(f[i])/binomial(n-i+degree(f[i]),degree(f[i]));
print(i-1,Tf[i-1],TEAftermEvaluationsEstimate(Tf[0],n,i-1,deg),tf[i-1],tf[i-1]*(1+degree(f[0])/(n-i+1)),1.0*Tf[i-1]/Tf[i],(1+1.0*degree(f[0])/(n-i)));
od:################# Example 19 (Table 5) in the papern:=7;
deg:=13;
r:=1.0*n/(n+deg):
Tf:=7752:
X:=[seq(x[i],i=1..n)];
p:=2^31-1;
f:=randpoly(X,degree=deg,terms=Tf,coeffs = rand(1..1000)):
Tf:=numterms(f);
Roll:=rand(1..p):
sd:=binomial(n+deg,deg):
tF:=evalf(Tf/sd);
alpha:=Roll();for i from 0 to deg-1 do
cf[i]:=coeftayl(f,x[n]=alpha,i) mod p; ## The actual i^th Taylor coefficient
ntcf[i]:=numterms(cf[i]); ## The actual number of terms of the i^th Taylor Coefficient
etf[i]:=1.0*ntcf[i]/binomial(n+deg-i,deg-i); ## The actual density ratio of the i^th Taylor coeff.
Fentcf[i]:=FirstApproxEstimate(tF,n,deg-i); ## The Estimation of numterms of the i^th Taylor coeff.
BoundOnTf[i]:=tF*binomial(n+deg-i,n); ## Bound on esti. of numterms of the i^th Taylor coeff.(6.15)
if i=0 then tfder[i]:=tF;
else
tfder[i]:=evalf(numterms(diff(f,(x[n]$(i))))/binomial(n+deg-i,n)); ## The actual dens ratio of the i^th der.
fi:
BoundOntf[i]:=tF*(n+deg-i)/n; ## Bound on esti. of dens. ratio of the i^th Taylor coeff.(6.15)
print(i,ntcf[i],Fentcf[i],BoundOnTf[i],tfder[i],etf[i],BoundOntf[i]);
od: