Euler Library

Let \( C_a(x) \) and \( C_b(x) \) be the characteristic polynomials of \( a_n \) and \( b_n \) respectively. Then we claim

\[ C_s(x) = \frac{C_a(x)C_b(x)}{hcf\left(C_a(x), C_b(x)\right)} \]

We've selected \( C_s(x) \) so that it has the same "solution equation" as \( s_n \) minus the constant coefficients. Now consider:

Where \( c_i \) are the components of the IC vector, and \( \lambda_k \) form the solution equation of \( s_n \):

\[ s_n = \lambda_1^n + n\lambda_1^n + \cdots + \lambda_n^n \]

Since the matrix in nonsingular (it is a confluent Vandermonde matrix), its columnspace is maximal and hence there exists an IC vector which satisfies the above.

import sympy as sp
from sympy import Matrix, Poly

# The base cases
ics1 = [1, 1]
ics2 = [1, 2, 5]

# The recurrences we want to add, defined here for clarity
def a_n(n, ics=ics1):
    if n <= len(ics):
        return ics[n - 1]
    else:
        return a_n(n - 1) + 2*a_n(n - 2)

def b_n(n, ics=ics2):
    if n <= len(ics):
        return ics[n - 1]
    else:
        return 4*b_n(n - 1) - 5*b_n(n - 2) + 2*b_n(n - 3)

x = sp.Symbol("x")
# x_n = x_{n - 1} + 2*x_{n - 2}
char_poly1 = x**2 - x - 2 
# x_n = 4*x_{n - 1} - 5*x{n - 2} + 2*x{n - 3}
char_poly2 = x**3 - 4*x**2 + 5*x - 2

p = Poly(
    sp.simplify(char_poly1 * char_poly2 / sp.gcd(char_poly1, char_poly2))
)
coeffs = p.all_coeffs()
print(p)

# We create the companion matrix of our polynomial
mat = Matrix.diag([1]*(len(coeffs) - 1))
mat.row_del(-1)
mat = mat.row_insert(0, Matrix([[-c for c in coeffs[1:]]]))
print(mat)

# Finally we show n applications of our matrix gives what we expect
ics = Matrix([[a_n(n) + b_n(n)] for n in range(1, len(coeffs))][::-1])
arr = list(ics)[::-1] + [(mat**n * ics)[0] for n in range(1, 10)]
print(arr)
print([a_n(n) + b_n(n) for n in range(1, 14)])

See also https://math.stackexchange.com/questions/1348838/sum-and-product-of-linear-recurrences.

Suppose \( b_n \) is a linear recurrence with characteristic polynomial \( f(x) = x^m + ... + c \), and a set of base cases \( b_0 ... b_{m - 1} \). Then it is well known:

\[ b_n = \sum_i f_i(n)\alpha_i^n \] Where \( f_i \) are some polynomials and \( \alpha_i \) are the roots of \( f \). Suppose we want to find the nth term of this linear recurrence, consider the polynomial \( x^n \) we can use the Euclidean algorithm to write:

Where \( R(x) \) is of degree \( < m \). Substituting \( x=\alpha_i \), we see that \( \alpha_i^n = Q(\alpha_i)*0 + R(\alpha_i) \). Let's assume first of all that \( f(x) \) has no repeated roots, so that \( b_n = \sum_i c_i\alpha_i^n \) for some constants \( c_i \). Then summing over the expression for each \( \alpha_i \) (multiplied by \( c_i \)), we obtain:

\[ \sum_i c_i\alpha_i^n = \sum_i c_i R(\alpha_i) = r_{m - 1}(c_1 \alpha_1^{m - 1} + c_2 \alpha_2^{m - 1} + ...) + r_{m - 2}(c_1 \alpha_1^{m - 2} + c_2 \alpha_2^{m - 2} + ...) + ... + r_0(c_1 + c_2 + ...) \]

Where \( r_i \) are the coefficients of \( R(x) \). Equivalently:

\[ b_n = r_{m - 1}b_{m - 1} + r_{m - 2}b_{m - 2} + ... + r_0b_0 \]

And so we have written \( b_n \) as a linear combination of its base cases.

Now, for the case that \( f \) has repeated roots, and \( f_i(n) \) is not constant. First we'll expand our earlier closed form for \( b_n \) as to obtain \( \sum_{i, k} c_in^ka_i^n \). Now instead of summing over multiples of (1), we first apply the operator \( x\frac{dx}{dx} \) (ie differentiate and multiply through by \( x \)) \( k \) times:

\[ n^kx^n = \left[Q_1(x)f(x) + Q_2(x)f'(x) + ... + Q_k(x) f^{(k)}(x)\right] + r_{m - 1}(m - 1)^kx^{m - 1} + r_{m - 2}(m - 2)^kx^{m - 2} + ... + 0 \]

Since \( \alpha_i \) is a root of mulitplicity \( k \) in \( f(x) \), the bracketed expression is equal to zero when evaluated at \( x=\alpha_i \), and so we can mulitply by \( c_i \) and sum over all \( i, k \) and achieve the same result as before.

Alternate explanation: https://codeforces.com/blog/entry/61306

And so our problem reduces to finding \( x^N \pmod{f(x)} \), where \( f(x) \) is the characteristic polynomial of our recurrence.

\( a + bi \in \mathbb{Z}[i] \) is prime iff one of the components is zero and the other is an associate of a (normal) prime of the form \( 3 \pmod 4 \), or its norm is (a normal) prime. This follows from the fact that \( p \) can be written as the sum of two squares iff -1 is a quadratic residue modulo \( p \). Note by Fermat's Two Square Theorem we can always write \( p \equiv 1 \pmod 4 \) as a sum of two squares.

It follows we can find the prime factorisation of \( n \in \mathbb{Z} \) by taking its normal factorisation and further decomposing primes of the form \( 1 \pmod 4 \). Solving \( p = (a + bi)(a - bi) = a^2 + b^2 \) can be done efficiently by first solving for \( x^2 \equiv -1 \pmod p \) (using Tonelli Shanks), we then find \( a + bi = gcd(p, \ x + i) \).

Consider nonnegative solutions to \( x_1 + \cdots + x_k = n \) where \( n \) is known. Using stars and bars, we can see the total numbers of n-tuple solutions is \( \binom{n + k - 1}{k - 1} \).

A fast method exists for \( a_1x_1 + \cdots + a_kx_k = n \) also. Let \( f(n) \) be the number of solutions to the equation for \( n \), and \( F(x) = \sum_{n = 0}^{\infty} f(n)x^n \) then we see:

Which we can rearrange to obtain: \[ (1 - x^{a_1})(1 - x^{a_2})\cdots(1 - x^{a_k})F(x) = 1 \]

Expanding, we can derive a linear recurrence for \( f(n) \). For example, suppose we want solutions to :

\[ 3x + 4y + 5z + w + 9v = n \]

We obtain:

\[ F(x)\left(-x^{22} + x^{21} + x^{19} - x^{16} - x^{15} + x^{13} - x^9 + x^7 + x^6 - x^3 - x + 1\right) = 1 \]

From which we can derive the recurrence:

\[ f(n) = f(n - 1) + f(n - 3) - f(n - 6) - f(n - 7) + f(n - 9) - f(n - 13) + f(n - 15) + f(n - 16) - f(n - 19) - f(n - 21) + f(n - 22) \]

The base cases \( f(0) \cdots f(k) \) can be obtained by manually expanding \( F(x) \), but only including powers \( \le k \).

• https://www.maths.dur.ac.uk/users/herbert.gangl/LIT/ww1chap2.pdf

Finding a maximal matching can be done in \( O(nm) \) time where \( n,m \) denote the number of vertices and edges respectively.

igraph is quite fast (though see also https://cp-algorithms.com/graph/kuhn_maximum_bipartite_matching.html):

import igraph as ig

edges = [(0, 6), (1, 7), (1, 8), (2, 5), (2, 7), (2, 8), (3, 6), (3, 7), (3, 8), (4, 7)]

vs = max(e + 1 for e, _ in edges)
G = ig.Graph.Bipartite(
    [False]*vs + [True]*vs,
    edges
)
matching = G.maximum_bipartite_matching()
# Matching(<igraph.Graph object at 0x7fb71fd1c950>,[6, 7, 5, 8, -1, 2, 0, 1, 3, -1],types=[False, False, False, False, False, True, True, True, True, True])