Next: QuestionsUp: Maths

# Euler Library

## 1 Linear Recurrence Closure

If $$a_n$$ and $$b_n$$ are linear recurrences, then:

1. $$s_n = a_n + b_n$$ is a linear recurrence
2. $$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(*(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) for n in range(1, 10)]
print(arr)
print([a_n(n) + b_n(n) for n in range(1, 14)])


## 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:

1. $$P^n_{i, j}$$ gives the probability of being in state $$j$$ starting at state $$i$$ after $$n$$ transitions
2. The ith entry of $$N1_n$$ gives the expected number of steps to absorption when starting from state $$i$$
3. $$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()
I = identity_matrix(ring, n)
N_inv = (I - Q)
ones_vector = vector(*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]


## 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$$

Created: 2023-10-23 Mon 18:48

Validate