Next: QuestionsUp: Maths

Euler Library

Table of Contents

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([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 \), 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 Factorisation in Quadratic Integer Rings

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

3.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.

3.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.

4 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 \]

5 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.

6 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()[0]
    I = identity_matrix(ring, n)
    N_inv = (I - Q)
    ones_vector = vector([1]*n)
    return N_inv.solve_right(ones_vector)[i]

7 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.

8 Diophantine Equations

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

8.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]

9 Nim

The Sprague-Grundy Theorem states that every impartial game is equivalent to a game of Nim with a single heap.

10 Partisan Game Theory

We may represent a partisan game as \( G = \{ G^L \left| G^R \} \) where the \( G^L, G^R \) represent the set of games a player can move to. We may assign a number to \( G \) according to the move advantage for either player; though only if there exist no \( x \in G^L, y \in G^R \) such that \( y < x \). A positive value represents an advantage for the left player (\( L \)), and negative an advantage for the right (\( R \)) (with zero representing a win for whoever moves second, ie no moves). Examples (for \( n \ge 0 \)):

\begin{align*} \{ \ | \ \} &= 0 \text{ (no moves for either player)}\\ \{ n | \ \} &= n + 1\text{ (L may move to game n - so +1, but R has no moves)}\\ \{ \ | -n \} &= -n - 1\text{ (R may move to game n, but L has no moves)} \end{align*}

Note there exist games for which one player is guarenteed a win, but has less than a one move advantage, for example:

\[ \{ -1 \left| 0 \} \]

The Simplicity Theorem states that \( \{ a | b \} \) is equal to the first value \( a < x < b \) which occurs first in a special construction, in particular all finite games will be dyadic.

11 Transforms

12 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 \)

Author: root

Created: 2024-03-23 Sat 11:44

Validate