######################################################################
##Dresner.txt: Save this file as  Dresner.txt                        #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read  Dresner.txt<Enter>                           #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Sept. 2018

print(`VERY PRELIMINARY`):
 
print(`Created:  Sept. 2018`):
print(` This is Dresner.txt `):
print(`It is one of the packages that accompany the article `):
print(`A tale of Magnets and Logic`):
print(`by  Doron Zeilberger and Noam Zeilberger`):

print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Extract1, GenPol, OneStep, TruncF `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:a, GuessDEz, GuessDEzx, Seqa, VerifyNoam `):
 print(` `):



elif nops([args])=1 and op(1,[args])=a then
print(`a(n): the n-th term of A267827 using Dresner's differential equation. Try:`):
print(`a(10); `):

elif nops([args])=1 and op(1,[args])=Extract1 then
print(`Extract1(POL,kList): extracts all the coeffs. of a POL in kList. Try:`):
print(`Extract1(a*x+b*y,[x,y]); `):

elif nops([args])=1 and op(1,[args])=GenPol then
print(`GenPol(kList,a,deg): The generic polynomial in kList of`):
print(`degree deg,  using the indexed variable a, followed by the set`):
print(`GenPol([x,y],A,2);`):



elif nops([args])=1 and op(1,[args])=GuessDEz then
print(`GuessDEz(P,z,F,Fz,k): inputs a polynomial P in the variable z, the variable z, symbols F and Fz, and an inteter k. Gueses a non-linear`):
print(`first-order differential equation of degree k in (z,F,Fz) that annihilated P. Try:`):
print(`GuessDEz(subs(x=0),TruncF(x,z,40)),z,F, Fz,4);`):

elif nops([args])=1 and op(1,[args])=GuessDEzx then
print(`GuessDEzx(P,z,x,F,Fz,Fx,k): inputs a polynomial P in the variable z, the variable z, symbols F and Fz, and an integer k. Gueses a non-linear`):
print(`first-order differential equation of degree k in (x,z,F,Fz,Fx) that annihilated P. Try:`):
print(`GuessDEzx(TruncF(x,z,30),z,x,F, Fz,Fx,2);`):

elif nops([args])=1 and op(1,[args])=Seqa then
print(`Seqa(N): the first N terms of A267827 using Dresner's differential equation. Try:`):
print(`a(10); `):

elif nops([args])=1 and op(1,[args])=TruncF then
print(`TruncF(x,z,N): The truncation of F(x,z) up to z^N. Try:`):
print(`TruncF(x,z,4);`):

elif nops([args])=1 and op(1,[args])=VerifyNoam then
print(`checks that Noam implies Dresner for the first N terms. Try:`):
print(`VerifyNoam(30);`):

else
print(`There is no ezra for`,args):
fi:
 
end:

##From Z
#Extract1(POL,kList): extracts all the coeffs. of a POL in kList
Extract1:=proc(POL,kList) local k1,kList1,i:
k1:=kList[nops(kList)]:
kList1:=[op(1..nops(kList)-1,kList)]:
if nops(kList)=1 then
 RETURN({seq(coeff(POL,k1,i),i=0..degree(POL,k1))}):

else
 RETURN({seq(
op(Extract1(coeff(POL,k1,i),kList1)), i=0..degree(POL,k1))}):
fi:

end:


#IV1(d,n): all the vectors of non-negative integres of length d whose
#sum is n
IV1:=proc(d,n)
local gu,i,gu1,i1:
if d=0 then
 if n=0 then
   RETURN({[]}):
 else
    RETURN({}):
fi:
fi:
 
gu:={}:

for i from 0 to n do
 gu1:=IV1(d-1,n-i):
 gu:=gu union {seq([op(gu1[i1]),i],i1=1..nops(gu1))}:
od:
gu:
end:

#IV(d,n): all the vectors of non-negative integres of length d whose
#sum is <=n
IV:=proc(d,n) local i: {seq(op(IV1(d,i)),i=0..n)}: end:

#GenPol(kList,a,deg): The generic polynomial in kList of
#degree deg,  using the indexed variable a, followed by the set
#of coeffs.

GenPol:=proc(kList,a,deg)
local gu,i,i1:
gu:=IV(nops(kList),deg):
convert([seq(a[i]*convert([seq(kList[i1]^gu[i][i1],i1=1..nops(kList))],`*`),
i=1..nops(gu))],`+`),{seq(a[i],i=1..nops(gu))}:

end:
##End From Z

#Dresner A267827
ez:=proc(): print(`  OneStep(F,x,z) , Seqax(N,x), SeqaxO(N,x), VerifyDE(N), GuessDE(L,z,d),  GuessDEg(L,z,d), ax(n,x) `):
print(`VerifyNoam(N), PN(F,Fz,Fx,Fxz,z), bx(n,x) `):
 end:

a:=proc(n) local k:
option remember:

if n=1 then
 2:
else
(6*n-2)*a(n-1)+add((6*k+2)*a(k)*a(n-1-k),k=1..n-2):
fi:
end:

Seqa:=proc(N) local i:
[1,seq(a(i),i=1..N-1)]:
end:


OneStep:=proc(F,x,z)
expand(x+z*(F-subs(x=0,F))^2+z*diff(F,x)):
end:

#SeqaxO(N,x): The first N terms of the odd coefficients, in z,  in F(x,z)
SeqaxO:=proc(N,x) local F,i,j,z:
F:=z:
for i from 2 to 2*N+1 do
F:=OneStep(F,x,z):
F:=add(coeff(F,z,j)*z^j,j=0..i):
od:

[seq(coeff(F,z,2*i-1),i=1..N)]:
end:



#SeqaxOld(N,x): The first N terms in the coefficients of z in F(x,z)
Seqax:=proc(N,x) local F,i,j,z:
F:=z:
for i from 2 to N+3 do
F:=OneStep(F,x,z):
F:=add(coeff(F,z,j)*z^j,j=0..i):
od:

[seq(coeff(F,z,i),i=1..N)]:
end:


VerifyDE:=proc(N) local gu,z,f,i,lu:
gu:=Seqa(N):

f:=add(gu[i+1]*z^i,i=0..N-1):
lu:=expand(6*z^2*f*diff(f,z)+2*z*f^2-f+1):

{seq(coeff(lu,z,i),i=1..N-1)}:
end:





ax:=proc(n,x) local k:
option remember:


if n=0 then
 x:
else
expand(add( (ax(k,x)-subs(x=0,ax(k,x)))*(ax(n-1-k,x)-subs(x=0,ax(n-k-1,x))),k=0..n-1)+diff(ax(n-1,x),x)):
fi:


end:


#Seqax(N,x): The first N terms in the coefficients of z in F(x,z)
Seqax:=proc(N,x) local i:
[seq(ax(i,x),i=0..N-1)]:
end:

#VerifyNoam(N): checks that Noam implies Dresner for the first N terms
VerifyNoam:=proc(N) local gu,i,x:
gu:=subs(x=0,Seqax(2*N,x)):
evalb([seq(gu[2*i],i=1..N)]=Seqa(N)):
end:



#PN(F,Fz,Fx,Fxz,z)
PN:=proc(F,Fz,Fx,Fxz,z) local x,gu1,gu2,gu3,y,yz:
gu1:=y-z+z*y^2-3*z^2*y*yz:
gu2:=expand(F-x-z*(F-y)^2-z*Fx):
gu3:=expand(Fz-(F-y)^2-2*z*(F-y)*(Fz-yz)-Fx-z*Fxz):
factor(subs(x=0,eliminate({gu1,gu2,gu3},{y,yz}))):
end:


#GuessDE(L,z,d): inputs a list L (starting at 0), and guesses a non-linear differential equation
#of the form f(z)=R(z)+P(z)*f(z)^2+ Q(z)*f(z)*f'(z), where R(z), P(z), Q(z) are polynomials in z of degree <=d
#satisfied by f(z)=add(L]i]*z^i,i=0..d).  It returns the pair [R(z),P(z),Q(z)]
#Try:
#GuessDEnew(Seqa(20),z,2);
GuessDEnew:=proc(L,z,d) local f,i,P,Q,R,a,b,c,eq,var,lu:

if nops(L)<=3*d+6 then
 print(`L must be at least`, 3*d+7, `long `):
 RETURN(FAIL):
fi:


f:=add(L[i+1]*z^i,i=0..nops(L)-1):

P:=add(a[i]*z^i,i=0..d):
Q:=add(b[i]*z^i,i=0..d):
R:=add(c[i]*z^i,i=0..d):
var:={seq(a[i],i=0..d), seq(b[i],i=0..d) ,  seq(c[i],i=0..d) }:

lu:=expand(Q*diff(f,z)*f+P*f^2+R-f):

eq:={seq(coeff(lu,z,i),i=0..nops(L)-1)}:

var:=solve(eq,var):

if var=NULL then
  RETURN(FAIL):
fi:

[subs(var,R),subs(var,P),subs(var,Q)]:

end:

#GuessDEold(L,z,d): inputs a list L (starting at 0), and guesses a non-linear differential equation
#of the form f(z)=R(z)+P(z)*f(z)^2+ Q(z)*f(z)*f'(z), where R(z), P(z), Q(z) are polynomials in z of degree <=d
#satisfied by f(z)=add(L]i]*z^i,i=0..d).  It returns the pair [R(z),P(z),Q(z)]
#Try:
#GuessDEold(Seqa(20),z,2);
GuessDEold:=proc(L,z,d) local f,i,P,Q,R,a,b,c,eq,var,lu:

if nops(L)<=3*d+6 then
 print(`L must be at least`, 3*d+7, `long `):
 RETURN(FAIL):
fi:


f:=add(L[i+1]*z^i,i=0..nops(L)-1):

P:=add(a[i]*z^i,i=0..d):
Q:=add(b[i]*z^i,i=0..d):
R:=add(c[i]*z^i,i=0..d):
var:={seq(a[i],i=0..d), seq(b[i],i=0..d) ,  seq(c[i],i=0..d) }:

lu:=expand(Q*diff(f,z)*f+P*f^2+R-f):

eq:={seq(coeff(lu,z,i),i=0..nops(L)-1)}:

var:=solve(eq,var):

if var=NULL then
  RETURN(FAIL):
fi:

[subs(var,R),subs(var,P),subs(var,Q)]:

end:




b:=proc(n) local k:
option remember:
if n=0 then
 0:
elif n=1 then
1:
else
-add(b(k)*b(n-1-k),k=0..n-1)+3*add((k+1)*b(k+1)*b(n-2-k),k=0..n-2):
fi:
end:

bx:=proc(n,x) local k:
option remember:


if n=0 then
 x:
else
expand(add( (ax(k,x)-b(k))*(ax(n-1-k,x)-b(n-1-k) ),k=0..n-1)+diff(bx(n-1,x),x)):
fi:


end:



#TruncF(x,z,N): The truncation of F(x,z) up to z^N. Try:
#TruncF(x,z,4);
TruncF:=proc(x,z,N) local i:
add(ax(i,x)*z^i,i=0..N):
end:

#GuessDEz(P,z,F,Fz,k): inputs a polynomial P in the variable z, the variable z, symbols F and Fz, and an inteter k. Gueses a non-linear
#first-order differential equation of degree k in (z,F,Fz) that annihilated P. Try:
#GuessDEz(subs(x=0),TruncF(x,z,20)),z,F, Fz,4);
GuessDEz:=proc(P,z,F,Fz,k) local gu,mu,var,eq,var1,a,i,var11,var2:
gu:=GenPol([z,F,Fz],a,k);
var:=gu[2]:
gu:=gu[1]:

mu:=expand(subs({F=P,Fz=diff(P,z)},gu)):
eq:={seq(coeff(mu,z,i),i=0..degree(P,z)-k)}:

if nops(eq)-nops(var)<10 then
 print(`Make the first argument bigger.`):
 RETURN(FAIL):
fi:


var1:=solve(eq,var):


gu:=subs(var1,gu):


if gu=0 then
 RETURN(FAIL):
fi:

var2:={}:

for var11 in var1 do
 if op(1,var11)=op(2,var11) then
   var2:=var2 union {op(1,var11)}:
 fi:
od:


if nops(var2)>1 then
 print(`There is more than one solution`):
 RETURN({seq(coeff(gu,var11,1),var11 in var2)}):
fi:

if nops(var2)=1 then
 RETURN(coeff(gu,var2[1],1)):
else
 RETURN(FAIL):
fi:

end:




#GuessDEzx(P,z,x,F,Fz,Fx,k): inputs a polynomial P in the variable z, the variable z, symbols F and Fz, and an integer k. Gueses a non-linear
#first-order differential equation of degree k in (x,z,F,Fz,Fx) that annihilated P. Try:
#GuessDEzx(TruncF(x,z,30),z,x,F, Fz,Fx,2);
GuessDEzx:=proc(P,z,x,F,Fz,Fx,k) local gu,mu,var,eq,var1,a,i,var11,var2,j:
gu:=GenPol([z,x,F,Fz,Fx],a,k);
var:=gu[2]:
gu:=gu[1]:

mu:=expand(subs({F=P,Fz=diff(P,z),Fx=diff(P,x)},gu)):
eq:={seq(seq(coeff(coeff(mu,z,i),x,j),i=0..degree(P,z)-k),j=0..degree(P,x)-k)}:

if nops(eq)-nops(var)<10 then
 print(`Make the first argument bigger.`):
 RETURN(FAIL):
fi:


var1:=solve(eq,var):


gu:=subs(var1,gu):


if gu=0 then
 RETURN(FAIL):
fi:

var2:={}:

for var11 in var1 do
 if op(1,var11)=op(2,var11) then
   var2:=var2 union {op(1,var11)}:
 fi:
od:


if nops(var2)>1 then
 print(`There is more than one solution`):
 RETURN({seq(coeff(gu,var11,1),var11 in var2)}):
fi:

if nops(var2)=1 then
 RETURN(coeff(gu,var2[1],1)):
else
 RETURN(FAIL):
fi:

end: