Before Reading...

1. Sorry for my poor English skill.
2. The propose of this article is "EASY" introduction to Kitamasa, so strict proofs can be omitted.

What is Kitamasa Method?

Kitamasa Method can compute N-th term of linear recurrence ($$$A_N = c_1A_{N-1} + c_2A_{N-2} + \cdots + c_KA_{N-K} = \sum_{i=1}^{K} c_iA_{N-i}$$$) in $$$O(K^2 \log N)$$$ or $$$O(K \log K \log N)$$$ when you manipulate polynomial with FFT.

Let's transform the equation

The goal of Kitamasa is find $$$d_i$$$, such that $$$A_N = \sum_{i=1}^{K} c_iA_{N-i} = \sum_{i=0}^{K-1} d_iA_i$$$. In other word, we can compute N-th term using first K elements instead of before K elements. If we can get $$$d_i$$$ in $$$T(N, K)$$$, we can compute N-th term in $$$O(T(N, K) + K)$$$ time.

We can find $$$d_i$$$ just using $$$A_i \leftarrow (c_1A_{i-1} + c_2A_{i-2} + \cdots)$$$ . For example, let's think about $$$A_N = 2A_{N-1} + A_{N-2}$$$. ($$$c_1 = 2, c_2 = 1$$$).

• $$$A_5 = 2A_4 + A_3$$$
• $$$A_5 = 2(2A_3 + A_2) + A_3 = 5A_3 + 2A_2$$$
• $$$A_5 = 5(2A_2 + A_1) + 2A_2 = 12A_2 + 5A_1$$$
• $$$A_5 = 12(2A_1 + A_0) + 5A_1 = 29A_1 + 12A_0$$$
• Finally, when $$$N = 5$$$ then $$$d_0 = 12, d_1 = 29$$$.

Actually, this process is subtracting $$$c(A_x - c_1A_{x-1} - c_2A_{x-2} - \cdots) = 0$$$ from both sides. In this example, we subtract $$$c(A_x - 2A_{x-1} - A_{x-2}) = 0$$$.

• $$$A_5 = 2A_4 + A_3$$$, Let's subtract $$$2(A_4 - 2A_3 - A_2) = 0$$$
• $$$A_5 = 5A_3 + 2A_2$$$, Let's subtract $$$5(A_3 - 2A_2 - A_1) = 0$$$
• $$$A_5 = 12A_2 + 5A_1$$$, Let's subtract $$$12(A_2 - 2A_1 - A_0) = 0$$$
• Finally, $$$A_5 = 29A_1 + 12A_0$$$

It is the same process as finding $$$x^N \mod f(x)$$$ when $$$f(x) = x^K - c_1x^{K-1} - c_2x^{K-2} - \cdots - c_Kx^0 = x^K - \sum_{i=1}^{K}c_ix^{K-i}$$$. For example, $$$x^5 = (x^2-2x-1)(x^3+2x^2+5x+12) + (29x + 12)$$$.

From now, we just calculate $$$x^N \mod f(x)$$$.

$$$x^N$$$?

We can't handle $$$x^N$$$ when $$$N$$$ is very large ($$$N \geq 10^9$$$). But we can calculate $$$x^N$$$ from $$$x^1, x^2, x^4, x^8, \cdots$$$ and it need $$$O(\log N)$$$ multiply operations.

Because the goal of this algorithm is get $$$x^N \mod f(x)$$$, we calculate it from $$$x^1 \mod f(x), x^2 \mod f(x), x^4 \mod f(x), \cdots$$$. We should multply/divide polynomial $$$O(\log N)$$$ time, and degree of each polynomial is $$$K-1$$$.

Let's see the pseudocode of Kitamasa:

using ll = long long;
using poly = vector<ll>;
ll kitamasa(poly c, poly a, ll n){
    poly d = {1}; // result
    poly xn = {0, 1}; // xn = x^1, x^2, x^4, ...
    poly f(c.size()+1); // f(x) = x^K - \sum c_{i}x^{i}
    f.back() = 1;
    for(int i=0; i<c.size(); i++) f[i] = -c[i];
    while(n){
        if(n & 1) d = Div(Mul(d, xn), f);
        n >>= 1; xn = Div(Mul(xn, xn), f);
    }

    ll ret = 0;
    for(int i=0; i<a.size(); i++) ret += a[i] * d[i];
    return ret;
}

$$$O(K^2 \log N)$$$ Kitamasa

When we calculate multiply/divide polynomial naively, time complexity of Kitamasa is $$$O(K^2 \log N)$$$.

int Mod(ll x){
    return (x %= MOD) < 0 ? x + MOD : x;
}

poly Mul(const poly &a, const poly &b){
    poly ret(a.size() + b.size() - 1);
    for(int i=0; i<a.size(); i++) for(int j=0; j<b.size(); j++) {
        ret[i+j] = (ret[i+j] + a[i] * b[j]) % MOD;
    }
    return ret;
}

poly Div(const poly &a, const poly &b){
    poly ret(all(a));
    for(int i=ret.size()-1; i>=b.size()-1; i--) for(int j=0; j<b.size(); j++) {
        ret[i+j-b.size()+1] = Mod(ret[i+j-b.size()+1] - ret[i] * b[j]);
    }
    ret.resize(b.size()-1);
    return ret;
}

/// A_{n} = \sum c_{i}A_{n-i} = \sum d_{i}A_{i}
ll kitamasa(poly c, poly a, ll n){
    poly d = {1}; // result
    poly xn = {0, 1}; // xn = x^1, x^2, x^4, ...
    poly f(c.size()+1); // f(x) = x^K - \sum c_{i}x^{i}
    f.back() = 1;
    for(int i=0; i<c.size(); i++) f[i] = Mod(-c[i]);
    while(n){
        if(n & 1) d = Div(Mul(d, xn), f);
        n >>= 1; xn = Div(Mul(xn, xn), f);
    }

    ll ret = 0;
    for(int i=0; i<a.size(); i++) ret = Mod(ret + a[i] * d[i]);
    return ret;
}

int main(){
    // calculate N-th Fibonacci number
    poly rec = {1, 1};
    poly dp = {0, 1};
    ll N; cin >> N;
    cout << kitamasa(rec, dp, N);
}

$$$O(K \log K \log N)$$$ Kitamasa

With FFT, we can multiply/divide polynomial in $$$O(K \log K)$$$. (https://cp-algorithms.com/algebra/polynomial.html) So time complexity of Kitamasa is $$$O(K \log K \log N)$$$.

This is my solution for Codechef RNG : https://github.com/justiceHui/AlgorithmImplement/blob/master/Math/Fast-Kitamasa.cpp