BM+线性递推学习笔记

qwaszx · 2020-10-20 18:15:47 · 个人记录

BM

BM 用于构造出给定数列的最短递推式，~~常用来打表~~

设给出的数列为 \{a_0,a_1,\cdots\}，我们需要构造一个长度最短的线性递推系数 \{f_1,f_2,\cdots f_m\} 使得对于任何 n\geq m 都有

a_n=\sum_{i=1}^ma_{n-i}f_i

使用增量法构造. 一开始的时候，设第 0 个版本的递推系数为 f^{(0)}=\{\}，即空递推系数，其长度为 m_0=0.

假设现在的递推系数的版本为 k，加入一个元素 a_n，我们需要检查递推关系是否合法. 设

\Delta_k=a_n-\sum_{i=1}^{m_k}a_{n-i}f_i^{(k)}

如果 \Delta_k=0，那么无事发生，继续加入后面的元素. 否则，设 fail_k=n 表示第 k 个版本的递推系数在第 n 个位置首次不合法. 我们需要用前面的版本构造一个新的合法的版本 f^{(k+1)}. 如果当前 k=0，那么我们直接取 f^{(1)}=\{0,0,\cdots,0\}（n+1 个 0）；否则，任取一个版本 k'<k，令 t=\dfrac{\Delta_k}{\Delta_{k'}}，我们构造一个递推系数

g=\{0,0,\cdots,t,-tf^{(k')}_1,\cdots,-tf^{(k')}_{m_{k'}}\}

即 fail_k-fail_{k'}-1 个 0 拼上 \dfrac{\Delta_k}{\Delta_{k'}} 再拼上 -\dfrac{\Delta_k}{\Delta_{k'}} 倍的 f^{(k')}. 我们可以发现这个 g 满足

\sum_{i}g_ia_{n'-i}=\Delta_k[n'=n]

于是我们只要让 f^{(k+1)}=f^{(k)}+g（对位相加，不足补 0）即可.

为了让递推系数最短，我们需要取让 g 长度最小的 k'，也就是选择 m_{k'}-fail_{k'} 最小的 k'. 空间可以滚动数组做到 O(n)，时间为 O(n^2).

线性递推

求出递推式之后需要快速计算第 n 项的 a. 假设现在已经有了递推式

a_n=\sum_{i=1}^ma_{n-i}f_i

考虑递推式的特征多项式

F(r)=r^m-f_1r^{m-1}-\cdots-f_m

我们知道 a_n 最终一定可以表为 a 的前 m 项的线性组合:

a_n=\sum_{i=0}^{m-1}h_ia_i

我们有结论

H(z)=z^n\bmod F(z)

我们有不同的方式来得出这个结论，常见的一种是通过线性代数的方法. 考虑转移矩阵

A=\left[\begin{matrix}f_1&f_2&f_3&\cdots\\1&0&0&\cdots\\0&1&0&\cdots\\\vdots&\vdots&\vdots&\ddots\end{matrix}\right]

其特征多项式为

F(\lambda)=|A-I\lambda|=\left|\begin{matrix}f_1-\lambda &f_2&f_3&\cdots\\1&-\lambda &0&\cdots\\0&1&-\lambda &\cdots\\\vdots&\vdots&\vdots&\ddots\end{matrix}\right|

通过展开第一行来计算就得到

\begin{aligned} F(\lambda)&=(f_1-\lambda)(-\lambda)^{m-1}-f_2(-\lambda)^{m-2}+\cdots\\ &=(-1)^{m}(\lambda^m-f_1\lambda^{m-1}-\cdots-f_m) \end{aligned}

Cayley-Hamilton 定理指出

F(A)=0

不过我线代水平8太行，不是很想证(

于是我们计算出

H(z)=z^n\bmod f(z)

这样我们就得到

a_n=(A^nA_0)_0=\sum_{i=0}^{m-1}h_i(A^iA_0)_0=\sum_{i=0}^{m-1}h_ia_i

实现的时候由于对 f(z) 取模和对 -f(z) 取模没有区别，所以可以忽略掉 (-1)^m. 然后由于乘法代价很大，可以使用四毛子方法来优化快速幂的乘法次数. 具体地，设一个阈值 B，我们预处理出 0 到 2^B-1 的幂，然后每倍增 B 次乘上预处理的幂来补齐余数. 这样需要 \log n+O(2^B+\dfrac{\log n}{B}) 次乘法，取 2^B=\Theta(\dfrac{\log n}{\log\log n}) 即可做到 \log n+O(\dfrac{\log n}{\log\log n}) 次乘法. 不过实际效果还没测试过，大概要等我闲的没事的时候(. 多项式取模的时候可以预处理模数的逆的点值，这样就只需要一次求逆. 当然暴力取模 m^2\log n 也是跑得很快的(

代码是暴力取模+BM，带多项式取模的快速幂靠着半在线卷积的速度完全不优化就进luogu第二页了(

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
using namespace std;
const int mod=998244353;
int n,m;
vector<int>a,f;
int qmod1(int x){return x>=mod?x-mod:x;}
int qmod2(int x){return x+(x>>31&mod);}
int qpower(int a,int b)
{
    int ans=1;for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ans=1ll*ans*a%mod;
    return ans;
}
void print(const vector<int>&a)
{
    for(int i=0;i<a.size();i++)cout<<a[i]<<" ";puts("");
}
vector<int>BM(const vector<int>&a)
{
    vector<int>ans,tmp,nw;int delta,tdelta,fail,tfail;
    for(int i=0;i<a.size();i++)
    {
        delta=a[i];
        for(int j=0;j<ans.size();j++)
            delta=qmod2(delta-1ll*ans[j]*a[i-j-1]%mod);
        if(!delta)continue;
        fail=i;
        if(!ans.size()){ans.resize(i+1),tmp.resize(i);tfail=i,tdelta=delta;continue;}
        nw=ans;
        int t=1ll*delta*qpower(tdelta,mod-2)%mod;
        if(ans.size()<i-tfail+tmp.size())ans.resize(i-tfail+tmp.size());
        ans[i-tfail-1]=qmod1(ans[i-tfail-1]+t);
        for(int j=0;j<tmp.size();j++)ans[i-tfail+j]=qmod2(ans[i-tfail+j]-1ll*t*tmp[j]%mod);
        if(fail-nw.size()>tfail-tmp.size()){tdelta=delta,tfail=fail,tmp=nw;}
    }
    return ans;
}
void Mod(vector<int>&a,const vector<int>&b)
{
    for(int i=a.size()-1;i>=b.size();i--)
        for(int j=0;j<b.size();j++)
            a[i-j-1]=(a[i-j-1]+1ll*a[i]*b[j])%mod;
    a.resize(b.size());
}
vector<int>operator *(const vector<int>&a,const vector<int>&b)
{
    vector<int>c;c.resize(a.size()+b.size()-1);
    for(int i=0;i<a.size();i++)
        for(int j=0;j<b.size();j++)
            c[i+j]=(c[i+j]+1ll*a[i]*b[j])%mod;
    return c;
}
vector<int>calc(int n,const vector<int>&m)
{
    vector<int>ans;ans.push_back(1);vector<int>a;a.push_back(0),a.push_back(1);
    for(;n;n>>=1,a=a*a,Mod(a,m))if(n&1)ans=ans*a,Mod(ans,m);
    return ans;
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=0,x;i<n;i++)scanf("%d",&x),a.push_back(x);
    f=BM(a);
    for(int i=0;i<f.size();i++)printf("%d ",f[i]);puts("");
    f=calc(m,f);int ans=0;
    for(int i=0;i<f.size();i++)ans=(ans+1ll*f[i]*a[i])%mod;
    cout<<ans<<endl;
}

求一段点值

求出 a_N,a_{N+1},\cdots ,a_{N+k-1}

我们已经可以求出 a_N 用 a_0,\cdots,a_{m-1} 线性表示的结果:

a_N=\sum_{i=0}^{m-1}h_ia_i

根据前面的推导，显然也有

a_{N+t}=\sum_{i=t}^{t+m-1}h_{i-t}a_i

于是做一次减法卷积即可.

2020.10.21 upd

~~怎么我昨天写完今天就有新科技了啊~~

~~2020年，在众多毒瘤的努力下，多项式取模成为时代的眼泪~~

~~300iq:well-known in Russia~~

EItxdy

一个更加快速且简单的线性递推方法

EI 他们好像还没研究完，所以这个东西能不能加特技还不清楚

我们知道线性递推的生成函数可以写为 \dfrac{P(z)}{Q(z)}，其中

Q(z)=1-\sum_{i=1}^df_iz^i

也就是特征多项式的系数翻转.

我们需要求的就是 [z^n]\dfrac{P(z)}{Q(z)}. 我们在分子分母上同乘一个 Q(-z) 就得到 [z^n]\dfrac{P(z)Q(-z)}{Q(z)Q(-z)}. 现在注意分母，设它为 V(z)=Q(z)Q(-z)，那么我们注意到 V(z)=V(-z)，这就是说 [z^{2k+1}]V(z)=0，于是这给出分解

\frac{P(z)}{Q(z)}=\frac{E(z^2)}{U(z^2)}+z\frac{O(z^2)}{U(z^2)},U(z^2)=Q(z)Q(-z),E(z^2)+zO(z^2)=P(z)Q(-z)

那么我们可以递归计算答案:

[z^n]\frac{P(z)}{Q(z)}=\begin{cases}[z^{n/2}]\dfrac{E(z)}{U(z)}&n\bmod 2=0\\\\ [z^{(n-1)/2}]\dfrac{O(z)}{U(z)}&n\bmod 2=1 \end{cases}

朴素实现（上下都做长为 d 的多项式乘法）即可做到 (2\log n+O(1))\bold{M}(d) 复杂度，这已经不比原来的快速幂+多项式取模做法劣了.

然而对于可以做 FFT 的情况下我们可以继续优化常数. 以下用 \widehat{F}^{(2n)}_i 或 \bold{DFT}[F]_i 表示对 F 做长为 2n 的 \bold{DFT} 的第 i 项的值. 上下文明确的时候上标会省略. 下面 n 是大于 d 的最小的 2 的幂，以下不说明的话默认 \bold{DFT} 长度为 2n.

用 \bold{DFT}[Q(z)] 计算 \bold{DFT}[Q(-z)]

我们有
\bold{DFT}[Q(-z)]_i=\sum_{k}q_k(-1)^k\omega_{2n}^{ki}=\sum_{k}q_k\omega_{2n}^{(i+n)k}=\bold{DFT}[Q(z)]_{i\oplus n}
这样就只要做 4 次长为 2n 的 FFT，时间优化到 \dfrac{4}{3}M(n).
只提取奇/偶项

如果我们要提取 A 的偶次幂项系数，也就是在求 B(z)=\dfrac{A(z)+A(-z)}{2}=A_{even}(z^2) 的系数，那么很显然有 \widehat{B}_k=\dfrac{\widehat{A}_k+\widehat{A}_{k\oplus n}}{2}，并且可以知道 \widehat{B} 的前 n 项就是 A_{even} 做长为 n 的 \bold{DFT} 的结果：
(\widehat{A}_{even}^{(n)})_i=\sum_{k}A_{2k}\omega_n^{ki}=\sum_{k}B_k\omega_{2n}^{ki}=\widehat{B}^{(2n)}_i
对奇数的情况也差不多，即计算 B(z)=\dfrac{A(z)-A(-z)}{2}=zA_{odd}(z^2) 的系数，那么这时就有 \widehat{B}=\dfrac{\widehat{A}_k-\widehat{A}_{k\oplus n}}{2}，并且可以导出 \widehat{A}_{odd} 和 \widehat{B} 之间的关系：
(\widehat{A}_{odd}^{(n)})_i=\sum_{k}A_{2k+1}\omega_n^{ki}=\sum_{k}B_k\omega_{2n}^{(k-1)i}=\omega_{2n}^{-i}\widehat{B}^{(2n)}_i
因此做 \bold{IDFT} 的时候可以把长度缩短为 n，这样需要做 2 次长为 2n 的 FFT 和两次长为 n 的 FFT，相当于一次 M(n)
不做 \bold{IDFT}

我们来考虑只维护P,Q 做 \bold{DFT} 的结果. 我们在前一次递归中其实已经算出了 P 和 Q 的长为 n 的 \bold{DFT} 的结果，我们可以利用它得到长为 2n 的 \bold{DFT} 的结果:
\widehat{A}^{(2n)}_{2i}=\sum_{k}A_k\omega_{2n}^{2ik}=\sum_{k}A_k\omega_n^{ik}=\widehat{A}^{(n)}_i \widehat{A}_{2i+1}^{(2n)}=\sum_kA_k\omega_{2n}^{(2i+1)k}=\sum_kA_k\omega_{2n}^k\omega_{n}^i
于是把 \widehat{A}^{(n)} \bold{IDFT} 回去之后点乘 \omega_{2n}^k 再 \bold{DFT} 就得到了 \widehat{A}^{(2n)} 的奇数项，使用了两次长为 n 的 FFT. 显然这一步不能继续优化复杂度了，不然我们就会得到一个低于 \Theta(n\log n) 的 \bold{DFT} 方法. 于是我们就得到一个每次递归只需要 4 次长为 n 的 FFT 的方法，相当于 \dfrac{2}{3}M(n).

跑得飞快，抢到了暂时的最优解(

void doublen(int a[],int b[],int l)
{
    copy(a,a+l,b);
    Poly::IDFT(b,l);
    for(int i=0;i<l;i++)b[i]=1ll*b[i]*Poly::w[l+i]%mod;
    Poly::DFT(b,l);
    for(int i=l-1;i>=0;i--)b[i<<1|1]=b[i];
    for(int i=0;i<l;i++)b[i<<1]=a[i];
}
int p[N],q[N];
int main()
{
//    freopen("data.txt","r",stdin);
    int n=getin(),K=getin();q[0]=1;
    for(int i=1;i<=K;i++)q[i]=qmod2(-getin());
    for(int i=0;i<K;i++)p[i]=qmod2(getin());
    Poly::mul(q,p,p,K);
    int l=1;while(l<=K)l<<=1;
//    fill(p+K,p+(l<<1),0),fill(q+K+1,q+(l<<1),0);
    int isfirst=1;
    static int tmpp[N],tmpq[N];    
    static int iw[N];iw[0]=1,iw[1]=qpower(3,mod-1-((mod-1)>>t));
    for(int i=2;i<l;i++)iw[i]=1ll*iw[i-1]*iw[1]%mod;//我嫌麻烦直接在这要把单位根的逆算出来了(
    while(n>=K)
    {
        if(isfirst)//注意第一次递归的时候要特判
        {
            isfirst=0;
            copy(p,p+K,tmpp),fill(tmpp+K,tmpp+(l<<1),0);
            copy(q,q+K+1,tmpq),fill(tmpq+K+1,tmpq+(l<<1),0);
            Poly::DFT(tmpp,l<<1),Poly::DFT(tmpq,l<<1);
        }
        else
        {
            doublen(p,tmpp,l),doublen(q,tmpq,l);
        }
        if(~n&1)
        {
            int iv2=qpower(2,mod-2);
            for(int i=0;i<l;i++)p[i]=(1ll*tmpp[i]*tmpq[i^l]+1ll*tmpp[i^l]*tmpq[i])%mod*iv2%mod;
        }
        else
        {
            int iv2=qpower(2,mod-2);
            for(int i=0;i<l;i++)p[i]=(1ll*tmpp[i]*tmpq[i^l]+1ll*(mod-tmpp[i^l])*tmpq[i])%mod*iv2%mod;
            for(int i=0;i<l;i++)p[i]=1ll*p[i]*iw[i]%mod;
        }
        for(int i=0;i<l;i++)q[i]=1ll*tmpq[i]*tmpq[i^l]%mod;
        n>>=1;
    }
    if(!isfirst){Poly::IDFT(p,l),Poly::IDFT(q,l);}//如果发生过递归那么p,q里存的是DFT的值，需要做一次IDFT
    Poly::Inv(q,q,K),Poly::mul(p,q,p,K);//n<K后直接做除法计算
    cout<<p[n]<<endl;
}

计算 z^n\bmod F(z)

我们还是来考虑多项式带余除法的式子:

z^n=A(z)F(z)+H(z)

翻转系数

1=A_{rev}(z)Q(z)+z^{n-d+1}H_{rev}(z)

两边用 Q(z) 除

\frac{1}{Q(z)}=A_{rev}(z)+z^{n-d+1}\frac{H_{rev}(z)}{Q(z)}

由于 A_{rev}(z) 的次数为 n-d，所以我们就知道

\mathcal{F}_{n,d}(\frac{1}{Q(z)})\equiv \frac{H_{rev}(z)}{Q(z)}\pmod {z^d}

其中 \mathcal{F}_{n,d}(\sum_{i\geq 0}c_iz^i)=\sum_{0\leq i<d}c_{n-d+1+i}z^i，也就是取幂级数的 n-d+1 到 n 项.

由于 H_{rev}(z) 的次数不超过 d，所以就有

H_{rev}(z)=Q(z)\mathcal{F}_{n,d}(\frac{1}{Q(z)})\bmod z^d

于是关键在于计算 \mathcal{F}_{n,d}(\dfrac{1}{Q(z)}). 我们还是使用前面线性递推的技巧:

\begin{aligned} \mathcal{F}_{n,d}(\frac{1}{Q(z)})&=\mathcal{F}_{n,d}\left(\frac{Q(-z)}{Q(z)Q(-z)}\right)\\ &=\mathcal{F}_{n,d}\left(Q(-z)z^{n-2d+1}\mathcal{F}_{n,2d}\left(\frac{1}{V(z^2)}\right)\right)\\ &=\mathcal{F}_{2d-1,d}\left(Q(-z)\mathcal{F}_{n,2d}\left(\frac{1}{V(z^2)}\right)\right) \end{aligned}

其中 V(z^2)=Q(z)Q(-z)，这里的想法是把 Q(-z) 拆出来做乘法，由于 Q(-z) 的次数只有 d，所以剩下的东西次数低于 n-2d+1 的部分没有贡献. 我们容易观察到 \mathcal{F}_{n,2d}\left(\dfrac{1}{V(z^2)}\right) 是能够递归的，令 S(z)=\mathcal{F}_{\lfloor n/2\rfloor,d}\left(\dfrac{1}{V(z)}\right)，那么:

\mathcal{F}_{n,2d}\left(\frac{1}{V(z^2)}\right)=\begin{cases}zS(z^2)&n\bmod 2=0\\ S(z^2)&n\bmod 2=1\end{cases}

不过这玩意的空间好像是 \Theta(d\log n) 的，有点大( 前面的常数优化可以继续类似地使用，不过我现在既不想写也不想分析乘法次数(

BM+线性递推 学习笔记

BM

线性递推

求一段点值

2020.10.21 upd

计算 z^n\bmod F(z)

BM+线性递推学习笔记