常系数齐次递推模板

NTT版，跑的贼慢_(:з」∠)_

#include <vector>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
using namespace std;

typedef long long ll;

int read() {
    char ch;
    for(ch=getchar();ch<'0'||ch>'9';ch=getchar());
    int x=ch-'0';
    for(ch=getchar();ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
    return x;
}

const int N=1e6+5,Mo=998244353;

void inc(ll &x,ll y) {
	x+=y;
	x-=x>=Mo?Mo:0;
	x+=x<0?Mo:0;
}

ll pwr(ll x,ll y) {
	ll z=1;
	for(;y;y>>=1,x=x*x%Mo)
		if (y&1) z=z*x%Mo;
	return z;
}

int n,k,K;

ll W[2][N];

void init(int b){
	for(int i=1;i<(1<<b);i<<=1){
		ll wn=pwr(3,(Mo-1)/(i<<1));
		for(int j=0;j<i;++j) W[1][i+j]=(j?wn*W[1][i+j-1]%Mo:1);
		wn=pwr(3,Mo-1-(Mo-1)/(i<<1));
		for(int j=0;j<i;++j) W[0][i+j]=(j?wn*W[0][i+j-1]%Mo:1);
	}
}

void DFT(ll *a,int len,int flag) {
	if (flag==-1) flag=0;
	for(int i=0,j=0;i<len;++i){
		if(i<j) swap(a[i],a[j]);
		for(int k=len>>1;(j^=k)<k;k>>=1);
	}
	for(int i=1;i<len;i<<=1)
		for(int j=0;j<len;j+=(i<<1))
			for(int k=0;k<i;++k) {
				ll x=a[j+k],y=a[j+k+i]*W[flag][i+k]%Mo;
				a[j+k]=(x+y)%Mo;
				a[j+k+i]=(x-y)%Mo;
			}
	ll inv=pwr(len,Mo-2);
	if (!flag) for(int i=0;i<len;i++) a[i]=a[i]*inv%Mo;
}

ll f[N],g[N],c[N];

void get_Inv(ll *a,ll *b,int n) {
	if (n==1) {b[0]=pwr(a[0],Mo-2);return;}
	get_Inv(a,b,(n+1)>>1);
	int len=1;
	for(;len<n<<1;len<<=1);
	fo(i,0,n-1) c[i]=a[i];fo(i,n,len-1) c[i]=0;
	fo(i,(n+1)>>1,len-1) b[i]=0;
	DFT(c,len,1);DFT(b,len,1);
	fo(i,0,len-1) b[i]=(2*b[i]-b[i]*b[i]%Mo*c[i])%Mo;
	DFT(b,len,-1);
	fo(i,n,len-1) b[i]=0;
}

ll ta[N],tb[N],tc[N],inv[N];

ll F[N],G[N],H[N],B[N],C[N],E[N];

ll da[N],db[N];

void pre(int n,int m) {
	int t=n-m+1,len;
	for(len=1;len<t<<1;len<<=1);
	fo(i,0,m) tb[i]=H[m-i];fo(i,m+1,t) tb[i]=0;
	get_Inv(tb,inv,t);
	for(len=1;len<n<<1;len<<=1);
	fo(i,t,len-1) inv[i]=0;
	DFT(inv,len,1);
}

void get_Mod(ll *a,ll *b,ll *d,ll *r,int n,int m) {
	int t=n-m+1,len;
	for(len=1;len<n<<1;len<<=1);
	fo(i,0,n) ta[i]=a[n-i];fo(i,n+1,len-1) ta[i]=0;
	DFT(ta,len,1);
	fo(i,0,len-1) tb[i]=ta[i]*inv[i]%Mo;
	DFT(tb,len,-1);
	fo(i,0,t-1) d[i]=tb[t-i-1];
	fo(i,t,len-1) d[i]=0;
	fo(i,0,len-1) ta[i]=tb[i]=0;
	fo(i,0,m) ta[i]=b[i];
	fo(i,0,t-1) tb[i]=d[i];
	DFT(ta,len,1);DFT(tb,len,1);
	fo(i,0,len-1) tc[i]=ta[i]*tb[i]%Mo;
	DFT(tc,len,-1);
	fo(i,0,m-1) r[i]=(a[i]-tc[i])%Mo;
	fo(i,m,len-1) r[i]=0;
}

void mult(ll *a,ll *b,ll *c,int len) {
	fo(i,0,len-1) da[i]=a[i],db[i]=b[i];
	DFT(da,len,1);DFT(db,len,1);
	fo(i,0,len-1) da[i]=da[i]*db[i]%Mo;
	DFT(da,len,-1);
	get_Mod(da,H,E,c,len-1,K);
}

int main() {
	init(17);
	scanf("%d%d",&n,&k);
	ll ans=1;
	fo(i,1,k) ans=ans*(n-i+1)%Mo;
	fo(i,1,k) ans=ans*pwr(i,Mo-2)%Mo;
	G[0]=1;
	fo(i,1,k-1) fd(j,k*k,i) inc(G[j],Mo-G[j-i]);
	fd(j,k*k,1) inc(G[j],Mo-G[j-1]);
	fd(j,k*k,1) inc(G[j],Mo-G[j-1]);
	K=k*(k-1)/2+2;
	int len=1;for(;len<K*2;len<<=1);
	get_Inv(G,F,len);
	fo(j,0,K) G[j]=Mo-G[j];G[0]=0;
	n-=k*(k-1)/2;
	if (n<len) {
		(ans-=F[n])%=Mo;
		printf("%lld\n",(ans+Mo)%Mo);
		return 0;
	}
	int m=n-K+1;fo(i,0,K-1) H[i]=Mo-G[K-i];H[K]=1;
	B[1]=1;C[0]=1;

	pre(len-1,K);

	for(;m;m>>=1) {
		if (m&1) mult(C,B,C,len);
		mult(B,B,B,len);
	}
	fo(i,0,K-1) (ans-=C[i]*F[i+K-1]%Mo)%=Mo;
	printf("%lld\n",(ans+Mo)%Mo);
	return 0;
}