Added the ability to evaluate linear recurrences without obtaining closed form expressions

sympy/discrete/__init__.py

@@ -3,8 +3,9 @@
 Transforms - fft, ifft, ntt, intt, fwht, ifwht, fzt, ifzt, fmt, ifmt
 Convolution - convolution, convolution_fft, convolution_ntt, convolution_fwht,
     convolution_subset, covering_product, intersecting_product
-Recurrence evaluation - reval_hcc
+Recurrence - linrec
 """
 
 from .transforms import (fft, ifft, ntt, intt, fwht, ifwht)
 from .convolution import convolution
+from .recurrence import linrec

sympy/discrete/recurrence.py

@@ -0,0 +1,107 @@
+"""
+Recurrence
+"""
+from __future__ import print_function, division
+
+from sympy.core import S, Symbol, sympify
+from sympy.core.compatibility import as_int, range, iterable
+
+def linrec(coeffs, init, n):
+    """
+    Linear recurrence evaluation of homogeneous type, having coefficients
+    independent of the variable on which the recurrence is defined.
+
+    Parameters
+    ==========
+
+    coeffs : iterable
+        coefficients of recurrence
+    init : iterable
+        initial values of recurrence
+    n : Integer
+        point of evaluation of recurrence
+
+    Explanation
+    ===========
+
+    Let, y(n) be the recurrence of given type, c be the sequence
+    of coefficients, b be the sequence of intial/base values of the
+    recurrence and k (equal to len(c)) be the order of recurrence.
+    Then,
+
+    if n < k,
+        y(n) = b[n]
+    otherwise,
+        y(n) = c[0]*y(n - 1) + c[1]*y(n - 2) + ... + c[k - 1]*y(n - k)
+
+    Examples
+    ========
+
+    >>> from sympy import linrec
+    >>> from sympy.abc import x, y, z
+
+    >>> linrec(coeffs=[1, 1], init=[0, 1], n=10)
+    55
+    >>> linrec(coeffs=[1, 1], init=[x, y], n=10)
+    34*x + 55*y
+    >>> linrec(coeffs=[x, y], init=[0, 1], n=5)
+    x**2*y + x*(x**3 + 2*x*y) + y**2
+    >>> linrec(coeffs=[1, 2, 3, 0, 0, 4], init=[x, y, z], n=16)
+    13576*x + 5676*y + 2356*z
+
+    """
+
+    if not coeffs:
+        return 0
+
+    if not iterable(coeffs):
+        raise TypeError("Expected a sequence of coefficients for"
+                        " the recurrence")
+
+    if not iterable(init):
+        raise TypeError("Expected a sequence of values for the initialization"
+                        " of the recurrence")
+
+    n = as_int(n)
+    if n < 0:
+        raise ValueError("Point of evaluation of recurrence must be a "
+                        "non-negative integer")
+
+    c = [sympify(arg) for arg in coeffs]
+    b = [sympify(arg) for arg in init]
+    k = len(c)
+
+    if len(b) > k:
+        raise TypeError("Count of initial values should not exceed the "
+                        "order of the recurrence")
+    else:
+        b += [S.Zero]*(k - len(b)) # remaining initial values default to zero
+
+    def _square_and_reduce(u, offset):
+        # squares (u[0] + u[1] * x + u[2] * x**2 + ... + u[k - 1] * x**(k - 1))
+        # (multiplies by x if offset is 1) and reduces the above result of length
+        # upto 2*k to k using the characterstic equation of the recurrence,
+        # which is, x**k = c[0] * x**(k - 1) + c[1] * x**(k - 2) + ... + c[k - 1]
+
+        w = [S.Zero]*(2*len(u) - 1 + offset)
+        for i, p in enumerate(u):
+            for j, q in enumerate(u):
+                w[offset + i + j] += p*q
+
+        for j in range(len(w) - 1, k - 1, -1):
+            for i in range(k):
+                w[j - i - 1] += w[j]*c[i]
+
+        return w[:k]
+
+    def _final_coeffs(n):
+        # computes the final coefficient list - `cf` corresponding to the
+        # point at which recurrence is to be evalauted - `n`, such that,
+        # y(n) = cf[0]*y(k - 1) + cf[1]*y(k - 2) + ... + cf[k - 1]*y(0)
+
+        if n < k:
+            return [S.Zero]*n + [S.One] + [S.Zero]*(k - n - 1)
+        else:
+            return _square_and_reduce(_final_coeffs(n // 2), offset=n%2)
+
+    return b[n] if n < k else sum(u*v for u, v in zip(_final_coeffs(n), b))

sympy/discrete/tests/test_recurrence.py

@@ -0,0 +1,59 @@
+from __future__ import print_function, division
+
+from sympy import sqrt, Rational, fibonacci
+from sympy.core import S, symbols, I
+from sympy.core.compatibility import range
+from sympy.utilities.pytest import raises
+from sympy.discrete.recurrence import linrec
+
+def test_linrec():
+    assert linrec(coeffs=[1, 1], init=[1, 1], n=20) == 10946
+    assert linrec(coeffs=[1, 2, 3, 4, 5], init=[1, 1, 0, 2], n=10) == 1040
+    assert linrec(coeffs=[0, 0, 11, 13], init=[23, 27], n=25) == 59628567384
+    assert linrec(coeffs=[0, 0, 1, 1, 2], init=[1, 5, 3], n=15) == 165
+    assert linrec(coeffs=[11, 13, 15, 17], init=[1, 2, 3, 4], n=70) == \
+        56889923441670659718376223533331214868804815612050381493741233489928913241
+    assert linrec(coeffs=[0]*55 + [1, 1, 2, 3], init=[0]*50 + [1, 2, 3], n=4000) == \
+        702633573874937994980598979769135096432444135301118916539
+
+    assert linrec(coeffs=[11, 13, 15, 17], init=[1, 2, 3, 4], n=10**4)
+    assert linrec(coeffs=[11, 13, 15, 17], init=[1, 2, 3, 4], n=10**5)
+
+    assert all(linrec(coeffs=[1, 1], init=[0, 1], n=n) == fibonacci(n)
+                                                    for n in range(95, 115))
+
+    assert all(linrec(coeffs=[1, 1], init=[1, 1], n=n) == fibonacci(n + 1)
+                                                    for n in range(595, 615))
+
+    a = [S(1)/2, S(3)/4, S(5)/6, 7, S(8)/9, S(3)/5]
+    b = [1, 2, 8, S(5)/7, S(3)/7, S(2)/9, 6]
+    x, y, z = symbols('x y z')
+
+    assert linrec(coeffs=a[:5], init=b[:4], n=80) == \
+        Rational(1726244235456268979436592226626304376013002142588105090705187189,
+            1960143456748895967474334873705475211264)
+
+    assert linrec(coeffs=a[:4], init=b[:4], n=50) == \
+        Rational(368949940033050147080268092104304441, 504857282956046106624)
+
+    assert linrec(coeffs=a[3:], init=b[:3], n=35) == \
+        Rational(97409272177295731943657945116791049305244422833125109,
+            814315512679031689453125)
+
+    assert linrec(coeffs=[0]*60 + [S(2)/3, S(4)/5], init=b, n=3000) == \
+        26777668739896791448594650497024/S(48084516708184142230517578125)
+
+    raises(TypeError, lambda: linrec(coeffs=[11, 13, 15, 17], init=[1, 2, 3, 4, 5], n=1))
+    raises(TypeError, lambda: linrec(coeffs=a[:4], init=b[:5], n=10000))
+    raises(ValueError, lambda: linrec(coeffs=a[:4], init=b[:4], n=-10000))
+    raises(TypeError, lambda: linrec(x, b, n=10000))
+    raises(TypeError, lambda: linrec(a, y, n=10000))
+
+    assert linrec(coeffs=[x, y, z], init=[1, 1, 1], n=4) == \
+        x**2  + x*y + x*z + y + z
+    assert linrec(coeffs=[1, 2, 1], init=[x, y, z], n=20) == \
+        269542*x + 664575*y + 578949*z
+    assert linrec(coeffs=[0, 3, 1, 2], init=[x, y], n=30) == \
+        58516436*x + 56372788*y
+    assert linrec(coeffs=[0]*50 + [1, 2, 3], init=[x, y, z], n=1000) == \
+        11477135884896*x + 25999077948732*y + 41975630244216*z