bzoj3992 & 洛谷3321 [SDOI2015]序列统计 NTT+快速幂优化dp

题目链接：
bzoj3992
洛咕3321

先考虑一个dp：
令$dp[i][j]$表示$i$个数，乘积为$j$的方案数。令$a[i]$表示集合$S$中的第$i$个数。
那么$dp[i][j\ast k\%m]=\sum dp[i-1][j]\ast v[k]\%mod$，其中$v[k]$表示是否存在$i$使得$a[i]=k$。
程序就不放了，只有10分，放一个快速幂优化的版本，珂以水到60：

#include<stdio.h>
#include<cstring>
#include<algorithm>
#define re register int
#define mod 1004535809
using namespace std;
typedef long long ll;
int read() {
	re x=0,f=1;
	char ch=getchar();
	while(ch<'0' || ch>'9') {
		if(ch=='-')	f=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9') {
		x=10*x+ch-'0';
		ch=getchar();
	}
	return x*f;
}
const int Size=10005;
int n,m,x,s,a[Size],v[Size];
void mul(int *a,int *b,int *c) {
	for(re i=1; i<m; i++) {
		for(re j=1; j<m; j++) {
			re now=i*j%m;
			c[now]=((ll)c[now]+(ll)a[i]*b[j])%mod;
		}
	}
}
int tmp[Size],dp[Size];
int main() {
	n=read();
	m=read();
	x=read();
	s=read();
	for(re i=1; i<=s; i++) {
		a[i]=read();
		v[a[i]]=1;
	}
	dp[1]=1;
	while(n) {
		if(n&1) {
			mul(dp,v,tmp);
			memcpy(dp,tmp,sizeof(tmp));
			memset(tmp,0,sizeof(tmp));
		}
		mul(v,v,tmp);
		memcpy(v,tmp,sizeof(tmp));
		memset(tmp,0,sizeof(tmp));
		n>>=1;
	}
	printf("%d",dp[x%m]);
	return 0;
}

$m$是质数，考虑$m$的原根$g$。根据原根的定义，对于0~m-1的任意值$x$，都有一个$k$使得$g^k=x(mod$ $m)$。
所以对于所有$a[i]$，令$g^k=a[i]$，把$a[i]$换成$k$。那么$a[i]\ast a[j]$就变成$a[i]+a[j]$（指数相加）。
所以dp方程变为$dp[i][(j+k)\%m]=\sum dp[i-1][j]\ast v[k]$
然后类似上面快速幂优化一下就珂以了qwq

毒瘤代码

#include<stdio.h>
#include<cstring>
#include<algorithm>
#include<math.h>
#define re register int
#define mod 1004535809
#define G 3
#define invG 334845270
#define inv2 502267905
using namespace std;
typedef long long ll;
int read() {
	re x=0,f=1;
	char ch=getchar();
	while(ch<'0' || ch>'9') {
		if(ch=='-')	f=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9') {
		x=10*x+ch-'0';
		ch=getchar();
	}
	return x*f;
}
int fpow(int b,int p,int k) {
	int ans=1;
	while(p) {
		if(p&1)	ans=(ll)ans*b%k;
		b=(ll)b*b%k;
		p>>=1;
	}
	return ans;
}
int fpow(int b,int p) {
	return fpow(b,p,mod);
}
inline int inv(int x) {
	return fpow(x,mod-2);
}
int getG(int x) {
	int sx=sqrt(x);
	for(re i=2; i<=x; i++) {
		bool fail=false;
		for(re j=2; j<sx; j++) {
			if(fpow(i,(x-1)/j,x)==1) {
				fail=true;
				break;
			}
		}
		if(!fail) {
			return i;
		}
	}
}
const int Size=30005;
int R[Size];
void NTT(int *a,int n,int type) {
	for(re i=0; i<n; i++) {
		if(i<R[i]) {
			swap(a[i],a[R[i]]);
		}
	}
	for(re i=1; i<n; i<<=1) {
		ll gn=fpow(type==1?G:invG,(mod-1)/(i<<1));
		for(re j=0; j<n; j+=i<<1) {
			ll g=1;
			for(re k=0; k<i; k++) {
				ll x=a[j+k],y=g*a[i+j+k]%mod;
				a[j+k]=(x+y)%mod;
				a[i+j+k]=(x-y+mod)%mod;
				g=g*gn%mod;
			}
		}
	}
	if(type==-1) {
		int iv=inv(n);
		for(re i=0; i<n; i++) {
			a[i]=(ll)a[i]*iv%mod;
		}
	}
}
int A[Size],B[Size],ans[Size];
void Inv(int *a,int *b,int n) {
	b[0]=inv(a[0]);
	int len;
	for(len=1; len<(n<<1); len<<=1) {
		int lim=len<<1;
		for(re i=0; i<len; i++) {
			A[i]=a[i];
			B[i]=b[i];
		}
		for(re i=0; i<lim; i++) {
			R[i]=(R[i>>1]>>1)|((i&1)*len);
		}
		NTT(A,lim,1);
		NTT(B,lim,1);
		for(re i=0; i<lim; i++) {
			b[i]=(ll)B[i]*(2ll-(ll)A[i]*B[i]%mod+mod)%mod;
		}
		NTT(b,lim,-1);
		for(re i=len; i<lim; i++) {
			b[i]=0;
		}
	}
	for(re i=0; i<len; i++)	A[i]=B[i]=0;
	for(re i=n; i<len; i++)	b[i]=0;
}
int tmp[Size];
void mul(int *a,int *b,int *c,int n) {
	int now=1,L=0;
	while(now<(n<<1)) {
		now<<=1;
		L++;
	}
	for(re i=0; i<n; i++) {
		A[i]=a[i];
		B[i]=b[i];
	}
	for(re i=0; i<now; i++) {
		R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
	}
	NTT(A,now,1);
	NTT(B,now,1);
	for(re i=0; i<now; i++) {
		tmp[i]=(ll)A[i]*B[i]%mod;
	}
	NTT(tmp,now,-1);
	for(re i=0; i<n; i++) {
		c[i]=tmp[i];
	}
	for(re i=n; i<now; i++) {
		c[i%n]=(c[i%n]+tmp[i])%mod;
	}
	for(re i=0; i<now; i++) {
		tmp[i]=A[i]=B[i]=0;
	}
}
int n,m,x,s,a[Size],b[Size],g[Size],dp[Size];
int main() {
//	freopen("WA.txt","w",stdout);
	n=read();
	m=read();
	x=read();
	s=read();
	int now=1,gn=getG(m);
	g[0]=-1;
	for(re i=1; i<m-1; i++) {
		g[now=now*gn%m]=i;
	}
	for(re i=1; i<=s; i++) {
		a[i]=read();
		b[g[a[i]]]++;
	}
	dp[0]=1;
	while(n) {
		if(n&1) {
			mul(dp,b,dp,m-1);
		}
		mul(b,b,b,m-1);
		n>>=1;
	}
	printf("%d",dp[g[x]]);
	return 0;
}
/*
819644598 71 64 5
0 3 6 7 10
*/