Help:=proc(): print(` Findd(q,r,n), FindDegf(p,q,r,n) `): print(`OneStep(L,n), GD(p0,a,n), Diff1(S,n) , Gosper(p0,a,n) `): print(`CheckGosper(S,n), GosperWP(p0,a,n,Para), WZ(F,a,n,k)`): end: #Findd(q,r,n): given polynomials q and r of n #finds out, if it is possible to have a non-negative #integer such that f(n):=n^d+c1*n^(d-1)+ ... can have #q(n+1)*f(n)-r(n)*f(n-1) have a vanishing leading term Findd:=proc(q,r,n) local A,B,F,L,d,c1: A:=degree(q,n): B:=degree(r,n): if A<>B then RETURN(FAIL): fi: if coeff(q,n,A)<>coeff(r,n,A) then RETURN(FAIL): fi: F:=expand(subs(n=n+1,q)*(1+c1/n)-r*(1+(c1-d)/n)): L:=coeff(F,n,degree(F,n)): d:=solve(L,d): if type(d,integer) and d>=0 then RETURN(d): else RETURN(FAIL): fi: end: #FindDegf(p,q,r,n): finds the degree of f(n) in the #functional equation: q(n+1)*f(n)-r(n)*f(n-1)=p(n) FindDegf:=proc(p,q,r,n) local d0,d1: d0:=Findd(q,r,n): d1:=max(degree(p,n)-degree(r,n),degree(p,n)-degree(q,n)): if d0=FAIL then RETURN(d1): else RETURN( max(d0,d1)): fi: end: #OneStep(L,n): Given a Gosper triple L=[p,q,r] (polynomials #in n) finds the next one obtained by finding g such that #such that g(n):=gcd(q(n),r(n+j)) <>1 #or returs the same thing OneStep:=proc(L,n) local S1,S2,a,b,j,p1,q1,r1,p,q,r,g,i: p:=L[1]: q:=L[2]: r:=L[3]: S1:={solve(q,n)}: S2:={solve(r,n)}: for a in S1 do for b in S2 do j:=b-a: if type(j,integer) and j>=0 then g:=gcd(q,subs(n=n+j,r)): r1:=normal(r/subs(n=n-j,g)): q1:=normal(q/g): p1:=normal(p*mul(subs(n=n-i,g),i=0..j-1)): RETURN([p1,q1,r1]): fi: od: od: L: end: #GD(p0,a,n): inputs a poly p0 and #hypergeometric sequence a(n) #and outputs three polynomials p(n), q(n), r(n) #such that p0(n)*a(n)/a(n-1)/p0(n-1)=(p(n)/p(n-1))*q(n)/r(n) #and gcd(q(n),r(n+j))=1 for j=0..infinity GD:=proc(p0,a,n) local R,p,q,r,newL,oldL: R:=normal(simplify(expand(simplify(a/subs(n=n-1,a))))); p:=p0: q:=normal(numer(R)): r:=normal(denom(R)): if not (type(q,polynom) and type(r,polynom)) then ERROR(`Not Hypergeometric, Stupid! (unless Maple is)`): fi: oldL:=[p,q,r]: newL:=OneStep(oldL,n) : while oldL<>newL do oldL:=newL: newL:=OneStep(newL,n) : od: newL: end: #Diff1(S,n): Given a hypergeometric sequence S(n) #finds the hypergeometric sequence S(n)-S(n-1) Diff1:=proc(S,n) local a,p0: a:=normal(1-simplify(expand(subs(n=n-1,S)/S))): p0:=numer(a): p0,S*normal(a/p0): end: #Gosper(p0,a,n): inputs a polynomial p and a hyper. sequence a of n #and outputs S(n) also hyper s.t. S(n)-S(n-1)=p(n)*a(n) #or FAIL Gosper:=proc(p0,a,n) local f,d,F,c,Mike,eq,p,q,r,var,i: Mike:=GD(p0,a,n): p:=Mike[1]: q:=Mike[2]; r:=Mike[3]; d:=FindDegf(p,q,r,n): f:=add(c[i]*n^i,i=0..d): var:={ seq(c[i],i=0..d)}: F:=expand(subs(n=n+1,q)*f-r*subs(n=n-1,f)-p): eq:={seq(coeff(F,n,i),i=0..degree(F,n))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: f:=subs(var,f): normal(p0*a*f/p*subs(n=n+1,q)); end: CheckGosper:=proc(S,n): normal(Gosper(Diff1(S,n),n)-S): end: #OneStepG(L,n): Given a Gosper triple L=[p,q,r] (polynomials #in n) finds the next one obtained by finding g such that #such that g(n):=gcd(q(n),r(n+j)) <>1 #or returs the same thing OneStepG:=proc(L,n) local S1,S2,a,b,j,p1,q1,r1,p,q,r,g,i,R,S: p:=L[1]: q:=L[2]: r:=L[3]: R:=resultant(q,subs(n=n+j,r),n): S:={solve(R,j)}: for j in S do if type(j,integer) and j>=0 then g:=gcd(q,subs(n=n+j,r)): r1:=normal(r/subs(n=n-j,g)): q1:=normal(q/g): p1:=normal(p*mul(subs(n=n-i,g),i=0..j-1)): RETURN([p1,q1,r1]): fi: od: [p,q,r]: end: #GDG(p0,a,n): inputs a poly p0 and #hypergeometric sequence a(n) #and outputs three polynomials p(n), q(n), r(n) #such that p0(n)*a(n)/a(n-1)/p0(n-1)=(p(n)/p(n-1))*q(n)/r(n) #and gcd(q(n),r(n+j))=1 for j=0..infinity GDG:=proc(p0,a,n) local R,p,q,r,newL,oldL: R:=normal(simplify(expand(simplify(a/subs(n=n-1,a))))); p:=p0: q:=normal(numer(R)): r:=normal(denom(R)): if not (type(q,polynom) and type(r,polynom)) then ERROR(`Not Hypergeometric, Stupid! (unless Maple is)`): fi: oldL:=[p,q,r]: newL:=OneStepG(oldL,n) : while oldL<>newL do oldL:=newL: newL:=OneStepG(newL,n) : od: newL: end: #GosperG(p0,a,n): inputs a polynomial p and a hyper. sequence a of n #and outputs S(n) also hyper s.t. S(n)-S(n-1)=p(n)*a(n) #or FAIL GosperG:=proc(p0,a,n) local f,d,F,c,Mike,eq,p,q,r,var,i: Mike:=GDG(p0,a,n): p:=Mike[1]: q:=Mike[2]; r:=Mike[3]; d:=FindDegf(p,q,r,n): f:=add(c[i]*n^i,i=0..d): var:={ seq(c[i],i=0..d)}: F:=expand(subs(n=n+1,q)*f-r*subs(n=n-1,f)-p): eq:={seq(coeff(F,n,i),i=0..degree(F,n))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: f:=subs(var,f): normal(p0*a*f/p*subs(n=n+1,q)); end: CheckGosperG:=proc(S,n): normal(GosperG(Diff1(S,n),n)-S): end: #GosperWP(p0,a,n,Para): inputs a polynomial p and a hyper. sequence a of n #and outputs S(n) and a choice of the parameters in the #set of parametrs Para that will this come true: #also hyper s.t. S(n)-S(n-1)=p(n)*a(n) #or FAIL #For example GosperW(n+c0,n!,n,{c0}); should give #n!, {c0=0} GosperWP:=proc(p0,a,n,Para) local f,d,F,c,Mike,eq,p,q,r,var,i,p1: Mike:=GD(p0,a,n): p:=Mike[1]: q:=Mike[2]; r:=Mike[3]; d:=FindDegf(p,q,r,n): f:=add(c[i]*n^i,i=0..d): var:={ seq(c[i],i=0..d)}: F:=expand(subs(n=n+1,q)*f-r*subs(n=n-1,f)-p): eq:={seq(coeff(F,n,i),i=0..degree(F,n))}: var:=solve(eq,var union Para): if var=NULL then RETURN(FAIL): fi: f:=subs(var,f): normal(p0*a*f/p*subs(n=n+1,q)),{seq(p1=subs(var,p1),p1 in Para)}; end: #WZ(F,a,n,k): inputs a closed-form F(n,k) (product of #factorials of linear forms in n and k, and a hyper. #a(n), and outputs the WZ-mate G(n,k) WZ:=proc(F,a,n,k) local F1, G,Justin: F1:=F/a: Justin:=normal(1-expand(subs(n=n-1,F1)/F1)): Gosper(numer(Justin),F1/denom(Justin),k): end: