Hello everyone!
In CP, linear recurrences with constant coefficients are a standard topic. Such problems are typically solved by building a transition matrix and using binary exponentiation to find the $$$N$$$-th term in $$$O(\log N)$$$ time (or by applying the Bostan-Mori algorithm for larger orders).
However, the approach must change when the coefficients are not constant and instead depend on $$$n$$$.
What is a holonomic sequence?
A sequence $$$a_n$$$ is called holonomic (or P-recursive) if it satisfies a linear recurrence relation where the coefficients are polynomials in $$$n$$$.
Formally, there exist polynomials $$$P_0(x), P_1(x), \dots, P_k(x)$$$ (where $$$P_k \neq 0$$$) such that for all $$$n$$$:
The integer $$$k$$$ is the order of the recurrence, and the maximum degree among the polynomials $$$P_i$$$ is the degree of the recurrence.
Classic Examples:
- Factorials: $$$1 \cdot a_n - n \cdot a_{n-1} = 0$$$. (Order 1, Degree 1).
- Catalan numbers: $$$(n+2) \cdot C_{n+1} - (4n+2) \cdot C_n = 0$$$. (Order 1, Degree 1).
- Derangements: $$$1 \cdot D_n - (n-1) \cdot D_{n-1} - (n-1) \cdot D_{n-2} = 0$$$. (Order 2, Degree 1).
Applications in CP
A key property of holonomic sequences is their closure under various operations. If sequences $$$A$$$ and $$$B$$$ are holonomic, then the following sequences are also holonomic:
- $$$A + B$$$
- $$$A \cdot B$$$
- Prefix sums: $$$\sum_{i=0}^n A_i$$$
This implies that for complex sums such as $$$\sum_{k=0}^n \binom{n}{k}^2 \binom{n+k}{k}^2$$$ (the Apery numbers), a polynomial recurrence is mathematically guaranteed to exist, as binomial coefficients and factorials are themselves holonomic.
Analytically, these recurrences can be found using Zeilberger's Algorithm (Creative Telescoping), which algebraically simplifies the summation.
However, implementing Zeilberger's algorithm during a contest is highly impractical. Since it is theoretically guaranteed that a recurrence exists, one can bypass the analytical derivation and instead guess the recurrence programmatically.
Automated Guessing with Gaussian Elimination
Consider a scenario where the $$$N$$$-th term of a sequence must be found in $$$O(N)$$$ time (for $$$N \le 10^7$$$), but the available direct calculation or dynamic programming approach takes $$$O(N^2)$$$. Linear algebra can be used to automatically deduce the underlying recurrence.
Assume the recurrence has a small order $$$R$$$ (e.g., $$$R=2$$$) and polynomial degree $$$D$$$ (e.g., $$$D=1$$$). The equation takes the following form:
Given that the values of $$$a_n, a_{n-1}, a_{n-2}$$$, and $$$n$$$ are known, the only unknowns are the constants $$$c_{i,j}$$$. There are exactly $$$U = (R+1) \times (D+1)$$$ unknown constants. In this example, $$$U = 3 \times 2 = 6$$$.
Reference Implementation (C++)
Below is a reference solution for the problem 506E (Mr. Kitayuta's Gift) utilizing this algorithm.
Problem 506E: Given a string $$$S$$$ ($$$|S| \le 200$$$), calculate the number of valid ways to insert exactly $$$N$$$ characters ($$$N \le 10^9$$$) such that the resulting string is a palindrome. All calculations should be modulo $$$10007$$$.
How my solution works:
- Naive DP: We write a simple $$$O(|S|^2 \cdot K)$$$ DP to compute the answers for small lengths.
- Parity Split (Crucial): Even and odd length palindromes behave like two independent sequences. Mixing them inflates the recurrence order to $$$\sim 600$$$, causing a Time Limit Exceeded (TLE) during matrix exponentiation. By exclusively computing terms for the target parity $$$(|S| + N) \pmod 2$$$, we halve the recurrence order to $$$\le 305$$$.
- Automated Guessing: We generate $$$K = 800$$$ terms and feed them into the Gaussian elimination algorithm with polynomial degree $$$D=0$$$ (which effectively finds a constant-coefficient recurrence).
- Optimized Matrix Exponentiation: We construct the transition matrix and raise it to the power of $$$N/2$$$. Since $$$MOD = 10007$$$, we can safely omit the modulo operation in the innermost loop of the matrix multiplication, allowing the solution to comfortably run in $$$\sim 0.3s$$$.
If there are any questions regarding this method or its implementation, please feel free to leave a comment or send a direct message.
