## Ex - Fibonacci: Revisited Editorial by en_translator

This problem features algorithms of computing the \(N\)-th term of a sequence defined by a linear recurrence relation. As usual, we first explain related topics before the actual editorial.

### Digit DP with carries

We first introduce a DP (Dynamic Programming) strongly related to computing the \(N\)-th term of a linear recurrence relation.

Consider the following problem:

We have \(a_1\)-yen, \(a_2\)-yen, \(\dots\) \(,a_K\)-yen coins. How many ways are there to pay \(N\) yen? \((\text{mod }998244353)\)

\(K \leq 100, S = \sum_i a_i \leq 1000, N \leq 10^{18}\)

This problem can be solved fast with a digit DP with a carry. (It may be as difficult as a problem for an orange coder. See also the pseudocode below.)

Let \(c_i\) be the number of \(a_i\)-coins to use. We determine \(c_1, \dots, c_K\) from their least significant bits.
Focusing on only one bit, we only have to consider whether or not to use a coin, so a sub-sum DP like \(dp[x] \gets dp[x] + dp[x-a]\) will do.

Once we’ve determine \((x-1)\) least significant bits of \(c_1, \dots, c_K\), what we haven’t determined is the \(x\)-th and greater bits, so the sum of prices \(\text{mod }2^x\) is a constant; thus, if the \((x-1)\)-th bit does not match, such a state can be discarded.

After all, after we compute the transition for the \((x-1)\)-th bit, we can eliminate the entries whose \((x-1)\)-th bits do not match from the DP table to halve the table size. We maintain the remaining parts as a carry and process them in the next bit. In every stage, the value is halved, so the DP table size is bounded by \(S + S/2 + S/4 + \dots \lt 2S\), so it works in a total of \(\mathrm{O}(K S \log N)\) time.

```
carry <- (zero-padding array of length 2S + 1)
carry[0] <- 1
while M != 0 do
for i = 1 to K do
for k = 2S - a_i to 0 do
carry[k + a_i] <- carry[k + a_i] + carry[k]
end for
end for
for k = 0 to 2S + 1:
p = 2k + (N mod 2)
if p < 2S + 1 then
carry[k] <- carry[p]
else
carry[k] <- 0
end if
end for
M <- floor(M / 2)
return carry[0]
```

Here are some exercises to check your understanding.

## Exercises (Beware spoilers of ARC, CF Div.1)

### Interpretation using digit DP

This chapter is largely based on Ryuhei Mori’s article (Japanese). (Translator’s note: for international reader, an article in CodeForces, and the original paper, may help understanding.)

The digit DP we used to solve the example problem can be interpreted as operations against a formal power series. We describe it in detail.

First, let us express the object of counting in this example problem with a generating function.

The generating function corresponding to using zero or more \(a_i\)-yen coins is a formal power series in the form of

\[ 1 + x^{a_1} + x^{2 a_1} + \dots = \frac{1}{1 - x^{a_1}}. \]

Thus, what we want to evaluate is

\[ \lbrack x^N \rbrack \frac{1}{1 - x^{a_1}} \cdot \frac{1}{1 - x^{a_2}} \cdot \dots \cdot \frac{1}{1 - x^{a_K}}. \]

Here, since

\[ \begin{aligned} \frac{1}{1-x} &= \frac{1 + x}{1 - x^2} \\ &= \frac{(1 + x)(1 + x^2)}{1 - x^4} \\ &= \frac{(1 + x)(1 + x^2)(1 + x^4)}{1 - x^8} \\ &\vdots \\ &= \prod_{b \geq 0} (1 + x^{2^b}), \end{aligned} \]

using

\[ F(x) = (1 + x^{a_1})(1 + x^{a_2}) \dots (1 + x^{a_K}), \]

the desired value is

\[ \lbrack x^N \rbrack \prod_{k \geq 0} F(x^{2^k}), \]

which is expressed as a product of \(F(x^{2^k})\).

As you can see, the objects that can be counted with a digit DP can sometimes be written in the form of (a function of \(x^{q^0}\)) \(\times\) (a function of \(x^{q^1}\)) \(\times\) (a function of \(x^{q^2}\)) \(\times \dots\), where \(q\) is a base. We consider how to find \(\lbrack x^N \rbrack\) of this expression fast. (We now assume \(q=2\).)

First, we call a formal power series \(B(x)\) a **2-digital formal power series** if we have

\[B(x) = \prod_{k \geq 0} Q_{k}(x^{2^k})\]

for some formal power series \(Q_0(x), Q_1(x), \dots\).

Next, we consider an operation of extracting the odd-order or even-order terms of a generating function to obtain a new generating function. We define \(\mathcal{S}_i\) (\(i \in \lbrace 0, 1 \rbrace\)) by

\[\mathcal{S}_i \left(\sum_{n \geq 0} a_n x^n \right) = \sum_{m \geq 0} a_{2m + i} x^m.\]

Consider applying \(\mathcal{S}_i\) to a 2-digital formal power series. Then the following fact holds.

Suppose that 2-digital formal power series \(B_0(x), B_1(x), \dots\) can be expressed as

\[B_n(x) = \prod_{k \geq 0} Q_{n + k}(x^{2^k})\]

with formal power series \(Q_0(x), Q_1(x), \dots\). Let \(A(x)\) be a formal power series. Then, for any non-negative integers \(N\) and \(n\),

\[\lbrack x^N \rbrack A(x) B_n(x) = \lbrack x^{\lfloor N/2 \rfloor} \rbrack \mathcal{S}_{N \bmod 2}(A(x) Q_n(x)) B_{n+1}(x).\]

(Proof) Applying \(\mathcal{S}_i\) to the product of the formal power series \(A(x)\) and \(B_n(x)\) yields

\[ \begin{aligned} \mathcal{S}_i \left( A(x) B_n(x) \right) &= \mathcal{S}_i(A(x) Q_n(x)) \prod_{k \geq 1} Q_{n+k}(x^{2^{k-1}}) \\ &= \mathcal{S}_i(A(x) Q_n(x)) B_{n+1}(x). \end{aligned} \]

Also, for a formal power series \(C(x)\),

\[ \lbrack x^N \rbrack C(x) = \lbrack x^{\lfloor N/2 \rfloor} \rbrack \mathcal{S}_{N \bmod 2}(C(x)) \]

holds. Assigning \(A(x) B_n(x)\) into \(C(x)\) above ends up in

\[ \begin{aligned} \lbrack x^N \rbrack A(x) B_n(x) &= \lbrack x^{\lfloor N/2 \rfloor} \rbrack \mathcal{S}_{N \bmod 2}(A(x) B_n(x)) \\ &= \lbrack x^{\lfloor N/2 \rfloor} \rbrack \mathcal{S}_{N \bmod 2}(A(x) Q_n(x)) B_{n+1}(x). \end{aligned} \]

(Q.E.D.)

In the equation we proved above, the left hand side evaluates \([x^N]\) and the right evaluates \([x^{\lfloor N / 2 \rfloor}]\). Thus, we can obtain \(\lbrack x^N \rbrack B_0(x)\) recursively as described by the following pseudocode. (`S(i, F)`

in the code denotes \(\mathcal{S}_i(F)\).)

```
A <- 1
k <- 0
while N > 0 do
A <- S(N mod 2, A Q_k)
k <- k + 1
N <- floor(N / 2)
end while
return [x^0] A
```

The complexity is dominated by the convolutions of \(A\) and \(Q_k(x)\). Thus, if we can express \(B_n(x)\) as a product of low-degree 2-digital expressions, this algorithm works fast enough.

In the example problem, what we wanted to evaluate was

\[ \lbrack x^N \rbrack \prod_{k \geq 0} F(x^{2^k}). \]

Thus, by assigning \(Q_n \gets F(x)\) to the algorithm above, we can derive the solution of using generating functions. Since the degree of \(A\) is always at most \(S\), each convolution costs \(\mathrm{O}(S \log S)\) time. Thus, the algorithm works in a total of \(\mathrm{O}(S \log S \log N)\) time.

Also, if you stare at this solution, you will see that what we compute is identical to the algorithm of digit DP. (The 4-th line is equivalent to the sub-sum DP + halving the DP table part.)

### Application to linear recurrence relation

Such an algorithm, which is a generalization of digit DP using generating functions, can also be applied to compute the \(N\)-th term of a sequence whose general term is expressed as a linear recurrence relation.

If the general term \(a_n\) of a sequence \(a_0, a_1, a_2, \dots\) can be represented using constants \(c_1, \dots, c_d\) as

\[a_n = \sum_{i=1}^d c_i a_{n-i} (n \geq d),\]

such a sequence is called a **linear recurrence sequence**.

While the \(N\)-th term of linear recurrence sequence can be obtained in a total of \(\mathrm{O}(d^3 \log N)\) time using fast matrix exponentiation, we describe an algorithm of complexity \(\mathrm{O}(d \log d \log N)\) using generating functions.

We can assert that the generating function of a linear recurrence sequence can be expressed as a rational function (a fraction whose numerator and denominator are both polynomials). Let \(A(x) = \sum_{n \geq 0} a_n x^n\) be the generating function of \(a_n\); then, in the range of \(n \geq d\), we have

\[ \begin{aligned} \lbrack x^n \rbrack \left(1 - \sum_{i = 1}^d c_i x^i\right)A(x) &= a_n - \sum_{i=1}^d c_i a_{n - i} \\ &= 0. \end{aligned} \]

Using

\[Q(x) = 1 - \sum_{i = 1}^d c_i x^i,\]

we can see that \(A(x)Q(x)\) is a polynomial with degree at most \(d-1\). Letting

\[ P(x) = \left(\sum_{n=0}^{d-1} a_n x^n \right) Q(x) \bmod{x^d}, \]

We have \(P(x) \equiv A(x) Q(x) \pmod{x^d}\), so \(A(x)\), \(P(x)\), and \(Q(x)\) have relations of

\[A(x) = \frac{P(x)}{Q(x)}.\]

Therefore, the \(N\)-th term of the linear recurrence sequence can be obtained as the \(N\)-th coefficient of the rational expression \(\frac{P(x)}{Q(x)}\) \((\deg(P) \lt \deg(Q))\). We introduce an algorithm called Botsan-Mori algorithm that evaluates this value.

Dividing the elements \(a_0, a_1, a_2, \dots\) of a linear recurrence sequence into even- and odd-indexed elements yields:

- a sequence consisting of even-indexed elements: \(a_0, a_2, a_4, \dots\), and
- a sequence consisting of odd-indexed elements: \(a_1, a_3, a_5, \dots\).

In fact, these sequences are both linear recurrence sequences of order at most \(d\). (We check this fact later.)

Additionally, let \(A(x)\) be the generating function of \((a_n)\); then the generating functions of the two sequences are expressed as \(\mathcal{S}_0(A(x))\) and \(\mathcal{S}_1(A(x))\).

Therefore, if we can evaluate \(\mathcal{S}_i(A(x))\) from \(A(x)\) fast, we can use the relation

\[\lbrack x^N \rbrack A(x) = \lbrack x^{\lfloor N/2 \rfloor} \rbrack \mathcal{S}_{N \bmod 2} (A(x))\]

to repeatedly halve the order of the desired term.

\(\mathcal{S}_i\left(\frac{P(x)}{Q(x)}\right)\) can be computed as follows. Multiplying \(Q(-x)\) to the numerator and denominator,

\[\frac{P(x)}{Q(x)} = \frac{P(x)Q(-x)}{Q(x)Q(-x)}.\]

Since the denominator \(Q(x)Q(-x)\) is an even polynomial, we can use a polynomial \(T(x)\) to get

\[T(x^2) = Q(x)Q(-x).\]

Using this \(T\), we can see that

\[ \begin{aligned} \mathcal{S}_i\left( \frac{P(x)}{Q(x)} \right) &= \mathcal{S}_i\left( \frac{P(x)Q(-x)}{T(x^2)} \right) \\ &= \frac{\mathcal{S}_i (P(x) Q(-x))}{T(x)}. \end{aligned} \]

Thus, the rational expression corresponding to \(\mathcal{S}_i(A(x))\) can be obtained by convolving \(P(x)\) and \(Q(-x)\), and \(Q(x)\) and \(Q(-x)\), in a total of \(\mathrm{O}(d \log d)\) time.

Therefore, the \(N\)-th coefficient of a rational expression \(\frac{P(x)}{Q(x)}\) can be obtained by the following algorithm in a total of \(\mathrm{O}(d \log d \log N)\) time (where \(\deg(Q) = d\)):

```
input : P(x), Q(x), N
output : [x^N] P(x) / Q(x)
while N > 0 do
P(x) <- S(N mod 2, P(x) Q(-x))
Q(x) <- S(0, Q(x) Q(-x))
N <- floor(N / 2)
end while
return P(0) / Q(0)
```

In fact, this algorithm can be considered as a variant of the aforementioned algorithm. That is, defining \((Q_n(x))_{n \geq 0}\) as a sequence of polynomials by \(Q_0(x) = Q(x)\) and \(Q_{n+1}(x^2) = Q_n(x) Q_n(-x)\),

\[ \begin{aligned} \frac{P(x)}{Q(x)} &= \frac{P(x)Q_0(-x)}{Q_0(x)Q_0(-x)} \\ &= P(x)Q_0(-x) \times \frac{1}{Q_1(x^2)} \\ &= P(x)Q_0(-x) \times \frac{Q_1(-x^2)}{Q_1(x^2)Q_1(-x^2)} \\ &= P(x)Q_0(-x) \times Q_1(-x^2) \times \frac{1}{Q_2(x^4)} \\ &\vdots \\ &= P(x) Q_0(-x) \times \prod_{n \geq 1} Q_n(-x^{2^n}), \end{aligned} \]

so \(\frac{P(x)}{Q(x)}\) is a 2-digital formal power series. Moreover, directly applying the algorithm in the last chapter is equivalent to the algorithm above.

### Solution of the problem

Now we describe the original problem. Since \(K\)-Fibonacci sequence is a linear recurrence sequence, we can express the original problem in terms of generating functions as follows.

- Given a rational expression \(\frac{P(x)}{Q(x)}\) where \(\deg(P) \lt \deg(Q) = K\), find the sum of \(\lbrack x^m \rbrack \frac{P(x)}{Q(x)}\) over all \(m\) such that \(m \text{ AND } N = m\).

We describe two ways to solve.

One way is to modify the Botsan-Mori algorithm.

In short, the Botsan-Mori algorithm performs the following operation in ascending order of \(i\). (If you start reading here: \(\mathcal{S}_i\) extract the even/odd-order terms from a generating function to get a new generating function; that is, \(\mathcal{S}_i \left(\sum_{n \geq 0} a_n x^n \right) = \sum_{m \geq 0} a_{2m + i} x^m\).)

- If the \(i\)-th bit of \(N\) is \(0\), then \(A(x) \gets \mathcal{S}_0(A(x))\).
- If the \(i\)-th bit of \(N\) is \(1\), then \(A(x) \gets \mathcal{S}_0(A(x))\).

The “\(m\) such that \(m \text{ AND } N = m\)” part in the original problem cam be rephrased as “if the \(i\)-th bit of \(N\) is \(1\), and the \(i\)-th bit of \(m\) is \(0\) or \(1\); if the \(i\)-th bit of \(N\) is \(0\), then the \(i\)-th bit is \(0\).” Thus, we can modify the conditional branch as follows:

- If the \(i\)-th bit of \(N\) is \(0\), then \(A(x) \gets \mathcal{S}_0(A(x))\).
- If the \(i\)-th bit of \(N\) is \(1\), then \(A(x) \gets \mathcal{S}_0(A(x))\).

All that left is to do the same recursion as Botsan-Mori.

Another way is to use the algorithm introduce in “Interpretation using digit DP.” The desired value can be expressed using the generating function \(F(x)\) of the \(K\)-Fibonacci as

\[\lbrack x^N \rbrack F(x) \left(\sum_{m \text{ AND } N = m} x^m\right).\]

Since \(F(x)\) is a rational expression, it is 2-digital. The summation is also 2-digital because

\[\lbrack x^N \rbrack F(x) \left(\sum_{m \text{ AND } N = m} x^m\right).\]

Thus, the algorithm in “Interpretation using digit DP” works by using appropriate expressions for \(Q_n(x)\).

The two solutions boil down to essentially the same algorithm, as you will see by implementing. The complexity is \(\mathrm{O}(K \log K \log N)\), which is fast enough.

- Sample code (C++)

```
#include <iostream>
#include <vector>
using namespace std;
#include "atcoder/convolution.hpp"
#include "atcoder/modint.hpp"
using mint = atcoder::modint998244353;
int main() {
long long K, N;
cin >> K >> N;
vector<mint> Q(K + 1, -1);
Q[0] = 1;
vector<mint> P = atcoder::convolution(Q, vector<mint>(K, 1));
P.resize(K);
for (; N; N >>= 1) {
vector<mint> Q2{Q};
for (int i = 1; i <= K; i += 2) Q2[i] = -Q2[i];
P = atcoder::convolution(P, Q2);
Q = atcoder::convolution(Q, Q2);
vector<mint> S(K), T(K + 1);
for (int i = 0; i < K; i++) S[i] = P[i * 2] + (N % 2 ? P[i * 2 + 1] : 0);
for (int i = 0; i < K + 1; i++) T[i] = Q[i * 2];
P = S, Q = T;
}
cout << (P[0] * (Q[0].inv())).val() << "\n";
}
```

posted:

last update: