1"""
2Recurrences
3"""
4
5from sympy.core import S, sympify
6from sympy.core.compatibility import as_int, iterable
7
8def linrec(coeffs, init, n):
9    r"""
10    Evaluation of univariate linear recurrences of homogeneous type
11    having coefficients independent of the recurrence variable.
12
13    Parameters
14    ==========
15
16    coeffs : iterable
17        Coefficients of the recurrence
18    init : iterable
19        Initial values of the recurrence
20    n : Integer
21        Point of evaluation for the recurrence
22
23    Notes
24    =====
25
26    Let `y(n)` be the recurrence of given type, ``c`` be the sequence
27    of coefficients, ``b`` be the sequence of initial/base values of the
28    recurrence and ``k`` (equal to ``len(c)``) be the order of recurrence.
29    Then,
30
31    .. math :: y(n) = \begin{cases} b_n & 0 \le n < k \\
32        c_0 y(n-1) + c_1 y(n-2) + \cdots + c_{k-1} y(n-k) & n \ge k
33        \end{cases}
34
35    Let `x_0, x_1, \ldots, x_n` be a sequence and consider the transformation
36    that maps each polynomial `f(x)` to `T(f(x))` where each power `x^i` is
37    replaced by the corresponding value `x_i`. The sequence is then a solution
38    of the recurrence if and only if `T(x^i p(x)) = 0` for each `i \ge 0` where
39    `p(x) = x^k - c_0 x^(k-1) - \cdots - c_{k-1}` is the characteristic
40    polynomial.
41
42    Then `T(f(x)p(x)) = 0` for each polynomial `f(x)` (as it is a linear
43    combination of powers `x^i`). Now, if `x^n` is congruent to
44    `g(x) = a_0 x^0 + a_1 x^1 + \cdots + a_{k-1} x^{k-1}` modulo `p(x)`, then
45    `T(x^n) = x_n` is equal to
46    `T(g(x)) = a_0 x_0 + a_1 x_1 + \cdots + a_{k-1} x_{k-1}`.
47
48    Computation of `x^n`,
49    given `x^k = c_0 x^{k-1} + c_1 x^{k-2} + \cdots + c_{k-1}`
50    is performed using exponentiation by squaring (refer to [1]_) with
51    an additional reduction step performed to retain only first `k` powers
52    of `x` in the representation of `x^n`.
53
54    Examples
55    ========
56
57    >>> from sympy.discrete.recurrences import linrec
58    >>> from sympy.abc import x, y, z
59
60    >>> linrec(coeffs=[1, 1], init=[0, 1], n=10)
61    55
62
63    >>> linrec(coeffs=[1, 1], init=[x, y], n=10)
64    34*x + 55*y
65
66    >>> linrec(coeffs=[x, y], init=[0, 1], n=5)
67    x**2*y + x*(x**3 + 2*x*y) + y**2
68
69    >>> linrec(coeffs=[1, 2, 3, 0, 0, 4], init=[x, y, z], n=16)
70    13576*x + 5676*y + 2356*z
71
72    References
73    ==========
74
75    .. [1] https://en.wikipedia.org/wiki/Exponentiation_by_squaring
76    .. [2] https://en.wikipedia.org/w/index.php?title=Modular_exponentiation&section=6#Matrices
77
78    See Also
79    ========
80
81    sympy.polys.agca.extensions.ExtensionElement.__pow__
82
83    """
84
85    if not coeffs:
86        return S.Zero
87
88    if not iterable(coeffs):
89        raise TypeError("Expected a sequence of coefficients for"
90                        " the recurrence")
91
92    if not iterable(init):
93        raise TypeError("Expected a sequence of values for the initialization"
94                        " of the recurrence")
95
96    n = as_int(n)
97    if n < 0:
98        raise ValueError("Point of evaluation of recurrence must be a "
99                        "non-negative integer")
100
101    c = [sympify(arg) for arg in coeffs]
102    b = [sympify(arg) for arg in init]
103    k = len(c)
104
105    if len(b) > k:
106        raise TypeError("Count of initial values should not exceed the "
107                        "order of the recurrence")
108    else:
109        b += [S.Zero]*(k - len(b)) # remaining initial values default to zero
110
111    if n < k:
112        return b[n]
113    terms = [u*v for u, v in zip(linrec_coeffs(c, n), b)]
114    return sum(terms[:-1], terms[-1])
115
116
117def linrec_coeffs(c, n):
118    r"""
119    Compute the coefficients of n'th term in linear recursion
120    sequence defined by c.
121
122    `x^k = c_0 x^{k-1} + c_1 x^{k-2} + \cdots + c_{k-1}`.
123
124    It computes the coefficients by using binary exponentiation.
125    This function is used by `linrec` and `_eval_pow_by_cayley`.
126
127    Parameters
128    ==========
129
130    c = coefficients of the divisor polynomial
131    n = exponent of x, so dividend is x^n
132
133    """
134
135    k = len(c)
136
137    def _square_and_reduce(u, offset):
138        # squares `(u_0 + u_1 x + u_2 x^2 + \cdots + u_{k-1} x^k)` (and
139        # multiplies by `x` if offset is 1) and reduces the above result of
140        # length upto `2k` to `k` using the characteristic equation of the
141        # recurrence given by, `x^k = c_0 x^{k-1} + c_1 x^{k-2} + \cdots + c_{k-1}`
142
143        w = [S.Zero]*(2*len(u) - 1 + offset)
144        for i, p in enumerate(u):
145            for j, q in enumerate(u):
146                w[offset + i + j] += p*q
147
148        for j in range(len(w) - 1, k - 1, -1):
149            for i in range(k):
150                w[j - i - 1] += w[j]*c[i]
151
152        return w[:k]
153
154    def _final_coeffs(n):
155        # computes the final coefficient list - `cf` corresponding to the
156        # point at which recurrence is to be evalauted - `n`, such that,
157        # `y(n) = cf_0 y(k-1) + cf_1 y(k-2) + \cdots + cf_{k-1} y(0)`
158
159        if n < k:
160            return [S.Zero]*n + [S.One] + [S.Zero]*(k - n - 1)
161        else:
162            return _square_and_reduce(_final_coeffs(n // 2), n % 2)
163
164    return _final_coeffs(n)
165