[51nod1538]一道难题

本文介绍了一种解决特定数学问题的算法,该算法利用多项式取模优化递归过程,实现了高效计算。通过反转系数变换和多项式求逆等手段,解决了高阶常系数线性递推问题。

先观察一下题目给出的式子:对所有满足$\begin{align*}\sum\limits_{i=1}^na_ib_i=m\end{align*}$的$b_{1\cdots n}$,计算$\begin{align*}\binom{\sum\limits_{i=1}^nb_i}{b_1\cdots b_n}\end{align*}$的和,注意这不是组合数而是多项式系数

先证一个微小的式子:$\begin{align*}\binom n{n_1\cdots n_t}=\sum\limits_{i=1}^t\binom{n-1}{n_1\cdots n_i-1\cdots n_t}\end{align*}$

先说明$\begin{align*}\binom n{n_1\cdots n_t}\end{align*}$计数了(把$n$个不同的球放入$t$个不同的盒子,使得第$i$个盒子有$n_i$个球)的方案数,我们先从$n$个球里选$n_1$个球放入第一个盒子,再从$n-n_1$个球中选出$n_2$个球放入第二个盒子,以此类推,这导出$\begin{align*}\binom n{n_1}\binom{n-n_1}{n_2}\cdots\binom{n-n_1-\cdots-n_{t-1}}{t_n}=\dfrac{n!}{\prod\limits_{i=1}^t{n_i!}}=\binom n{n_1\cdots n_t}\end{align*}$(左边是组合数,右边是多项式系数)

我们还可以这样放球:随便选一个球,硬点这个球放入第$i$个盒子,再对剩下$n-1$个球进行分配,所以$\begin{align*}\binom n{n_1\cdots n_t}=\sum\limits_{i=1}^t\binom{n-1}{n_1\cdots n_i-1\cdots n_t}\end{align*}$

回到这个题目,$b_i$减少$1$对应着$m$减少$a_i$,设答案为$f_m$,则我们有$\begin{align*}f_m=\sum\limits_{i=1}^nf_{m-a_i}\end{align*}$,然后它就变成了常系数线性递推的模板题

但是我们发现这题的递推阶数可以到$23333$,暴力取模不可行,所以我们来写一个多项式取模

多项式除法:给定$n$次多项式$A(x)$和$m$次多项式$B(x)$,要求$n-m$次多项式$D(x)$和$m-1$次多项式$R(x)$使得$A(x)=B(x)D(x)+R(x)$

先定义一个反转系数变换:对$n$次多项式$P(x)$,$P^R(x)=x^nP\left(\dfrac1x\right)$,相当于把这个多项式的系数reverse一下

然后把原式中的$x$全部用$\dfrac1x$代替,等式两边同乘$x^n$,看看会发生什么

$\begin{align*}x^nA\left(\dfrac1x\right)&=x^nB\left(\dfrac1x\right)D\left(\dfrac1x\right)+x^nR\left(\dfrac1x\right)\\A^R(x)&=B^R(x)D^R(x)+x^{n-m+1}R^R(x)\\A^R(x)&\equiv B^R(x)D^R(x)(\text{mod }x^{n-m+1})\end{align*}$

然后你发现余数被削除了,所以$D(x)=\left(\dfrac{A^R(x)}{B^R(x)}\%x^{n-m+1}\right)^R$,直接用多项式求逆就可以计算了

算出$D(x)$之后就可以直接算$R(x)$了,总时间复杂度$O(n\log_2n)$

然后就做完了,一开始T爆,在本机跑了差不多20s才跑得过极限数据

后来加了两个优化:①从高位到低位算快速幂,这样可以每次只做一次多项式取模(平方),乘$x$再取模可以直接$O(n)$处理②预处理出特征多项式反转后求逆结果的点值,因为全程都在对它取模

还是过不了,以为是辣鸡卡常题,但猛然发现自己把模数声明成变量了,改成常量后3s出解然后就过了(捂脸

经过大量常数优化的代码简直不能看,但我也懒得改了==

#include<stdio.h>
#include<string.h>
const int mod=104857601,lol=23333;
typedef long long ll;
int mul(int a,int b){return a*(ll)b%mod;}
int ad(int a,int b){return(a+b)%mod;}
int de(int a,int b){return(a-b)%mod;}
void swap(int&a,int&b){
	int c=a;
	a=b;
	b=c;
}
int max(int a,int b){return a>b?a:b;}
int pow(int a,int b){
	int s=1;
	while(b){
		if(b&1)s=mul(s,a);
		a=mul(a,a);
		b>>=1;
	}
	return s;
}
int rev[65536],N,iN;
void pre(int n){
	int i,k;
	for(N=1,k=0;N<n;N<<=1)k++;
	for(i=0;i<N;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
	iN=pow(N,mod-2);
}
void ntt(int*a,int on){
	int i,j,k,t,w,wn;
	for(i=0;i<N;i++){
		if(i<rev[i])swap(a[i],a[rev[i]]);
	}
	for(i=2;i<=N;i<<=1){
		wn=pow(3,(on==1)?(mod-1)/i:(mod-1-(mod-1)/i));
		for(j=0;j<N;j+=i){
			w=1;
			for(k=0;k<i>>1;k++){
				t=mul(w,a[i/2+j+k]);
				a[i/2+j+k]=de(a[j+k],t);
				a[j+k]=ad(a[j+k],t);
				w=mul(w,wn);
			}
		}
	}
	if(on==-1){
		for(i=0;i<N;i++)a[i]=mul(a[i],iN);
	}
}
int t0[65536];
void getinv(int*a,int*b,int n){
	if(n==1){
		b[0]=pow(a[0],mod-2);
		return;
	}
	int i;
	getinv(a,b,n>>1);
	pre(n<<1);
	memset(t0,0,sizeof(t0));
	memcpy(t0,a,(n+1)<<2);
	ntt(t0,1);
	ntt(b,1);
	for(i=0;i<N;i++)b[i]=(2ll-b[i]*(ll)t0[i])%mod*b[i]%mod;
	ntt(b,-1);
	memset(b+n,0,(N-n)<<2);
}
void reverse(int*a,int n){
	for(int i=0;i<=n>>1;i++)swap(a[i],a[n-i]);
}
void add(int*a,int n,int*b,int m,int*c,int&k){
	k=max(n,m);
	for(int i=0;i<=k;i++)c[i]=ad(a[i],b[i]);
	while(k!=0&&c[k]==0)k--;
}
void dec(int*a,int n,int*b,int m,int*c,int&k){
	k=max(n,m);
	for(int i=0;i<=k;i++)c[i]=de(a[i],b[i]);
	while(k!=0&&c[k]==0)k--;
}
int ta[65536],tb[65536];
void mul(int*a,int n,int*b,int m,int*c,int&k){
	int i;
	k=n+m;
	pre(k+2);
	memset(ta,0,N<<2);
	memset(tb,0,N<<2);
	memcpy(ta,a,(n+1)<<2);
	memcpy(tb,b,(m+1)<<2);
	ntt(ta,1);
	ntt(tb,1);
	for(i=0;i<N;i++)c[i]=mul(ta[i],tb[i]);
	ntt(c,-1);
	memset(c+k+1,0,(N-k-1)<<2);
}
int t1[65536],r1[65536],r2[65536];
void div(int*a,int n,int*b,int m,int*c,int&k){
	if(n<m){
		k=0;
		return;
	}
	int i,rn;
	memset(tb,0,sizeof(tb));
	memcpy(ta,a,(n+1)<<2);
	memcpy(tb,b,(m+1)<<2);
	rn=32768;
	pre(65536);
	reverse(ta,n);
	memset(ta+rn,0,(N-rn)<<2);
	ntt(ta,1);
	for(i=0;i<N;i++)c[i]=mul(ta[i],r2[i]);
	ntt(c,-1);
	k=n-m;
	reverse(c,k);
	memset(c+k+1,0,(N-k-1)<<2);
	while(k!=0&&c[k]==0)k--;
}
void modulo(int*a,int n,int*b,int m,int*c,int&k){
	if(n<m){
		k=n;
		memcpy(c,a,(k+1)<<2);
		return;
	}
	div(a,n,b,m,c,k);
	mul(c,k,b,m,c,k);
	dec(a,n,c,k,c,k);
}
int t2[65536];
void sqr(int*a,int&n,int*c,int k){
	int i;
	pre((n+1)<<1);
	memset(ta,0,N<<2);
	memcpy(ta,a,(n+1)<<2);
	ntt(ta,1);
	for(i=0;i<N;i++)t2[i]=mul(ta[i],ta[i]);
	ntt(t2,-1);
	memset(a,0,N<<2);
	modulo(t2,n<<1,c,k,a,n);
}
int d[30010],f[30010],ch[65536],one[65536];
int work(ll m){
	int n1,n2,M,i,l;
	for(M=lol;d[M]==0;M--);
	n1=M;
	ch[M]=1;
	for(i=1;i<=M;i++)ch[i-1]=-d[M+1-i];
	for(i=0;i<=M;i++)r1[i]=ch[M-i];
	getinv(r1,r2,32768);
	ntt(r2,1);
	n2=0;
	one[0]=1;
	m+=M-1;
	for(l=0;1ll<<l<=m;l++);
	l--;
	while(~l){
		if((m>>l)&1){
			for(i=n2;i>=0;i--)one[i+1]=one[i];
			one[0]=0;
			n2++;
			if(n2>n1){
				for(i=n2-1;i>0;i--)one[i]=de(one[i],mul(ch[i-1],one[n2]));
				one[n2]=0;
				while(n2!=0&&one[n2]==0)n2--;
			}
			if(n2==n1){
				for(i=n2-1;i>=0;i--)one[i]=de(one[i],mul(ch[i],one[n2]));
				one[n2]=0;
				while(n2!=0&&one[n2]==0)n2--;
			}
		}
		if(l)sqr(one,n2,ch,n1);
		l--;
	}
	return(one[M-1]+mod)%mod;
}
int main(){
	int n,i,a,A,B;
	ll m;
	scanf("%d%lld%d%d%d",&n,&m,&a,&A,&B);
	d[a]++;
	for(i=1;i<n;i++){
		a=(a*A+B)%lol+1;
		d[a]++;
	}
	printf("%d",work(m));
}

转载于:https://www.cnblogs.com/jefflyy/p/8743909.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值