bzoj3992-序列统计

给出\(n,m,x,S\)，其中\(S\subseteq [0,m)\)，问有多少个长度为\(n\)的数列\(a\)使得\(a_i\in S\)，并且数列中所有元素的乘积mod \(m\)为\(x\)。

答案对\(1004535809=479\times 2^{21}+1\)取模。\(m\)为质数。\(n\le 10^9,m< 8000\)。

分析

在实数意义下，一个经典的处理乘法的方法是通过取对数转化成加法，例如在计算最小乘积生成树时，我们对边权取对数，转化成一般的最小生成树问题。在模意义下也有类似的巧妙方法，利用模数的原根。

可以证明\(1,2,4,p,2p,p^n\)都有原根，而且只有这些数有原根（\(p\)表示奇素数）。原根\(g\)满足\(g^x\equiv 1\)成立当且仅当\(x=p-1\)。这就是说，用\(g\)的幂可以在模意义下表示出\([1,p)\)的所有数。也就是说，如果模数有原根，那么模意义下的乘法可以转化成原根的指数的加法。

令\(x=g^b\)，那么问题就变成了，在数集中，选出\(n\)个数，一个数可以选多次，有多少种方法能够让选出来的数的和为\(b\)。

运用生成函数\(f\)，那么问题就变成了\(f^n(x)\)，注意这里其实是一个\(\mod m-1\)的循环卷积，最后累加回来即可。\(n\)比较大用快速幂。

这题的关键在于利用原根把乘法转化成加法。

代码

写这个题主要是想试一下任意模数fft的做法，所以谢了两份代码，第一份是ntt的，第二个是拆模数的fft。

ntt

#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long giant;
const int q=479<<21|1;
const int g=3;
int read() {
	int x=0,f=1;
	char c=getchar();
	for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
	for (;isdigit(c);c=getchar()) x=x*10+c-'0';
	return x*f;
}
const int maxj=16;
const int maxn=1<<maxj|1;
int a[maxn],b[maxn],c[maxn],mg,M,xj,rev[maxn],n,m,d[maxn];
int Wn[maxj][2];
int sub(int x,int y) {
	int ret=((giant)x-y)%q;
	if (ret<0) ret+=q;
	return ret;
}
int Plus(int x,int y) {
	return ((giant)x+y)%q;
}
int multi(int x,int y) {
	return (giant)x*y%q;
}
int mi(int x,int y) {
	int ret=1;
	for (;y;y>>=1,x=multi(x,x)) if (y&1) ret=multi(ret,x);
	return ret;
}
int inv(int x) {
	return mi(x,q-2);
}
void ntt(int a[],int op=0) {
	for (int i=0;i<M;++i) if (i<rev[i]) swap(a[i],a[rev[i]]);
	for (int the=1,i=2;i<=M;i<<=1,++the) {
		for (int j=0;j<M;j+=i) {
			for (int k=j,w=1;k<j+i/2;++k,w=multi(w,Wn[the][op])) {
				int u=a[k],v=multi(a[k+i/2],w);
				a[k]=Plus(u,v),a[k+i/2]=sub(u,v);
			}
		}
	}
	if (op) {
		int iv=inv(M);
		for (int i=0;i<M;++i) a[i]=multi(a[i],iv);
	}
}
void self(int a[]) {
	ntt(a);
	for (int i=0;i<M;++i) a[i]=multi(a[i],a[i]);
	ntt(a,1);
}
void mul(int c[],int a[]) {
	ntt(a);
	ntt(c);
	for (int i=0;i<M;++i) c[i]=multi(c[i],a[i]);
	ntt(a,1);
	ntt(c,1);
}
void back(int a[]) {
	memset(d,0,sizeof d);
	for (int i=0;i<M;++i) d[i%(m-1)]=Plus(d[i%(m-1)],a[i]);
	memcpy(a,d,(sizeof d[0])*M);
}
void ans(int a[],int len,int y) {
	memset(c,0,sizeof c);
	c[0]=1;
	for (xj=0,M=1;M<=(len<<1);++xj,M<<=1);
	for (int i=0;i<M;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(xj-1));
	while (y) {
		if (y&1) {
			mul(c,a);
			back(c);
		}
		y>>=1;
		self(a);
		back(a);
	}
}
int getg(int m) {
	for (int i=2;i<m;++i) {
		int tmp=i;
		bool flag=true;
		for (int j=1;j<m;++j,tmp=tmp*i%m) if (tmp==1 && j!=m-1) {
			flag=false;
			break;
		}
		if (flag) return i;
	}	
	return 0;
}
int S[maxn];
int main() {
#ifndef ONLINE_JUDGE
	freopen("test.in","r",stdin);
	freopen("my.out","w",stdout);
#endif
	Wn[0][0]=Wn[0][1]=1;
	for (int i=1;i<maxj;++i) Wn[i][1]=inv(Wn[i][0]=mi(g,(q-1)>>i));
	n=read(),m=read();
	int x=read(),s=read();
	mg=getg(m);
	for (int tmp=1,i=0;i<m-1;++i,tmp=tmp*mg%m) b[tmp]=i;
	for (int i=1;i<=s;++i) S[i]=read();
	for (int i=1;i<=s;++i) if (S[i]) ++a[b[S[i]]];
	ans(a,m,n);
	printf("%d\n",c[b[x]]);
	return 0;
}

fft

任意模数fft的做法是，设模数为\(p\)，把多项式的系数拆成\(a\sqrt p+b\)，分成两个多项式\((A,B)\)。

假设\(X=A\sqrt p+B,Y=C\sqrt p+D\)，那么进行分步计算：

\[X\times Y=A\times C*(\sqrt p*\sqrt p) +(B\times C+A\times D)*\sqrt p+B\times D \mod p \]

分开取模，这样就不会受到精度问题影响了。要开long double。

#include<cstdio>
#include<cctype>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long double lb;
int read() {
	int x=0,f=1;
	char c=getchar();
	for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
	for (;isdigit(c);c=getchar()) x=x*10+c-'0';
	return x*f;
}
typedef long long giant;
const int q=479<<21|1;
const int sq=31695;
const int sqb=(giant)sq*sq%q;
const lb pi=acos(-1.0);
struct C {
	lb r,i;
	C (lb r=0,lb i=0):r(r),i(i) {}
};
C operator * (C a,C b) {return C(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
C operator + (C a,C b) {return C(a.r+b.r,a.i+b.i);}
C operator - (C a,C b) {return C(a.r-b.r,a.i-b.i);}
C operator / (C a,lb x) {return C(a.r/x,a.i/x);}
const int maxj=16;
const int maxn=1<<maxj|1;
int M,a[maxn],c[maxn],xj,rev[maxn],b[maxn];
struct Array {
	C a[maxn];
	C& operator [] (int x) {return a[x];}
	void back(int m) {
		for (int i=m-1;i<M;++i) a[i%(m-1)]=a[i%(m-1)]+a[i],a[i]=0;
	}
	void init(int b[]) {
		for (int i=0;i<M;++i) a[i]=b[i];
	}
	void one() {
		for (int i=0;i<M;++i) a[i]=0;
		a[0]=1;
	}
	void fft(int op=1) {
		for (int i=0;i<M;++i) if (i<rev[i]) swap(a[i],a[rev[i]]);
		for (int i=2;i<=M;i<<=1) {
			C wn(cos(2*pi/i),op*sin(2*pi/i));
			for (int j=0;j<M;j+=i) {
				C w(1);
				for (int k=j;k<j+i/2;++k,w=w*wn) {
					C u=a[k],v=a[k+i/2]*w;
					a[k]=u+v,a[k+i/2]=u-v;
				}
			}
		}
		if (op==-1) {
			for (int i=0;i<M;++i) a[i]=a[i]/M;
		}
	}
} A,B,D,E,F,G,H,TA,TC;
Array operator * (Array &a,Array &b) {
	Array ret;
	a.fft(),b.fft();
	for (int i=0;i<M;++i) ret[i]=a[i]*b[i];
	a.fft(-1),b.fft(-1),ret.fft(-1);
	return ret;
}
Array operator + (Array &a,Array &b) {
	Array ret;
	for (int i=0;i<M;++i) ret[i]=a[i]+b[i];
	return ret;
}
void operator *= (Array &a,Array &b) {
	a.fft(),b.fft();
	for (int i=0;i<M;++i) a[i]=a[i]*b[i];
	a.fft(-1),b.fft(-1);
}
void split(Array &a,Array &c,Array &d) {
	for (int i=0;i<M;++i) {
		giant x=(giant)(a[i].r+0.5);
		giant y=x%sq;
		giant z=x/sq;
		c[i]=z;
		d[i]=y;
	}
}
int multi(giant x,giant y) {
	return (giant)x*y%q;
}
int Plus(giant x,giant y) {
	return ((giant)x+y)%q;
}
Array mul(Array &a,Array &b) {
	Array ret;
	split(a,A,B);
	split(b,D,E);
	F=A*D;
	G=A*E;
	H=B*D;
	G=G+H;
	H=B*E;
	for (int i=0;i<M;++i) {
		giant x=(giant)(F[i].r+0.5);
		giant y=(giant)(G[i].r+0.5);
		giant z=(giant)(H[i].r+0.5);
		x=multi(x,sqb);
		y=multi(y,sq);
		x=Plus(x,Plus(y,z));
		ret[i]=x;
	}
	return ret;
}
void ans(int a[],int m,int y) {
	for (xj=0,M=1;M<=(m<<1);M<<=1,++xj);
	for (int i=0;i<M;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(xj-1));
	TA.init(a);
	TC.one();
	while (y) {
		if (y&1) TC=mul(TC,TA),TC.back(m);
		y>>=1;
		TA=mul(TA,TA),TA.back(m);
	}
	for (int i=0;i<M;++i) c[i]=((giant)(TC[i].r+0.5))%q;
}
int getm(int m) {
	for (int i=2;i<m;++i) {
		bool flag=true;
		for (int tmp=i,j=1;j<m;++j,tmp=tmp*i%m) if (tmp==1 && j!=m-1) {
			flag=false;
			break;
		}
		if (flag) return i;
	}
	return 0;
}
int main() {
#ifndef ONLINE_JUDGE
	freopen("test.in","r",stdin);
	freopen("my.out","w",stdout);
#endif
	int n=read(),m=read(),x=read(),s=read();
	int mg=getm(m);
	for (int tmp=1,i=0;i<m-1;++i,tmp=tmp*mg%m) b[tmp]=i;
	for (int i=1,the;i<=s;++i) if ((the=read())!=0) ++a[b[the]];
	ans(a,m,n);
	printf("%d\n",c[b[x]]);
	return 0;
}