Help:=proc(): print(` GD(a,n) `): end:

#GD(a,n): inputs a hypergeometric sequence a(n)
#and outputs three polynomials p(n), q(n), r(n)
#such that a(n)/a(n-1)=(p(n)/p(n-1))*q(n)/r(n)
#and gcd(q(n),r(n+j))=1 for j=0..infinity
#infinity:=1000
GD:=proc(a,n) local p,q,r,j,R,g,i:

R:=simplify(a/subs(n=n-1,a));

p:=1:
q:=numer(R):
r:=denom(R):

if not (type(q,polynom) and type(r,polynom)) then
 ERROR(`Not Hypergeometric, Stupid! (unless Maple is)`):
fi:

for j from 1 to 1000 do
 g:=gcd(q,subs(n=n+j,r)):
  if g<>1 then
   q:=q/g:
   r:=r/subs(n=n-j,g):
   p:=p*mul(subs(n=n-i,g),i=0..j-1):
  fi:
od:

[p,q,r]:

end:

#Gosper(a,n): inputs a hyper. sequence a of n
#and outputs S(n) also hyper s.t. S(n)-S(n-1)=a(n)
#or FAIL
Gosper:=proc(a,n) local f,d,F,c,Mike,eq:

Mike:=GD(a,n):

p:=Mike[1]: q:=Mike[2]; r:=Mike[3];

#find d in a more precise way
d:=max(degree(p,n)-degree(q,n),degree(p,n)-degree(r,n))+10:

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):

a*f/p*subs(n=n+1,q);

end: