Kitamasa Algorithm

masonxiong · 2026-01-28 20:06:23 · 算法·理论

Kitamasa Algorithm

上次联考集训的时候学了这个，这次直接用这个打穿了联考，但是因为忘了原理被拷打了。遂写此文加深一下印象。

简介

Kitamasa Algorithm 主要用于解决常系数齐次线性递推数列的远端求值问题。换句话说，Kitamasa 用于解决下面这个问题：

已知 0-indexed 数列 \{a_n\} 满足下面这个递推关系：

a_n=\sum_{i=1}^ka_{n-i}\cdot c_{k-i}

给定 \{c_k\} 和 \{a_n\} 的前 k 项，求 a_n。数据范围满足 1\le k\le 2\cdot 10^3,1\le n\le 2\cdot 10^9。

一种广为人知的解法是矩阵快速幂，但是这种解法的时间复杂度是 \mathbf O(k^3\log n) 的，不够优秀。使用 Kitamasa 可以做到 \mathbf O(\operatorname{M}(k)\log n)，其中 \operatorname{M}(k) 是两个 k 次多项式相乘的时间复杂度。

• 朴素乘是 \mathbf O(k^2\log n) 的。
• 如果可能的话，使用快速数论变换可以做到 \mathbf O(k\log k\log n)。

相较于与之实现相同功能且时间复杂度一致的 Bostan-Mori 和 Fiduccia，Kitamasa 更容易理解，且实现更为简单，最重要的：Kitamasa 不依赖减法操作。下面我们将详细介绍 Kitamasa 的原理，并给出一个实现。

原理

对于任何一个 k 阶的常系数齐次线性递推数列，确定了其前 k 项的值，就确定了整个数列。换句话说，对于满足 a_n=\sum\limits_{i=1}^ka_{n-i}\cdot c_{k-i} 的数列 \{a_n\}，其中的每一项 a_n 都可以唯一表示为 \sum\limits_{i=0}^{k-1}a_i\cdot b_i。

于是我们称数列 \{b_k\} 为 a_n 这一项的表示数列，记为 f_n。那么首先我们知道 a_k 这一项的表示数列 f_k 就是数列 \{c_k\}，其次我们还知道 i\in[0,k-1],a_i 的表示数列 f_i 满足 f_{i,j}=[i=j]。

下面是 Kitamasa 的核心思路：

• 求 a_n 这一项的值，相当于求 a_n 这一项的表示数列 f_n。
• 如果我们能够通过 f_m 和 f_p 求出 f_{m+p}=\operatorname{merge}(f_m,f_p)，再加上 f_{[0,k-1]} 均已知，我们就可以使用快速幂，调用 \mathbf O(\log n) 次 \operatorname{merge} 求出 f_n。

现在问题转化为实现 \operatorname{merge}。重述一遍问题：数列 f_n 满足 a_n=\sum\limits_{i=0}^{k-1} a_i\cdot f_{n,i}。下面我们要做的是：给定 f_m,f_p，求 f_{m+p}。首先有两个等式：

感性理解上是：给定了这个数列的递推关系，那么它前后的关系应该是不变的。严格证明略去。

根据第一个等式，为了求 f_{m+p} 我们需要求 a_{m+i},i\in[0,k-1] 的表示数列 f_{m+i},i\in[0,k-1]，这共有 k 项，其中 i=0 的时候是 f_m，这是我们已知的；于是我们只需要掌握通过 f_m 求出 f_{m+1}，剩下的都可以递推求出。

解释：这是根据第二个等式将 p=1 带入的结果。

解释：将求和式 i=k-1 一项提出来。

解释：将 a_k 按照定义展开。

解释：提出左求和式的 i=0 一项。

解释：合并两个求和式。现在 a_{m+1} 的表示数列 f_{m+1} 满足：

于是你掌握了通过 f_m 递推 \mathbf O(k) 求出 f_{m+1}。连续递推，你可以 \mathbf O(k^2) 求出所有 f_{m+i},i\in[0,k-1]。根据前面的关系 a_{m+p}=\sum\limits_{i=0}^{k-1}a_{m+i}\cdot f_{p,i}，你又知道 f_p，你就得到了 f_{m+p}。具体来说：

解释：按照表示数列的定义展开 a_{m+i} 和 a_{m+p}。

解释：交换求和符号。

解释：对位拆开。以上我们 \mathbf O(k^2) 实现了 \operatorname{merge}。使用快速幂思想，我们就可以用 \mathbf O(k^2\log n) 的时间复杂度实现求 f_n，那也就求出了 a_n。下面给出一个参考实现：

auto kitamasa(std::vector<std::uint32_t> Ak, std::vector<std::uint32_t> Ck, std::size_t n) -> std::uint32_t {
    auto k = Ak.size();

    auto merge = [&](auto Fm, auto Fp) {
        std::vector<std::uint32_t> Fmp(k, 0), Fmi(Fm);
        for (std::size_t i = 0; i != k; ++i) {
            for (std::size_t j = 0; j != k; ++j)
                Fmp[j] = Fmp[j] + Fmi[j] * Fp[i];
            auto FmiBack = Fmi[k - 1];
            for (std::size_t j = k - 1; ~j; --j)
                Fmi[j] = (j ? Ck[j] * FmiBack + Fmi[j - 1] : Ck[0] * FmiBack);
        }
        return Fmp;
    };

    std::vector<std::uint32_t> Fm(k, 0), Fn(k, 0);
    for (Fm[1] = Fn[0] = 1; n; (n & 1) && (Fn = merge(Fn, Fm), 1), Fm = merge(Fm, Fm), n >>= 1);

    std::uint32_t An = 0;
    for (std::size_t i = 0; i != k; ++i)
        An = An + Ak[i] * Fn[i];
    return An;
}

推广

Kitamasa 和矩阵快速幂一样，并不要求序列项一定是“数”。更严格地说：令 \mathbb R 为一个交换幺环（或更一般地，交换半环），\mathbb M 为一个 \mathbb R-模（或 \mathbb R-半模）。系数 c_0,\dots,c_{k-1}\in\mathbb R，初值 a_0,\dots,a_{k-1}\in\mathbb M，并定义递推

其中“\oplus”是 \mathbb M 上的加法，“\otimes”表示 \mathbb R 对 \mathbb M 的标量乘法（若 \mathbb M=\mathbb R，则就是环乘法）。Kitamasa 的整个过程只会用到：

• 在 \mathbb R 中对系数做有限次的加法与乘法（用来计算表示系数 f_{n,i}\in\mathbb R）；
• 在 \mathbb M 中做有限次的加法与 \mathbb R-标量乘法（用来计算 \bigoplus\limits_{i=0}^{k-1} a_i\otimes f_{n,i}）。

因此，只要上述代数结构满足加法乘法的结合律、交换律以及分配律，并具有单位元，就可以在该结构上使用 Kitamasa 实现远端求值；并不依赖于“数值大小”“大小比较”或“除法”。

优化

我们现在 \operatorname{merge} 的时间复杂度是 \mathbf O(k^2) 的。如果 \mathbb R 上有快速卷积算法，并且能够快速求出加法逆元，我们就可以用 \mathbf O(\operatorname M(k)) 的时间复杂度实现 \operatorname{merge}。

例如当 \mathbb R=\mathbb{F}_{998244353} 时就可以使用快速数论变换将 \operatorname{merge} 的时间复杂度优化到 \mathbf O(k\log k)，总体时间复杂度为 \mathbf O(k\log k\log n)，这也是常系数齐次线性递推数列的远端求值问题的最优时间复杂度。

因此，我们要找到一种方法把 \operatorname{merge} 改写为卷积的形式。首先为避免混淆，我们再次强调两层运算：

根据我们推广的定义，表示数列 f_n 的意义就变成了 a_n=\bigoplus\limits_{i=0}^{k-1} a_i\otimes f_{n,i}。我们已经在前面推出了从 f_m 得到 f_{m+1} 的递推，在推广中其形式如下：

注意，花括号右侧的 +、\cdot 都在 \mathbb R 中。现在我们把它抽象成一个映射 \operatorname T，其定义如下：

根据这个定义就有 f_{n+1}=\operatorname T(f_n)，从而推出一个非常重要的等式：f_{n+t}=\operatorname T^t(f_n)，其中 \operatorname T^t 可以理解为将 \operatorname T 连续作用 t 次。直觉上，\operatorname T 就是“把要表示的那一项往后推一位（n\mapsto n+1）”在系数空间里的对应操作。

根据前面的等式 a_{m+p}=\bigoplus\limits_{i=0}^{k-1} a_{m+i}\otimes f_{p,i}，我们把每个 a_{m+i} 用初值展开，代入并交换求和（能够交换求和是因为 \mathbb M 的 \oplus 与 \mathbb R 的标量乘法满足分配律）：

而 a_{m+p} 的表示数列满足 a_{m+p}=\bigoplus\limits_{j=0}^{k-1} a_j\otimes f_{m+p,j}。因为表示数列是唯一的，我们对位拆开即可得到 f_{m+p,j}=\sum\limits_{i=0}^{k-1} f_{m+i,j}\cdot f_{p,i}。把向量写法合起来，就是 f_{m+p}=\sum\limits_{i=0}^{k-1} f_{p,i}\cdot f_{m+i}。再用 f_{m+i}=\operatorname T^i(f_m)，我们把 \operatorname{merge} 写成了下面这个形式：

这个式子告诉我们：固定一个 p，向量 f_p 决定了一个映射 \operatorname A_p=\sum\limits_{i=0}^{k-1}f_{p,i}\operatorname T^i。也就是说，上面这个式子表达了 f_{m+p}=\operatorname A_p(f_m)，即“把下标加上 p（从 m 变到 m+p）”这件事，在系数空间里就是叠上一个映射 \operatorname A_p。

现在上面这个式子的用处就出来了，如果我们还要继续加，比如再加一个 q， 一方面按定义：

另一方面，根据前面式子对应的加法规律，我们也应当有 f_{m+p+q}=\operatorname A_{p+q}(f_m)。由于上式对任意 f_m 都成立，这迫使映射本身满足下面这个关系：

因此，\operatorname{merge} 描述的问题“如何由 f_p,f_q 求 f_{p+q}”，等价于“如何计算两个算子 \operatorname A_p,\operatorname A_q 的复合 \operatorname A_q\circ\operatorname A_p，并把结果再写回 \sum\limits_{s=0}^{k-1}r_s\operatorname T^s 的形式”。

先做第一步，我们按照定义大力复合 \operatorname A_p,\operatorname A_q：

注意到一件事情：复合结果中 \operatorname T^s 一项的系数是 \sum\limits_{i+j=s}f_{q,j}f_{p,i}，这是非常标准的卷积形式！确切地说，这是加法卷积形式，也就是我们常说的多项式乘法。

现在 s 的范围是 [0,2k-2]，我们需要做第二步，把 s 的范围缩小到 [0,k-1] 并表示 \operatorname A_q\circ\operatorname A_p。

首先，a_n 的递推公式的推广版本为 a_{n+k}=\bigoplus\limits_{j=0}^{k-1} a_{n+j}\otimes c_j。把两边都用初值展开并对位比较，会得到对任意 n，都有 f_{n+k}=\sum\limits_{j=0}^{k-1} c_j\cdot f_{n+j}。而左边 f_{n+k}=\operatorname T^k(f_n)，右边 \sum c_j f_{n+j}=\sum c_j \operatorname T^j(f_n)，因此上式等价于：

然后，我们尝试将所有 s\ge k 的 \operatorname T^s 全部用 t\in[0,k-1],\operatorname T^t 的式子转化掉。不妨令 s=k+t，由上面这个式子可知：

我们发现 t+j=s-k+j，而 j\in[0,k-1]，所以 t+j<s，我们成功地将 \operatorname T^s 一项用比它小的若干项表示出来，那我们就消掉了 \operatorname T^s！

于是第二步我们也完成了，设第一步得到 \operatorname T^s 的系数是 r_s，我们只需要从大到小（为了不重复产生高次项）枚举 s\in[k,2k-2],j\in[0,k-1] 执行 r_{s-k+j}\xleftarrow+r_s\cdot c_j 即可。

以上卷积和消高次项朴素实现仍然是 \mathbf O(k^2) 的，但是这个形式就很有优化前途了。一种参考实现如下：

auto kitamasa(std::vector<std::uint32_t> Ak, std::vector<std::uint32_t> Ck, std::size_t n) -> std::uint32_t {
    auto k = Ak.size();

    auto merge = [&](auto Fm, auto Fp) {
        std::vector<std::uint32_t> Fmp(k * 2 - 1, 0);
        for (std::size_t i = 0; i != k; ++i)
            for (std::size_t j = 0; j != k; ++j)
                Fmp[i + j] = Fmp[i + j] + Fm[i] * Fp[j];
        for (std::size_t s = k * 2 - 2; s >= k; --s)
            for (std::size_t j = 0; j != k; ++j)
                Fmp[s - k + j] = Fmp[s - k + j] + Fmp[s] * Ck[j];
        Fmp.erase(Fmp.begin() + std::ptrdiff_t(k), Fmp.end());
        return Fmp;
    };

    std::vector<std::uint32_t> Fm(k, 0), Fn(k, 0);
    for (Fm[1] = Fn[0] = 1; n; (n & 1) && (Fn = merge(Fn, Fm), 1), Fm = merge(Fm, Fm), n >>= 1);

    std::uint32_t An = 0;
    for (std::size_t i = 0; i != k; ++i)
        An = An + Ak[i] * Fn[i];
    return An;
}

加法卷积非常好优化：直接使用快速数论变换即可。下面解释如何优化“消高次项”。首先我们在推导里有等式 \operatorname T^k=\sum\limits_{j=0}^{k-1}c_j\operatorname T^j，它告诉我们：出现 \operatorname T^k 时，可以用右边的线性组合替换。

这件事和“多项式取余”完全同构：令形式变量 x 扮演 \operatorname T，考虑关系 x^k=\sum\limits_{j=0}^{k-1}c_jx^j，把 x^k 移到右边（注意这一步需要减法），得到一个多项式 P(x)=x^k-\sum\limits_{j=0}^{k-1}c_jx^j。

此时，“把卷积得到的 G(x)=\sum\limits_{s=0}^{2k-2}r_sx^s 消成次数 <k”这件事，等价于求 R(x)=G(x)\bmod P(x) 的系数。

感性理解一下：对 P(x) 取模，相当于减去若干倍的 P(x)。减掉 r_tx^t 倍的 P(x)，就相当于消掉了 G(x) 中 x^{t+k} 的一项。我们的目标是把 G(x) 所有不低于 x^k 的项全部消掉，就是将 G(x) 消到次数不超过 P(x)，这就是取模的定义：持续减去一个值，直到不超过这个值为止。

我们现在举一个例子，比如斐波那契数列：a_{n+2}=a_{n+1}+a_n，可以把它写成算符的形式：\Phi^2 a_n=\Phi a_{n}+a_n，经过移项并提出 a_n 可以获得 a_n(\Phi^2-\Phi-1)=0。

因此可以得到 \Phi^2-\Phi-1=0，但是大家不能把 \Phi 看成一个数字，他本身是一个符号，这个等式的意义是在等式的左边乘以一个任何在数列中的数，算出的结果必须为 0。

同时，它也提醒我们 \Phi^2=\Phi+1，即 \Phi^2 代表的操作与 \Phi+1 代表的操作等价。

我们再看一个例子，比如说我们想要计算斐波那契数列的第 4 项，即：

读者可以自行验证正确性。

所以我们可以拿着这个等式，甚至将这个“符号”当成一个未知数（当然它并不是一种未知数）进行大部分的代数运算。

现在考虑一般的情况，比如说一个递推数列可以写成以下的形式：

这启示我们可以将一般的递推数列转化为关于 \Phi 的等式。

Kitamasa 解决的问题是已知 a_n=\sum\limits_{i=0}^kb_ia_i（b_i 是已经算出来的系数）以及 a_{m}=\sum\limits_{i=0}^kd_ia_i（d_i 是已经算出来的系数）时，计算 a_{n+m}=\sum\limits_{i=0}^{k} w_ia_i 的 w_i 的系数是多少。

仍然考虑使用算符代换，分别将前面两个式子写成 \Phi^n=\sum\limits_{i=0}^kb_i\Phi^i 和 \Phi^m=\sum\limits_{i=0}^kd_i\Phi^i。

那么接下来一步便是魔法：

发现目标序列就是将序列 b_0,b_1,\cdots,b_k 以及 d_0,d_1,\cdots,d_k 进行卷积的结果，直接证明了 kitamasa 的第一个重要结论。

将卷积计算出来之后，设 a_{n+m}=\sum\limits_{i=0}^{2k}v_i\Phi^ia_0，对于 i > k 的 \Phi^i 可以直接代入方程 (*)：\Phi^{k+1} = \sum\limits_{i=0}^k c_i \Phi^i，达到降次效果。朴素实现仅需要从后往前模拟代入过程即可。

不难发现这个实际上与多项式取模是同构的。直接说明了 Kitamasa 的第二个重要结论。