######################################################################
## DET: Save this file as DET to use it,                             #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read DET: <Enter>                                 #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  zeilberg@math.rutgers.edu.                                       # 
######################################################################

with(linalg):

with(SolveTools);

print(`First Written: June 2005: tested for Maple 10 `):
print(`Version of Feb. 17, 2006:  `):
print():
print(`This is DET, A Maple package`):
print(`accompanying Doron Zeilberger's  article: `):
print(`The Holonomic Ansatz II: Automatic Determinant Evaluation.`):
print(`It guesses and then proves determinant-evaluations of`):
print(` holonomic matrices.`): 
print():
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/DET .`):
print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`):
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(` DET: A Maple package for automatically guesing and then proving Holnomic `):
print(`Determinant Evaluations. The supporting procedures are`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(` The supporting  procedures are:  CheckZug, CheckZugI `):
print(` Dave, DaveH, DaveHI, DaveI, DetRat,  DetRatDirect, DetRatI, `):
print(`DetRatSeq, DetRatSeqI, Findrec, findrec , OpeToCF `):
print(` PtorOpeR, SeqFromRec, Zug , ZugI, ZugP  `):
else
 ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then

print(` DET: A Maple package for automatically guesing and then proving Holnomic `):
print(`Determinant Evaluations. The MAIN procedures are`):

print(`  DaveH, DaveHP, DetRatRec, DetRatRecI, DetRatRecP, DetRatRecPI `):
print(`  Rproof, RproofP, SRproof, SRproofI, SRproofP `):


elif nargs=1 and args[1]=CheckZug then
print(`CheckZug(a,m,n,N,Z,K): given an expression a of m and n`):
print(`checks whether the winning pair Z agrees up to K as describing`):
print(`the determiant sequence of the matrix a(m,n)`):
print(`For example, try: Z:=Zug(1/(m+n+1),m,n,N,20);`):
print(`CheckZug(1/(m+n+1),m,n,N,Z,100);`):

elif nargs=1 and args[1]=CheckZugI then
print(`CheckZugI(a,m,n,N,Z,K): given an expression a of m and n`):
print(`checks whether the winning pair Z agrees up to K as describing`):
print(`the determiant sequence`):
print(`the determiant sequence of the matrix delta(m,n)+a(m,n)`):
print(`For example, try: Z:=ZugI(binomial(m+n,n),m,n,N,20);`):
print(`CheckZugI(binomial(m+n,n),m,n,N,Z,30);`):

elif nargs=1 and args[1]=DaveH then
print(`DaveH(a,m,n,N,Max): Given an expression a (of  m and n)`):
print(`guesses a recurrence operator P(m,n,M) `):
print(`annihilating A(m,n):= the cofactor of the (m,n) entry `):
print(`of the (m+1) by (m+1) matrix a(i,j), 0<=i,j<=m`):
print(`Max (>20) is the largest row you are willing to try`):
print(`Here N is the shift in n`):
print(`For example, try:`):
print(`DaveH(1/(m+n+1),m,n,N,20);`):

elif nargs=1 and args[1]=DaveHI then
print(`DaveHI(a,m,n,N,Max): Given an expression a (of  m and n)`):
print(`guesses a recurrence operator P(m,n,M) `):
print(`annihilating A(m,n):= the cofactor of the (m,n) entry `):
print(`of the (m+1) by (m+1) matrix (delta(i,j)+a(i,j)), 0<=i,j<=n`):
print(`Max (>20) is the largest row you are willing to try`):
print(`Here N is the shift in n`):
print(`For example, try:`):
print(`DaveHI(binomial(m+n,m),m,n,N,20);`):


elif nargs=1 and args[1]=DaveHP then
print(`DaveHP(a,m,n,N,Max,p,S1,S2): Given an expression a (of  m and n)`):
pritnt(`that also depends on a parameter p`):
print(`guesses a recurrence operator P(m,n,M) `):
print(`annihilating A(m,n):= the cofactor of the (m,n) entry `):
print(`of the (m+1) by (m+1) matrix a(i,j), 0<=i,j<=mn`):
print(`Max (>20) is the largest row you are willing to try`):
print(`Here N is the shift in n`):
print(`S1 and S2 are the starting and ending trial points for p`):
print(`For example, try:`):
print(`DaveHP(1/(m+n+p),m,n,N,20,p,3,13);`):

elif nargs=1 and args[1]=DetRat then
print(`DetRat(a,m,n,m0): det(a(i,j), 0<=i,j<=m0)/det(a(i,j), 0<=i,j<=m0-1)`):
print(`where a is an expression in m and n`):
print(`For example, try: DetRat(1/(m+n+1),m,n,10);`):

elif nargs=1 and args[1]=DetRatDirect then
print(`DetRatDirect(a,m,n,m0): det(a(i,j), 0<=i,j<=m0)/det(a(i,j), 0<=i,j<=m0-1)`):
print(`where a is an expression in m and n, computed directly`):
print(`For example, try: DetRatDirect(1/(m+n+1),m,n,10);`):

elif nargs=1 and args[1]=DetRatI then
print(`DetRatI(a,m,n,m0): det(b(i,j), 0<=i,j<=m0)/det(b(i,j), 0<=i,j<=m0-1)`):
print(`where b(i,j)=a(i,j)+delta(i,j), and delta(i,j) is the identity matrix`):
hprint(` (Kroneker delta), and  a is an expression in m and n`):
print(`For example, try: DetRatI(-binomial(20+m+n,n),m,n,10);`):

elif nargs=1 and args[1]=DetRatRec1 then
print(`DetRatRec1(a,m,n,N,Max) : Conjectures a First-Order recurrence `):
print(` , of degree (in n) <=Max, for`):
print(`det(a(i,j), 0<=i,j<=n)/det(a(i,j), 0<=i,j<=n-1) `):
print(`If there is no such recurrence, then it returns FAIL.`):
print(`For example, try:`):
print(`DetRatRec1(1/(m+n+1),m,n,N,10);`):

elif nargs=1 and args[1]=DetRatRec then
print(`DetRatRec(a,m,n,N,Max) : Conjectures a linear recurrence `):
print(`operator of the form ope(n,N), where N is the shift in n and`):
print(` with ORDER+ DEGREE <=Max, (not necessairly first-order) annihilating`):
print(`det(a(i,j), 0<=i,j<=n)/det(a(i,j), 0<=i,j<=n-1) `):
print(`If there is no such recurrence, then it returns FAIL.`):
print(`For example, try:`):
print(`DetRatRec(1/(m+n+1),m,n,N,10);`):

elif nargs=1 and args[1]=DetRatRecI then
print(`DetRatRecI(a,m,n,N,Max) : Conjectures a linear recurrence `):
print(`operator of the form ope(n,N), where N is the shift in n and`):
print(` , ORDER+ DEGREE <=Max, (not necessairly first-order) annihilating`):
print(`det(delta(i,j)+a(i,j), 0<=i,j<=n)/det(delta(i,j)+a(i,j), 0<=i,j<=n-1) `):
print(`If there is no such recurrence, then it returns FAIL.`):
print(`For example, try:`):
print(`DetRatRecI(binomial(m+n,m),m,n,N,10);`):

elif nargs=1 and args[1]=DetRatRecP then
print(`DetRatRecP(a,m,n,N,Max,p,S1,S2): Like DetRatRec (q.v.) but the `):
print(` expression `):
print(`a(i,j) also depends on a parameter p, and it tries by interpolating`):
print(`from DetRatRec with p=S1..S2`):
print(`For example, try:`):
print(`DetRatRecP(1/(m+n+p),m,n,N,7,p,1,20);`):

elif nargs=1 and args[1]=DetRatRecPI then
print(`DetRatRecPI(a,m,n,N,Max,p,S1,S2): Like DetRatRecP (q.v.) but the `):
print(` expression `):
print(`a(i,j) also depends on a parameter p, and it tries by interpolating`):
print(`from DetRatRec with p=S1..S2`):
print(`For example, try:`):
print(`DetRatRecPI(binomial(p+m+n,n),m,n,N,10,p,1,20);`):

elif nargs=1 and args[1]=DetRatSeq then
print(`DetRatSeq(a,m,n,n0): the sequence of DetRat for i from 1 to n0`):
print(`For example, try: DetRatSeq(1/(m+n+1),m,n,10);`):

elif nargs=1 and args[1]=DetRatSeqI then
print(`DetRatSeqI(a,m,n,n0): the sequence of DetRatI for i from 1 to n0`):
print(`For example, try: DetRatSeqI(binomial(m+n,n),m,n,10);`):

elif nargs=1 and args[1]=Dave then
print(`Dave(a,m,n,n0): Given an expression a of m and n`):
print(`outputs the cofactors of the last row`):
print(`(normalized with the LAST being 1) of the`):
print(`(n+1) by (n+1) matrix a(i,j) (0<=i,j<=n)`):
print(`For example, try:`):
print(`Dave(1/(m+n+1),m,n,10);`):

elif nargs=1 and args[1]=DaveI then
print(`DaveI(a,m,n,n0): Given an expression a of m and n`):
print(`outputs the cofactors of the last row`):
print(`(normalized with the LAST being 1) of the`):
print(`(n+1) by (n+1) matrix (delta(i,j)+a(i,j) (0<=i,j<=n)`):
print(`For example, try:`):
print(`DaveI(1/(m+n+1),m,n,10);`):

elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the INTEGER-sequence f of degree DEGREE and order ORDER.`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f of INTEGERS tries to find a linear recurrence equation with`):
print(`poly coeffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):


elif nargs=1 and args[1]=OpeToCF then
print(`OpeToCF(ope,n,N): converts ope(n,N) of first order to Closed-Form with f(0)=1`):
print(`For example, try: opeToCF(N-n-1,n,N);`):


elif nargs=1 and args[1]=PtorOpeR then
print(`PtorOpeR(ope,n,N,R): solving the recurrence ope(n,N)x(n)=0, x(0)=1`):
print(`in terms of rising factorial R(n,k) (R is a symbol)`):
print(`For example try PtorOpeR((n+1)*N-1,n,N,R);`):

elif nargs=1 and args[1]=Rproof then
print(`Rproof(a,m,n,N,Max,R): Given an expression a of m and n`):
print(`outputs the determinant ratios`):
print(` [det(a(i,j), 0<=i,j<=n)]/[det(a(i,j), 0<=i,j<=n-1)]`):
print(`if it happens to be closed-form`):
print(` and its rigorous, WZ-style, proof, if available.`):
print(` Max is the trial paramter, and R is the symbol `):
print(`for the rising factorial`):
print(`For example, try:`):
print(`Rproof(1/(m+n+1),m,n,N,30,R);`):

elif nargs=1 and args[1]=RproofP then
print(`RproofP(a,m,n,N,Max,R,p,S1,S2): Given an expression a of m and n`):
print(`that also depends on a paramter p`):
print(`outputs the determinant ratios`):
print(` [det(a(i,j), 0<=i,j<=n)]/[det(a(i,j), 0<=i,j<=n-1)]`):
print(`if it happens to be closed-form`):
print(` and its rigorous, WZ-style, proof, if available.`):
print(`It tries the range p=S1..S2 `):
print(`For example, try:`):
print(`RproofP(1/(m+n+p),m,n,N,20,R,p,1,10);`):

elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try: SeqFromRec(N-n-1,n,N,[1],10);`):

elif nargs=1 and args[1]=SRproof then
print(`SRproof(a,m,n,N,Max,K): Given an expression a of m and n`):
print(`outputs the determinant ratios`):
print(` [det(a(i,j), 0<=i,j<=n)]/[det(a(i,j), 0<=i,j<=n-1)]`):
print(` and its semi-rigorous proof.`):
print(`It checks that the "proof" indeed is correct up to K`):
print(`For example, try:`):
print(`SRproof(1/(m+n+1),m,n,N,30,60);`):


elif nargs=1 and args[1]=SRproofI then
print(`SRproofI(a,m,n,N,Max,K): Given an expression a of m and n`):
print(`outputs the determinant ratios `):
print(` [det(delta(i,j)+a(i,j), 0<=i,j<=n)]/[det(delta(i,j)+ a(i,j)), 0<=i,j<=n-1)]`):
print(`and its semi-rigorous proof.`):
print(`It checks that the "proof" indeed is correct up to K`):
print(`For example, try:`):
print(`SRproofI(binomial(m+n,n),m,n,N,30,60);`):

elif nargs=1 and args[1]=SRproofP then
print(`SRproofP(a,m,n,N,Max,p,S1,S2,K): Given an expression a of m and n`):
print(`that also depends on a parameter p,`):
print(`outputs the determinant ratios`):
print(` [det(a(i,j), 0<=i,j<=m)]/[det(a(i,j), 0<=i,j<=m-1)]`):
print(` and its semi-rigorous proof.`):
print(`It checks that the "proof" indeed is correct up to K`):
print(`For example, try:`):
print(`SRproofP(1/(m+n+p),m,n,N,30,p,1,10,20);`):


elif nargs=1 and args[1]=Zug then
print(`Zug(a,m,n,N,Max): The pair [conj(operator) for ratio, conj. horiz. `):
print(`recurrence for cofactors]`):
print(`For the matrix a(m,n) `):
print(`e.g. Zug(1/(m+n+1),m,n,N,20);`):

elif nargs=1 and args[1]=ZugI then
print(`ZugI(a,m,n,N,Max): The pair [conj(operator) for ratio, conj. horiz. `):
print(`recurrence for cofactors] for the matrix I+a(m,n)`):
print(`e.g. ZugI(binomial(m+n,n),m,n,N,20);`):

elif nargs=1 and args[1]=ZugP then
print(`ZugP(a,m,n,N,Max,p,S1,S2): The pair [conj(operator) for ratio, conj. horiz. `):
print(`recurrence for cofactors]`):
print(`For the matrix a(m,n) `):
print(` where a is closed-form in m and n that also depends on p`):
print(`e.g. Zug(1/(m+n+p),m,n,N,20,p,1,10);`):

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

fi:
end:



###From EKHAD
#yafe just translates from operator notation to everyday notation

solve1:=solve:
yafe:=proc(ope,N,n,SUM)
local gu,i:
gu:=0:
 
for i from 0 to degree(ope,N) do
 gu:=gu+coeff(ope,N,i)*subs(n=n+i,SUM):
od:
 
gu:
end:
 
 
 
yafec:=proc(ope,D,x,INTEGRAND)
local gu,i:
 
 gu:=coeff(ope,D,0)*INTEGRAND:
 
 
for i from 1 to degree(ope,D) do
 gu:=gu+coeff(ope,D,i)*diff(INTEGRAND,x$i):
od:
 
gu:
end:
 
  
simplify1:=proc(bitu,k,a)
local gu,gu1,i,khez,sp:
sp:=1:
gu:=bitu:
 
if type(gu,`*`) then
 
 for i from 1 to nops(gu) do
  gu1:=op(i,gu):
 
 if type(gu1,`^`) and type(op(2,gu1), integer) then
 
     khez:=op(2,gu1):
     gu1:=op(1,gu1):
     sp:=sp*simplify(subs(k=k+a,gu1)/gu1)^khez:
  
   else 
     sp:=sp*simplify(subs(k=k+a,gu1)/gu1):
 
   fi:
 
od:
 
elif type(gu,`^`) and type(op(2,gu), integer) then
   khez:=op(2,gu):
   gu1:=op(1,gu):
   sp:=sp*simplify(subs(k=k+a,gu1)/gu1)^khez:  
 
 
else
 
 sp:=simplify(subs(k=k+a,gu)/gu):
 
fi:
 
sp:
 
end:
 
 
rf:=proc(a,b):
 (a+b-1)!/(a-1)!:
end:
 
 
 
RootOf1:=proc(f,x)
local kv,kvi,i:
kv:=[solve(f=0,x)]:
kvi:={}:
for i from 1 to nops(kv) do
 
  if type(kv[i],integer) and kv[i]>0 then
   kvi:=kvi union {kv[i]}
  fi:
 
od:
 
RETURN(kvi):
 
end:
 
pashet:=proc(p,N)
local i,gu1,gu,p1:
p1:=normal(p):
gu1:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu/gu1):
end:
 
 
 
hafel:=proc(POL,SUMMAND,ope,n,N)
local i,FAC,degg,bi,rat,POL1,SUMMAND1,zub:
degg:=degree(ope,N):
FAC:=0:
 
for i from 0 to degg do
 bi:=coeff(ope,N,i):
 rat:=simplify1(SUMMAND,n,i):
 FAC:=FAC+bi*subs(n=n+i,POL)*rat:
 FAC:=normal(FAC):
od:
 
POL1:=numer(FAC):
zub:=denom(FAC):
SUMMAND1:=SUMMAND/zub:
RETURN(POL1,SUMMAND1,zub):
end:
 
ctold:=proc(SUMMAND,ORDER,k,n,N)
local gu,i,ope,POL1,SUMMAND1,yakhas,P,Q,R,j,res1,kv,hakhi,g,
A1, B1,  P1, SDZ, SDZ1,
degg, eq, fu, gugu, i1, ia1, k1, kvutsa, l1, l2, meka1, mekb1, mumu, 
va1, va2,ope1,denFAC,ope1a:
 
if nargs<>5 then
ERROR(`Syntax: ct(SUMMAND,ORDER,summation_variable,auxiliary_var, Shift_n)`):
fi:
 
 
ope:=0:
 
for i from 0 to ORDER do
 ope:=ope+bpc[i]*N^i:
od: 
gu:=hafel(1,convert(SUMMAND,factorial),ope,n,N):
POL1:=gu[1]:
SUMMAND1:=gu[2]:
denFAC:=gu[3]:
yakhas:=simplify1(SUMMAND1,k,-1):
yakhas:=normal(1/yakhas):
 
print(`yakhas is`,yakhas):
P1:=1:
Q:=numer(yakhas):
R:=denom(yakhas):
res1:=resultant(Q,subs(k=k+j,R),k):
print(`res1 is`,res1):
kv:=RootOf1(res1,j):
 
 
 print(`kv is`,kv):
while nops(kv) > 0 do
 hakhi:=max(op(kv)):
 g:=gcd(Q,subs(k=k+hakhi,R)):
 Q:=normal(Q/g):
 R:=normal(R/subs(k=k-hakhi,g)):
 P1:=P1*product(subs(k=k-i1,g),i1=0..hakhi-1):
 res1:=resultant(Q,subs(k=k+j,R),k):
 kv:=RootOf1(res1,j):
od:
 
P:=POL1*P1:
 
A1:=subs(k=k+1,Q)+R:
A1:=expand(A1):
B1:=subs(k=k+1,Q)-R:
B1:=expand(B1):
l1:=degree(A1,k):
if B1=0 then
l2:=-100:
else
l2:=degree(B1,k):
fi:
meka1:=coeff(A1,k,l1):
mekb1:=coeff(B1,k,l2):
if l1<=l2 
then
k1:=degree(P,k)-l2:
elif l2=l1-1 
then
 mumu:= (-2)*mekb1/meka1:
  if type(mumu,integer) 
  then
   k1:=max(mumu, degree(P,k)-l1+1):
  else
   k1:=degree(P,k)-l1+1:
  fi:
elif l2 < l1-1
then
 k1:= max( 0, degree(P,k)-l1+1 ):
fi:
fu:=0:
 
 
print(`k1 is`,k1):
if k1 < 0 then
RETURN(0):
fi:
if k1 >= 0 then
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*k^ia1:
od:
gugu:=subs(k=k+1,Q)*fu-R*subs(k=k-1,fu)-P:
gugu:=expand(gugu):
degg:=degree(gugu,k):
 
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,k,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2:
va1:=solve1(eq,va1):
kvutsa:={va1}:
fu:=subs(va1,fu):
ope:=subs(va1,ope):
fi:
 
 
 
if ope=0 or kvutsa={} or fu=0 then
RETURN(0):
fi:
gu:={}:
for i1 from 0 to k1 do
gu:=gu union {apc[i1]=1}:
od:
for i1 from 0 to ORDER do
gu:=gu union {bpc[i1]=1}:
od:
fu:=subs(gu,fu):
ope:=subs(gu,ope):
 
ope:=pashet(ope,N):
ope1:=ope*denom(ope):
 
SDZ:=denom(ope)*subs(k=k+1,Q)*fu/P1/denFAC :
SDZ1:=subs(k=k-1,SDZ)*simplify1(convert(SUMMAND,factorial),k,-1):
SDZ1:=simplify(SDZ1):
 
ope1a:=0:
for i from 0 to degree(ope1,N) do
ope1a:=ope1a+sort(coeff(ope1,N,i)*N^i):
od:
 
RETURN(ope1a,SDZ1):
end:
 
 
 
 
ct:=proc(SUMMAND,ORDER,k,n,N)
local gu,i,ope,POL1,SUMMAND1,yakhas,P,Q,R,j,res1,kv,hakhi,g,
A1, B1,  P1, SDZ, SDZ1,
degg, eq, fu, gugu, i1, ia1, k1, kvutsa, l1, l2, meka1, mekb1, mumu, 
va1, va2,ope1,denFAC,ope1a:
 
if nargs<>5 then
ERROR(`Syntax: ct(SUMMAND,ORDER,summation_variable,auxiliary_var, Shift_n)`):
fi:
 
if tovu(convert(SUMMAND,factorial),k,n)=0 then
ERROR(`The summand can be separated into f(`,k,`)g(`,n,`)`): 
fi:
 
ope:=0:
 
for i from 0 to ORDER do
 ope:=ope+bpc[i]*N^i:
od:
 
gu:=hafel(1,convert(SUMMAND,factorial),ope,n,N):
POL1:=gu[1]:
SUMMAND1:=gu[2]:
denFAC:=gu[3]:
yakhas:=simplify1(SUMMAND1,k,-1):
yakhas:=normal(1/yakhas):
 
P1:=1:
Q:=numer(yakhas):
R:=denom(yakhas):
 
 
res1:=findgQR(Q,R,k,100):
 
while res1[2]<>0 do
 g:=res1[1]:
 hakhi:=res1[2]:
 Q:=normal(Q/g):
 R:=normal(R/subs(k=k-hakhi,g)):
 P1:=P1*product(subs(k=k-i1,g),i1=0..hakhi-1):
 res1:=findgQR(Q,R,k,100):
od:
 
P:=POL1*P1:
 
A1:=subs(k=k+1,Q)+R:
A1:=expand(A1):
B1:=subs(k=k+1,Q)-R:
B1:=expand(B1):
l1:=degree(A1,k):
if B1=0 then
l2:=-100:
else
l2:=degree(B1,k):
fi:
meka1:=coeff(A1,k,l1):
mekb1:=coeff(B1,k,l2):
if l1<=l2 
then
k1:=degree(P,k)-l2:
elif l2=l1-1 
then
 mumu:= (-2)*mekb1/meka1:
  if type(mumu,integer) 
  then
   k1:=max(mumu, degree(P,k)-l1+1):
  else
   k1:=degree(P,k)-l1+1:
  fi:
elif l2 < l1-1
then
 k1:= max( 0, degree(P,k)-l1+1 ):
fi:
fu:=0:
 
 
if k1 < 0 then
RETURN(0):
fi:
if k1 >= 0 then
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*k^ia1:
od:
gugu:=subs(k=k+1,Q)*fu-R*subs(k=k-1,fu)-P:
gugu:=expand(gugu):
degg:=degree(gugu,k):
 
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,k,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2:
va1:=solve1(eq,va1):
kvutsa:={va1}:

fu:=subs(va1,fu):
ope:=subs(va1,ope):
fi:
 
 
 
if ope=0 or kvutsa={} or fu=0 then
RETURN(0):
fi:
gu:={}:
for i1 from 0 to k1 do
gu:=gu union {apc[i1]=1}:
od:
for i1 from 0 to ORDER do
gu:=gu union {bpc[i1]=1}:
od:
fu:=subs(gu,fu):
ope:=subs(gu,ope):
 
ope:=pashet(ope,N):
ope1:=ope*denom(ope):
 
SDZ:=denom(ope)*subs(k=k+1,Q)*fu/P1/denFAC :
SDZ1:=subs(k=k-1,SDZ)*simplify1(convert(SUMMAND,factorial),k,-1):
SDZ1:=simplify(SDZ1):
 
ope1a:=0:
for i from 0 to degree(ope1,N) do
ope1a:=ope1a+sort(coeff(ope1,N,i)*N^i):
od:
 
RETURN(ope1a,SDZ1):
end:
 
 
 
cttry:=proc(SUMMAND,ORDER,k,n,resh,N)
local gu,i,ope,POL1,SUMMAND1,yakhas,P,Q,R,j,res1,kv,hakhi,g,
A1, B1,  P1, 
degg, eq, fu, gugu, i1, ia1, k1, kvutsa, l1, l2, meka1, mekb1, mumu, 
va1, va2,eqg:
 
print(`SUMMAND is`, SUMMAND):
if nargs<>6 then
ERROR(`The syntax is cttry(term,ORDER,sum_variable,aux_var,para_list,Shift)`):
fi:
 
if tovu(convert(SUMMAND,factorial),k,n)=0 then
ERROR(`The summand can be separated into f(`,k,`)g(`,n,`)`): 
fi:
 
ope:=0:
 
for i from 0 to ORDER do
 ope:=ope+bpc[i]*N^i:
od:
 
gu:=hafel(1,convert(SUMMAND,factorial),ope,n,N):
POL1:=gu[1]:
SUMMAND1:=gu[2]:
yakhas:=simplify1(SUMMAND1,k,-1):
yakhas:=normal(1/yakhas):
 
P1:=1:
Q:=numer(yakhas):
R:=denom(yakhas):
 
res1:=findgQR(Q,R,k,100):
 
while res1[2]<>0 do
 g:=res1[1]:
 hakhi:=res1[2]:
 Q:=normal(Q/g):
 R:=normal(R/subs(k=k-hakhi,g)):
 P1:=P1*product(subs(k=k-i1,g),i1=0..hakhi-1):
 res1:=findgQR(Q,R,k,100):
od:
 
P:=POL1*P1:
 
A1:=subs(k=k+1,Q)+R:
A1:=expand(A1):
B1:=subs(k=k+1,Q)-R:
B1:=expand(B1):
l1:=degree(A1,k):
if B1=0 then
l2:=-100:
else
l2:=degree(B1,k):
fi:
meka1:=coeff(A1,k,l1):
mekb1:=coeff(B1,k,l2):
if l1<=l2 
then
k1:=degree(P,k)-l2:
elif l2=l1-1 
then
 mumu:= (-2)*mekb1/meka1:
  if type(mumu,integer) 
  then
   k1:=max(mumu, degree(P,k)-l1+1):
  else
   k1:=degree(P,k)-l1+1:
  fi:
elif l2 < l1-1
then
 k1:= max( 0, degree(P,k)-l1+1 ):
fi:
fu:=0:
 
 
if k1 < 0 then
RETURN(0):
fi:
if k1 >= 0 then
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*k^ia1:
od:
gugu:=subs(k=k+1,Q)*fu-R*subs(k=k-1,fu)-P:
gugu:=expand(gugu):
degg:=degree(gugu,k):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,k,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2:
eqg:=subs(n=1/2,eq):
 
for i from 1 to nops(resh) do
 eqg:=subs(op(i,resh)=i^2+3,eqg):
od:
 
va1:=solve1(eq,va1):
kvutsa:={va1}:
fu:=subs(va1,fu):
ope:=subs(va1,ope):
fi:
 
 
 
if ope=0 or kvutsa={} or fu=0 then
RETURN(0):
 else
  RETURN(1):
fi:
 
end:
 
 
paper:=proc(SUMMAND,k,n,N,ope1,SDZ1,NAME,REF)
local SHEM,IDENTITY,RECURRENCE:
 
if degree(ope1,N)=1 then
SHEM:=IDENTITY:
else
SHEM:=RECURRENCE:
fi:
lprint(``):
lprint(`A PROOF OF THE`,NAME,SHEM):
lprint(``):
lprint(`By Shalosh B. Ekhad, Temple University, ekhad@math.temple.edu`):
lprint(``):
lprint(`I will give a short proof of the following result(`,REF,`).`):
lprint(``):
if degree(ope1,N)=1 then
lprint(`(Note that since the recurrence below is first order, this`):
lprint(`means that the sum`, SUM(n), `has closed form,and it is`):
lprint(`easily seen to be equivalent.)`):
lprint(``):
fi:
lprint(`Theorem:Let`, F(n,k), `be given by`):
lprint(``):
print(SUMMAND):
lprint(``):
lprint(`and let`, SUM(n),`be the sum of`, F(n,k),` with`):
lprint(`respect to`, k,`.`):
lprint(``):
lprint(SUM(n),` satisfies the following linear recurrence equation`):
lprint(``):
print(yafe(ope1,N,n,SUM(n))):
print(`=0.`):
lprint(``):
lprint(`PROOF: We cleverly construct`, G(n,k),`:=`):
lprint(``):
print(SDZ1*SUMMAND):
lprint(`with the motive that`):
lprint(``):
print(yafe(ope1,N,n,F(n,k))):
lprint(`=`,G(n,k+1)-G(n,k), `(check!)`):
lprint(``):
lprint(`and the theorem follows upon summing with respect to`, k,`.QED.`):
end:
 
 
 
paper3:=proc(SUMMAND,k,n,N,ope1,SDZ1):
 
lprint(``):
lprint(`A PROOF OF A RECURRENCE`):
lprint(``):
lprint(`By Shalosh B. Ekhad, Temple University, ekhad@math.temple.edu`):
lprint(``):
 
lprint(`Theorem:Let`, F(n,k), `be given by`):
lprint(``):
print(SUMMAND):
lprint(``):
lprint(`and let`, SUM(n),`be the sum of`, F(n,k),` with`):
lprint(`respect to`, k,`.`):
lprint(``):
lprint(SUM(n),` satisfies the following linear recurrence equation`):
lprint(``):
print(yafe(ope1,N,n,SUM(n))):
print(`=0.`):
lprint(``):
lprint(`PROOF: We cleverly construct`, G(n,k),`:=`):
lprint(``):
print(SDZ1*SUMMAND):
lprint(`with the motive that`):
lprint(``):
print(yafe(ope1,N,n,F(n,k))):
lprint(`=`,G(n,k+1)-G(n,k), `(check!)`):
lprint(``):
lprint(`and the theorem follows upon summing with respect to`, k,`.QED.`):
end:
 
 
paperc:=proc(INTEGRAND,y,x,D,ope1,SDZ1):
lprint(``):
lprint(``):
lprint(`A PROOF OF A DIFFERENTIAL EQUATION SATISFIED BY AN INTEGRAL`):
lprint(``):
lprint(`By Shalosh B. Ekhad, Temple University, ekhad@math.temple.edu`):
lprint(``):
lprint(`I will give a short proof of the following result.`):
lprint(``):
 
lprint(`Theorem:Let`, F(x,y), `be given by`):
lprint(``):
print(INTEGRAND):
lprint(``):
lprint(`and let`, INTEGRAL(x),`be the integral of`, F(x,y),` with`):
lprint(`respect to`, y,`.`):
lprint(``):
lprint(INTEGRAL(x),` satisfies the following linear differential equation`):
lprint(``):
print(yafec(ope1,D,x,INTEGRAL(x))):
print(`=0.`):
lprint(``):
lprint(`PROOF: We cleverly construct`, G(x,y),`:=`):
lprint(``):
print(SDZ1*INTEGRAND):
lprint(`with the motive that`):
lprint(``):
print(yafec(ope1,D,x,F(x,y))):
lprint(`=`,diff(G(x,y),y), `(check!)`):
lprint(``):
lprint(`and the theorem follows upon integrating with respect to`, y,`.`):
end:
 
paperd:=proc(INTEGRAND,x,n,N,ope1,SDZ1):
lprint(``):
lprint(`A PROOF OF A LINEAR RECURRENCE SATISFIED BY AN INTEGRAL`):
lprint(``):
lprint(`By Shalosh B. Ekhad, Temple University, ekhad@math.temple.edu`):
lprint(``):
lprint(`I will give a short proof of the following result.`):
lprint(``):
 
lprint(`Theorem:Let`, F(n,x), `be given by`):
lprint(``):
print(INTEGRAND):
lprint(``):
lprint(`and let`, INTEGRAL(n),`be the integral of`, F(n,x),` with`):
lprint(`respect to`, x,`.`):
lprint(``):
lprint(INTEGRAL(n),` satisfies the following linear recurrence equation`):
lprint(``):
print(yafe(ope1,N,n,INTEGRAL(n))):
print(`=0.`):
lprint(``):
lprint(`PROOF: We cleverly construct`, G(n,x),`:=`):
lprint(``):
print(SDZ1*INTEGRAND):
lprint(`with the motive that`):
lprint(``):
print(yafe(ope1,N,n,F(n,x))):
lprint(`=`,diff(G(n,x),x)):
lprint(``):
lprint(`and the theorem follows upon integrating with respect to`, x,`.QED.`):
end:
 
 
paperc:=proc(INTEGRAND,y,x,D,ope1,SDZ1):
lprint(`A PROOF OF A DIFFERENTIAL EQUATION SATISFIED BY AN INTEGRAL`):
lprint(``):
lprint(`By Shalosh B. Ekhad, Temple University, ekhad@math.temple.edu`):
lprint(``):
lprint(`I will give a short proof of the following result.`):
lprint(``):
 
lprint(`Theorem:Let`, F(x,y), `be given by`):
lprint(``):
print(INTEGRAND):
lprint(``):
lprint(`and let`, INTEGRAL(x),`be the integral of`, F(x,y),` with`):
lprint(`respect to`, y,`.`):
lprint(``):
lprint(INTEGRAL(x),` satisfies the following linear differential equation`):
lprint(``):
print(yafec(ope1,D,x,INTEGRAL(x))):
print(`=0.`):
lprint(``):
lprint(`PROOF: We cleverly construct`, G(x,y),`:=`):
lprint(``):
print(SDZ1*INTEGRAND):
lprint(`with the motive that`):
lprint(``):
print(yafec(ope1,D,x,F(x,y))):
lprint(`=`,diff(G(x,y),y), `(check!)`):
lprint(``):
lprint(`and the theorem follows upon integrating with respect to`, y,`.QED.`):
end:
 
 
 
zeil4:=proc(SUMMAND,k,n,N)
local ORDER,MAXORDER,gu,gu1,resh:
 
 
 MAXORDER:=6:
 
  resh:=[]:
 
 
if simplify1(SUMMAND,n,1)=1 then
ERROR(`Summand does not depend on`,n,` hence the sum is identically constant`):
fi:
for ORDER from 0 to MAXORDER do
# gu1:=cttry(SUMMAND,ORDER,k,n,resh,N):
# if gu1=1 then
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 then
     RETURN(gu):
   fi:
 #fi:
od:
 
ERROR(`No recurrence of order<=`,MAXORDER,`was found`):
 
end:
 
 
 
 
zeil5:=proc(SUMMAND,k,n,N,MAXORDER)
local ORDER,gu,gu1,resh:
 
resh:=[]:
 
if simplify1(SUMMAND,n,1)=1 then
ERROR(`Summand does not depend on`,n,` hence the sum is identically constant`):
fi:
for ORDER from 0 to MAXORDER do
 #gu1:=cttry(SUMMAND,ORDER,k,n,resh,N):
 #if gu1=1 then
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 then
     RETURN(gu):
   fi:
 #fi:
od:
 
ERROR(`No recurrence of order<=`,MAXORDER,`was found`):
 
end:
 
 
 
 
zeil6:=proc(SUMMAND,k,n,N,MAXORDER,resh)
local ORDER,gu,gu1:
 
 
if simplify1(SUMMAND,n,1)=1 then
ERROR(`Summand does not depend on`,n,` hence the sum is identically constant`):
fi:
for ORDER from 0 to MAXORDER do
 #gu1:=cttry(SUMMAND,ORDER,k,n,resh,N):
 #if gu1=1 then
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 then
     RETURN(gu):
   fi:
 #fi:
od:
 
ERROR(`No recurrence of order<=`,MAXORDER,`was found`):
 
end:
 
 
zeil:=proc():
 
if nargs=4 then
 zeil4(args):
elif nargs=5 then
 zeil5(args):
elif nargs=6 then
 zeil6(args):
else
  print(`zeil(SUMMAND,k,n,N) or zeil(SUMMAND,k,n,N,MAXORDER) or`):
  ERROR(`zeil(SUMMAND,k,n,N,MAXORDER,param_list)`):
fi:
 
 
 
end:
 
 
 
 
 
zeillim:=proc(SUMMAND,k,n,N,alpha,beta)
local ope,CERT,lu,k1,i,gu,lu1,lu2:
 
gu:=Zeillim(SUMMAND,k,n,N,alpha+1,beta+1):
ope:=gu[1]:
CERT:=gu[2]:
lu:=gu[3]: 
lu1:=subs(k=alpha,SUMMAND)+subs(k=n-beta,SUMMAND):
lu1:=simplify(lu1):
lu1:=normal(lu1):
lu2:=0:
for i from 0 to degree(ope,N) do
   lu2:=lu2+coeff(ope,N,i)*simplify(subs(n=n+i,lu1)):
 od:
 lu2:=normal(lu2):
ope,CERT,normal(expand(normal(simplify(expand(normal(lu+lu2)))))):
 
end:
 
Zeillim:=proc(SUMMAND,k,n,N,alpha,beta)
local ope,CERT,lu,k1,i,gu:
 
gu:=zeil(SUMMAND,k,n,N):
ope:=gu[1]:
CERT:=gu[2]:
 
lu:=simplify(subs(k=n-beta+1,CERT)*subs(k=n-beta+1,SUMMAND))
-simplify(subs(k=alpha,CERT)*subs(k=alpha,SUMMAND)):
 
for i from 0 to degree(ope,N) do
 for k1 from 1 to i do
   lu:=lu+coeff(ope,N,i)*simplify(subs({n=n+i,k=n-beta+k1},SUMMAND)):
 od:
od:
 
gu,normal(expand(lu)):
 
end:
 
 
 
zeilpap3:=proc(SUMMAND,k,n)
local  SDZ1, gu, ope1,N:
 
gu:=zeil4(SUMMAND,k,n,N):
ope1:=gu[1]:
SDZ1:=gu[2]:
paper3(SUMMAND,k,n,N,ope1,SDZ1):
end:
 
 
 
 
zeilpap5:=proc(SUMMAND,k,n,NAME,REF)
local SDZ1, gu, ope1,N:
 
gu:=zeil4(SUMMAND,k,n,N):
ope1:=gu[1]:
SDZ1:=gu[2]:
paper(SUMMAND,k,n,N,ope1,SDZ1,NAME,REF):
end:
 
 zeilpap:=proc()
 
if nargs=5 then
zeilpap5(args):
elif nargs=3 then
zeilpap3(args):
else
 ERROR(`zeilpap(SUMMAND,k,n,NAME,REF) or  zeilpap(SUMMAND,k,n)`):
fi:
end:
 
 
 
 
AZpapc:=proc(INTEGRAND,y,x)
local D,SDZ1,gu,ope1:
gu:=AZc(INTEGRAND,y,x,D):
ope1:=gu[1]:
SDZ1:=gu[2]:
paperc(INTEGRAND,y,x,D,ope1,SDZ1):
end:
 
 
 
 
AZpapd:=proc(INTEGRAND,x,n)
local D,SDZ1, gu, ope1:
gu:=AZd(INTEGRAND,x,n,D):
ope1:=gu[1]:
SDZ1:=gu[2]:
paperd(INTEGRAND,x,n,D,ope1,SDZ1):
end:
 
 
 
 
 
goremd:=proc(N,ORDER)
local i, gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*N^i:
od:
RETURN(gu):
end:
 
goremc:=proc(D,ORDER)
local i,gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*D^i:
od:
RETURN(gu):
end:
 
gzor:=proc(f,x,n)
local i,gu:
gu:=f:
for i from 1 to n do
gu:=diff(gu,x):
od:
RETURN(gu):
end:
 
gzor1:=proc(a,D,x)
local i,gu:
gu:=0:
for i from 0 to degree(a,D) do
gu:=gu+diff(coeff(a,D,i),x)*D^i+coeff(a,D,i)*D^(i+1):
od:
end:
 
pashetd:=proc(p,N)
local gu,p1,mekh,i:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu,mekh):
end:
 
pashetc:=proc(p,D)
local gu,p1,i,mekh:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,D) do
gu:=gu+factor(coeff(p1,D,i))*D^i:
od:
RETURN(gu,mekh):
end:
 
 
AZd:= proc(A,y,n,N)
local ORDER,gu,KAMA:
KAMA:=20:
 
for ORDER from 0 to KAMA do
gu:=duisd(A,ORDER,y,n,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
0:
end:
 
 
AZc:=proc(A,y,x,D)
local ORDER,gu,KAMA:
KAMA:=20:
 
for ORDER from 0 to KAMA do
gu:=duisc(A,ORDER,y,x,D):
if gu<>0 then
 RETURN(gu):
fi:
od:
0:
end:
 
 
duisd:= proc(A,ORDER,y,n,N)
local gorem, conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,shad,KAMA,i1:
KAMA:=10:
gorem:=goremd(N,ORDER):
 
 
conj:=gorem:
yakhas:=0:
 
for i1 from 0 to degree(conj,N) do
 
yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve1(eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):
 
 
shad:=pashetd(gorem,N):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:
 
RETURN(shad[1],S):
end:
 
 
duisc:= proc(A,ORDER,y,x,D)
local KAMA,gorem,conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,i,shad:
 
KAMA:=10:
gorem:=goremc(D,ORDER):
 
conj:=gorem:
yakhas:=0:
 
for i from 0 to degree(conj,D) do
yakhas:=yakhas+normal(coeff(conj,D,i)*gzor(A,x,i)/A):
yakhas:=normal(yakhas):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve1(eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):
shad:=pashetc(gorem,D):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:
shad[1],S:
end:
 
 
 
bdokcertc:=proc(A,y,x,D,ope,S)
local gu,i:
gu:=0:
 
for i from 0 to degree(ope,D) do
gu:=gu+coeff(ope,D,i)*simplify(gzor(A,x,i)/A):
gu:=normal(gu):
od:
 
gu:=gu/simplify(diff(S*A,y)/A):
normal(gu);
end:
 
 
bdokcertd:=proc(A,y,n,N,ope,S)
local gu,i:
gu:=0:
 
for i from 0 to degree(ope,N) do
gu:=gu+coeff(ope,N,i)*simplify( subs(n=n+i,A)/A):
gu:=normal(gu):
od:
 
gu:=gu/simplify(diff(S*A,y)/A):
normal(gu);
end:
 
 
 
bdokcto:=proc(SUMMAND1,ORDER,k,n,N)
local mu,gu,i,G1,ope,lu,SUMMAND:
 
SUMMAND:=convert(SUMMAND1,factorial):
mu:=ct(SUMMAND,ORDER,k,n,N):
 
if mu=0 then
 RETURN(0):
fi:
 
ope:=mu[1]:
G1:=mu[2]*SUMMAND:
gu:=0:
 
for i from 0 to degree(ope,N) do
gu:=gu+coeff(ope,N,i)*simplify( subs(n=n+i,SUMMAND)/SUMMAND):
gu:=normal(gu):
od:
 
lu:=simplify(subs(k=k+1,G1)/SUMMAND)-mu[2]:
lu:=normal(lu):
normal(gu/lu);
end:
 
 
 
 
bdokct:=proc(SUMMAND1,ORDER,k,n,N)
local mu,gu,i,G1,ope,lu,SUMMAND:
 
SUMMAND:=convert(SUMMAND1,factorial):
mu:=ct(SUMMAND,ORDER,k,n,N):
 
if mu=0 then
 RETURN(0):
fi:
 
ope:=mu[1]:
G1:=mu[2]*SUMMAND:
gu:=0:
 
for i from 0 to degree(ope,N) do
gu:=gu+coeff(ope,N,i)*simplify1(SUMMAND,n,i):
gu:=normal(gu):
od:
 
lu:=subs(k=k+1,mu[2])*simplify1(SUMMAND,k,1)-mu[2]:
lu:=normal(lu):
normal(gu/lu);
end:
 
 
        
celine:=proc(f,ii,jj)
local m,j,l,i,tt,yy,zz,rr,xx,ss,zy,
 b, iii, k, n, r, zzz:
 
m:='m':j:='j':l:='l':i:='i':
iii:=ii*(jj+1)+jj:
xx:=seq(b[i],i=0..iii):i:='i':
ss:=numer(simplify(expand(sum(sum(b[i*(jj+1)+j]*f(n-j,k-i)/f(n,k), 
            i=0..ii),j=0..jj)))):
tt:=coeffs(collect(ss,k),k):
yy:=solve1({tt},{xx}):
zz:=simplify(sum(sum(b[l*(jj+1)+m]*F(n-m,k-l),l=0..ii),m=0..jj)): 
i:='i':j:='j':lprint(` ` ):lprint(` `):lprint(`The full  recurrence is`):
lprint(` `):
zz:=numer(simplify(subs(yy,zz))):
for i from 0 to ii do for j from 0 to jj do                
zz:=collect(zz,F(n-j,k-i)) od od:
lprint(zz,`==0`);zzz:=zz:
l:='l':j:='j':m:='m':
r:=subs(yy,simplify(expand(sum(sum(sum(b[l*(jj+1)+m],l=j+1..ii)*simplify(expand(f(n-m,k-j)/f(n,k))),j=0..ii),m=0..jj)))): 
 
rr:=simplify(factor(simplify(r))):lprint(` ` ):lprint(` `):
lprint(`The telescoped form is`);lprint(` `);
zy:=factor(subs(yy,sum(simplify(sum(b[l*(jj+1)+m],l=0..ii))*F(n-m,k),m=0..jj))):
lprint(zy,`==G(n,k)-G(n,k-1)`); 
 
lprint(` `);
lprint(`where G(n,k)=R(n,k)*F(n,k) and the rational function  R(n,k) is `):
lprint(` `);
rr;
end:
 
findgQR:=proc(Q,R,k,L)
local j,g:
 
for j from 1 to L do
 g:=gcd(Q,subs(k=k+j,R)):
 if degree(g,k)>0 then
  RETURN(g,j):
 fi:
od:
1,0:
end:
 
  
 
 
tovu:=proc(SU,k,n)
local gu:
 
gu:=simplify1(simplify1(SU,n,1),k,1):
 
if gu=1 then
RETURN(0):
else
RETURN(1):
fi:
end:
 
 
 
 
 
###End From EKHAD

####################Begin Findrec




#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f1,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a,K,f,i1:
option remember:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f1) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:

K:=lcm(seq(denom(f1[i1]),i1=1..nops(f1))):
f:=[seq(f1[i1]*K,i1=1..nops(f1))]:



if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


#print(`There is hope for  a recurrence of order`, ORDER, `and degree`, DEGREE):

ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=Linear(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

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

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(-1):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 

#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f1,n,N,MaxC)
local DEGREE, ORDER,ope,L,K,f,ope1,DEGREE1,ORDER1,i1:
K:=lcm(seq(denom(f1[i1]),i1=1..nops(f1))):

if K=0 then
print(`Contains 0 s`):
 RETURN(FAIL):
fi:

f:=[seq(f1[i1]*K,i1=1..nops(f1))]:


for L from 0 to MaxC do
for ORDER from 1 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):

if ope<>FAIL and ope<>-1 then
  RETURN(ope):
fi:
 if ope=-1 then
       ORDER1:=ORDER-1:
 
   for i1 from 2 while (2+DEGREE+i1)*(ORDER)+4<=nops(f) do

      DEGREE1:=DEGREE+i1:
      ope1:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE1,ORDER1,n,N):
       if ope1<>FAIL then
             RETURN(ope1):
       fi:
    od:

 fi:
 od:
od:
FAIL:

end:


#findrecK(f,DEGREE,ORDER): Is there a unique recurrence if degree DEGREE and order ORDER
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrecK([seq(i,i=1..10)],0,2);
findrecK:=proc(f,DEGREE,ORDER)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,n,N,a:
option remember:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
print(`Insufficient data for a recurrence of order`, ORDER, ` degree `, DEGREE):
RETURN(false):
fi:


ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=Linear(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(false):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)=1 then
 RETURN(true):
else
 RETURN(false):
fi:

end:

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:



#findrecR(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of rational numbers degree DEGREE and order ORDER
#For example, try: findrecR([seq(i,i=1..10)],0,2,n,N);
findrecR:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:



if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=Linear(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

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

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 

#FindrecR(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try FindrecR([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
FindrecR:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 1 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrecR([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:


#Findrec1(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs. of FIRST ORDER
#of maximum DEGREE+1<=MaxC
#e.g. try Findrec1([1$20],n,N,1);
Findrec1:=proc(f1,n,N,MaxC)
local DEGREE, ORDER,ope,K,f,ope1,DEGREE1,ORDER1,i1:
K:=lcm(seq(denom(f1[i1]),i1=1..nops(f1))):

if K=0 then
print(`Contains 0 s`):
 RETURN(FAIL):
fi:

f:=[seq(f1[i1]*K,i1=1..nops(f1))]:

for ORDER from 1 to 1 do
for DEGREE from 1 to MaxC do
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):

if ope<>FAIL and ope<>-1 then
  RETURN(ope):
fi:
 if ope=-1 then
       ORDER1:=ORDER-1:
 
   for i1 from 2 while (2+DEGREE+i1)*(ORDER)+4<=nops(f) do

      DEGREE1:=DEGREE+i1:
      ope1:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE1,ORDER1,n,N):
       if ope1<>FAIL then
             RETURN(ope1):
       fi:
    od:

 fi:
 od:
od:

FAIL:

end:

##########end Findrec


######Begin One variable guessing of rational functions



#GenRat1(x,dxt,dxb,a): a generic (monic top)
#rational function of degree of top dxt 
#and degree of bottom dxb
#with coeffs.parametrized by a[i] followed by the set
#of coeffs.
GenRat1:=proc(x,dxt,dxb,a) local i:

(add(a[i]*x^i,i=0..dxt))/
add(a[i+dxt+1]*x^i,i=0..dxb)
,{seq(a[i],i=0..dxt+dxb+1)}:
end:




#GR1PolRat(S1,N,x,M): Guesses a poly in N rational in x for
#the data set S1 of the form {[i,POL(N,Rat(M))]}:
GR1PolRat:=proc(S1,N,x,M)
local gu,S1a,Deg,lu,i,i1:
option remember:
if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:

Deg:=max(seq(degree(S1[i1][2],N),i1=1..nops(S1))):
gu:=0:
for i from 0 to Deg do
S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,i)],i1=1..nops(S1))}:
 lu:=GR1Rat(S1a,M,x):
  if lu=FAIL then
   RETURN(FAIL):
  else
  gu:=gu+factor(normal(lu))*N^i:
  fi:
 od:
gu:
end:


#GR1Rat(S1,N,x): Guesses a rational function in N rational in x for
#the data set S1 of the form {[i,RAT(N)]}:
GR1Rat:=proc(S1,N,x)
local Deg1,Deg2,i,i1,S1Top,S1Bot,R1,c,Top,Bot:

if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:

Deg1:=max(seq(degree(numer(S1[i1][2]),N),i1=1..nops(S1))):

Deg2:=max(seq(degree(denom(S1[i1][2]),N),i1=1..nops(S1))):

S1Top:={}:
S1Bot:={}:

for i from 1 to nops(S1) do
R1:=S1[i][2]:
 if degree(numer(R1),N)=Deg1 and degree(denom(R1),N)=Deg2  then
 c:=coeff(numer(R1),N,Deg1):
 S1Top:=S1Top union {[S1[i][1],numer(R1)/c]}:
 S1Bot:=S1Bot union {[S1[i][1],denom(R1)/c]}:
 fi:
od:
Top:=GR1Pol(S1Top,N,x):
Bot:=GR1Pol(S1Bot,N,x):

 if Top<>FAIL and Bot<>FAIL then
  RETURN(Top/Bot):
 else
  RETURN(FAIL):
 fi:

end:






#GR1d1d2(S1,x,d1,d2): guesses a rational function in x of degree d1/d2
#that fits the data in DS1
GR1d1d2:=proc(S1,x,d1,d2)
local eq,var,P,a,var1,i:
P:=GenRat1(x,d1,d2,a):
var:=P[2]:
P:=P[1]:

if nops(S1)<=d1+d2+3 then
# print(`Insufficient data`):
 RETURN(FAIL):
fi:

eq:={seq(numer(normal(subs(x=S1[i][1],P)-S1[i][2])),i=1..d1+d2+3)}:

var1:=Linear(eq,var):
if expand(subs(var1,numer(P)))=0 then
  RETURN(FAIL):
fi:
P:=subs(var1,P):
if Bdok1(S1,P,x) then 
   RETURN(factor(normal(P))):
 else
 RETURN(FAIL):
fi:

end:

#GR1(S1,x): guesses a polynomial from the data-set S1
GR1:=proc(S1,x) local d1,d2,gu,L,i1:
if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:

for L from 0 to nops(S1)-4 do
 for d1 from 0 to L do
   d2:=L-d1:
 gu:=GR1d1d2(S1,x,d1,d2):
  if gu<>FAIL then
    RETURN(gu):
  fi:
 od:
od:

FAIL:
end:

#GR1start(S1,x,D1,D2): guesses a polynomial from the data-set S1
#starting with degree D1/D2
GR1start:=proc(S1,x,D1,D2) local d1,d2,gu,i1:
if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:
 for d1 from D1 to nops(S1)-3-D2 do
  for d2 from D2 to nops(S1)-3-d1 do
 gu:=GR1d1d2(S1,x,d1,d2):
  if gu<>FAIL then
    RETURN(gu):
  fi:
 od:
od:

FAIL:
end:





#GR1Pol(S1,N,x): Guesses a poly in N rational in x for
#the data set S1 of the form {[i,POL(N)]}:
GR1Pol:=proc(S1,N,x)
local gu,S1a,Deg,lu,i,i1,d1,d2,lDeg:
if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:

Deg:=max(seq(degree(S1[i1][2],N),i1=1..nops(S1))):
lDeg:=min(seq(ldegree(S1[i1][2],N),i1=1..nops(S1))):
gu:=0:

S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,lDeg)],i1=1..nops(S1))}:
 lu:=GR1(S1a,x):
  if lu=FAIL  then
   RETURN(FAIL):
  else
  gu:=gu+factor(normal(lu)*N^lDeg):
  fi:


d1:=degree(numer(lu),x):
d2:=degree(denom(lu),x):
for i from lDeg+1 to Deg do
S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,i)],i1=1..nops(S1))}:
 lu:=GR1start(S1a,x,d1,d2):
  if lu=FAIL then
   RETURN(FAIL):
  else
  gu:=gu+factor(normal(lu))*N^i:
  fi:
 od:

gu:
end:




#Bdok1(DS1,P,x): checks that P in the single variable x  indeed fits the data set DS1
Bdok1:=proc(DaS,P,x) local i:

for i from 1 to nops(DaS) do
if subs(x=DaS[i][1],denom(P))<>0 and
normal(subs(x=DaS[i][1],P)-DaS[i][2])<>0 then

  RETURN(false):
fi:
od:
true:
end:

#############End guessing of rational functions of oner variable

Kr:=proc(i,j): if i=j then 1 else 0:fi:end:

#DetRat(a,m,n,m0): det(a(i,j), 0<=i,j<=n0)/det(a(i,j), 0<=i,j<=n0-1)
#For example, try:DetRat(1/(m+n+1),10);
DetRat:=proc(a,m,n,m0) local gu,j:
gu:=Dave(a,m,n,m0):
if gu=FAIL then
 RETURN(FAIL):
else
add(gu[j+1]*eval(subs({m=m0,n=j},a)),j=0..m0):
fi:
end:

#DetRatI(a,m,n,m0): det(b(i,j), 0<=i,j<=n0)/det(b(i,j), 0<=i,j<=n0-1)
#where (b(i,j)=delta(i,j)+a(i,j))
#For example, try:DetRatI(-binomial(20+m+n,n),10);
DetRatI:=proc(a,m,n,m0) local gu,j:
gu:=DaveI(a,m,n,m0):
if gu=FAIL then
 RETURN(FAIL):
else
add(gu[j+1]*eval(subs({m=m0,n=j},a)),j=0..m0)+
gu[m0+1]:
fi:
end:

#DetRatSeq(a,m,n,n0): the sequence of DetRat for i from 1 to n0
DetRatSeq:=proc(a,m,n,n0) local i,lu,lu1:
lu:=[]:
for i from 1 to n0 do
lu1:=DetRat(a,m,n,i):
if lu1=FAIL then
 RETURN(FAIL):
else
 lu:=[op(lu),lu1]:
fi:
od:
end:

#DetRatSeqI(a,m,n,n0): the sequence of DetRatI for i from 1 to n0
DetRatSeqI:=proc(a,m,n,n0) local i,lu,lu1:
lu:=[]:
for i from 1 to n0 do
lu1:=DetRatI(a,m,n,i):
if lu1=FAIL then
 RETURN(FAIL):
else
 lu:=[op(lu),lu1]:
fi:
od:
end:


#Dave(a,m,n,m0): Given an expression  a of 
#m and n outputs the cofactors of the last row
#(normalized with the first being 1) of the
#(n+1) by (n+1) matrix a(i,j) (0<=i,j<=n)
#For example, try:
#Dave(1/(m+n+1),m,n,10);
Dave:=proc(a,m,n,n0) local i1,j1,x,var,var1,eq,ind,gu:
option remember:
eq:=eval({seq(add(subs({m=i1,n=j1},a)*x[j1],j1=0..n0),i1=0..n0-1)}):

var:={seq(x[j1],j1=0..n0)}:


var1:=Linear(eq,var):

ind:={}:
for i1 from 1 to nops(var1) do
if op(1,op(i1,var1))=op(2,op(i1,var1)) then
 ind:=ind union {op(1,op(i1,var1))}:
fi:
od:

if nops(ind)<>1 then
 RETURN(FAIL):
fi:

gu:=subs(ind=1,[seq(subs(var1,x[j1]),j1=0..n0)]):
if gu[nops(gu)]<>0 then
 RETURN([seq(gu[i1]/gu[nops(gu)],i1=1..nops(gu))]):
else
RETURN(FAIL):
fi:

end:





Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:



DetRatRecS:=proc(a,m,n,N,Max) local lu:
lu:=DetRatSeq(a,m,n,trunc((Max+5)*(Max+1)/4)+7):
if lu=FAIL then
 RETURN(FAIL):
fi:
Findrec(lu,n,N,Max):
end:
 

DetRatRecSI:=proc(a,m,n,N,Max) local lu:
lu:=DetRatSeqI(a,m,n,trunc((Max+5)*(Max+1)/4)+7):
if lu=FAIL then
 RETURN(FAIL):
fi:
Findrec(lu,n,N,Max):
end:
 

#DetRatRec(a,m,n,N,Max) : Conjectures a recurrence for 
#det(a(i,j), 0<=i,j<=n)/det(a(i,j), 0<=i,j<=n-1) 
#For example, try:
#DetRatRec(1/(m+n+1),m,n,N,10);
DetRatRec:=proc(a,m,n,N,Max) local lu,i:
for i from 1 to Max do
 lu:=DetRatRecS(a,m,n,N,i) :
if lu<>FAIL then
  RETURN(lu):
fi:
od:
FAIL:
end:

#DetRatRecI(a,m,n,N,Max) : Conjectures a recurrence for 
#det(a(i,j)+delta(i,j), 0<=i,j<=n)/det(a(i,j), 0<=i,j<=n-1) 
#For example, try:
#DetRatRecI(binomial(m+n,m),m,n,N,10);
DetRatRecI:=proc(a,m,n,N,Max) local lu,i:
for i from 1 to Max do
 lu:=DetRatRecSI(a,m,n,N,i) :
if lu<>FAIL then
  RETURN(lu):
fi:
od:
FAIL:
end:

#DetRatRec1(a,m,n,N,Max) : Conjectures a recurrence for 
#det(a(i,j), 0<=i,j<=n)/det(a(i,j), 0<=i,j<=n-1) 
#For example, try:
#DetRatRec1(1/(m+n+1),m,n,N,10);
DetRatRec1:=proc(a,m,n,N,Max) local lu:
lu:=DetRatSeq(a,m,n,(Max+1)*2+6):
if lu=FAIL then
 RETURN(FAIL):
fi:
Findrec1(lu,n,N,Max):
end:
 


##Trying a new DetRatRecP
DetRatRecP:=proc(a,m,n,N,Max,p,S1,S2) local i1,gu:

gu:={seq([i1,DetRatRec(subs(p=i1,a),m,n,N,Max)],i1=S1..S2)};
GR1Pol(gu,N,p):
end:


##End Trying a new DetRatRecP
#DetRatRecP(a,n,N,Max,p,S1,S2): Like DetRatRec (q.v.) but the function
#a(i,j) is now an expression also depends on a parameter p, and it tries by interpolating
#from DetRatRec with p=S1..S2
#For example, try:
#DetRatRecP(1/(m+n+p),n,N,10,p,1,20);

DetRatRecPOld:=proc(a,m,n,N,Max,p,S1,S2) local i1,gu,gu1,gu2,gu3,Max1,lu:

gu1:=DetRatRec(subs(p=S1,a),m,n,N,Max);
gu2:=DetRatRec(subs(p=S1+1,a),m,n,N,Max);
gu3:=DetRatRec(subs(p=S1+2,a),m,n,N,Max);
Max1:=max(
degree(numer(gu1),n)+degree(numer(gu1),N),
degree(numer(gu2),n)+degree(numer(gu2),N),
degree(numer(gu3),n)+degree(numer(gu3),N)):

gu:=[[S1,gu1],[S1+1,gu2],[S1+2,gu3]]:

for i1 from S1+3 to min(S2,S1+8) do
gu1:=DetRatRecS(subs(p=i1,a),m,n,N,Max1);

 if gu1=FAIL then
  RETURN(FAIL):
 else
  gu:=[op(gu),[i1,gu1]]:
 fi:
od:


lu:=GR1Pol(gu,N,p):

if lu<>FAIL then 
 RETURN(lu):
fi:


for i1 from S1+8 to S2 do
gu1:=DetRatRecS(subs(p=i1,a),m,n,N,Max1);

 if gu1=FAIL then
  RETURN(FAIL):
 else
  gu:=[op(gu),[i1,gu1]]:
 fi:

lu:=GR1Pol(gu,N,p):

if lu<>FAIL
 then RETURN(lu):
fi:

od:

FAIL: 
end:

#DetRatRecPI(a,n,N,Max,p,S1,S2): Like DetRatRecI (q.v.) but the function
#a(i,j) is now an expression also depends on a parameter p, and it tries by interpolating
#from DetRatRec with p=S1..S2
#For example, try:
#DetRatRecPI(1/(m+n+p),n,N,10,p,1,20);

DetRatRecPI:=proc(a,m,n,N,Max,p,S1,S2) local i1,gu,gu1,gu2,gu3,lu:

gu1:=DetRatRecI(subs(p=S1,a),m,n,N,Max);
gu2:=DetRatRecI(subs(p=S1+1,a),m,n,N,Max);
gu3:=DetRatRecI(subs(p=S1+2,a),m,n,N,Max);

gu:=[[S1,gu1],[S1+1,gu2],[S1+2,gu3]]:

for i1 from S1+3 to min(S2,S1+8) do
gu1:=DetRatRecSI(subs(p=i1,a),m,n,N,Max);

 if gu1=FAIL then
  RETURN(FAIL):
 else
  gu:=[op(gu),[i1,gu1]]:
 fi:
od:


lu:=GR1Pol(gu,N,p):

if lu<>FAIL then 
 RETURN(lu):
fi:


for i1 from S1+8 to S2 do
gu1:=DetRatRecSI(subs(p=i1,a),m,n,N,Max);

 if gu1=FAIL then
  RETURN(FAIL):
 else
  gu:=[op(gu),[i1,gu1]]:
 fi:

lu:=GR1Pol(gu,N,p):

if lu<>FAIL
 then RETURN(lu):
fi:

od:

FAIL: 
end:

 
#DaveI(a,m,n,m0): Given an expression  a of 
#m and n outputs the cofactors of the last row
#(normalized with the first being 1) of the
#(n+1) by (n+1) matrix (delta(i,j)+a(i,j)) (0<=i,j<=n)
#For example, try:
#DaveI(1/(m+n+1),m,n,10);
DaveI:=proc(a,m,n,n0) local i1,j1,x,var,var1,eq,ind,gu:
option remember:
eq:=eval({seq(
add(subs({m=i1,n=j1},a)*x[j1],j1=0..i1-1)+
add((1+subs({m=i1,n=j1},a))*x[j1],j1=i1..i1)+
add(subs({m=i1,n=j1},a)*x[j1],j1=i1+1..n0),
i1=0..n0-1)}):

var:={seq(x[j1],j1=0..n0)}:


var1:=Linear(eq,var):

ind:={}:
for i1 from 1 to nops(var1) do
if op(1,op(i1,var1))=op(2,op(i1,var1)) then
 ind:=ind union {op(1,op(i1,var1))}:
fi:
od:

if nops(ind)<>1 then
 RETURN(FAIL):
fi:

gu:=subs(ind=1,[seq(subs(var1,x[j1]),j1=0..n0)]):
if gu[nops(gu)]<>0 then
 RETURN([seq(gu[i1]/gu[nops(gu)],i1=1..nops(gu))]):
else
RETURN(FAIL):
fi:

end:



#DaveH1(a,m,n,M,DEG,ORD,Max): Given an expression of two non-negative
#arguments guesses a recurrence operator OPE(n,N) of degree DEG in n
#and order (i.e. degree in N) ORD 
#annihilating A(m,n):= the cofactor of the (n,m) entry 
#of the (n+1) by (n+1) matrix a(i,j), 0<=i,j<=n
#Here N is the shift in n
#Max is the farthest you are willing to explore
#For example, try:
#DaveH1(1/(m+n+1),m,n,M,3,1,10);
DaveH1:=proc(a,m,n,M,DEG,ORD,Max) local gu,gu0,ope,m0,
Ope,m1,m2,m3:

gu:={}:

m1:=(1+DEG)*(1+ORD)+5+ORD+1:

for m0 from m1 to m1+10 do

gu0:=Dave(a,m,n,m0):


gu0:=[op(2..nops(gu0),gu0)]:

ope:=findrec(gu0,DEG,ORD,m,M):
  if ope<>FAIL and ope<>-1 then
   gu:=gu union {[m0,ope]}:
  fi:
 od:


#RETURN(gu):
if nops(gu)<5 then

 RETURN(gu,FAIL):
fi:

Ope:=GR1PolRat(gu,M,n,m):


if Ope<>FAIL then

      m2:=m0+1:
         for m3 from m2 to m2+5 do

          if DaveFromOpe(Ope,m,n,M,m3)<>FAIL and Dave(a,m,n,m3)<>DaveFromOpe(Ope,n,m,M,m3) then
                 print(Ope, ` is not quite right for`, m3):
             fi:         
        od:

RETURN(Ope):

fi:


for m0 from m1+11 to m1+Max do
gu0:=Dave(a,m,n,m0):
gu0:=[op(2..nops(gu0),gu0)]:
ope:=findrec(gu0,DEG,ORD,m,M):
  if ope<>FAIL and ope<>-1 then
   gu:=gu union {[m0,ope]}:
    Ope:=GR1PolRat(gu,M,n,m):
     if Ope<>FAIL then

      m2:=m0+1:
         for m3 from m2 to m2+5 do
          if Dave(a,m,n,m3)<>DaveFromOpe(Ope,m,n,M,m3) then
                  print(Ope, `failed for`, m3):
                 RETURN(FAIL):
          fi:         
         od:
      RETURN(Ope):

     fi:
    fi:
 od:

FAIL:
end:





#FindDegOrdDaveH(a,m,n,Max): Given a function a of two non-negative
#arguments finds the Degree and Order
#annihilating A(m,n):= the cofactor of the (n,m) entry 
#of the (n+1) by (n+1) matrix a(i,j), 0<=i,j<=n
#Max (>20) is the largest row you are willing to try
#Here N is the shift in n
#For example, try:
#FindDegOrdDaveH(1/(m+n+1),m,n,20);
FindDegOrdDaveH:=proc(a,m,n,Max) local ope,m0,gu0,gu0A,gu0B,gu0C,L,DEG,ORD,opeA,
opeB,opeC,M:

ope:=FAIL:
for m0 from 20 to Max do
gu0:=Dave(a,m,n,m0):
if gu0=FAIL then
  RETURN(FAIL):
fi:
gu0:=[op(2..nops(gu0),gu0)]:
L:=2*trunc(sqrt(nops(gu0)-4))-4:
ope:=Findrec(gu0,m,M,L):

  if ope<>FAIL then
    DEG:=degree(numer(ope),m):
    ORD:=degree(numer(ope),M):
    if findrec(gu0,DEG,ORD,m,M)<>FAIL  and  findrecK(gu0,DEG,ORD) then
      gu0A:=Dave(a,m,n,m0+1):
      gu0A:=[op(2..nops(gu0A),gu0A)]:
      opeA:=findrecK(gu0A,DEG,ORD):

      gu0B:=Dave(a,m,n,m0+2):
      gu0B:=[op(2..nops(gu0B),gu0B)]:
      opeB:=findrecK(gu0B,DEG,ORD):

      gu0C:=Dave(a,m,n,m0+3):
      gu0C:=[op(2..nops(gu0B),gu0C)]:
      opeC:=findrecK(gu0C,DEG,ORD):

         if opeA and opeB and opeC then
          RETURN(DEG,ORD):
         fi:
    fi:
 fi:
od:

FAIL:
end:






#DaveH(a,m,n,N,Max): Given an expression a (of  m and n)
#guesses a recurrence operator P(m,n,M) 
#annihilating A(m,n):= the cofactor of the (m,n) entry 
#of the (n+1) by (n+1) matrix a(i,j), 0<=i,j<=n
#Max (>20) is the largest row you are willing to try
#Here N is the shift in n
#For example, try:
#DaveH(1/(m+n+1),m,n,N,20);
DaveH:=proc(a,m,n,N,Max) local gu,DEG,ORD:
gu:=FindDegOrdDaveH(a,m,n,Max):

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

DEG:=gu[1]: ORD:=gu[2]:
subs({m=n,n=m},DaveH1(a,m,n,N,DEG,ORD,Max)):
end:


#DaveFromOpe(Ope,m,n,N,m0): Finds Dave(a,n0) from Ope(m,n,M) 
#For example try: 
#Ope:=DaveH(1/(m+n+1),m,n,N,30):
#DaveFromOpe(Ope,m,n,N,10);
DaveFromOpe:=proc(Ope1,m,n,N,m0)
local T,i,i1,Ope,ku:
Ope:=numer(Ope1):

T[m0]:=1:

for i from m0+1 to m0+degree(Ope,N) do
 T[i]:=0:
od:

for i from m0-1 by -1 to 0 do
ku:= subs({n=i,m=m0},coeff(Ope,N,0)):
 if ku=0 then
  RETURN(FAIL):
 fi:

 T[i]:= -add(subs({n=i,m=m0},coeff(Ope,N,i1))*T[i+i1],i1=1..degree(Ope,N))/
ku
:
od:

[seq(T[i1],i1=0..m0)]:
end:


#DaveHP(a,m,n,N,Max,p,S1,S2): Like DaveH (q.v.) but 
#a is now an expression also depends also on a parameter p, and it tries by interpolating
#from DetRatRec with p=S1..S2
#For example, try:
#DaveHP(1/(m+n+p),m,n,N,10,p,1,20);
DaveHP:=proc(a,m,n,N,Max,p,S1,S2) local i1,gu,gu1,gu2,gu3,lu,DEG,ORD:

if S2-S1<0 then
 print(`S2-S1 should be >=0`):
 RETURN(FAIL):
fi:

gu1:=[FindDegOrdDaveH(subs(p=S1,a),m,n,Max)]:
gu2:=[FindDegOrdDaveH(subs(p=S1+1,a),m,n,Max)]:
gu3:=[FindDegOrdDaveH(subs(p=S1+2,a),m,n,Max)]:

if gu1=[FAIL] or gu2=[FAIL] or gu3=[FAIL] then
  RETURN(FAIL):
fi:

 if nops({gu1,gu2,gu3})<>1 then
 print(`Can't agree on order and/or degree`):
 RETURN(FAIL):
fi:


DEG:=gu1[1]: ORD:=gu1[2]:
gu:=[]:

for i1 from S1 to min(S2,S1+8) do
gu1:=subs({m=n,n=m},DaveH1(subs(p=i1,a),m,n,N,DEG,ORD,Max)):

 if gu1=FAIL then
  RETURN(FAIL):
 else
  gu:=[op(gu),[i1,gu1]]:
 fi:
od:


lu:=GR1Pol(gu,N,p):

if lu<>FAIL then 
 RETURN(lu):
fi:


for i1 from  min(S2,S1+8)+1 to S2 do
gu1:=subs({m=n,n=m},DaveH1(subs(p=i1,a),m,n,N,DEG,ORD,Max)):

 if gu1=FAIL then
  RETURN(FAIL):
 else
  gu:=[op(gu),[i1,gu1]]:
 fi:

lu:=GR1Pol(gu,N,p):

if lu<>FAIL then 
 RETURN(lu):
fi:


od:



FAIL: 
end:



#Zug(a,m,n,N,Max): 
#The pair [conj(operator) for ratio, conj. horiz. recurrence for cofactors]
#e.g. Zug(1/(m+n+1),m,n,N,20);
Zug:=proc(a,m,n,N,Max) local lu1,lu2:
lu1:=DetRatRec(a,m,n,N,Max):
if lu1=FAIL then
  RETURN(FAIL):
fi:
lu2:=DaveH(a,m,n,N,Max):
if lu2=FAIL then
 print(`The conjectured operator for the ratios is`, lu1):
 RETURN(FAIL):
fi:
[lu1,lu2]:
end:

#ZugP(a,m,n,N,Max,p,S1,S2): 
#The pair [conj(operator) for ratio, conj. horiz. recurrence for cofactors]
#where a also depends on   a parameter p that is tried from S1 to S2
#e.g. ZugP(1/(m+n+p),m,n,N,20,p,1,10);
ZugP:=proc(a,m,n,N,Max,p,S1,S2) local lu1,lu2:
lu1:=DetRatRecP(a,m,n,N,Max,p,S1,S2):
if lu1=FAIL then
  RETURN(FAIL):
fi:
lu2:=DaveHP(a,m,n,N,Max,p,S1,S2):
if lu2=FAIL then
 print(`The conjectured operator for the ratios is`, lu1):
 RETURN(FAIL):
fi:
[lu1,lu2]:
end:


#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:



#CheckZug(a,m,n,N,Z,K): checks whether the winning pair Z agrees up to K
#For example, try: Z:=Zug(1/(m+n+1),m,n,N,20);
#CheckZug(1/(m+n+1),m,n,N,Z,100);
CheckZug:=proc(a,m,n,N,Z,K) local ope1,ope2,Ini,m0,S,m1,lu,i1:
ope1:=Z[1]:
ope2:=Z[2]:
Ini:=[seq(DetRat(a,m,n,m0), m0=1..degree(ope1,N))]:
S:=SeqFromRec(ope1,n,N,Ini,K);

for m0 from 1 to K do
 lu:=DaveFromOpe(ope2,m,n,N,m0):

for m1 from 1 to m0-1 do
 if normal(add(lu[i1+1]*eval(subs({m=m1,n=i1},a)),i1=0..m0))<>0 then

  print(`The innner product of the conjectured factors for row number`, m0, `with row number`, m1 ):
  print(` of the matrix is not 0, as it should be`):
  RETURN(false):
 fi:
od:

if normal(add(lu[i1+1]*eval(subs({m=m0,n=i1},a)),i1=0..m0))<>normal(S[m0]) then
 print(`The innner product of the conjectured cofactors for row number`, m0, ` with that row of the matrix`):
 print(`is not as conjectured` ):
 RETURN(false):
fi:
od:
true:
end:




##begin DaveHI
#DaveHI1(a,m,n,M,DEG,ORD,Max): 
#Like DaveH1 but for I+A. For example, try:
#DaveHI1(binomial(m+n,n),m,n,M,3,1,10);
DaveHI1:=proc(a,m,n,M,DEG,ORD,Max) local gu,gu0,ope,m0,
Ope,m1,m2,m3:

gu:={}:

m1:=(1+DEG)*(1+ORD)+5+ORD+1:

for m0 from m1 to m1+10 do

gu0:=DaveI(a,m,n,m0):


gu0:=[op(2..nops(gu0),gu0)]:

ope:=findrec(gu0,DEG,ORD,m,M):
  if ope<>FAIL and ope<>-1 then
   gu:=gu union {[m0,ope]}:
  fi:
 od:

if nops(gu)<5 then

 RETURN(gu,FAIL):
fi:

Ope:=GR1PolRat(gu,M,n,m):


if Ope<>FAIL then

      m2:=m0+1:
         for m3 from m2 to m2+5 do

          if DaveFromOpe(Ope,m,n,M,m3)<>FAIL and DaveI(a,m,n,m3)<>DaveFromOpe(Ope,n,m,M,m3) then
                 print(Ope, ` is not quite right for`, m3):
             fi:         
        od:

RETURN(Ope):

fi:


for m0 from m1+11 to m1+Max do
gu0:=DaveI(a,m,n,m0):
gu0:=[op(2..nops(gu0),gu0)]:
ope:=findrec(gu0,DEG,ORD,m,M):
  if ope<>FAIL and ope<>-1 then
   gu:=gu union {[m0,ope]}:
    Ope:=GR1PolRat(gu,M,n,m):
     if Ope<>FAIL then

      m2:=m0+1:
         for m3 from m2 to m2+5 do
          if DaveI(a,m,n,m3)<>DaveFromOpe(Ope,m,n,M,m3) then
                  print(Ope, `failed for`, m3):
                 RETURN(FAIL):
          fi:         
         od:
      RETURN(Ope):

     fi:
    fi:
 od:

FAIL:

end:





#FindDegOrdDaveHI(a,m,n,Max): Given an expression of two non-negative
#arguments finds the Degree and Order
#annihilating A(m,n):= the cofactor of the (n,m) entry 
#of the (n+1) by (n+1) matrix (delta(i,j)+a(i,j)), 0<=i,j<=n
#Max (>20) is the largest row you are willing to try
#Here N is the shift in n
#For example, try:
#FindDegOrdDaveHI(1/(m+n+1),m,n,20);
FindDegOrdDaveHI:=proc(a,m,n,Max) local ope,m0,gu0,gu0A,gu0B,gu0C,L,DEG,ORD,opeA,
opeB,opeC,M:

ope:=FAIL:
for m0 from 20 to Max do
gu0:=DaveI(a,m,n,m0):
if gu0=FAIL then
  RETURN(FAIL):
fi:
gu0:=[op(2..nops(gu0),gu0)]:
L:=2*trunc(sqrt(nops(gu0)-4))-4:
ope:=Findrec(gu0,m,M,L):


  if ope<>FAIL then
    DEG:=degree(numer(ope),m):
    ORD:=degree(numer(ope),M):
    if findrec(gu0,DEG,ORD,m,M)<>FAIL  and  findrecK(gu0,DEG,ORD) then
      gu0A:=DaveI(a,m,n,m0+1):
      gu0A:=[op(2..nops(gu0A),gu0A)]:
      opeA:=findrecK(gu0A,DEG,ORD):

      gu0B:=DaveI(a,m,n,m0+2):
      gu0B:=[op(2..nops(gu0B),gu0B)]:
      opeB:=findrecK(gu0B,DEG,ORD):

      gu0C:=DaveI(a,m,n,m0+3):
      gu0C:=[op(2..nops(gu0B),gu0C)]:
      opeC:=findrecK(gu0C,DEG,ORD):


         if opeA and opeB and opeC then

          RETURN(DEG,ORD):
         fi:
    fi:
 fi:
od:

FAIL:
end:



#DaveHI(a,m,n,N,Max): Given an expression a (of  m and n)
#guesses a recurrence operator P(m,n,M) 
#annihilating A(m,n):= the cofactor of the (m,n) entry 
#of the (m+1) by (m+1) matrix (delta(i,j)+a(i,j)), 0<=i,j<=n
#Max (>20) is the largest row you are willing to try
#Here N is the shift in n
#For example, try:
#DaveHI(binomial(m+n,m),m,n,N,20);
DaveHI:=proc(a,m,n,N,Max) local gu,DEG,ORD:
gu:=FindDegOrdDaveHI(a,m,n,Max):


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

DEG:=gu[1]: ORD:=gu[2]:
subs({m=n,n=m},DaveHI1(a,m,n,N,DEG,ORD,Max)):
end:

#ZugI(a,m,n,N,Max): 
#The pair [conj(operator) for ratio, conj. horiz. recurrence for cofactors]
#e.g. ZugI(binomial(m+n,n),m,n,N,20);
ZugI:=proc(a,m,n,N,Max) local lu1,lu2:
lu1:=DetRatRecI(a,m,n,N,Max):
if lu1=FAIL then
  RETURN(FAIL):
fi:
lu2:=DaveHI(a,m,n,N,Max):
if lu2=FAIL then
 print(`The conjectured operator for the ratios is`, lu1):
 RETURN(FAIL):
fi:
[lu1,lu2]:
end:

#CheckZugI(a,m,n,N,Z,K): checks whether the winning pair Z agrees up to K
#For example, try: Z:=Zug(1/(m+n+1),m,n,N,20);
#CheckZugI(binomial(m+n,n),m,n,N,Z,100);
CheckZugI:=proc(a,m,n,N,Z,K) local ope1,ope2,Ini,m0,S,m1,lu,i1:
ope1:=Z[1]:
ope2:=Z[2]:

Ini:=[seq(DetRatI(a,m,n,m0), m0=1..degree(ope1,N))]:
S:=SeqFromRec(ope1,n,N,Ini,K);

for m0 from 1 to K do
 lu:=DaveFromOpe(ope2,m,n,N,m0):

for m1 from 1 to m0-1 do
 if add(lu[i1+1]*eval(subs({m=m1,n=i1},a)),i1=0..m1-1)+
add(lu[i1+1]*(1+eval(subs({m=m1,n=i1},a))),i1=m1..m1)+
add(lu[i1+1]*eval(subs({m=m1,n=i1},a)),i1=m1+1..m0)
<>0 then

  print(`The innner product of the conjectured factors for row number`, m0, `with row number`, m1 ):
  print(` of the matrix is not 0, as it should be`):
  RETURN(false):
 fi:
od:

if add(lu[i1+1]*eval(subs({m=m0,n=i1},a)),i1=0..m0)+
add(lu[i1+1],i1=m0..m0)<>S[m0] then
 print(`The innner product of the conjectured cofactors for row number`, m0, ` with that row of the matrix`):
 print(`is not as conjectured` ):
 RETURN(false):
fi:
od:
true:
end:

###end DaveHI


#SRproof(a,m,n,N,Max,K): Given an expression a of m and n
#outputs the determinant ratios and its semi-rigorous proof
#It checks that the "proof" indeed is correct up to K
#For example, try:
#SRproof(1/(m+n+1),m,n,N,30,60);
SRproof:=proc(a,m,n,N,Max,K) local Z,i,j,ku:
Z:=Zug(a,m,n,N,Max):
if Z=FAIL then
 RETURN(FAIL):
fi:

print(`Theorem: Let A(m) be the (m+1) by (m+1) matrix`):
print(` whose (i,j) entry is a(i,j):=`, subs({m=i,n=j},a)):
print(` (0<=i,j<=m) . `):
print(` Let  b(m) (m=0,1,2, ...) be the sequene annihilated by the`):
print(` recurrence operator `, subs({n=m,N=M},Z[1])):
ku:=degree(Z[1],N):
if ku=1 then
 print(`with the initial condition b(1)= `, DetRat(a,m,n,1) , ` . `):
else
print(`with the initial values from 1 to`, ku , `being `):
print([seq(DetRat(a,m,n,i),i=1..ku)], ` . ` ):
fi:
print(` Then [det A(m)]/[det A(m-1)] = b(m) .`):
print(``):
print(`Semi-Rigorous Proof: `):
print(`Let c(m,n) be the unique vector of length m+1`):
print(`(c(m,0), ..., c(m,m)) `):
print(`satisfying c(m,m)=1, c(m,n)=0, for n>m, and annihilated by`):
print(`the operator`, Z[2], ` . `):
print(`I claim that `):
print(`(1) Sum(c(k,n)*a(m,n),n=0..m)=0 for 0<=k<m and `):
print(`(2) Sum(c(m,n)*a(m,n),n=0..m)=b(m) . `):
print(`These two statements imply the theorem since (1)`):
print(` is the defining property for the cofactors of the m-th row`):
print(`up to normalization, and with the normalization such that  c(n,n)=1, `):
print(` which is tantamount to dividing the cofactors of the m-th row by the cofactor `):
print(` of the (m,m) entry, which is det(A(m-1)),`):
print(`the left of (2) is det[A(m)]/det[A(m-1)].`):
print(`But both c(m,n) and b(m) are holonomic functions,`):
print(` hence (1) and (2) are routinely provable within the Holonomic Paradigm.`):
print(` Since we know for sure that a proof exists (if it is true),`):
print(` and we know how to find it, we don't feel that it is necessary`):
print(` to actually bother the computer with the time- and memory- consuming chore of finding it.`):
print( ` Instead,  we will be content to verify`):
print(` (1) for 0<=k < m <=`, K, `,`):
print(` and (2) for 0<=m <=`, K, ` . `):
print(`For this range these are indeed`):
print(CheckZug(a,m,n,N,Z,K) , ` QED(sr) . `  ):
end:

#DetRatDirect(a,m,n,n0): det(a(i,j), 0<=i,j<=n0)/det(a(i,j), 0<=i,j<=n0-1) Directly
#here a is an expression in m and n
#For example, try:DetRatDirect(1/(m+n+1),m,n,10);
DetRatDirect:=proc(a,m,n,n0) local A,B,i,j:
if n0=0 then
 RETURN(eval(subs({m=0,n=0},a))):
fi:
A:=matrix(n0,n0):
B:=matrix(n0+1,n0+1):

for i from 0 to n0-1 do
for j from 0 to n0-1 do
 A[i+1,j+1]:=subs({m=i,n=j},a):
 B[i+1,j+1]:=subs({m=i,n=j},a):
od:
od:

for i from n0 to n0 do
for j from 0 to n0 do
 B[i+1,j+1]:=subs({m=i,n=j},a):
od:
od:


for i from 0 to n0-1 do
for j from n0 to n0 do
 B[i+1,j+1]:=subs({m=i,n=j},a):
od:
od:


if det(A)<>0 then
det(B)/det(A):
else
FAIL:
fi:
end:




#SRproofI(a,m,n,N,Max,K): Given an expression a of m and n
#outputs the determinant ratios of det(a(i,j)+delta(i,j)) and its semi-rigorous proof
#It checks that the "proof" indeed is correct up to K
#For example, try:
#SRproofI(binomial(m+n),m,n,N,30,60);
SRproofI:=proc(a,m,n,N,Max,K) local Z,i,j,ku:
Z:=ZugI(a,m,n,N,Max):
if Z=FAIL then
 RETURN(FAIL):
fi:

print(`Theorem: Let A(m) be the (m+1) by (m+1) matrix`):
print(` whose (i,j) entry is a(i,j):=`, subs({m=i,n=j},a)):
print(` (0<=i,j<=m) , and let I(m) be the (m+1) by (m+1) identity matrix `):
print(` Let  b(m) (m=1,2, ...) be the sequene annihilated by the`):
print(` recurrence operator `, subs({N=M,n=m},Z[1])):
ku:=degree(Z[1],N):
if ku=1 then
 print(`with the initial condition b(1)= `, DetRat(a,m,n,1) , ` . `):
else
print(`with the initial values from 1 to`, ku , `being `):
print([seq(DetRatI(a,m,n,i),i=1..ku)], ` . ` ):
fi:

print(` Then [det (I(m)+A(m))]/[det (I(m-1)+A(m-1))] = b(m) .`):
print(``):
print(`Semi-Rigorous Proof: `):
print(`Let c(m,n) be the unique vector of length m+1`):
print(`(c(m,0), ..., c(m,m)) `):
print(`satisfying c(m,m)=1, c(m,n)=0, for n>m, and annihilated by`):
print(`the operator`, Z[2], ` . `):
print(`I claim that `):
print(`(1) Sum(c(k,n)*a(m,n),n=0..m)=0 for 0<=k<m and `):
print(`(2) Sum(c(m,n)*a(m,n),n=0..m)=b(m) . `):
print(`These two statements imply the theorem since (1)`):
print(` is the defining property for the cofactors of the m-th row`):
print(`up to normalization, and with the normalization such that  c(m,m)=1, `):
print(` which is tantamount to dividing the cofactors of the m-th row by the cofactor `):
print(` of the (m,m) entry, which is det(A(m-1)),`):
print(`the left of (2) is det[A(m)]/det[A(m-1)].`):
print(`But both c(m,n) and b(m) are holonomic functions,`):
print(` hence (1) and (2) are routinely provable within the Holonomic Paradigm.`):
print(` Since we know for sure that a proof exists (if it is true),`):
print(` and we know how to find it, we don't feel that it is necessary`):
print(` to actually bother the computer with the time- and memory- consuming chore of finding it.`):
print( ` Instead,  we will be content to verify`):
print(` (1) for 0<=k < m <=`, K, `,`):
print(` and (2) for 0<=m <=`, K, ` . `):
print(`For this range these are indeed`):
print(CheckZugI(a,m,n,N,Z,K) , ` QED(sr) . `  ):
end:


#OpeToCF(ope,n,N): converts ope(n,N) of first order to Closed-Form with f(0)=1
#For example, try: opeToCF(N-n-1,n,N);
OpeToCF:=proc(ope,n,N) local g:
if degree(ope,N)<>1 then
 RETURN(FAIL):
fi:
rsolve({coeff(ope,N,1)*g(n+1)+coeff(ope,N,0)*g(n)=0,g(0)=1},g(n)):
end:


#OpeToCF1(ope,m,n,N): Solves ope(m,n,N)f(n)=0, f(n)=1
#For example, try: opeToCF(N-(m-n),m,n,N);
OpeToCF1:=proc(ope,m,n,N) local g,i:
if degree(ope,N)<>1 then
 RETURN(FAIL):
fi:
subs(i=m-n,rsolve({subs(n=m-i-1,coeff(ope,N,0))*g(i+1)+
subs(n=m-i-1,coeff(ope,N,1))*g(i)=0,g(0)=1},g(i))):
end:




#Rproof(a,m,n,N,Max,R): Given an expression a of m and n
#outputs the determinant ratios and its rigorous proof
#by WZ, in the lucky case where both the ratio and
#the cofactors are closed-form
#R is the symbol for the rising factorial R(a,n)=a(a+1)*...*(a+n-1)
#For example, try:
#Rproof(1/(m+n+1),m,n,N,30,R);
Rproof:=proc(a,m,n,N,Max,R) local Z,i,j,b,c,lu1,lu2,M,k,ope,c1,n0,b2,c2:

Z:=Zug(a,m,n,N,Max):


if Z=FAIL then
 RETURN(FAIL):
fi:

if degree(Z[1],N)<>1 then
 print(`The sequence of determinant ratios does not obey a first-order recurrence`):
 print(`and is (probably) not closed-form`):
 RETURN(FAIL):
fi:

if degree(Z[2],N)<>1 then
 print(`The cofactors double-sequence does not obey a first-order recurrence`):
 print(`and is (probably) not closed-form`):
 RETURN(FAIL):
fi:

b:=PtorOpeR(Z[1],n,N,R):
if b=FAIL then
 RETURN(Z):
fi:

b:=normal(simplify(expand(DetRat(a,m,n,1))))*b:



if not evalb([seq(eval(subs(R=rf1,subs(n=n0,b))),n0=1..10)]=[seq(DetRatDirect(a,m,n,n0),n0=1..10)]) then
print(`Something is wrong`):
RETURN(FAIL):
fi:


c:=PtorOpeR(Z[2],n,N,R):

if c=FAIL then
 print(`Couldn't solve`, Z[2]):
 RETURN(Z):
fi:

c:=normal(c/subs(n=m,c)):

c2:=[seq([seq(eval(subs({m=m0,n=n0},subs(R=rf1,c))),n0=0..m0)],m0=1..10)]:
if not evalb(c2=[seq(Dave(a,m,n,m0),m0=1..10)])then
print(`Something is wrong`):
RETURN(FAIL):
fi:


c1:=normal(OpeToCF1(Z[2],m,n,N)):

print(`Notation:  R(a,k) denotes the rising factorial a(a+1)...(a+k-1)`):
print(`Theorem: Let A(m) be the (m+1) by (m+1) matrix`):
print(` whose (i,j) entry is a(i,j):=`, subs({m=i,n=j},a)):
print(` (0<=i,j<=m) . `):
print(`Let b(m):= `):
print(subs(n=m,b), ` . `):
print(` Then [det A(m)]/[det A(m-1)] = b(m)`):
print(` and hence det [A(m)]=b(1)b(2)...b(m) .`):
print(``):
print(`Rigorous Proof: `):

print(`Let c(m,n) be`):
print(c, ` . ` ):
print(`Note that c(m,m)=1, c(m,n)=0, for n>m .`):
print(`I claim that `):
print(`(1) Sum(c(k,n)*a(m,n),n=0..m)=0 for 0<=k<m and `):
print(`(2) Sum(c(m,n)*a(m,n),n=0..m)=b(m) for m>=0 . `):
print(`These two statements imply the theorem since (1)`):
print(` is the defining property for the cofactors of the m-th row`):
print(`up to normalization, and with the normalization such that  c(m,m)=1, `):
print(` which is tantamount to dividing the cofactors of the m-th row by the cofactor `):
print(` of the (m,m) entry, which is det(A(m-1)),`):
print(`the left of (2) is det[A(m)]/det[A(m-1)].`):

lu1:=zeil(a*subs(m=k,c1),n,m,M):
print(`According to EKHAD, the left side of (1) is annihilated by the operator`):
print(lu1[1]):
print(`the certificate being`, lu1[2], ` . `):
print(`This implies that indeed the sum is zero for k<m .`):
lu2:=zeil(a*c1,n,m,M):
print(`Also according to EKHAD, the sum on the left of (2) is annihilated by`):

ope:=Yafe(lu2[1],M)[2]:
print(ope):
if ope=subs({n=m,N=M},Z[1]) then
print(`which is the same operator annihilating`, subs(n=m,b)):
 else
  print(`Something is wrong`):
 RETURN(FAIL):
fi:

print(`the certificate being`):
print(lu2[2] , ` QED .`):

end:


#RproofPold(a,m,n,N,Max,p,S1,S2): Given an expression a of m and n
#that also depends on a parameter p
#(it tries for p=S1..S2)
#outputs the determinant ratios and its rigorous proof
#by WZ, in the lucky case where both the ratio and
#the cofactors are closed-form
#For example, try:
#RproofPold(1/(m+n+p),m,n,N,30,p,1,10);
RproofPold:=proc(a,m,n,N,Max,p,S1,S2) local Z,i,j,b,c,lu1,lu2,M,k,ope:

Z:=ZugP(a,m,n,N,Max,p,S1,S2):


if Z=FAIL then
 RETURN(FAIL):
fi:

if degree(Z[1],N)<>1 then
 print(`The sequence of determinant ratios does not obey a first-order recurrence`):
 print(`and is (probably) not closed-form`):
 RETURN(FAIL):
fi:

if degree(Z[2],N)<>1 then
 print(`The cofactors double-sequence does not obey a first-order recurrence`):
 print(`and is (probably) not closed-form`):
 RETURN(FAIL):
fi:

b:=normal(DetRat(a,m,n,1))*OpeToCF(Z[1],n,N):
c:=normal(OpeToCF1(Z[2],m,n,N)):

b:=convert(b,factorial):
b:=normal(expand(b)):
c:=convert(c,factorial):
print(`Theorem: Let A(m)=(a(i,j)) be the (m+1) by (m+1) matrix`):
print(` whose (i,j) entry is a(i,j):= `, subs({m=i,n=j},a)):
print(` (0<=i,j<=m) . `):
print(` Then [det A(m)]/[det A(m-1)] =`):
print(subs(n=m,b) , ` . `):
print(` and hence det [A(m)]=b(1)b(2)...b(m) .`):
print(``):
print(`Rigorous Proof: `):

print(`Let c(m,n) be`):
print(c,  ` . `):
print(`Note that c(m,m)=1, c(m,n)=0, for n>m .`):
print(`I claim that `):
print(`(1) Sum(c(k,n)*a(m,n),n=0..m)=0 for 0<=k<m and `):
print(`(2) Sum(c(m,n)*a(m,n),n=0..m)=b(m) . `):
print(`These two statements imply the theorem since (1)`):
print(` is the defining property for the cofactors of the m-th row`):
print(`up to normalization, and with the normalization such that  c(m,m)=1, `):
print(` which is tantamount to dividing the cofactors of the m-th row by the cofactor `):
print(` of the (m,m) entry, which is det(A(m-1)),`):
print(`the left of (2) is det[A(m)]/det[A(m-1)].`):

lu1:=zeil(a*subs(m=k,c),n,m,M):
print(`According to EKHAD, the left side of (1) is annihilated by the operator`):
print(lu1[1], ` , `):
print(`the certificate being`, lu1[2] , ` . `):
print(`This implies that indeed the sum is zero for k<m .`):
lu2:=zeil(c*subs(m=k,c),n,m,M):
print(`Also according to EKHAD, the sum on the left of (2) is annihilated by`):
lu2:=zeil(a*c,n,m,M):
ope:=Yafe(lu2[1],M)[2]:
print(ope):
if ope=subs({n=m,N=M},Z[1]) then
print(`which is the same operator annihilating`, subs(n=m,b), ` , `):
 else
  print(`Something is wrong`):
 RETURN(FAIL):
fi:

print(`the certificate being`):
print(lu2[2], ` QED . `):

end:






#SRproofP(a,m,n,N,Max,p,S1,S2,K): Given an expression a of m and n
#outputs the determinant ratios and its semi-rigorous proof
#It checks that the "proof" indeed is correct up to K
#For example, try:
#SRproofP(1/(m+n+p),m,n,N,30,p,1,10,60);
SRproofP:=proc(a,m,n,N,Max,p,S1,S2,K) local Z,i,j,ku,M:
Z:=ZugP(a,m,n,N,Max,p,S1,S2):
if Z=FAIL then
 RETURN(FAIL):
fi:

print(`Theorem: Let A(m) be the (m+1) by (m+1) matrix`):
print(` whose (i,j) entry is a(i,j):=`, subs({m=i,n=j},a)):
print(` (0<=i,j<=m) . `):
print(` Let  b(m) (m=0,1,2, ...) be the sequene annihilated by the`):
print(` recurrence operator `, subs({n=m,N=M},Z[1])):
ku:=degree(Z[1],N):
if ku=1 then
 print(`with the initial condition b(1)= `, normal(simplify(expand(DetRat(a,m,n,1)))) , ` . `):
else
print(`with the initial values from 1 to`, ku , `being `):
print([seq(DetRat(a,m,n,i),i=1..ku)], ` . ` ):
fi:
print(` Then [det A(m)]/[det A(m-1)] = b(m) .`):
print(``):
print(`Semi-Rigorous Proof: `):
print(`Let c(m,n) be the unique vector of length m+1`):
print(`(c(m,0), ..., c(m,m)) `):
print(`satisfying c(m,m)=1, c(m,n)=0, for n>m, and annihilated by`):
print(`the operator`, Z[2], ` . `):
print(`I claim that `):
print(`(1) Sum(c(k,n)*a(m,n),n=0..m)=0 for 0<=k<m and `):
print(`(2) Sum(c(m,n)*a(m,n),n=0..m)=b(m) . `):
print(`These two statements imply the theorem since (1)`):
print(` is the defining property for the cofactors of the m-th row`):
print(`up to normalization, and with the normalization such that  c(n,n)=1, `):
print(` which is tantamount to dividing the cofactors of the m-th row by the cofactor `):
print(` of the (m,m) entry, which is det(A(m-1)),`):
print(`the left of (2) is det[A(m)]/det[A(m-1)].`):
print(`But both c(m,n) and b(m) are holonomic functions,`):
print(` hence (1) and (2) are routinely provable within the Holonomic Paradigm.`):
print(` Since we know for sure that a proof exists (if it is true),`):
print(` and we know how to find it, we don't feel that it is necessary`):
print(` to actually bother the computer with the time- and memory- consuming chore of finding it.`):
print( ` Instead,  we will be content to verify`):
print(` (1) for 0<=k < m <=`, K, `,`):
print(` and (2) for 0<=m <=`, K, ` . `):
print(`For this range these are indeed`):
print(CheckZug(a,m,n,N,Z,K) , ` QED(sr) . `  ):
end:


#PtorOpeR(ope,n,N,R): solving the recurrence ope(n,N)x(n)=0, x(0)=1
#in terms of rising factorial R(n,k) (R is a symbol)
#For example try PtorOpeR((n+1)*N-1,n,N,R);
PtorOpeR:=proc(ope,n,N,R)
local rat,c,t1,b1,t1c,b1c,i,ku:

if degree(ope,N)<>1 then
 RETURN(FAIL):
fi:

ku:=1:
rat:=normal(-coeff(ope,N,0)/coeff(ope,N,1)):
t1:=expand(numer(rat)):
b1:=expand(denom(rat)):

t1c:=coeff(t1,n,degree(t1,n)):
b1c:=coeff(b1,n,degree(b1,n)):

if b1c=0 or t1c=0 then
 RETURN(FAIL):
fi:

c:=t1c/b1c:

ku:=c:
t1:=expand(t1/t1c):
b1:=expand(b1/b1c):


t1:=factor(t1):
b1:=factor(b1):

t1:=[solve(t1,n)]:
b1:= [solve(b1,n)]:

if convert([seq(-b1[i],i=1..nops(b1))],`*`)=0 then
 RETURN(FAIL):
fi:
ku:=ku*convert([seq(-t1[i],i=1..nops(t1))],`*`)/
convert([seq(-b1[i],i=1..nops(b1))],`*`):

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

c^n* convert([seq(R(-t1[i],n),i=1..nops(t1))],`*`)/
convert([seq(R(-b1[i],n),i=1..nops(b1))],`*`)/ku:

end:





#RproofP(a,m,n,N,Max,R,p,S1,S2): Given an expression a of m and n
#outputs the determinant ratios and its rigorous proof
#by WZ, in the lucky case where both the ratio and
#the cofactors are closed-form
#R is the symbol for the rising factorial R(a,n)=a(a+1)*...*(a+n-1)
#For example, try:
#RproofP(1/(m+n+1),m,n,N,30,R,p,S1,S2);
RproofP:=proc(a,m,n,N,Max,R,p,S1,S2) local Z,i,j,b,c,lu1,lu2,M,k,ope,c1:

Z:=ZugP(a,m,n,N,Max,p,S1,S2):


if Z=FAIL then
 RETURN(FAIL):
fi:

if degree(Z[1],N)<>1 then
 print(`The sequence of determinant ratios does not obey a first-order recurrence`):
 print(`and is (probably) not closed-form`):
 RETURN(FAIL):
fi:

if degree(Z[2],N)<>1 then
 print(`The cofactors double-sequence does not obey a first-order recurrence`):
 print(`and is (probably) not closed-form`):
 RETURN(FAIL):
fi:

b:=PtorOpeR(Z[1],n,N,R):
if b=FAIL then
 print(`Couldn't solve`, Z[1]):
 RETURN(Z):
fi:

b:=normal(normal(expand(convert(simplify(DetRat(a,m,n,1)),factorial)))*b):


c:=PtorOpeR(Z[2],n,N,R):
if c=FAIL then
 print(`Couldn't solve`, Z[2]):
 RETURN(Z):
fi:

c:=normal(c/subs(n=m,c)):

c1:=normal(OpeToCF1(Z[2],m,n,N)):

print(`Notation:  R(a,k) denotes the rising factorial a(a+1)...(a+k-1)`):
print(`Theorem: Let A(m) be the (m+1) by (m+1) matrix`):
print(` whose (i,j) entry is a(i,j):=`, subs({m=i,n=j},a)):
print(` (0<=i,j<=m) . `):
print(`Let b(m):= `):
print(subs(n=m,b), ` . `):
print(` Then [det A(m)]/[det A(m-1)] = b(m)`):
print(` and hence det [A(m)]=b(1)b(2)...b(m) .`):
print(``):
print(`Rigorous Proof: `):

print(`Let c(m,n) be`):
print(c, ` . ` ):
print(`Note that c(m,m)=1, c(m,n)=0, for n>m .`):
print(`I claim that `):
print(`(1) Sum(c(k,n)*a(m,n),n=0..m)=0 for 0<=k<m and `):
print(`(2) Sum(c(m,n)*a(m,n),n=0..m)=b(m) for m>=0 . `):
print(`These two statements imply the theorem since (1)`):
print(` is the defining property for the cofactors of the m-th row`):
print(`up to normalization, and with the normalization such that  c(m,m)=1, `):
print(` which is tantamount to dividing the cofactors of the m-th row by the cofactor `):
print(` of the (m,m) entry, which is det(A(m-1)),`):
print(`the left of (2) is det[A(m)]/det[A(m-1)].`):

lu1:=zeil(a*subs(m=k,c1),n,m,M):
print(`According to EKHAD, the left side of (1) is annihilated by the operator`):
print(lu1[1]):
print(`the certificate being`, lu1[2], ` . `):
print(`This implies that indeed the sum is zero for k<m .`):
lu2:=zeil(a*c1,n,m,M):
print(`Also according to EKHAD, the sum on the left of (2) is annihilated by`):

ope:=Yafe(lu2[1],M)[2]:
print(ope):
if ope=subs({n=m,N=M},Z[1]) then
print(`which is the same operator annihilating`, subs(n=m,b)):
 else
  print(`Something is wrong`):
 RETURN(FAIL):
fi:

print(`the certificate being`):
print(lu2[2] , ` QED .`):

end:

rf1:=proc(a,b) local i:
 mul(a+i,i=0..b-1):
end: