Official

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: