gosper.py

from sympy import *
from functools import reduce
import operator
x, t = symbols('x t')

def gosper_form(u_1, u_2):
    res = resultant(u_1, u_2.subs(x, x + t), x)
    roots = poly(res, t).all_roots()
    int_roots = [rt for rt in roots if rt.is_Integer and rt > 0]
    if int_roots == []:
        return (1, u_1, u_2.subs(x, x - 1))
    else:
        i = max(int_roots)
        g = gcd(u_1, u_2.subs(x, x + i))
        p = reduce(operator.mul, [g.subs(x, x-j) for j in range(1, i + 1)], 1)
        (pp, q, r) = gosper_form(u_1 / g, u_2 / g.subs(x, x - i))
        return (p * pp, q, r)

(p, q, r) = gosper_form(x, x - 2)
print((p, q, r))
print(simplify(poly(p, x).subs(x, x + 1) / p * q / poly(r, x).subs(x, x + 1)))