Finalized implementation of Gosper's algorithm

In [1]: from sympy.concrete.gosper import gosper_sum

In [2]: gosper_sum(binomial(2*k, k)/4**k, (k, 0, n))
Out[2]:
 -n
4  ⋅(1 + 2⋅n)⋅Binomial(2⋅n, n)

In [3]: gosper_sum(_, (n, 0, n))
Out[3]:
 -n
4  ⋅(1 + 2⋅n)⋅(3 + 2⋅n)⋅Binomial(2⋅n, n)
────────────────────────────────────────
                   3

In [4]: gosper_sum(_, (n, 0, n))
Out[4]:
 -n
4  ⋅(1 + 2⋅n)⋅(3 + 2⋅n)⋅(5 + 2⋅n)⋅Binomial(2⋅n, n)
──────────────────────────────────────────────────
                        15

sympy/concrete/__init__.py

@@ -1,3 +1,2 @@
 from products import product, Product
 from summations import summation, Sum
-from gosper import normal, gosper

sympy/concrete/gosper.py

@@ -1,100 +1,200 @@
-from sympy.core.singleton import S
-from sympy.core.symbol import Dummy
-from sympy.core.mul import Mul
-from sympy.core import sympify
+"""Gosper's algorithm for hypergeometric summation. """
 
-from sympy.polys import gcd, quo, roots, resultant
+from sympy.core import S, Dummy, sympify, symbols
+from sympy.polys import Poly, parallel_poly_from_expr, factor
+from sympy.solvers import solve
+from sympy.simplify import hypersimp
 
-def normal(f, g, n=None):
-    """Given relatively prime univariate polynomials 'f' and 'g',
-       rewrite their quotient to a normal form defined as follows:
+def gosper_normal(f, g, n, polys=True):
+    r"""
+    Compute the Gosper's normal form of ``f`` and ``g``.
 
-                       f(n)       A(n) C(n+1)
-                       ----  =  Z -----------
-                       g(n)       B(n)  C(n)
+    Given relatively prime univariate polynomials ``f`` and ``g``,
+    rewrite their quotient to a normal form defined as follows::
 
-       where Z is arbitrary constant and A, B, C are monic
-       polynomials in 'n' with following properties:
+    .. math::
 
-           (1) gcd(A(n), B(n+h)) = 1 for all 'h' in N
-           (2) gcd(B(n), C(n+1)) = 1
-           (3) gcd(A(n), C(n)) = 1
+        \frac{f(n)}{g(n)} = Z \cdot \frac{A(n) C(n+1)}{B(n) C(n)}
 
-       This normal form, or rational factorization in other words,
-       is crucial step in Gosper's algorithm and in difference
-       equations solving. It can be also used to decide if two
-       hypergeometric are similar or not.
+    where ``Z`` is an arbitrary constant and ``A``, ``B``, ``C`` are
+    monic polynomials in ``n`` with the following properties:
 
-       This procedure will return a tuple containing elements
-       of this factorization in the form (Z*A, B, C). For example:
+    1. $\gcd(A(n), B(n+h)) = 1 \forall h \in \mathbb{N}$
+    2. $\gcd(B(n), C(n+1)) = 1$
+    3. $\gcd(A(n), C(n)) = 1$
 
-       >>> from sympy import Symbol, normal
-       >>> n = Symbol('n', integer=True)
+    This normal form, or rational factorization in other words, is a
+    crucial step in Gosper's algorithm and in solving of difference
+    equations. It can be also used to decide if two hypergeometric
+    terms are similar or not.
 
-       >>> normal(4*n+5, 2*(4*n+1)*(2*n+3), n)
-       (1/4, 3/2 + n, 1/4 + n)
+    This procedure will return a tuple containing elements of this
+    factorization in the form ``(Z*A, B, C)``.
 
-    """
-    f, g = map(sympify, (f, g))
+    **Examples**
 
-    p = f.as_poly(n, field=True)
-    q = g.as_poly(n, field=True)
+    >>> from sympy.concrete.gosper import gosper_normal
+    >>> from sympy.abc import n
 
-    a, p = p.LC(), p.monic()
-    b, q = q.LC(), q.monic()
+    >>> gosper_normal(4*n+5, 2*(4*n+1)*(2*n+3), n, polys=False)
+    (1/4, 3/2 + n, 1/4 + n)
 
-    A = p.as_expr()
-    B = q.as_expr()
+    """
+    (p, q), opt = parallel_poly_from_expr((f, g), n, field=True)
 
-    C, Z = S.One, a / b
+    a, A = p.LC(), p.monic()
+    b, B = q.LC(), q.monic()
 
+    C, Z = A.one, a/b
     h = Dummy('h')
 
-    res = resultant(A, B.subs(n, n+h), n)
+    D = Poly(n + h, n, h, domain=opt.domain)
+
+    R = A.resultant(B.compose(D))
+    roots = set(R.ground_roots().keys())
+
+    for r in set(roots):
+        if not r.is_Integer or r < 0:
+            roots.remove(r)
+
+    for i in sorted(roots):
+        d = A.gcd(B.shift(+i))
+
+        A = A.quo(d)
+        B = B.quo(d.shift(-i))
+
+        for j in xrange(1, i+1):
+            C *= d.shift(-j)
+
+    A = A.mul_ground(Z)
+
+    if not polys:
+        A = A.as_expr()
+        B = B.as_expr()
+        C = C.as_expr()
+
+    return A, B, C
+
+def gosper_term(f, n):
+    """
+    Compute Gosper's hypergeometric term for ``f``.
+
+    Suppose ``f`` is a hypergeometric term such that::
+
+    .. math::
+
+        s_n = \sum_{k=0}^{n-1} f_k
+
+    and $f_k$ doesn't depend on $n$. Returns a hypergeometric
+    term $g_n$ such that $g_{n+1} - g_n = f_n$.
+
+    **Examples**
+
+    >>> from sympy.concrete.gosper import gosper_term
+    >>> from sympy.functions import factorial
+    >>> from sympy.abc import n
+
+    >>> gosper_term((4*n + 1)*factorial(n)/factorial(2*n + 1), n)
+    (-1/2 - n)/(1/4 + n)
+
+    """
+    r = hypersimp(f, n)
+
+    if r is None:
+        return None    # 'f' is *not* a hypergeometric term
+
+    p, q = r.as_numer_denom()
 
-    nni_roots = roots(res, h, filter='Z',
-        predicate=lambda r: r >= 0).keys()
+    A, B, C = gosper_normal(p, q, n)
+    B = B.shift(-1)
 
-    if not nni_roots:
-        return (f, g, S.One)
+    N = S(A.degree())
+    M = S(B.degree())
+    K = S(C.degree())
+
+    if (N != M) or (A.LC() != B.LC()):
+        D = set([K - max(N, M)])
+    elif not N:
+        D = set([K - N + 1, S(0)])
     else:
-        for i in sorted(nni_roots):
-            d = gcd(A, B.subs(n, n+i), n)
+        D = set([K - N + 1, (B.nth(N - 1) - A.nth(N - 1))/A.LC()])
 
-            A = quo(A, d, n)
-            B = quo(B, d.subs(n, n-i), n)
+    for d in set(D):
+        if not d.is_Integer or d < 0:
+            D.remove(d)
 
-            C *= Mul(*[ d.subs(n, n-j) for j in xrange(1, i+1) ])
+    if not D:
+        return None    # 'f(n)' is *not* Gosper-summable
 
-        return (Z*A, B, C)
+    d = max(D)
 
-def gosper(term, k, a, n):
-    from sympy.solvers import rsolve_poly
+    coeffs = symbols('c:%s' % (d + 1), cls=Dummy)
+    domain = A.get_domain().inject(*coeffs)
 
-    if not term:
-        return None
+    x = Poly(coeffs, n, domain=domain)
+    H = A*x.shift(1) - B*x - C
+
+    solution = solve(H.coeffs(), coeffs)
+
+    if solution is None:
+        return None    # 'f(n)' is *not* Gosper-summable
+
+    for coeff in coeffs:
+        if coeff not in solution:
+            solution[coeff] = S(0)
+
+    x = x.as_expr().subs(solution)
+
+    if x is S.Zero:
+        return None    # 'f(n)' is *not* Gosper-summable
     else:
-        p, q = term.as_numer_denom()
-        A, B, C = normal(p, q, k)
+        return B*x/C
+
+def gosper_sum(f, k):
+    """
+    Gosper's hypergeometric summation algorithm.
+
+    Given a hypergeometric term ``f`` such that::
+
+    .. math::
 
-        B = B.subs(k, k-1)
+        s_n = \sum_{k=0}^{n-1} f_k
 
-        R = rsolve_poly([-B, A], C, k)
-        symbol = []
+    and $f(n)$ doesn't depend on $n$, returns $g_{n} - g(0)$ where
+    $g_{n+1} - g_n = f_n$, or ``None`` if $s_n$ can not be expressed
+    in closed form as a sum of hypergeometric terms.
 
-        if not (R is None  or  R is S.Zero):
-            if symbol != []:
-                symbol = symbol[0]
+    **Examples**
 
-                W = R.subs(symbol, S.Zero)
+    >>> from sympy.concrete.gosper import gosper_sum
+    >>> from sympy.functions import factorial
+    >>> from sympy.abc import n, k
 
-                if W is S.Zero:
-                    R = R.subs(symbol, S.One)
-                else:
-                    R = W
+    >>> gosper_sum((4*k + 1)*factorial(k)/factorial(2*k + 1), (k, 0, n))
+    (-n! + 2*(1 + 2*n)!)/(1 + 2*n)!
+
+    **References**
+
+    .. [1] Marko Petkovsek, Herbert S. Wilf, Doron Zeilberger, A = B,
+           AK Peters, Ltd., Wellesley, MA, USA, 1997, pp. 73--100
+
+    """
+    indefinite = False
+
+    if isinstance(k, tuple):
+        k, a, b = k
+    else:
+        indefinite = True
+
+    g = gosper_term(f, k)
+
+    if g is None:
+        return None
+
+    if indefinite:
+        result = f*g
+    else:
+        result = (f*(g+1)).subs(k, b) - (f*g).subs(k, a)
 
-            Z = B*R*term/C
-            return simplify(Z.subs(k, n+1) - Z.subs(k, a))
-        else:
-            return None
+    return factor(result)
 

sympy/concrete/tests/test_gosper.py

@@ -1,8 +1,61 @@
-from sympy import Symbol, normal
-from sympy.abc import n
+"""Tests for Gosper's algorithm for hypergeometric summation. """
 
-def test_normal():
-    assert normal(4*n+5, 2*(4*n+1)*(2*n+3), n)
+from sympy.concrete.gosper import gosper_normal, gosper_term, gosper_sum
+from sympy import S, Poly, symbols, factorial, binomial, gamma
+
+n, m, k, i, j = symbols('n,m,k,i,j', integer=True)
+
+def test_gosper_normal():
+    assert gosper_normal(4*n + 5, 2*(4*n + 1)*(2*n + 3), n) == \
+        (Poly(S(1)/4, n), Poly(n + S(3)/2), Poly(n + S(1)/4))
+
+def test_gosper_term():
+    assert gosper_term((4*k + 1)*factorial(k)/factorial(2*k + 1), k) == (-k - S(1)/2)/(k + S(1)/4)
+
+def test_gosper_sum():
+    assert gosper_sum(1, (k, 0, n)) == 1 + n
+    assert gosper_sum(k, (k, 0, n)) == n*(1 + n)/2
+    assert gosper_sum(k**2, (k, 0, n)) == n*(1 + n)*(1 + 2*n)/6
+    assert gosper_sum(k**3, (k, 0, n)) == n**2*(1 + n)**2/4
+
+    assert gosper_sum(2**k, (k, 0, n)) == 2*2**n - 1
+
+    assert gosper_sum(factorial(k), (k, 0, n)) is None
+    assert gosper_sum(binomial(n, k), (k, 0, n)) is None
+
+    assert gosper_sum(factorial(k)/k**2, (k, 0, n)) is None
+    assert gosper_sum((k - 3)*factorial(k), (k, 0, n)) is None
+
+    assert gosper_sum(k*factorial(k), k) == factorial(k)
+    assert gosper_sum(k*factorial(k), (k, 0, n)) == n*factorial(n) + factorial(n) - 1
+
+    assert gosper_sum((-1)**k*binomial(n, k), (k, 0, n)) == 0
+    assert gosper_sum((-1)**k*binomial(n, k), (k, 0, m)) == -(-1)**m*(m - n)*binomial(n, m)/n
+
+    assert gosper_sum((4*k + 1)*factorial(k)/factorial(2*k + 1), (k, 0, n)) == \
+        (2*factorial(2*n + 1) - factorial(n))/factorial(2*n + 1)
+
+def test_gosper_sum_indefinite():
+    assert gosper_sum(k, k) == k*(k - 1)/2
+    assert gosper_sum(k**2, k) == k*(k - 1)*(2*k - 1)/6
+
+    assert gosper_sum(1/(k*(k + 1)), k) == -1/k
+    assert gosper_sum(-(27*k**4 + 158*k**3 + 430*k**2 + 678*k + 445)*gamma(2*k + 4)/(3*(3*k + 7)*gamma(3*k + 6)), k) == \
+        (3*k + 5)*(k**2 + 2*k + 5)*gamma(2*k + 4)/gamma(3*k + 6)
+
+def test_gosper_sum_parametric():
+    assert gosper_sum(binomial(S(1)/2, m - j + 1)*binomial(S(1)/2, m + j), (j, 1, n)) == \
+        n*(1 + m - n)*(-1 + 2*m + 2*n)*binomial(S(1)/2, 1 + m - n)*binomial(S(1)/2, m + n)/(m*(1 + 2*m))
+
+def test_gosper_sum_iterated():
+    f1 = binomial(2*k, k)/4**k
+    f2 = (1 + 2*n)*binomial(2*n, n)/4**n
+    f3 = (1 + 2*n)*(3 + 2*n)*binomial(2*n, n)/(3*4**n)
+    f4 = (1 + 2*n)*(3 + 2*n)*(5 + 2*n)*binomial(2*n, n)/(15*4**n)
+    f5 = (1 + 2*n)*(3 + 2*n)*(5 + 2*n)*(7 + 2*n)*binomial(2*n, n)/(105*4**n)
+
+    assert gosper_sum(f1, (k, 0, n)) == f2
+    assert gosper_sum(f2, (n, 0, n)) == f3
+    assert gosper_sum(f3, (n, 0, n)) == f4
+    assert gosper_sum(f4, (n, 0, n)) == f5
 
-def test_gosper():
-    pass

sympy/functions/combinatorial/factorials.py

@@ -126,12 +126,8 @@ def eval(cls, n):
 
                     return C.Integer(result)
 
-        if n.is_integer:
-            if n.is_negative:
-                return S.Zero
-        else:
-            return C.gamma(n+1)
-
+        if n.is_negative:
+            return S.Zero
 
     @classmethod    # ?
     def _eval_rewrite_as_gamma(self, arg):
@@ -419,18 +415,13 @@ def eval(cls, r, k):
 
                         return result
 
-        if k.is_integer:
-            if k.is_negative:
-                return S.Zero
-        else:
-            return C.gamma(r+1)/(C.gamma(r-k+1)*C.gamma(k+1))
-
+        if k.is_negative:
+            return S.Zero
 
     def _eval_rewrite_as_gamma(self, r, k):
         return C.gamma(r+1) / (C.gamma(r-k+1)*C.gamma(k+1))
 
     def _eval_is_integer(self):
         return self.args[0].is_integer and self.args[1].is_integer
 
-
 binomial = Binomial

sympy/polys/domains/compositedomain.py

@@ -1,8 +1,16 @@
 """Implementation of :class:`CompositeDomain` class. """
 
 from sympy.polys.domains.domain import Domain
+from sympy.polys.polyerrors import GeneratorsError
 
 class CompositeDomain(Domain):
-    """Base class for composite domains, e.g. ZZ[x]. """
+    """Base class for composite domains, e.g. ZZ[x], ZZ(X). """
 
     is_Composite = True
+
+    def inject(self, *gens):
+        """Inject generators into this domain. """
+        if not (set(self.gens) & set(gens)):
+            return self.__class__(self.dom, *(self.gens + gens))
+        else:
+            raise GeneratorsError("common generators in %s and %s" % (self.gens, gens))