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


print(`First Written: Feb. 2024: tested for Maple 2020 `):
print(`Version  : Feb. 8, 2024:  `):
print():
print(`This is Jackmond.txt, A Maple package`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(`"The Jackson-Richmod 4CT constant is EXACTLY 10/27" `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Jackmond.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
 



ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: GR, GR1`):
print(``):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` Jackmond.txt: A Maple package for  `):

print(`The MAIN procedures are: JR, Paper `):
print(``):

elif nargs=1 and args[1]=GR then
print(`GR(L,n): inputs a list of numbers and guesses a rational function in n for them such that its value at i is L[i]. Try:`):
print(`GR([seq((i+1)/(i+2),i=1..20)],n);`):

elif nargs=1 and args[1]=GR1 then
print(`GR1(L,n,d1,d2): inputs a list of numbers and guesses a rational function in n, with numerator degree d1 and denominator degree d2, for them such that its value at i is L[i]. Try:`):
print(`GR1([seq((i+1)/(i+2),i=1..20)],n,1,1);`):


elif nargs=1 and args[1]=JR then
print(`JR(f,n,r,d): inputs an expression f in n, a positive integer r and a positive integer d`):
print(`guesses a rational function of degree d for [x^n] F(x)^r/f(n) where`):
print(`F(x)=Sum(f(i)*x^i,i=1..infinity). For example, to get the Lemma in the paper (the original case needed by Bruce Richmond and David Jackson, type:`):
print(`JR(2/(n+1)!*(4*n+1)!/(3*n+2)!,n,2,4);`):


elif nargs=1 and args[1]=Paper then
print(`Paper(f,n,Maxd): inputs an expression f of n, and guesses explicit expressions for`):
print(`Coefficient of x^n in Sum(f(i)*x^i,i=1..infinity)^r divided by f(n). Try:`):
print(`Paper(2/(n+1)!*(4*n+1)!/(3*n+2)!,n,10);`):

else
 print(`There is no such thing as`, args):

fi:


end:




#GR1(L,n,d1,d2): inputs a list of numbers and guesses a rational function in n, with numerator degree d1 and denominator degree d2, for them such that its value at i is L[i]. Try:
#GR1([seq((i+1)/(i+2),i=1..20)],n,1,1);
GR1:=proc(L,n,d1,d2) local a,b,eq,var,i,f:

if nops(L)<=(d1+2)*(d2+2)+3 then
 print(`Insufficient data`):
 RETURN(FAIL):
fi:
 

f:=(n^d1+add(a[i]*n^i,i=0..d1-1))/add(b[i]*n^i,i=0..d2):

eq:={seq(numer(subs(n=i,f)-L[i]),i=1..(d1+2)*(d2+2)+3)}:
var:={seq(a[i],i=0..d1-1),seq(b[i],i=0..d2)}:
var:=solve(eq,var):
if var=NULL then
 RETURN(FAIL):
fi:


f:=factor(subs(var,f)):
if {seq(subs(n=i,f)-L[i],i=(d1+2)*(d2+2)+3..nops(L))}<>{0} then
 print(f, `didn't work out`):
RETURN(FAIL):
fi:
factor(f):

end:

#GR(L,n): inputs a list of numbers and guesses a rational function in n for them such that its value at i is L[i]. Try:
#GR([seq((i+1)/(i+2),i=1..20)],n);
GR:=proc(L,n) local d1, d2,K,K1,f:
K:=trunc(sqrt(nops(L)-3))-2:


for K1 from 0 to K do
 for d1 from 0 to K1 do
 d2:=K-d1:

 f:=GR1(L,n,d1,d2):

 if f<>FAIL then
   RETURN(f):
 fi:
 od:
od:
FAIL:
end:





#JR(f,n,r,d): inputs an expression f in n, a positive integer r and a positive integer d
#guesses a rational function of degree d for [x^n] F(x)^r/f(n) where
#F(x)=Sum(f(i)*x^i,i=1..infinity):
#JR(2/(n+1)!*(4*n+1)!/(3*n+2)!,n,2,4);
JR:=proc(f,n,r,d) local F,i,L,L1,N,x,g2:
N:=(d+2)^2+5:
F:=add(eval(subs(n=i,f))*x^i,i=1..N):
L:=[seq(eval(subs(n=i,f)),i=1..N)]:
g2:=expand(F^r):
L1:=[seq(coeff(g2,x,i),i=1..N)]:
GR1([seq(L1[i]/L[i],i=1..N)],n,d,d):
end:



#Paper(f,n,Maxd): inputs an expression f of n, and guesses explicit expressions for
#Coefficient of x^n in Sum(f(i)*x^i,i=1..infinity)^r divided by f(n). Try:
#Paper(2/(n+1)!*(4*n+1)!/(3*n+2)!,n,10);

Paper:=proc(f,n,Maxd) local gu,x,r,r1,gu1,lu,A,B,t0:
t0:=time():
gu:=[]:
lu:=[]:
gu1:=JR(f,n,2,Maxd):

if gu1=FAIL then
 RETURN([]):
fi:
gu:=[gu1]:
lu:=[lcoeff(numer(gu1),n)/lcoeff(denom(gu1),n)]:
for r1 from 3 while gu1<>FAIL do
 gu1:=JR(f,n,r1,Maxd):
 if gu1<>FAIL then
  gu:=[op(gu),gu1]:
lu:=[op(lu),lcoeff(numer(gu1),n)/lcoeff(denom(gu1),n)]:
 fi:
od:



print(`Explicit Expressions for the coefficient of x^n in`, (Sum(f*x^n,n=1..infinity))^r, `divided by`, f, `for r from 2 to`, nops(gu)+1):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
print(`Theorem: `):
print(``):
print(`Let A[r](n) be the coefficient of x^n in `,  (Sum(f*x^n,n=1..infinity))^r, `divided by `, f):
print(``):
print(`Let B[r] be the limit of A[r][n] as n goes to infinity`):
print(`then we have the following`):
print(``):
for r1 from 2 to nops(gu)+1 do
print(``):
print( A[r1,n], ` =`  ):
print(``):
print(gu[r1-1]):
print(``):
print( B[r1], `= `):
print(``):
print(lu[r1-1]):
print(``):
print(`and in Maple notation`):
lprint(gu[r1-1]):
print(``):
lprint(lu[r1-1]):
od:

print(``):
print(`To sum up, here are the`, A[r,n], `for f from 2 to `, nops(gu)+1):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):


print(``):
print(`and here are the`, B[r], `for f from 2 to `, nops(gu)+1):
print(``):
print(lu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu):

print(``):

print(`-----------------------------------------`):

print(``):
print(`This ends this paper that took`, time()-t0, `seconds to generate`):

end: