专栏文章

Kitamasa Algorithm

mmasonxiong
算法·理论参与者 21已保存评论 23

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
22 条
当前快照
1 份
快照标识符
@mlia5win
此快照首次捕获于
2026/02/12 01:05
上周
此快照最后确认于
2026/02/19 01:13
11 小时前
查看原文

Kitamasa Algorithm

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

简介

Kitamasa Algorithm 主要用于解决常系数齐次线性递推数列的远端求值问题。换句话说，Kitamasa 用于解决下面这个问题：
已知 00-indexed 数列 {an}\{a_n\} 满足下面这个递推关系：
an=∑i=1kan−i⋅ck−ia_n=\sum_{i=1}^ka_{n-i}\cdot c_{k-i}
给定 {ck}\{c_k\} 和 {an}\{a_n\} 的前 kk 项，求 ana_n。数据范围满足 1≤k≤2⋅103,1≤n≤2⋅1091\le k\le 2\cdot 10^3,1\le n\le 2\cdot 10^9。
一种广为人知的解法是矩阵快速幂，但是这种解法的时间复杂度是 O(k3log⁡n)\mathbf O(k^3\log n) 的，不够优秀。使用 Kitamasa 可以做到 O(M⁡(k)log⁡n)\mathbf O(\operatorname{M}(k)\log n)，其中 M⁡(k)\operatorname{M}(k) 是两个 kk 次多项式相乘的时间复杂度。
  • 朴素乘是 O(k2log⁡n)\mathbf O(k^2\log n) 的。
  • 如果可能的话，使用快速数论变换可以做到 O(klog⁡klog⁡n)\mathbf O(k\log k\log n)。
相较于与之实现相同功能且时间复杂度一致的 Bostan-Mori 和 Fiduccia，Kitamasa 更容易理解，且实现更为简单，最重要的：Kitamasa 不依赖减法操作。下面我们将详细介绍 Kitamasa 的原理，并给出一个实现。

原理

对于任何一个 kk 阶的常系数齐次线性递推数列，确定了其前 kk 项的值，就确定了整个数列。换句话说，对于满足 an=∑i=1kan−i⋅ck−ia_n=\sum\limits_{i=1}^ka_{n-i}\cdot c_{k-i} 的数列 {an}\{a_n\}，其中的每一项 ana_n 都可以唯一表示为 ∑i=0k−1ai⋅bi\sum\limits_{i=0}^{k-1}a_i\cdot b_i。
于是我们称数列 {bk}\{b_k\} 为 ana_n 这一项的表示数列，记为 fnf_n。那么首先我们知道 aka_k 这一项的表示数列 fkf_k 就是数列 {ck}\{c_k\}，其次我们还知道 i∈[0,k−1],aii\in[0,k-1],a_i 的表示数列 fif_i 满足 fi,j=[i=j]f_{i,j}=[i=j]。
下面是 Kitamasa 的核心思路：
  • 求 ana_n 这一项的值，相当于求 ana_n 这一项的表示数列 fnf_n。
  • 如果我们能够通过 fmf_m 和 fpf_p 求出 fm+p=merge⁡(fm,fp)f_{m+p}=\operatorname{merge}(f_m,f_p)，再加上 f[0,k−1]f_{[0,k-1]} 均已知，我们就可以使用快速幂，调用 O(log⁡n)\mathbf O(\log n) 次 merge⁡\operatorname{merge} 求出 fnf_n。
现在问题转化为实现 merge⁡\operatorname{merge}。重述一遍问题：数列 fnf_n 满足 an=∑i=0k−1ai⋅fn,ia_n=\sum\limits_{i=0}^{k-1} a_i\cdot f_{n,i}。下面我们要做的是：给定 fm,fpf_m,f_p，求 fm+pf_{m+p}。首先有两个等式：
am+p=∑i=0k−1am+i⋅fp,ia_{m+p}=\sum_{i=0}^{k-1}a_{m+i}\cdot f_{p,i} am+p=∑i=0k−1ap+i⋅fm,ia_{m+p}=\sum_{i=0}^{k-1}a_{p+i}\cdot f_{m,i}
感性理解上是：给定了这个数列的递推关系，那么它前后的关系应该是不变的。严格证明略去。
根据第一个等式，为了求 fm+pf_{m+p} 我们需要求 am+i,i∈[0,k−1]a_{m+i},i\in[0,k-1] 的表示数列 fm+i,i∈[0,k−1]f_{m+i},i\in[0,k-1]，这共有 kk 项，其中 i=0i=0 的时候是 fmf_m，这是我们已知的；于是我们只需要掌握通过 fmf_m 求出 fm+1f_{m+1}，剩下的都可以递推求出。
am+1=∑i=0k−1a1+i⋅fm,ia_{m+1}=\sum_{i=0}^{k-1}a_{1+i}\cdot f_{m,i}
解释：这是根据第二个等式将 p=1p=1 带入的结果。
am+1=(ak⋅fm,k−1)+(∑i=0k−2a1+i⋅fm,i)a_{m+1}=\left(a_k\cdot f_{m,k-1}\right)+\left(\sum_{i=0}^{k-2}a_{1+i}\cdot f_{m,i}\right)
解释：将求和式 i=k−1i=k-1 一项提出来。
am+1=((∑i=0k−1ai⋅ci)⋅fm,k−1)+(∑i=0k−2a1+i⋅fm,i)a_{m+1}=\left(\left(\sum_{i=0}^{k-1}a_i\cdot c_i\right)\cdot f_{m,k-1}\right)+\left(\sum_{i=0}^{k-2}a_{1+i}\cdot f_{m,i}\right)
解释：将 aka_k 按照定义展开。
am+1=((a0⋅c0⋅fm,k−1)+(∑i=1k−1ai⋅ci⋅fm,k−1))+(∑i=0k−2a1+i⋅fm,i)a_{m+1}=\left(\left(a_0\cdot c_0\cdot f_{m,k-1}\right)+\left(\sum_{i=1}^{k-1}a_i\cdot c_i\cdot f_{m,k-1}\right)\right)+\left(\sum_{i=0}^{k-2}a_{1+i}\cdot f_{m,i}\right)
解释：提出左求和式的 i=0i=0 一项。
am+1=a0⋅(c0⋅fm,k−1)+∑i=1k−1ai⋅(ci⋅fm,k−1+fm,i−1)a_{m+1}=a_0\cdot(c_0\cdot f_{m,k-1})+\sum_{i=1}^{k-1}a_i\cdot(c_i\cdot f_{m,k-1}+f_{m,i-1})
解释：合并两个求和式。现在 am+1a_{m+1} 的表示数列 fm+1f_{m+1} 满足：
fm+1,i={ci⋅fm,k−1i=0ci⋅fm,k−1+fm,i−1i≠0f_{m+1,i}=\begin{cases}c_i\cdot f_{m,k-1}&i=0\\c_i\cdot f_{m,k-1}+f_{m,i-1}&i\ne 0\end{cases}
于是你掌握了通过 fmf_m 递推 O(k)\mathbf O(k) 求出 fm+1f_{m+1}。连续递推，你可以 O(k2)\mathbf O(k^2) 求出所有 fm+i,i∈[0,k−1]f_{m+i},i\in[0,k-1]。根据前面的关系 am+p=∑i=0k−1am+i⋅fp,ia_{m+p}=\sum\limits_{i=0}^{k-1}a_{m+i}\cdot f_{p,i}，你又知道 fpf_p，你就得到了 fm+pf_{m+p}。具体来说：
∑j=0k−1aj⋅fm+p,j=∑i=0k−1(∑j=0k−1aj⋅fm+i,j)⋅fp,i\sum_{j=0}^{k-1}a_j\cdot f_{m+p,j}=\sum_{i=0}^{k-1}\left(\sum_{j=0}^{k-1}a_j\cdot f_{m+i,j}\right)\cdot f_{p,i}
解释：按照表示数列的定义展开 am+ia_{m+i} 和 am+pa_{m+p}。
∑j=0k−1aj⋅fm+p,j=∑j=0k−1aj(∑i=0k−1fm+i,j⋅fp,i)\sum_{j=0}^{k-1}a_j\cdot f_{m+p,j}=\sum_{j=0}^{k-1}a_j\left(\sum_{i=0}^{k-1}f_{m+i,j}\cdot f_{p,i}\right)
解释：交换求和符号。
fm+p,j=∑i=0k−1fm+i,j⋅fp,if_{m+p,j}=\sum_{i=0}^{k-1}f_{m+i,j}\cdot f_{p,i}
解释：对位拆开。以上我们 O(k2)\mathbf O(k^2) 实现了 merge⁡\operatorname{merge}。使用快速幂思想，我们就可以用 O(k2log⁡n)\mathbf O(k^2\log n) 的时间复杂度实现求 fnf_n，那也就求出了 ana_n。下面给出一个参考实现：
CPP
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 和矩阵快速幂一样，并不要求序列项一定是“数”。更严格地说：令 R\mathbb R 为一个交换幺环（或更一般地，交换半环），M\mathbb M 为一个 R\mathbb R-模（或 R\mathbb R-半模）。系数 c0,…,ck−1∈Rc_0,\dots,c_{k-1}\in\mathbb R，初值 a0,…,ak−1∈Ma_0,\dots,a_{k-1}\in\mathbb M，并定义递推
an=⨁i=1kan−i⊗ck−i(n≥k),a_n=\bigoplus_{i=1}^k a_{n-i}\otimes c_{k-i}\quad(n\ge k),
其中“⊕\oplus”是 M\mathbb M 上的加法，“⊗\otimes”表示 R\mathbb R 对 M\mathbb M 的标量乘法（若 M=R\mathbb M=\mathbb R，则就是环乘法）。Kitamasa 的整个过程只会用到：
  • 在 R\mathbb R 中对系数做有限次的加法与乘法（用来计算表示系数 fn,i∈Rf_{n,i}\in\mathbb R）；
  • 在 M\mathbb M 中做有限次的加法与 R\mathbb R-标量乘法（用来计算 ⨁i=0k−1ai⊗fn,i\bigoplus\limits_{i=0}^{k-1} a_i\otimes f_{n,i}）。
因此，只要上述代数结构满足加法乘法的结合律、交换律以及分配律，并具有单位元，就可以在该结构上使用 Kitamasa 实现远端求值；并不依赖于“数值大小”“大小比较”或“除法”。

优化

我们现在 merge⁡\operatorname{merge} 的时间复杂度是 O(k2)\mathbf O(k^2) 的。如果 R\mathbb R 上有快速卷积算法，并且能够快速求出加法逆元，我们就可以用 O(M⁡(k))\mathbf O(\operatorname M(k)) 的时间复杂度实现 merge⁡\operatorname{merge}。
例如当 R=F998244353\mathbb R=\mathbb{F}_{998244353} 时就可以使用快速数论变换将 merge⁡\operatorname{merge} 的时间复杂度优化到 O(klog⁡k)\mathbf O(k\log k)，总体时间复杂度为 O(klog⁡klog⁡n)\mathbf O(k\log k\log n)，这也是常系数齐次线性递推数列的远端求值问题的最优时间复杂度。
因此，我们要找到一种方法把 merge⁡\operatorname{merge} 改写为卷积的形式。首先为避免混淆，我们再次强调两层运算：
  • ai∈Ma_i\in\mathbb M，在 M\mathbb M 中我们用 ⊕\oplus 表示加法，用 ⊗\otimes 表示 R\mathbb R 对 M\mathbb M 的标量乘法；
  • fn,i,ci∈Rf_{n,i},c_i\in\mathbb R，因此接下来推导卷积时出现的“++”与“⋅\cdot”都发生在 R\mathbb R 中。
根据我们推广的定义，表示数列 fnf_n 的意义就变成了 an=⨁i=0k−1ai⊗fn,ia_n=\bigoplus\limits_{i=0}^{k-1} a_i\otimes f_{n,i}。我们已经在前面推出了从 fmf_m 得到 fm+1f_{m+1} 的递推，在推广中其形式如下：
fm+1,i={c0⋅fm,k−1i=0fm,i−1+ci⋅fm,k−11≤i≤k−1f_{m+1,i}=\begin{cases}c_0\cdot f_{m,k-1}&i=0\\f_{m,i-1}+c_i\cdot f_{m,k-1}&1\le i\le k-1\end{cases}
注意，花括号右侧的 ++、⋅\cdot 都在 R\mathbb R 中。现在我们把它抽象成一个映射 T⁡\operatorname T，其定义如下：
T⁡(u)i={c0⋅uk−1i=0ui−1+ci⋅uk−11≤i≤k−1\operatorname T(u)_i=\begin{cases}c_0\cdot u_{k-1}&i=0\\u_{i-1}+c_i\cdot u_{k-1}&1\le i\le k-1\end{cases}
根据这个定义就有 fn+1=T⁡(fn)f_{n+1}=\operatorname T(f_n)，从而推出一个非常重要的等式：fn+t=T⁡t(fn)f_{n+t}=\operatorname T^t(f_n)，其中 T⁡t\operatorname T^t 可以理解为将 T⁡\operatorname T 连续作用 tt 次。直觉上，T⁡\operatorname T 就是“把要表示的那一项往后推一位（n↦n+1n\mapsto n+1）”在系数空间里的对应操作。
根据前面的等式 am+p=⨁i=0k−1am+i⊗fp,ia_{m+p}=\bigoplus\limits_{i=0}^{k-1} a_{m+i}\otimes f_{p,i}，我们把每个 am+ia_{m+i} 用初值展开，代入并交换求和（能够交换求和是因为 M\mathbb M 的 ⊕\oplus 与 R\mathbb R 的标量乘法满足分配律）：
am+p=⨁i=0k−1(⨁j=0k−1aj⊗fm+i,j)⊗fp,i=⨁j=0k−1aj⊗(∑i=0k−1fm+i,j⋅fp,i)\begin{aligned}a_{m+p}&=\bigoplus_{i=0}^{k-1}\left(\bigoplus_{j=0}^{k-1} a_j\otimes f_{m+i,j}\right)\otimes f_{p,i}\\&=\bigoplus_{j=0}^{k-1} a_j\otimes\left(\sum_{i=0}^{k-1} f_{m+i,j}\cdot f_{p,i}\right)\end{aligned}
而 am+pa_{m+p} 的表示数列满足 am+p=⨁j=0k−1aj⊗fm+p,ja_{m+p}=\bigoplus\limits_{j=0}^{k-1} a_j\otimes f_{m+p,j}。因为表示数列是唯一的，我们对位拆开即可得到 fm+p,j=∑i=0k−1fm+i,j⋅fp,if_{m+p,j}=\sum\limits_{i=0}^{k-1} f_{m+i,j}\cdot f_{p,i}。把向量写法合起来，就是 fm+p=∑i=0k−1fp,i⋅fm+if_{m+p}=\sum\limits_{i=0}^{k-1} f_{p,i}\cdot f_{m+i}。再用 fm+i=T⁡i(fm)f_{m+i}=\operatorname T^i(f_m)，我们把 merge⁡\operatorname{merge} 写成了下面这个形式：
fm+p=(∑i=0k−1fp,i T⁡i)(fm)f_{m+p}=\left(\sum_{i=0}^{k-1} f_{p,i}\,\operatorname T^i\right)(f_m)
这个式子告诉我们：固定一个 pp，向量 fpf_p 决定了一个映射 A⁡p=∑i=0k−1fp,iT⁡i\operatorname A_p=\sum\limits_{i=0}^{k-1}f_{p,i}\operatorname T^i。也就是说，上面这个式子表达了 fm+p=A⁡p(fm)f_{m+p}=\operatorname A_p(f_m)，即“把下标加上 pp（从 mm 变到 m+pm+p）”这件事，在系数空间里就是叠上一个映射 A⁡p\operatorname A_p。
现在上面这个式子的用处就出来了，如果我们还要继续加，比如再加一个 qq， 一方面按定义：
fm+p+q=A⁡q(fm+p)=A⁡q(A⁡p(fm))=(A⁡q∘A⁡p)(fm)f_{m+p+q}=\operatorname A_q(f_{m+p})=\operatorname A_q(\operatorname A_p(f_m))=(\operatorname A_q\circ\operatorname A_p)(f_m)
另一方面，根据前面式子对应的加法规律，我们也应当有 fm+p+q=A⁡p+q(fm)f_{m+p+q}=\operatorname A_{p+q}(f_m)。由于上式对任意 fmf_m 都成立，这迫使映射本身满足下面这个关系：
A⁡p+q=A⁡q∘A⁡p\operatorname A_{p+q}=\operatorname A_q\circ\operatorname A_p
因此，merge⁡\operatorname{merge} 描述的问题“如何由 fp,fqf_p,f_q 求 fp+qf_{p+q}”，等价于“如何计算两个算子 A⁡p,A⁡q\operatorname A_p,\operatorname A_q 的复合 A⁡q∘A⁡p\operatorname A_q\circ\operatorname A_p，并把结果再写回 ∑s=0k−1rsT⁡s\sum\limits_{s=0}^{k-1}r_s\operatorname T^s 的形式”。
先做第一步，我们按照定义大力复合 A⁡p,A⁡q\operatorname A_p,\operatorname A_q：
A⁡p=∑i=0k−1fp,iT⁡i,A⁡q=∑j=0k−1fq,jT⁡j\operatorname A_p=\sum_{i=0}^{k-1}f_{p,i}\operatorname T^i,\qquad\operatorname A_q=\sum_{j=0}^{k-1}f_{q,j}\operatorname T^j A⁡q∘A⁡p=∑i=0k−1∑j=0k−1fq,jfp,iT⁡i+j=∑s=02k−2(∑i+j=sfq,jfp,i)T⁡s\operatorname A_q\circ\operatorname A_p=\sum_{i=0}^{k-1}\sum_{j=0}^{k-1}f_{q,j}f_{p,i}\operatorname T^{i+j}=\sum_{s=0}^{2k-2}\left(\sum_{i+j=s}f_{q,j}f_{p,i}\right)\operatorname T^s
注意到一件事情：复合结果中 T⁡s\operatorname T^s 一项的系数是 ∑i+j=sfq,jfp,i\sum\limits_{i+j=s}f_{q,j}f_{p,i}，这是非常标准的卷积形式！确切地说，这是加法卷积形式，也就是我们常说的多项式乘法。
现在 ss 的范围是 [0,2k−2][0,2k-2]，我们需要做第二步，把 ss 的范围缩小到 [0,k−1][0,k-1] 并表示 A⁡q∘A⁡p\operatorname A_q\circ\operatorname A_p。
首先，ana_n 的递推公式的推广版本为 an+k=⨁j=0k−1an+j⊗cja_{n+k}=\bigoplus\limits_{j=0}^{k-1} a_{n+j}\otimes c_j。把两边都用初值展开并对位比较，会得到对任意 nn，都有 fn+k=∑j=0k−1cj⋅fn+jf_{n+k}=\sum\limits_{j=0}^{k-1} c_j\cdot f_{n+j}。而左边 fn+k=T⁡k(fn)f_{n+k}=\operatorname T^k(f_n)，右边 ∑cjfn+j=∑cjT⁡j(fn)\sum c_j f_{n+j}=\sum c_j \operatorname T^j(f_n)，因此上式等价于：
T⁡k=∑j=0k−1cjT⁡j\operatorname T^k=\sum_{j=0}^{k-1} c_j\operatorname T^j
然后，我们尝试将所有 s≥ks\ge k 的 T⁡s\operatorname T^s 全部用 t∈[0,k−1],T⁡tt\in[0,k-1],\operatorname T^t 的式子转化掉。不妨令 s=k+ts=k+t，由上面这个式子可知：
T⁡s=T⁡k+t=T⁡tT⁡k=T⁡t(∑j=0k−1cjT⁡j)=∑j=0k−1cjT⁡t+j\operatorname T^s=\operatorname T^{k+t}=\operatorname T^t\operatorname T^k=\operatorname T^t\left(\sum_{j=0}^{k-1}c_j\operatorname T^j\right)=\sum_{j=0}^{k-1}c_j\operatorname T^{t+j}
我们发现 t+j=s−k+jt+j=s-k+j，而 j∈[0,k−1]j\in[0,k-1]，所以 t+j<st+j<s，我们成功地将 T⁡s\operatorname T^s 一项用比它小的若干项表示出来，那我们就消掉了 T⁡s\operatorname T^s！
于是第二步我们也完成了，设第一步得到 T⁡s\operatorname T^s 的系数是 rsr_s，我们只需要从大到小（为了不重复产生高次项）枚举 s∈[k,2k−2],j∈[0,k−1]s\in[k,2k-2],j\in[0,k-1] 执行 rs−k+j←+rs⋅cjr_{s-k+j}\xleftarrow+r_s\cdot c_j 即可。
以上卷积和消高次项朴素实现仍然是 O(k2)\mathbf O(k^2) 的，但是这个形式就很有优化前途了。一种参考实现如下：
CPP
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;
}
加法卷积非常好优化：直接使用快速数论变换即可。下面解释如何优化“消高次项”。首先我们在推导里有等式 T⁡k=∑j=0k−1cjT⁡j\operatorname T^k=\sum\limits_{j=0}^{k-1}c_j\operatorname T^j，它告诉我们：出现 T⁡k\operatorname T^k 时，可以用右边的线性组合替换。
这件事和“多项式取余”完全同构：令形式变量 xx 扮演 T⁡\operatorname T，考虑关系 xk=∑j=0k−1cjxjx^k=\sum\limits_{j=0}^{k-1}c_jx^j，把 xkx^k 移到右边（注意这一步需要减法），得到一个多项式 P(x)=xk−∑j=0k−1cjxjP(x)=x^k-\sum\limits_{j=0}^{k-1}c_jx^j。
此时，“把卷积得到的 G(x)=∑s=02k−2rsxsG(x)=\sum\limits_{s=0}^{2k-2}r_sx^s 消成次数 <k<k”这件事，等价于求 R(x)=G(x) mod P(x)R(x)=G(x)\bmod P(x) 的系数。
感性理解一下：对 P(x)P(x) 取模，相当于减去若干倍的 P(x)P(x)。减掉 rtxtr_tx^t 倍的 P(x)P(x)，就相当于消掉了 G(x)G(x) 中 xt+kx^{t+k} 的一项。我们的目标是把 G(x)G(x) 所有不低于 xkx^k 的项全部消掉，就是将 G(x)G(x) 消到次数不超过 P(x)P(x)，这就是取模的定义：持续减去一个值，直到不超过这个值为止。
kk 次多项式取模也可以用快速数论变换优化到 O(klog⁡k)\mathbf O(k\log k)，这是多项式算法内容，这里不过多赘述，感兴趣的可以去【模板】多项式除法了解相关内容。
CPP
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> P(k + 1, 0);
        for (std::size_t i = 0; i != k; ++i)
            P[i] = (Ck[i] ? Modulus - Ck[i] : 0);
        P[k] = 1;
        return polymod(Fmp, polymul(Fm, Fp));
    };

    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 + std::uint64_t(Ak[i] * Fn[i])) % Modulus;
    return An;
}

算符

上面的推导过程足足有 11K11\text{K} 字符，真是复杂过头了。@reinforest 给出了一个利用算符的简洁推导过程。
对于一个递推数列，我们定义“位移算符” Φai=ai+1\Phi a_i=a_{i+1}，即将一项变成它的下一项。我们定义 Φnai=ai+n\Phi^na_i=a_{i+n}，即连续进行 nn 次位移操作的结果。
这个算符满足很多的性质，比如：
  • a0a_0 与 axa_x 的转化：ax=Φxa0a_x=\Phi^xa_0。
  • 第一结合律：Φp(Φqax)=Φp+qax\Phi^p(\Phi^q a_x)=\Phi^{p+q}a_x。
  • 第一类分配律：Φ(ax+ay)=Φax+Φay\Phi(a_x+a_y)=\Phi a_x+\Phi a_y。
  • 第二类分配律：Φpax+Φqax=ax(Φp+Φq)\Phi^pa_x+\Phi^q a_x=a_x(\Phi^p+\Phi^q)，这个理解成一种记号即可。
  • ⋯\cdots
我们现在举一个例子，比如斐波那契数列：an+2=an+1+ana_{n+2}=a_{n+1}+a_n，可以把它写成算符的形式：Φ2an=Φan+an\Phi^2 a_n=\Phi a_{n}+a_n，经过移项并提出 ana_n 可以获得 an(Φ2−Φ−1)=0a_n(\Phi^2-\Phi-1)=0。
因此可以得到 Φ2−Φ−1=0\Phi^2-\Phi-1=0，但是大家不能把 Φ\Phi 看成一个数字，他本身是一个符号，这个等式的意义是在等式的左边乘以一个任何在数列中的数，算出的结果必须为 00。
同时，它也提醒我们 Φ2=Φ+1\Phi^2=\Phi+1，即 Φ2\Phi^2 代表的操作与 Φ+1\Phi+1 代表的操作等价。
我们再看一个例子，比如说我们想要计算斐波那契数列的第 44 项，即：
a_4 &= \Phi^4a_0 \\ &= (\Phi+1)^2a_0 \\ &= (\Phi^2+2\Phi+1)a_0 \\ &= [(\Phi+1)+2\Phi+1]a_0 \\ &= [3\Phi+2]a_0 \\ &= 3a_1+2a_0 \end{aligned}$$ 读者可以自行验证正确性。 所以我们可以拿着这个等式，甚至将这个“符号”当成一个未知数（当然它并不是一种未知数）进行大部分的代数运算。 现在考虑一般的情况，比如说一个递推数列可以写成以下的形式： $$\begin{aligned} a_n &= c_0a_{n-k-1}+c_1a_{n-k}+c_2a_{n-k+1}+\cdots+c_{k}a_{n-1} \\ a_n &= \sum_{i=0}^k c_ia_{n-k-1+i} & (\text{写成求和形式}) \\ \Phi^{k+1} a_{n-k-1} & = \sum_{i=0}^k c_i \Phi^i a_{n-k-1} & (\text{使用算符代换}) \\ \Phi^{k+1} &= \sum_{i=0}^k c_i \Phi^i (*)& (\text{消去 $a_{n-k-1}$}) \end{aligned}$$ 这启示我们可以将一般的递推数列转化为关于 $\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$。 那么接下来一步便是魔法： $$\begin{aligned} a_{n+m}&=\Phi^{n+m}a_0=(\Phi^n)(\Phi^m)a_0 \\ &=\left(\sum\limits_{i=0}^kb_i\Phi^i\right)\left(\sum\limits_{i=0}^kd_i\Phi^i\right)a_0 \end{aligned}$$ 发现目标序列就是将序列 $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 的第二个重要结论。

评论

共 23 条评论，欢迎与作者交流。

正在加载评论...