# Euler Library

## Table of Contents

## 1 Linear Recurrence Closure

If \( a_n \) and \( b_n \) are linear recurrences, then:

- \( s_n = a_n + b_n \) is a linear recurrence
- \( p_n = a_nb_n \) is a linear recurrence

### 1.1 Proof (1) collapsible

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:

\begin{align*} \left( \begin{array}{cccc} \lambda_1^n & n\lambda_1^n & \cdots & \lambda_n^n \\ \lambda_1^{n - 1} & (n - 1)\lambda_1^{n-1} & \cdots & \lambda_n^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ \lambda_1 & \lambda_1 & \cdots & \lambda_n \end{array} \right) \left( \begin{array}{c} c_n \\ c_{n - 1} \\ \vdots \\ c_1 \end{array} \right) = \left( \begin{array}{c} s_n \\ s_{n - 1} \\ \vdots \\ s_1 \end{array} \right) \end{align*}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.

## 2 Kitamasa Method for Linear Recurrences

### 2.1 Polynomials

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:

\begin{equation} x^n = Q(x)f(x) + R(x) \end{equation}Where \( R(x) \) is of degree \( < m \). Substituting \( x=\alpha_i^n \), 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

#### 2.1.1 Example

Suppose \( b_n = 2^n + n = 2^n + n(1)^n \), and so has characteristic polynomial \( f(x) = (x - 1)^2(x - 2) = x^3 - 4x^2 + 5x - 2 \) , and so \( b_n = 4b_{n - 1} - 5b_{n - 2} + 2b_{n - 3} \), with base cases \( b_0 = 1, \ b_1 = 3, \ b_2 = 6 \). Suppose we want to find \( b_7 \). Then:

\[ x^7 = \left(x^4 + 4x^3 + 11x^2 + 26x + 57\right)f(x) + \left[120x^2 - 233x + 114\right] \]

And so we find that \( 120b_2 - 233b_1 + 114b_0 = 135 = 2^7 + 7 \).

### 2.2 ðŸ’– Matrices ðŸ’–

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

## 3 Power Modulo

If \( a \) and \( p \) are coprime then:

\[ a^x \pmod p = a^{x \pmod{\phi (p)}}\pmod p \]

We can also reduce \( \phi (p) \) to \( |\langle a \rangle | \) if necessary. More generally if no restrictions are placed on \( a \) and \( p \) then:

\[ a^x \pmod p = \begin{cases} a^{\left(x\pmod{\phi (p)}\right) + \phi(p)} \pmod p & \text{for } x \ge \phi(n) \\ a^x \pmod p & \text{for } x < \phi(p)\\ \end{cases} \]

### 3.1 Proof (first condition): collapsible

Let \( p_1^{e_1}p_2^{e_2}... \) be the prime factorisation of \( p \), and \( e_i' \) be the largest \( e \) s.t. \( p^e | a \). Consider

\[ a^{\phi(n)} \pmod{p_i^{e_i}} = (\lambda p_i^{e_i'})^{\phi(n)} \pmod{p_i^{e_i}} \]

Where \( \lambda \) is coprime to \( p \). Now if \( \phi(n) > e_i \) then

\[ (\lambda p_i^{e_i'})^{\phi(n)} \equiv 0 \pmod{p_i^{e_i}} \]

For all nonzero \( e_i' \), ie if \( a \) is \( p \)-free. We have:

\[ e_i < p^{e_i}(p - 1) = \phi(p^{e_i}) \le \phi(n) \]

Now consider the simultaneous congruences we have for \( a^{q\phi(n) + r} \), with \( q > 0 \) :

\[ a^{q\phi(n) + r} = \ m_i \pmod{p_i^{e_i}} \]

For \( p_i \) with \( e_i' > 0 \), we must have \( m_i = 0 \), and so given \( q > 0 \), and applying Euler's Theorem for \( p_i \not | a \), we arrive at the conclusion that the simultaneous congruences do not depend on \( q > 0 \). From here, applying The Chinese remainder theorem gets us our claim.

## 4 Factorisation in Quadratic Integer Rings

### 4.1 Gaussian Integers

\( 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) \).

#### 4.1.1 Proof collapsible

If \( p | x^2 + 1 \) then \( p | (x + i)(x - i) \), arguing for a contradiction we see we cannot have \( p|(x + i) \) or \( p | (x - i) \). Thus \( gcd(p, \ x + i) \) and \( gcd(p, \ x - i) \) are not units so we can find a proper divisor of \( p \) by computing either of them. Clearly \( a + bi = gcd(p, x \pm i) \) has nonzero real and imaginary parts (since \( p \) is prime in \( \mathbb{Z} \)), and so if \( (c + di)(a + bi) = p \) we must have \( c + di = a - bi \) since \( p \) is real. This completes the proof.

#### 4.1.2 Sum of Squares

This gives us a way to compute all the ways in which \( n \) can be written as the sum of squares: https://mathoverflow.net/a/319809/332110, and also implies if \( S(n) \) denotes the total number of ways an integer \( n = 2^\lambda p_1^{a_1}p_2^{a_2}...p_n^{a_n}q_1^{b_1}q_2^{b_2}...q_m^{b_m} \) (\( p_i \equiv 1 \pmod 4, q_i \equiv 3 \pmod 4 \)), can be written as the sum of two positive squares then:

\begin{equation*} S(n) = \begin{cases} \frac{1}{2}\left[\prod_{i=1}^m (2a_i + 1) - 1\right] & \text{if } b_i \text{ are all even}\\ 0 & \text{otherwise } \end{cases} \end{equation*}Which follows from the Brahmaguptaâ€“Fibonacci identity.

## 5 Euler Transform

Given we know the number of connected graph \( a_n \) satsifying some property, we can calculate the total number of graphs satsifying the property (\( b_n \)) using the Euler Transform:

\[ b_n = \frac{1}{n} \left( c_n + \sum^{n - 1}_{k=1} c_k b_{n - k} \right) \]

With

\[ c_n = \sum_{d|n} d a_d \]

## 6 Floor Sums

We can compute the sum: \[ D(n) = \sum_{k = 1}^n \left \lfloor \frac{n}{k} \right \rfloor \]

Also known as the Divisor Summatory Function, in \( O(\sqrt{n}) \) time:

def divisor_sum(n): u = int(math.sqrt(n)) # equivalent to math.isqrt(n) in python 3.8 return 2*sum(n//k for k in range(1, u + 1)) - u*u

Also known are fast methods to compute \( \sum_{x = 0}^n \left \lfloor \frac{ax + b}{c} \right \rfloor \), \( \sum_{x = 0}^n x \left \lfloor \frac{ax + b}{c} \right \rfloor \), \( \sum_{x = 0}^n \left \lfloor \frac{ax + b}{c} \right \rfloor^2 \) and even more generally \( \sum_{x = 0}^n x^{k_1} \left \lfloor \frac{ax + b}{c} \right \rfloor^{k_2} \), see this blog posts for more details.

## 7 Markov Chains

An absorbing Markov chain has the following form:

\begin{align*} P = \left( \begin{array}{cc} Q & R \\ 0 & I_n \\ \end{array} \right) \end{align*}Where \( Q \) is the matrix of transitions between non-absorbing states, and \( R \) the matrix between absorbing and non-absorbing states. We further define the fundamental matrix:

\[ N = \sum_{k = 0}^{\infty} Q^k = (I_t - Q)^{-1} \]

The following properties hold:

- \( P^n_{i, j} \) gives the probability of being in state \( j \) starting at state \( i \) after \( n \) transitions
- The ith entry of \( N1_n \) gives the expected number of steps to absorption when starting from state \( i \)
- \( P^k = \left(\begin{array}{cc}Q^k & (1 - Q^k)NR \\ 0 & I_n \\\end{array}\right)\)

If we want for example to compute the expected number of steps until absorption for some initial state, e.g. \( N1_n \), it is usually faster (and generally when using matrix-vector calculations) to skip matrix inversion and do the following:

\begin{align*} N*1_n = (I_n - Q)^{-1}*1_n \\ \Rightarrow (I_n - Q)^{-1}*1_n = v \\ \Rightarrow 1_n = (I_n - Q)v \end{align*}
Then for example `sage`

`(In - Q).solve_right(ones_vector)`

gives us \( v \):

def expected(i, Q, ring=QQ): n = Q.dimensions()[0] I = identity_matrix(ring, n) N_inv = (I - Q) ones_vector = vector([1]*n) return M.solve_right(ones_vector)[i]

## 8 Binomial Coefficients Modulo \( p \)

Lucas' Theorem states that if \( m_km_{k - 1}\cdots m_0 \) and \( n_kn_{k - 1}\cdots n_0 \) are the base \( p \) digit representations of \( m \) and \( n \) respectively, then:

\[ \binom{m}{n} \equiv \prod_{i = 0}^k \binom{m_i}{n_i} \pmod p \]

In base two, this implies that \( \binom{m}{n} \) is odd iff the bits of \( n \) are a subset of the bits of \( m \). There exist extensions to this theorem for prime powers (see https://web.archive.org/web/20170202003812/http://www.dms.umontreal.ca/~andrew/PDF/BinCoeff.pdf, https://blog.prabowodjonatan.id/posts/binomial-mod-pe/, https://blog.prabowodjonatan.id/posts/binomial-mod-pe/).

Kummer's Theorem states that \( v_p(\binom{m}{n}) \) is equal to the number of carries when \( m \) is added to \( n - m \) in base \( p \).

The number of entries in the mth row of Pascal's triangle that are not divisible by \( p \) is equal to the product over all digits \( d \) of \( m \) written in base \( p \) of \( 1 + d \). Similar formulations exist for the number of entries divisible by \( p \) in the first \( m \) rows of Pascal's triangle.

## 9 Diophantine Equations

### 9.1 Nonnegative Solutions to Linear Equations

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:

\begin{align*} F(x) &= (1 + x^{a_1} + x^{2a_1} + \cdots)(1 + x^{a_2} + x^{2a_2} + \cdots)\cdots(1 + x^{a_k} + x^{2a_k} + \cdots) \\ &= \frac{1}{1 - x^{a_1}}\frac{1}{1 - x^{a_2}}\cdots\frac{1}{1 - x^{a_k}} \end{align*}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 \).

#### 9.1.1 Python collapsible

from collections import defaultdict def recurrence(coeffs): # Usage: recurrence([3, 4, 5, 1, 9]) # return an array [a1, a2, a3...] s.t. f(n) = a1*f(n - 1) + a2*f(n - 2) + ... d = defaultdict(int, {0: 1}) for pw in coeffs: new_dct = d.copy() for k, v in d.items(): new_dct[k + pw] += -1*v d = new_dct mx = max(d) mul = -1*d[mx] // abs(d[mx]) return [d[i]*mul for i in range(mx)][::-1]

## 10 Transforms

## 11 Assorted

- \( f_{n + 1} = a_nf_n + \cdots + a_0f_{n - k} \) is totally periodic modulo \( n \) if \( a_0 \) is a unit in \( \mathbb{Z}_n \)