多项式全家桶 - 常系数齐次线性递推

常系数齐次线性递推

注意：这篇是一个不用很高线性代数基础的理解方法

（图解党狂喜）

问题描述

$f_0,f_1...f_{k-1}$ 已知，并已知递推系数 $a_1,a_2...a_k$

$f_n=\sum\limits_{i=1}^{k} f_{n-k}a_k$

求 $f_n$

$n\le 1145141919810,k\le 114514$

常规的矩阵快速幂过不去了。我们需要一个带 $\log$ 的做法。

换一个思路

考虑把转移看成一个带权 DAG 图：点 $u$ 和点 $u-i$ 有一个有向边，边权为 $a_i$。（$1\le i\le k$，后面如果没提到 $i$ 的取值范围就是这个范围）

然后我们来考虑点 $n$ 的状态把点 $[0,k-1]$ 算了多少遍。假设我们能求出来“多少遍”，那直接就可以求出 $f_n$ 了：每个 $f_{0...k-1}$，乘以对应的遍数，加起来，就行了。

然后我们考虑用拓扑排序求这个东西。设 $c$ 表示遍数数组。

我们的边界条件是 $n$：一遍，$c_n=1$

对于每个当前入度为 $0$ 的点 $u$，遍历它的出边 $v$，设边权为 $w$。那么 c[v]+=c[u]*w，显然。

然后这么做一遍就可以得到 $c_{0...k-1}$，对应的乘起来即可。复杂度 $O(nk)$

然后你会问：艹这nm不就是直接递推么

是的

但是我们反过来考虑了个寂寞 ... 吗？

魔法

考虑 $c$ 的生成函数，设为 $C(x)=\sum\limits_{i=0}^{\infin} x^ic_i$。

一开始我们知道这玩意就是 $x^n$ （在 $n$ 位置有一个 $n$，其余为 $0$）

我们令拓扑排序删边的时候，如果一个点的出度为 $0$ （即出边被删没了），且不是 $[0,k-1]$ 中的点，就把这个点也删除，同时标记它的 $c$ 为 $0$ （这样显然没问题，因为这个点已经不用了）

然后这样显然只会剩下 $c_{0..k-1}$ 有可能有值（其它都是 $0$）

考虑一次拓扑排序做了甚么：假设当前入度为 $0$ 的点是 $u$，那相当于：

• 标记 $c_{u-i}$ 加上 $c_u\times a_i$
• 然后标记 $c_u$ 为 $0$，即 $c_u$ 减去 $c_u$

一共操作了 $k+1$ 个数，并且操作数都是 $c_u$ 的倍数。

在生成函数上考虑，它相当于减去了 $c_u$ 倍的：

$x^u-\sum\limits_{i=1}^{k} x^{u-i}a_i$

这个东西的系数是 $\{1,-a_1,-a_2...-a_k\}$，然后后面都是 $0$。把后面的 $0$ 去掉，只保留这 $k+1$ 项系数，构成的多项式被称作这个递推式的 特征多项式。设它为 $F$

注：如果您有了解过斐波那契的通项是怎么推的，那您应该接触过这个“特征多项式”的定义

斐波那契的特征多项式为：$x^2-x-1$

然后每次我们就是把 $F$ 和当前 $C$ 的最高位（设为 $u$）对齐，然后减去 $F$ 的 $c_u$ 倍 —— 这是拓扑排序干的事情。

我们发现它的本质就是在求余数

回想一下小学竖式计算除法的过程，又因为 $F$ 的最高位为 $1$，所以减去 $F$ 的 $c_u$ 倍其实就是在做带余除法求余数

所以，$C=x^n\mod F$

然后我们只要求出 $x^n\mod F$ 的值，就可以知道遍数数组 $c_{0..k-1}$，对应的乘以初始值 $f_{0..k-1}$，就得到了 $f_n$

然后这个 $C$ 怎么算呢？$n$ 老大了

考虑多项式快速幂 —— 不用 $\ln/\exp$，而是普通的倍增快速幂，参考整数的写法，并把“模”的操作改成多项式取模就可以了。

#include <bits/stdc++.h>
using namespace std;
namespace Flandre_Scarlet
{
	#define N   2000006
	#define GG  3
	#define GI  332748118
	#define mod 998244353
	#define int long long
	#define F(i,l,r) for(int i=l;i<=r;++i)
	#define D(i,r,l) for(int i=r;i>=l;--i)
	#define Fs(i,l,r,c) for(int i=l;i<=r;c)
	#define Ds(i,r,l,c) for(int i=r;i>=l;c)
	#define MEM(x,a) memset(x,a,sizeof(x))
	#define FK(x) MEM(x,0)
	#define Tra(i,u) for(int i=G.st(u),v=G.to(i);~i;i=G.nx(i),v=G.to(i))
	#define p_b push_back
	#define sz(a) ((int)a.size())
	#define all(a) a.begin(),a.end()
	#define iter(a,p) (a.begin()+p)
	#define PUT(a,n) F(i,1,n) printf("%d ",a[i]); puts("");
	int I() {char c=getchar(); int x=0; int f=1; while(c<'0' or c>'9') f=(c=='-')?-1:1,c=getchar(); while(c>='0' and c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar(); return ((f==1)?x:-x);}
	template <typename T> void Rd(T& arg){arg=I();}
	template <typename T,typename...Types> void Rd(T& arg,Types&...args){arg=I(); Rd(args...);}
	void RA(int *p,int n) {F(i,1,n) *p=I(),++p;}
	
	int qpow(int a,int b,int m=mod) {int r=1; while(b) {if (b&1) r=r*a%m; a=a*a%m,b>>=1;} return r;}
	int pgg[30],pgi[30];
	void init()
	{
		F(i,0,23)
		{
			pgg[i]=qpow(GG,(mod-1)>>i);
			pgi[i]=qpow(GI,(mod-1)>>i);
		}
	}
	#define ad(x,y) ((x)+(y)>mod?(x)+(y)-mod:(x)+(y))
	#define dc(x,y) ((x)-(y)<0?(x)-(y)+mod:(x)-(y))
	namespace poly
	{
		int w[N],r[N];
		void NTT(int f[N],int lim,int type)
		{
			F(i,0,lim-1) if (i<r[i]) swap(f[i],f[r[i]]);
			for(int mid=1,pp=1;mid<lim;mid<<=1,++pp)
			{
				int Wn=(type==1?pgg[pp]:pgi[pp]);
				w[0]=1; F(i,1,mid-1) w[i]=w[i-1]*Wn%mod;
				for(int i=0;i<lim;i+=(mid<<1))
				{
					for(int j=0;j<mid;++j)
					{
						register int y=f[i|mid|j]*w[j]%mod;
						f[i|mid|j]=dc(f[i|j],y);
						f[i|j]=ad(f[i|j],y);
					}
				}
			}
			if (type==-1)
			{
				int ivv=qpow(lim,mod-2);
				F(i,0,lim-1) f[i]=f[i]*ivv%mod;
			}
		}
		#define REV F(i,0,lim-1) r[i]=((r[i>>1]>>1)|((i&1)?len:0));
		int pool[10][N];
		inline void Inv(int f[N],int g[N],int n) // 求逆, pool 0~2
		{
			int *iv=pool[0],*a=pool[1],*b=pool[2];
			F(i,0,4*n) iv[i]=a[i]=b[i]=0;
			iv[0]=qpow(f[0],mod-2);

			int len=1,lim=1;
			for(len=1;len<=(n<<1);len<<=1)
			{
				lim=len<<1;
				
				F(i,0,4*len) a[i]=b[i]=0;
				F(i,0,len-1) a[i]=f[i],b[i]=iv[i];
				REV; 
				NTT(a,lim,1); NTT(b,lim,1);
				F(i,0,lim-1) a[i]=(2*b[i]-a[i]*b[i]%mod*b[i]%mod+mod)%mod;
				NTT(a,lim,-1);
				F(i,0,len-1) iv[i]=a[i];
			}
			F(i,0,n-1) g[i]=iv[i];
			F(i,n,lim-1) g[i]=0;
		}
		inline void mul(int f[N],int g[N],int n,bool flag=1) // *=, pool 3~4
        // flag: 是否对n取膜
        {
            int len=1,lim=1; while(len<=(n<<1)) len=lim,lim<<=1;
            int *a=pool[3],*b=pool[4];
            F(i,0,lim-1) a[i]=b[i]=0;
            F(i,0,n-1) a[i]=f[i],b[i]=g[i];
            REV;
            NTT(a,lim,1); NTT(b,lim,1);
            F(i,0,lim-1) a[i]=a[i]*b[i]%mod; 
            NTT(a,lim,-1);
            F(i,0,2*n-1) f[i]=a[i]; F(i,2*n,lim) f[i]=0;
            if (flag) F(i,n,2*n-1) f[i]=0;
        }
        void Mod(int f1[N],int f2[N],int n,int m,int R[N]) // 多项式取模, pool 5~6
        // 这里只保留了余数, 商开到 pool 里而不是传参数修改了
        {
            int *a=pool[5],*b=pool[6],*Q=pool[7];
            F(i,0,4*n) a[i]=b[i]=0;
            F(i,0,n-1) a[i]=f1[n-1-i]; F(i,n-m+1,n-1) a[i]=0; // a=f1_r%(x^(n-m+1))
            F(i,0,m-1) b[i]=f2[m-1-i]; F(i,n-m+1,m-1) b[i]=0; 
            Inv(b,b,n-m+1);
            mul(a,b,n-m+1);
            F(i,0,n-m) Q[i]=a[n-m-i];

            F(i,0,4*n) a[i]=b[i]=0;
            F(i,0,n-1) a[i]=f1[i];
            F(i,0,m-1) b[i]=f2[i];
            int len=1,lim=1; while(len<=n) len=lim,lim<<=1; 
            REV;
            NTT(b,lim,1); NTT(Q,lim,1);
            F(i,0,lim-1) b[i]=b[i]*Q[i]%mod;
            NTT(b,lim,-1); NTT(Q,lim,-1);
            F(i,0,m-2) R[i]=(a[i]-b[i]%mod+mod)%mod; F(i,m-1,lim-1) R[i]=0;
        }
		void PowMod(int f[N],int p,int g[N],int n,int h[N]) // f^p%g, g有n项, 保存在h
		{
			int *res=pool[8];
			F(i,0,8*n) res[i]=0;
			res[0]=1; // res=1
			while(p)
			{
				if (p&1)
				{
					mul(res,f,n,0); // res*=f
					Mod(res,g,2*n,n,res); // res%=g
				}
				mul(f,f,n,0); // f=f*f
				Mod(f,g,2*n,n,f); // f%=g
				p>>=1;
			}
			F(i,0,n-2) h[i]=res[i]; F(i,n,8*n) h[i]=0;
		}
	}
	int n,k;
	int a[N],t[N];
    void Input()
    {
    	Rd(n,k);
    	F(i,1,k) t[i]=(I()%mod+mod)%mod;
    	F(i,0,k-1) a[i]=(I()%mod+mod)%mod;
    }
    int f[N],g[N],c[N];
    void Sakuya()
    {
    	init();

    	f[1]=1; // 这玩意就是 x
    	F(i,1,k) g[k-i]=(mod-t[i]); g[k]=1; // 特征多项式
    	// f^n%g -> c
    	poly::PowMod(f,n,g,k+1,c);
    	
    	int ans=0;
    	F(i,0,k-1) ans+=a[i]*c[i]%mod;
    	ans%=mod;
    	printf("%lld\n",ans);
    }
    void IsMyWife()
    {
        Input();
        Sakuya();
    }
}
#undef int //long long
int main()
{
    Flandre_Scarlet::IsMyWife();
    getchar();
    return 0;
}