######################################################################
#GuessMulPol.txt: Save this file as  GuessMulPol.txt                 #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read  `GuessMulPol.txt`<Enter>                     # 
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: 
 
print(`Created: Jan. 8, 2020`):
print(` This is GuessMulPol.txt `):
print(`by written by Doron Zeilberger`):
print(`It guesses multi-variable polynomial expressions for sequences in terms of given sequences. `):
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 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: ApplyPolSeq, Comps1, Comps, CompsG, Hmn, GenPol,GenPolG `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  GuessPol, GuessPol1, GuessPol1G `):
 print(` `):


elif nops([args])=1 and op(1,[args])=ApplyPolSeq then
print(`ApplyPolSeq(P,x,LL): Inputs a polynomial P of x[1], ..., x[k] (where k=nops(LL)) and a list LL of length k`):
print(`consisting of the values of k sequences LL[1], ..., LL[k] all of length N, outputs`):
print(`the sequence of length N whose n-th entry is P(LL[1][n], ..., LL[k][n]). Try:`):
print(`ApplyPolSeq(x[1]+x[2]+2*1]*x[2],[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]);`):


elif nops([args])=1 and op(1,[args])=Comps then
print(`Comps(k,r): the set of vectors of non-negative integers [v1,..., vk] that adds up to <=r. Try:`):
print(`Comps(3,4);`):

elif nops([args])=1 and op(1,[args])=Comps1 then
print(`Comps1(k,r): the set of vectors of non-negative integers [v1,..., vk] that adds up to r. Try:`):
print(`Comps1(3,4);`):

elif nops([args])=1 and op(1,[args])=CompsG then
print(`CompsG(A,k,r): The set of vectors of non-negative integers with k componets such that the first component is between 0 and A`):
print(`and the remaining ones have sum <=r. Try:`):
print(`CompsG(5,3,1);`):

elif nops([args])=1 and op(1,[args])=GenPol then
print(`GenPol(x,k,r,a): A generic polynomial of total degree r in x[1], ..., x[k] with coefficients indexed by a,`):
print(`followed by the set of coefficients. Try:`):
print(`GenPol(x,3,3,a);`):

elif nops([args])=1 and op(1,[args])=GenPolG then
print(`GenPolG(x,k,A,r,a): A generic polynomial of total degree r in x[2], ..., x[k] and degree A in x[1],with coefficients indexed by a,`):
print(`followed by the set of coefficients. Try:`):
print(`GenPolG(x,3,3,1,a);`):

elif nops([args])=1 and op(1,[args])=GuessPol then
print(`GuessPol(LL,L1,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and`):
print(`a positive integer r, outputs a polynomial P , if it exists, such that `):
print(`L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Otherwise it returns FAIL. Try:`):
print(` GuessPol([[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]], ApplyPolSeq(x[1]+x[2],x,[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]),x);`):

elif nops([args])=1 and op(1,[args])=GuessPol1 then
print(`GuessPol1(LL,L1,r,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and`):
print(`a positive integer r, outputs a polynomial P of total degree <=r , if it exists, such that `):
print(`L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Otherwise it returns FAIL. Try:`):
print(` GuessPol1([[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]], ApplyPolSeq(x[1]+x[2],x,[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]),1,x);`):

elif nops([args])=1 and op(1,[args])=GuessPol1G then
print(`GuessPol1G(LL,L1,A,r,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and`):
print(`a positive integers A and r, outputs a polynomial P of degree<=A in x[1] and total degree <=r in {x[2], ..., x[r]}, if it exists, such that `):
print(`L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Otherwise it returns FAIL. Try:`):
print(` GuessPol1G([[seq(i,i=1..20)], [seq(Hnm(i,1),i=1..20)]], ApplyPolSeq(x[1]^3+x[2],x,[[seq(i,i=1..20)], [seq(Hnm(i,1),i=1..20)]]),3,1,x);`):

elif nops([args])=1 and op(1,[args])=Hnm then
print(`Hnm(n,m): add(1/i^m,i=1..n). Try: `):
print(`Hnm(10,2);`):


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




ez:=proc(): 
print(`GuessPol1(LL,L1,r,x)`):
print(`GuessPol1([[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]], ApplyPolSeq(x[1]+x[2],x,[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]),1,x);`):
print(`GuessPol(LL,L1,x)`):
print(`GuessPol([[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]], ApplyPolSeq(x[1]+x[2],x,[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]),x);`):
end:


Hnm:=proc(n,m) local i:add(1/i^m,i=1..n):end:
#Comps1(k,r): the set of vectors of non-negative integers [v1,..., vk] that adds up to r. Try:
#Comps1(3,4);
Comps1:=proc(k,r) local gu,i,mu,mu1:
option  remember:

if k=0 then
 if r=0 then
   RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:

gu:={}:

for i from 0 to r do
 mu:=Comps1(k-1,r-i):
 gu:=gu union {seq([op(mu1),i],mu1 in mu)}:
od:
gu:
end:


#Comps(k,r): the set of vectors of non-negative integers [v1,..., vk] that adds up to <=r. Try:
#Comps(3,4);
Comps:=proc(k,r) local r1:
{seq(op(Comps1(k,r1)),r1=0..r)}:
end:

#CompsG(A,k,r): The set of vectors of non-negative integers with k componets such that the first component is between 0 and A
#and the remaining ones have sum <=r. Try:
#CompsG(5,3,1);
CompsG:=proc(A,k,r) local a,v:
{seq(seq([a,op(v)], v in Comps(k-1,r)),a=0..A)}:
end:



#GenPol(x,k,r,a): A generic polynomial of total degree r in x[1], ..., x[k] with coefficients indexed by a,
#followed by the set of coefficients. Try:
#GenPol(x,3,3,a);
GenPol:=proc(x,k,r,a) local gu,i,j:
gu:=Comps(k,r):

[add(a[i]*mul(x[j]^gu[i][j],j=1..k),i=1..nops(gu)),{seq(a[i],i=1..nops(gu))}]:
end:


#GenPolG(x,k,A,r,a): A generic polynomial of total degree r in x[2], ..., x[k] and degree A in x[1],with coefficients indexed by a,
#followed by the set of coefficients. Try:
#GenPolG(x,3,3,1,a);
GenPolG:=proc(x,k,A,r,a) local gu,i,j:
gu:=CompsG(A,k,r):
[add(a[i]*mul(x[j]^gu[i][j],j=1..k),i=1..nops(gu)),{seq(a[i],i=1..nops(gu))}]:
end:

#ApplyPolSeq(P,x,LL): Inputs a polynomial P of x[1], ..., x[k] (where k=nops(LL)) and a list LL of length k
#consisting of the values of k sequences LL[1], ..., LL[k] all of length N, outputs
#the sequence of length N whose n-th entry is P(LL[1][n], ..., LL[k][n]). Try
#ApplyPolSeq(x[1]+x[2]+2*1]*x[2],[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]);
ApplyPolSeq:=proc(P,x,LL) local k,i,n,N:

k:=nops(LL):

if not  (type(LL,list)  and nops(LL)>0 and {seq(type(LL[i],list),i=1..k)}={true} and nops({seq(nops(LL[i]),i=1..k)})=1) then
 print(LL, `must be a non-empty list of lists all of the same length `):
 RETURN(FAIL):
fi:

N:=nops(LL[1]):
[seq(subs({seq(x[i]=LL[i][n],i=1..k)},P),n=1..N)]:
end:

#GuessPol1(LL,L1,r,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and
#a positive integer r, outputs a polynomial P of total degree <=r , if it exists, such that 
#L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Otherwise it returns FAIL. Try:
# GuessPol1([[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]], ApplyPolSeq(x[1]+x[2],x,[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]),1,x);
GuessPol1:=proc(LL,L1,r,x) local i,P,eq,var,gu,a,var1,N,k,n:
k:=nops(LL):

if not  (type(LL,list)  and nops(LL)>0 and {seq(type(LL[i],list),i=1..k)}={true} and nops({seq(nops(LL[i]),i=1..k)})=1) then
 print(LL, `must be a non-empty list of lists all of the same length `):
 RETURN(FAIL):
fi:

N:=nops(LL[1]):

if not (type(L1,list) and nops(L1)=N )then
 print(L1, `must be a list of numbers of length`, N):
 RETURN(FAIL):
fi:

gu:=GenPol(x,k,r,a):

P:=gu[1]: var:=gu[2]:

if N-nops(var)<=5 then
 print(`The length of the sequences must be at least`, nops(var)+6 ):
 RETURN(FAIL):
fi:

eq:={seq(subs({seq(x[i]=LL[i][n],i=1..k)},P)-L1[n],n=1..N)}:
var1:=solve(eq,var):

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

subs(var1,P):


end:



#GuessPol(LL,L1,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and
#a positive integer r, outputs a polynomial P such that 
#L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Try:
# GuessPol([[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]], ApplyPolSeq(x[1]+x[2],x,[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]),x);
GuessPol:=proc(LL,L1,x) local r,gu,N,i,k:
k:=nops(LL):

if not  (type(LL,list)  and nops(LL)>0 and {seq(type(LL[i],list),i=1..k)}={true} and nops({seq(nops(LL[i]),i=1..k)})=1) then
 print(LL, `must be a non-empty list of lists all of the same length `):
 RETURN(FAIL):
fi:

N:=nops(LL[1]):

if not (type(L1,list) and nops(L1)=N )then
 print(L1, `must be a list of numbers of length`, N):
 RETURN(FAIL):
fi:


for r from 0 while nops(Comps(k,r))<=N-6 do
gu:=GuessPol1(LL,L1,r,x):
if gu<>FAIL then
 RETURN(gu):
fi:
od:

FAIL:
end:




#GuessPol1G(LL,L1,A,r,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and
#a positive integers A and r, outputs a polynomial P of degree<=A in x[1] and total degree <=r in {x[2], ..., x[r]}, if it exists, such that 
#L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Otherwise it returns FAIL. Try:
# GuessPol1G([[seq(i,i=1..20)], [seq(Hnm(i,1),i=1..20)]], ApplyPolSeq(x[1]^3+x[2],x,[[seq(i,i=1..20)], [seq(Hnm(i,1),i=1..20)]]),3,1,x);
GuessPol1G:=proc(LL,L1,A,r,x) local i,P,eq,var,gu,a,var1,N,k,n:
k:=nops(LL):

if not  (type(LL,list)  and nops(LL)>0 and {seq(type(LL[i],list),i=1..k)}={true} and nops({seq(nops(LL[i]),i=1..k)})=1) then
 print(LL, `must be a non-empty list of lists all of the same length `):
 RETURN(FAIL):
fi:

N:=nops(LL[1]):

if not (type(L1,list) and nops(L1)=N )then
 print(L1, `must be a list of numbers of length`, N):
 RETURN(FAIL):
fi:

gu:=GenPolG(x,k,A,r,a):

P:=gu[1]: var:=gu[2]:

if N-nops(var)<=5 then
 print(`The length of the sequences must be at least`, nops(var)+6 ):
 RETURN(FAIL):
fi:

eq:={seq(subs({seq(x[i]=LL[i][n],i=1..k)},P)-L1[n],n=1..N)}:
var1:=solve(eq,var):

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

subs(var1,P):


end: