22090: Gosper algorithm and Wilf-Zeilberger certificate

rwst · rwst · commit f089c1345497 · 2017-04-13T10:04:41.000+02:00
diff --git a/src/sage/libs/pynac/pynac.pxd b/src/sage/libs/pynac/pynac.pxd
@@ -150,7 +150,12 @@ cdef extern from "sage/libs/pynac/wrap.h":
     py_object_from_numeric(GEx e)     except +
 
     # Algorithms
-    GEx g_gcd "gcd"(GEx a, GEx b)   except +
+    GEx g_gcd "gcd"(GEx a, GEx b) except +
+    GEx g_gosper_term "gosper_term"(GEx the_ex, GEx n) except +
+    GEx g_gosper_sum_definite "gosper_sum_definite"(GEx the_ex,
+            GEx n, GEx a, GEx b, int* p) except +
+    GEx g_gosper_sum_indefinite "gosper_sum_indefinite"(GEx the_ex,
+            GEx n, int* p) except +
     GEx to_gamma(GEx expr)          except +
     GEx gamma_normalize(GEx expr)   except +
     GEx g_resultant "resultant"(GEx a, GEx b, GEx v) except +
diff --git a/src/sage/symbolic/expression.pyx b/src/sage/symbolic/expression.pyx
@@ -6654,6 +6654,153 @@ cdef class Expression(CommutativeRingElement):
             sig_off()
         return new_Expression_from_GEx(self._parent, x)
 
+    def gosper_sum(self, *args):
+        """
+        Return the summation of this hypergeometric expression using
+        Gosper's algorithm.
+
+        INPUT:
+
+        - a symbolic expression that may contain rational functions,
+          powers, factorials, gamma function terms, binomial
+          coefficients, and Pochhammer symbols that are rational-linear
+          in their arguments
+
+        - the main variable and, optionally, summation limits
+
+        EXAMPLES::
+
+            sage: a,b,k,m,n = var('a b k m n')
+            sage: SR(1).gosper_sum(n)
+            n
+            sage: SR(1).gosper_sum(n,5,8)
+            4
+            sage: n.gosper_sum(n)
+            1/2*(n - 1)*n
+            sage: n.gosper_sum(n,0,5)
+            15
+            sage: n.gosper_sum(n,0,m)
+            1/2*(m + 1)*m
+            sage: n.gosper_sum(n,a,b)
+            -1/2*(a + b)*(a - b - 1)
+
+        ::
+
+            sage: (factorial(m + n)/factorial(n)).gosper_sum(n)
+            n*factorial(m + n)/((m + 1)*factorial(n))
+            sage: (binomial(m + n, n)).gosper_sum(n)
+            n*binomial(m + n, n)/(m + 1)
+            sage: (binomial(m + n, n)).gosper_sum(n, 0, a)
+            (a + m + 1)*binomial(a + m, a)/(m + 1)
+            sage: (binomial(m + n, n)).gosper_sum(n, 0, 5)
+            1/120*(m + 6)*(m + 5)*(m + 4)*(m + 3)*(m + 2)
+            sage: (rising_factorial(a,n)/rising_factorial(b,n)).gosper_sum(n)
+            (b + n - 1)*gamma(a + n)*gamma(b)/((a - b + 1)*gamma(a)*gamma(b + n))
+            sage: factorial(n).gosper_term(n)
+            Traceback (most recent call last):
+            ...
+            ValueError: expression not Gosper-summable
+        """
+        cdef Expression s, a, b
+        cdef GEx x
+        cdef int i = 1
+        if len(args) > 1:
+            n, l1, l2 = args
+            s = self.coerce_in(n)
+            a = self.coerce_in(l1)
+            b = self.coerce_in(l2)
+            sig_on()
+            try:
+                x = g_gosper_sum_definite(self._gobj, s._gobj,
+                        a._gobj, b._gobj, &i)
+            finally:
+                sig_off()
+            if i == 0:
+                raise ValueError("expression not Gosper-summable")
+            return new_Expression_from_GEx(self._parent, x)
+        else:
+            s = self.coerce_in(args[0])
+            sig_on()
+            try:
+                x = g_gosper_sum_indefinite(self._gobj, s._gobj, &i)
+            finally:
+                sig_off()
+            if i == 0:
+                raise ValueError("expression not Gosper-summable")
+            return new_Expression_from_GEx(self._parent, x)
+
+    def gosper_term(self, n):
+        """
+        Return Gosper's hypergeometric term for ``self``.
+
+        Suppose ``f``=``self`` is a hypergeometric term such that:
+
+        .. math::
+            s_n = \sum_{k=0}^{n-1} f_k
+
+        and `f_k` doesn't depend on `n`. Return a hypergeometric
+        term `g_n` such that `g_{n+1} - g_n = f_n`.
+
+        EXAMPLES::
+
+            sage: _ = var('n')
+            sage: SR(1).gosper_term(n)
+            n
+            sage: n.gosper_term(n)
+            1/2*(n^2 - n)/n
+            sage: (n*factorial(n)).gosper_term(n)
+            1/n
+            sage: factorial(n).gosper_term(n)
+            Traceback (most recent call last):
+            ...
+            ValueError: expression not Gosper-summable
+        """
+        cdef Expression s
+        cdef int i = 1
+        s = self.coerce_in(n)
+        sig_on()
+        try:
+            x = g_gosper_term(self._gobj, s._gobj)
+        except ValueError:
+            raise ValueError("expression not Gosper-summable")
+        finally:
+            sig_off()
+        return new_Expression_from_GEx(self._parent, x)
+
+    def WZ_certificate(self, n, k):
+        r"""
+        Return the Wilf-Zeilberger certificate for this hypergeometric
+        summand in ``n``, ``k``.
+
+        To prove the identity `\sum_k F(n,k)=\textrm{const}` it suffices
+        to show that `F(n+1,k)-F(n,k)=G(n,k+1)-G(n,k),` with `G=RF` and
+        `R` the WZ certificate.
+
+        EXAMPLES:
+
+        To show that `\sum_k \binom{n}{k} = 2^n` do::
+
+            sage: _ = var('k n')
+            sage: F(n,k) = binomial(n,k) / 2^n
+            sage: c = F(n,k).WZ_certificate(n,k); c
+            1/2*k/(k - n - 1)
+            sage: G(n,k) = c * F(n,k); G
+            (n, k) |--> 1/2*k*binomial(n, k)/(2^n*(k - n - 1))
+            sage: (F(n+1,k) - F(n,k) - G(n,k+1) + G(n,k)).simplify_full()
+            0
+        """
+        cdef Expression s, f
+        cdef int i = 1
+        s = self.coerce_in(k)
+        f = self.subs(n == n+1) - self
+        sig_on()
+        try:
+            x = g_gosper_term(f._gobj, s._gobj)
+        finally:
+            sig_off()
+        g = new_Expression_from_GEx(self._parent, x)
+        return (f*g / self).to_gamma().gamma_normalize().simplify_full().factor()
+
     def lcm(self, b):
         """
         Return the lcm of self and b.
diff --git a/src/sage/tests/gosper-sum.py b/src/sage/tests/gosper-sum.py
@@ -0,0 +1,216 @@
+"""
+        References
+        ==========
+
+        .. [1] Marko Petkovsek, Herbert S. Wilf, Doron Zeilberger, A = B,
+               AK Peters, Ltd., Wellesley, MA, USA, 1997, pp. 73--100
+
+TESTS::
+
+    sage: _ = var('a b k m n r')
+    sage: SR(1).gosper_sum(k)
+    k
+    sage: SR(1).gosper_sum(k,0,n)
+    n + 1
+    sage: SR(1).gosper_sum(k, a, b)
+    -a + b + 1
+    sage: a.gosper_sum(k)
+    a*k
+    sage: a.gosper_sum(k,0,n)
+    a*(n + 1)
+    sage: a.gosper_sum(k,a,b)
+    -(a - b - 1)*a
+    sage: k.gosper_sum(k)
+    1/2*(k - 1)*k
+    sage: k.gosper_sum(k,0,n)
+    1/2*(n + 1)*n
+    sage: k.gosper_sum(k, a, b)
+    -1/2*(a + b)*(a - b - 1)
+    sage: k.gosper_sum(k, a+1, b)
+    -1/2*(a + b + 1)*(a - b)
+    sage: k.gosper_sum(k, a, b+1)
+    -1/2*(a + b + 1)*(a - b - 2)
+    sage: (k^3).gosper_sum(k)
+    1/4*(k - 1)^2*k^2
+    sage: (k^3).gosper_sum(k, a, b)
+    -1/4*(a^2 + b^2 - a + b)*(a + b)*(a - b - 1)
+
+    sage: (1/k).gosper_sum(k)
+    Traceback (most recent call last):
+    ...
+    ValueError: expression not Gosper-summable
+    sage: x = (1/k/(k+1)/(k+2)/(k+3)/(k+5)/(k+7)).gosper_sum(k,a,b)
+    sage: y = sum(1/k/(k+1)/(k+2)/(k+3)/(k+5)/(k+7) for k in range(5,500))
+    sage: assert x.subs(a==5,b==499) == y
+
+The following are from A==B, p.78 ff. Since the resulting expressions
+get more and more complicated with many correct representations that
+may differ depending on future capabilities, we check correctness by
+doing random summations::
+
+    sage: def check(ex, var, val1, val2from, val2to):
+    ....:     import random
+    ....:     symb = SR.var('symb')
+    ....:     s1 = ex.gosper_sum(var, val1, symb)
+    ....:     R = random.SystemRandom
+    ....:     val2 = R().randint(val2from, val2to)
+    ....:     s2 = sum(ex.subs(var==i) for i in range(val1, val2+1))
+    ....:     assert s1.subs(symb==val2) == s2
+
+    sage: def check_unsolvable(ex, *args):
+    ....:     try:
+    ....:         SR(ex).gosper_sum(*args)
+    ....:         raise AssertionError
+    ....:     except ValueError:
+    ....:         pass
+
+    sage: check((4*n+1) * factorial(n)/factorial(2*n+1), n, 0, 10, 20)
+    sage: check_unsolvable(factorial(k), k,0,n)
+    sage: (k * factorial(k)).gosper_sum(k,1,n)
+    n*factorial(n) + factorial(n) - 1
+    sage: (k * factorial(k)).gosper_sum(k)
+    factorial(k)
+    sage: check_unsolvable(binomial(n,k), k,0,n)
+    sage: ((-1)^k*binomial(n,k)).gosper_sum(k,0,n)
+    0
+    sage: ((-1)^k*binomial(n,k)).gosper_sum(k,0,a)
+    -(-1)^a*(a - n)*binomial(n, a)/n
+    sage: (binomial(1/2,a-k+1)*binomial(1/2,a+k)).gosper_sum(k,1,b)
+    (2*a + 2*b - 1)*(a - b + 1)*b*binomial(1/2, a + b)*binomial(1/2, a - b + 1)/((2*a + 1)*a)
+    sage: t = (binomial(2*k,k)/4^k).gosper_sum(k,0,n); t
+    (2*n + 1)*binomial(2*n, n)/4^n
+    sage: t = t.gosper_sum(n,0,n); t
+    1/3*(2*n + 3)*(2*n + 1)*binomial(2*n, n)/4^n
+    sage: t = t.gosper_sum(n,0,n); t
+    1/15*(2*n + 5)*(2*n + 3)*(2*n + 1)*binomial(2*n, n)/4^n
+    sage: t = t.gosper_sum(n,0,n); t
+    1/105*(2*n + 7)*(2*n + 5)*(2*n + 3)*(2*n + 1)*binomial(2*n, n)/4^n
+    sage: (binomial(2*k+2*a,2*a)*binomial(2*k,k)/binomial(k+a,a)/4^k).gosper_sum(k,0,n)
+    (2*a + 2*n + 1)*binomial(2*a + 2*n, 2*a)*binomial(2*n, n)/(4^n*(2*a + 1)*binomial(a + n, a))
+    sage: (4^k/binomial(2*k,k)).gosper_sum(k,0,n)
+    1/3*(2*4^n*n + 2*4^n + binomial(2*n, n))/binomial(2*n, n)
+
+    # The following are from A==B, 5.7 Exercises
+    sage: for k in range(1,5): (n^k).gosper_sum(n,0,m)
+    1/2*(m + 1)*m
+    1/6*(2*m + 1)*(m + 1)*m
+    1/4*(m + 1)^2*m^2
+    1/30*(3*m^2 + 3*m - 1)*(2*m + 1)*(m + 1)*m
+    sage: for k in range(1,4): (n^k*2^n).gosper_sum(n,0,m)
+    2*2^m*m - 2*2^m + 2
+    2*2^m*m^2 - 4*2^m*m + 6*2^m - 6
+    2*2^m*m^3 - 6*2^m*m^2 + 18*2^m*m - 26*2^m + 26
+    sage: (1 / (n^2 + sqrt(5)*n - 1)).gosper_sum(n,0,m) # known bug
+    ....:                            # TODO: algebraic solutions
+    sage: check((n^4 * 4^n / binomial(2*n, n)), n, 0, 10, 20)
+    sage: check((factorial(3*n) / (factorial(n) * factorial(n+1) * factorial(n+2) * 27^n)), n,0,10,20)
+    sage: (binomial(2*n, n)^2 / (n+1) / 4^(2*n)).gosper_sum(n,0,m)
+    (2*m + 1)^2*binomial(2*m, m)^2/(4^(2*m)*(m + 1))
+    sage: (((4*n-1) * binomial(2*n, n)^2) / (2*n-1)^2 / 4^(2*n)).gosper_sum(n,0,m)
+    -binomial(2*m, m)^2/4^(2*m)
+    sage: check(n * factorial(n-1/2)^2 / factorial(n+1)^2, n,0,10,20)
+
+    sage: (n^2 * a^n).gosper_sum(n,0,m)
+    (a^2*a^m*m^2 - 2*a*a^m*m^2 - 2*a*a^m*m + a^m*m^2 + a*a^m + 2*a^m*m + a^m - a - 1)*a/(a - 1)^3
+    sage: ((n - r/2)*binomial(r, n)).gosper_sum(n,0,m)
+    1/2*(m - r)*binomial(r, m)
+    sage: x = var('x')
+    sage: (factorial(n-1)^2 / factorial(n-x) / factorial(n+x)).gosper_sum(n,1,m)
+    (m^2*factorial(m - 1)^2*factorial(x + 1)*factorial(-x + 1) + x^2*factorial(m + x)*factorial(m - x) - factorial(m + x)*factorial(m - x))/(x^2*factorial(m + x)*factorial(m - x)*factorial(x + 1)*factorial(-x + 1))
+    sage: ((n*(n+a+b)*a^n*b^n)/factorial(n+a)/factorial(n+b)).gosper_sum(n,1,m).simplify_full()
+    -(a^(m + 1)*b^(m + 1)*factorial(a - 1)*factorial(b - 1) - factorial(a + m)*factorial(b + m))/(factorial(a + m)*factorial(a - 1)*factorial(b + m)*factorial(b - 1))
+
+    sage: check_unsolvable(1/n, n,1,m)
+    sage: check_unsolvable(1/n^2, n,1,m)
+    sage: check_unsolvable(1/n^3, n,1,m)
+    sage: ((6*n + 3) / (4*n^4 + 8*n^3 + 8*n^2 + 4*n + 3)).gosper_sum(n,1,m)
+    (m + 2)*m/(2*m^2 + 4*m + 3)
+    sage: (2^n * (n^2 - 2*n - 1)/(n^2 * (n+1)^2)).gosper_sum(n,1,m)
+    -2*(m^2 - 2^m + 2*m + 1)/(m + 1)^2
+    sage: ((4^n * n^2)/((n+2) * (n+1))).gosper_sum(n,1,m)
+    2/3*(2*4^m*m - 2*4^m + m + 2)/(m + 2)
+    sage: check_unsolvable(2^n / (n+1), n,0,m-1)
+    sage: check((4*(1-n) * (n^2-2*n-1) / n^2 / (n+1)^2 / (n-2)^2 / (n-3)^2), n, 4, 10, 20)
+    sage: check(((n^4-14*n^2-24*n-9) * 2^n / n^2 / (n+1)^2 / (n+2)^2 / (n+3)^2), n, 1, 10, 20)
+
+Exercises 3 (h), (i), (j) require symbolic product support so we leave
+them out for now.
+
+::
+
+    sage: _ = var('a b k m n r')
+    sage: check_unsolvable(binomial(2*n, n) * a^n, n)
+    sage: (binomial(2*n,n)*(1/4)^n).gosper_sum(n)
+    2*(1/4)^n*n*binomial(2*n, n)
+    sage: ((k-1) / factorial(k)).gosper_sum(k)
+    -k/factorial(k)
+    sage: F(n, k) = binomial(n, k) / 2^n
+    sage: check_unsolvable(F(n, k), k)
+    sage: _ = (F(n+1,k)-F(n,k)).gosper_term(k)
+    sage: F(n,k).WZ_certificate(n,k)
+    1/2*k/(k - n - 1)
+    sage: F(n, k) = binomial(n, k)^2 / binomial(2*n, n)
+    sage: check_unsolvable(F(n, k), k)
+    sage: _ =(F(n+1, k) - F(n, k)).gosper_term(k)
+    sage: F(n,k).WZ_certificate(n,k)
+    1/2*(2*k - 3*n - 3)*k^2/((k - n - 1)^2*(2*n + 1))
+    sage: F(n, k) = binomial(n,k) * factorial(n) / factorial(k) / factorial(a-k) / factorial(a+n)
+    sage: check_unsolvable(F(n, k), k)
+    sage: _ = (F(n+1, k) - F(n, k)).gosper_term(k)
+    sage: F(n,k).WZ_certificate(n,k)
+    k^2/((a + n + 1)*(k - n - 1))
+
+    sage: (1/n/(n+1)/(n+2)/(n+5)).gosper_sum(n)
+    1/720*(55*n^5 + 550*n^4 + 1925*n^3 + 2510*n^2 - 1728)/((n + 4)*(n + 3)*(n + 2)*(n + 1)*n)
+    sage: (1/n/(n+1)/(n+2)/(n+7)).gosper_sum(n)
+    1/1050*(91*n^7 + 1911*n^6 + 15925*n^5 + 66535*n^4 + 142534*n^3 + 132104*n^2 - 54000)/((n + 6)*(n + 5)*(n + 4)*(n + 3)*(n + 2)*(n + 1)*n)
+    sage: (1/n/(n+1)/(n+2)/(n+5)/(n+7)).gosper_sum(n)
+    1/10080*(133*n^7 + 2793*n^6 + 23275*n^5 + 97755*n^4 + 213472*n^3 + 206892*n^2 - 103680)/((n + 6)*(n + 5)*(n + 4)*(n + 3)*(n + 2)*(n + 1)*n)
+
+The following are from A=B, 7.2 WZ Proofs of the hypergeometric database::
+
+    sage: _ = var('a b c i j k m n r')
+    sage: F(n,k) = factorial(n+k)*factorial(b+k)*factorial(c-n-1)*factorial(c-b-1)/factorial(c+k)/factorial(n-1)/factorial(c-n-b-1)/factorial(k+1)/factorial(b-1)
+    sage: F(n,k).WZ_certificate(n,k)
+    -(c + k)*(k + 1)/((c - n - 1)*n)
+    sage: F(n,k)=(-1)^(n+k)*factorial(2*n+c-1)*factorial(n)*factorial(n+c-1)/factorial(2*n+c-1-k)/factorial(2*n-k)/factorial(c+k-1)/factorial(k)
+    sage: F(n,k).WZ_certificate(n,k)
+    1/2*(c^2 - 2*c*k + k^2 + 7*c*n - 6*k*n + 10*n^2 + 4*c - 3*k + 10*n + 2)*(c + k - 1)*k/((c - k + 2*n + 1)*(c - k + 2*n)*(k - 2*n - 1)*(k - 2*n - 2))
+    sage: F(n,k)=factorial(a+k-1)*factorial(b+k-1)*factorial(n)*factorial(n+c-b-a-k-1)*factorial(n+c-1)/factorial(k)/factorial(n-k)/factorial(k+c-1)/factorial(n+c-a-1)/factorial(n+c-b-1)
+    sage: F(n,k).WZ_certificate(n,k)
+    -(a + b - c + k - n)*(c + k - 1)*k/((a - c - n)*(b - c - n)*(k - n - 1))
+    sage: F(n,k)=(-1)^k*binomial(n+b,n+k)*binomial(n+c,c+k)*binomial(b+c,b+k)/factorial(n+b+c)*factorial(n)*factorial(b)*factorial(c)
+    sage: F(n,k).WZ_certificate(n,k)
+    1/2*(b + k)*(c + k)/((b + c + n + 1)*(k - n - 1))
+    sage: phi(t) = factorial(a+t-1)*factorial(b+t-1)/factorial(t)/factorial(-1/2+a+b+t)
+    sage: psi(t) = factorial(t+a+b-1/2)*factorial(t)*factorial(t+2*a+2*b-1)/factorial(t+2*a-1)/factorial(t+a+b-1)/factorial(t+2*b-1)
+    sage: F(n,k) = phi(k)*phi(n-k)*psi(n)
+    sage: F(n,k).WZ_certificate(n,k)
+    (2*a + 2*b + 2*k - 1)*(2*a + 2*b - 2*k + 3*n + 2)*(a - k + n)*(b - k + n)*k/((2*a + 2*b - 2*k + 2*n + 1)*(2*a + n)*(a + b + n)*(2*b + n)*(k - n - 1))
+
+The following are also from A=B, 7 The WZ Phenomenon::
+
+    sage: F(n,k) = factorial(n-i)*factorial(n-j)*factorial(i-1)*factorial(j-1)/factorial(n-1)/factorial(k-1)/factorial(n-i-j+k)/factorial(i-k)/factorial(j-k)
+    sage: F(n,k).WZ_certificate(n,k)
+    (k - 1)/n
+    sage: F(n,k) = binomial(3*n,k)/8^n
+    sage: F(n,k).WZ_certificate(n,k)
+    1/8*(4*k^2 - 30*k*n + 63*n^2 - 22*k + 93*n + 32)*k/((k - 3*n - 1)*(k - 3*n - 2)*(k - 3*n - 3))
+
+Example 7.5.1 gets a different but correct certificate (the certificate
+in the book fails the proof)::
+
+    sage: F(n,k) = 2^(k+1)*(k+1)*factorial(2*n-k-2)*factorial(n)/factorial(n-k-1)/factorial(2*n)
+    sage: c = F(n,k).WZ_certificate(n,k); c
+    -1/2*(k - 2*n + 1)*k/((k - n)*(2*n + 1))
+    sage: G(n,k) = c*F(n,k)
+    sage: t = (F(n+1,k) - F(n,k) - G(n,k+1) + G(n,k))
+    sage: t.simplify_full().is_trivial_zero()
+    True
+    sage: c = k/2/(-1+k-n)
+    sage: GG(n,k) = c*F(n,k)
+    sage: t = (F(n+1,k) - F(n,k) - GG(n,k+1) + GG(n,k))
+    sage: t.simplify_full().is_trivial_zero()
+    False
+"""
+
