# 
# 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;

printf("### Eqn (3) in the paper\n");
;
TEstimate:=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));








printf("### Theoretical estimate after m random evaluations Eqn(8)\n");
;

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 example
;
G:=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));





printf("### Approximation to the theoretical estimate Eqn (9)\n");
;

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 example
;
G:=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));


printf("###  Example 8, Table (1) in the paper\n");

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);


printf("###Example 9, Table (2)  in the Paper \n");
;
#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);





printf("### Example 11,  Table (3) in the paper (Sparse case, on Zippel's assumption)\n");
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:
printf("### Example 12,  Table (4) in the paper (Dense case on Zippel's assumption)\n");
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:
printf("### Example 19 (Table 5) in the paper\n");
n:=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: