[51nod1538]一道难题

news/2024/9/6 6:10:47

先观察一下题目给出的式子:对所有满足$\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


http://www.niftyadmin.cn/n/976768.html

相关文章

[phvia/firman] PHP多进程服务器模型中的惊群

[ 典型场景 ] 典型的多进程服务器模型是这样的&#xff0c;主进程绑定ip&#xff0c;监听port&#xff0c;fork几个子进程&#xff0c;子进程安装信号处理器&#xff0c;随后轮询资源描述符检查是否可读可写&#xff1b; 子进程的轮询又涉及到 IO复用&#xff0c;accept连接&am…

Mysql学习总结(42)——MySql常用脚本大全

备份 &#xff08;所有&#xff09; C:\Program Files\MySQL\MySQL Server 5.6\bin>mysqldump --no-defaults -hlocalhost -P3306 -uroot -p -R test > h:\test.sql 备份 &#xff08;结构&#xff09; C:\Program Files\MySQL\MySQL Server 5.6\bin>mysqldump --no-d…

使用VMware克隆Linux系统

最近在学习使用solr云技术&#xff0c;因为是用来学习操作&#xff0c;因此需要在一台虚拟机上&#xff0c;安装多台LinuxOS。 但是又想偷懒&#xff0c;不想每安装一个LinuxOS&#xff0c;就重新配置Linux环境&#xff0c;所以使用克隆&#xff0c;只需安装好一个模板LinuxOS就…

Learn Forge tutorial - 向导式Forge进阶教程

Autodesk Forge团队陆续的推出了很多学习资料和样例工程&#xff0c;这些资源在Autodesk Forge 学习简谈做了介绍。而部分资源由于API本身的变动&#xff0c;或者代码的变动&#xff0c;可能无法正常使用了&#xff0c;我们会逐个的进行清理&#xff0c;让大家获知最新的更新。…

java @param参数注解

注解&#xff0c;param是参数的解释。如/***param s 这里表示对s的文字说明&#xff0c;描述*/public void aa(String s){}一般java中表示注解&#xff0c;解释一个方法&#xff0c;类&#xff0c;属性的作用

适配器类(便利类)的由来:当你自己写的类中想用某个接口中个别方法的时候(注意:不是所有的方法),肿么办?...

有的时候需要将接口和抽象类配合起来使用&#xff0c;这样可以为开发者提供相当的便利性&#xff0c;开发者觉得哪个方便就选用哪个。这样的抽象类称为便利类。此时&#xff0c;便利类并不需要实现接口的所有方法&#xff0c;可以留给继承它的子类去实现它们。 抽象父类提供给子…

Apache kafka 简介

2019独角兽企业重金招聘Python工程师标准>>> kafka是一种高吞吐量的分布式发布订阅消息系统&#xff0c;它可以处理消费者规模的网站中的所有动作流数据 kafka名词解释 Broker&#xff1a;Kafka 集群包含一个或多个服务器&#xff0c;这种服务器被称为 broker。 Top…

华为招聘机试整理5:简单四则运算

华为招聘机试整理5&#xff1a;简单四则运算 题目&#xff1a;简单四则运算 问题描写叙述&#xff1a; 输入一个仅仅包括个位数字的简单四则运算表达式字符串&#xff0c;计算该表达式的值 注&#xff1a; 1、表达式仅仅含 , -, 乘, / 四则运算符&#xff0c;不含括号 2、表达式…