BZOJ 3992: [SDOI2015]序列统计（NTT+生成函数）

Posted On 2019年2月25日

3992: [SDOI2015]序列统计

Simplified Description

\(给定一个[0,m-1]范围内的数字集合S，从中选择n个数字（可重复）构成序列。给定x，求序列所有数字乘积%m后为x的序列方案数%1004535809。1<=n<=10^9，3<=m<=8000，m为素数\)

Solution

乘积为定值怎么办？化为对数即可转化为加法

对数是实数怎么办？模意义下用原根代替即可

因为取得数字不限，对于生成函数求\(n\)次方的\(X\)（\(x\)的对应原根对数）次项即可

注意NTT时对于\(>=m-1\)的次项应加回\(i \% (m-1)\)次项

#include <cstdio>
#include <iostream>
#include <cmath>
using namespace std;
inline int read()
{
	int f=1,k=0;char c=getchar();
	while (c<'0'||c>'9')c=='-'&&(f=-1),c=getchar();
	while (c>='0'&&c<='9')k=k*10+c-'0',c=getchar();
	return k*f;
}
const int maxn=80005<<3,P=1004535809;
int contain[maxn],mi[maxn],mp[maxn],M,a[maxn],ans[maxn];
inline int Plus(const int&a,const int&b)
{
	int tp=a+b;
	if (tp<0)tp+=P;
	if (tp>=P)tp-=P;
	return tp;
}
inline int Minus(const int&a,const int&b)
{
	int tp=a-b;
	if (tp<0)tp+=P;
	if (tp>=P)tp-=P;
	return tp;
}
inline int Pow(int x,int b,int P)
{
	int ans=1;
	for (;b;b>>=1,x=1ll*x*x%P)
		if (b&1)ans=1ll*ans*x%P;
	return ans;
}
struct Template_getg
{
	int prime[maxn],p[maxn],cnt;
	bool nprime[maxn];
	inline bool ok(const int x,const int M)
	{
		for (int i=1;i<=cnt;i++)
			if (Pow(x,(M-1)/p[i],M)%M==1)return 0;
		return 1;
	}
	inline int main(int m)
	{
		if (m==2)return 1;
		int tot=0,M=m;--m;
		for (int i=2;i<=m;i++)
		{
			if (!nprime[i])prime[++tot]=i;
			for (int j=1;i*prime[j]<=m&&j<=tot;j++)
			{
				nprime[i*prime[j]]=1;
				if (i%prime[j]==0)break;
			}
		}
		for (int i=1;i<=tot;i++)
			if (m%prime[i]==0)
			{
				p[++cnt]=prime[i];
				m/=prime[i];
				while (m%prime[i]==0)m/=prime[i];
			}
		for (int i=2;;i++)
			if (ok(i,M))return i;
		return -1;
	}
}getg;
struct Template_NTT
{
	int R[maxn];
	inline void pre(int&n)
	{
		int m=n<<1,L=0;--m;
		for (n=1;n<m;n<<=1)++L;
		for (int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
		//		cout << "N:"<<n<<endl;
	}
	inline void main(int*a,const int n,const int f)
	{
		for (int i=0;i<n;i++)
			if (i<R[i])swap(a[i],a[R[i]]);
		for (int i=1;i<n;i<<=1)
		{
			int wn=Pow(3,(P-1)/(i<<1),P);
			if (f==-1)wn=Pow(wn,P-2,P);
			for (int j=0;j<n;j+=(i<<1))
				for (int k=0,w=1;k<i;k++,w=1ll*w*wn%P)
				{
					int x=a[j+k],y=1ll*w*a[j+k+i]%P;
					a[j+k]=Plus(x,y);
					a[j+k+i]=Minus(x,y);
				}
		}
		if (f==-1)for (int i=0,inv=Pow(n,P-2,P);i<n;i++)a[i]=1ll*a[i]*inv%P;
	}
	int tmp[maxn];
	inline void mul(int*a,int*b,int n)
	{
		for (int i=0;i<n;i++)tmp[i]=b[i];
		main(a,n,1);
		main(tmp,n,1);
		for (int i=0;i<n;i++)a[i]=1ll*a[i]*tmp[i]%P;
		main(a,n,-1);
		for (int i=M-1;i<n;i++)if (a[i])a[i%(M-1)]=Plus(a[i],a[i%(M-1)]),a[i]=0;
	}
}NTT;
signed main()
{
	int n=read(),m=read(),x=read()%m,siz=read(),N=m;
	M=m;
	int g=getg.main(m);
	NTT.pre(N);
	mp[mi[0]=1]=0;
	for (int i=1;i<M-1;i++)mp[mi[i]=1ll*mi[i-1]*g%m]=i;
	for (int i=1,v;i<=siz;i++)(v=read())&&(contain[mp[v]]=1);
	for (int i=0;i<M;i++)a[i]=contain[i];
	ans[0]=1;
	for (;n;n>>=1,NTT.mul(a,a,N))
		if (n&1)
			NTT.mul(ans,a,N);
	printf("%d\n",ans[mp[x]]);
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

#include <cstdio>

#include <iostream>

#include <cmath>

using namespace std;

inline int read()

{

int f=1,k=0;char c=getchar();

while (c<'0'||c>'9')c=='-'&&(f=-1),c=getchar();

while (c>='0'&&c<='9')k=k*10+c-'0',c=getchar();

return k*f;

}

const int maxn=80005<<3,P=1004535809;

int contain[maxn],mi[maxn],mp[maxn],M,a[maxn],ans[maxn];

inline int Plus(const int&a,const int&b)

{

int tp=a+b;

if (tp<0)tp+=P;

if (tp>=P)tp-=P;

return tp;

}

inline int Minus(const int&a,const int&b)

{

int tp=a-b;

if (tp<0)tp+=P;

if (tp>=P)tp-=P;

return tp;

}

inline int Pow(int x,int b,int P)

{

int ans=1;

for (;b;b>>=1,x=1ll*x*x%P)

if (b&1)ans=1ll*ans*x%P;

return ans;

}

struct Template_getg

{

int prime[maxn],p[maxn],cnt;

bool nprime[maxn];

inline bool ok(const int x,const int M)

{

for (int i=1;i<=cnt;i++)

if (Pow(x,(M-1)/p[i],M)%M==1)return 0;

return 1;

}

inline int main(int m)

{

if (m==2)return 1;

int tot=0,M=m;--m;

for (int i=2;i<=m;i++)

{

if (!nprime[i])prime[++tot]=i;

for (int j=1;i*prime[j]<=m&&j<=tot;j++)

{

nprime[i*prime[j]]=1;

if (i%prime[j]==0)break;

}

for (int i=1;i<=tot;i++)

if (m%prime[i]==0)

{

p[++cnt]=prime[i];

m/=prime[i];

while (m%prime[i]==0)m/=prime[i];

}

for (int i=2;;i++)

if (ok(i,M))return i;

return -1;

}

}getg;

struct Template_NTT

{

int R[maxn];

inline void pre(int&n)

{

int m=n<<1,L=0;--m;

for (n=1;n<m;n<<=1)++L;

for (int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));

// cout << "N:"<<n<<endl;

}

inline void main(int*a,const int n,const int f)

{

for (int i=0;i<n;i++)

if (i<R[i])swap(a[i],a[R[i]]);

for (int i=1;i<n;i<<=1)

{

int wn=Pow(3,(P-1)/(i<<1),P);

if (f==-1)wn=Pow(wn,P-2,P);

for (int j=0;j<n;j+=(i<<1))

for (int k=0,w=1;k<i;k++,w=1ll*w*wn%P)

{

int x=a[j+k],y=1ll*w*a[j+k+i]%P;

a[j+k]=Plus(x,y);

a[j+k+i]=Minus(x,y);

}

if (f==-1)for (int i=0,inv=Pow(n,P-2,P);i<n;i++)a[i]=1ll*a[i]*inv%P;

}

int tmp[maxn];

inline void mul(int*a,int*b,int n)

{

for (int i=0;i<n;i++)tmp[i]=b[i];

main(a,n,1);

main(tmp,n,1);

for (int i=0;i<n;i++)a[i]=1ll*a[i]*tmp[i]%P;

main(a,n,-1);

for (int i=M-1;i<n;i++)if (a[i])a[i%(M-1)]=Plus(a[i],a[i%(M-1)]),a[i]=0;

}

}NTT;

signed main()

{

int n=read(),m=read(),x=read()%m,siz=read(),N=m;

M=m;

int g=getg.main(m);

NTT.pre(N);

mp[mi[0]=1]=0;

for (int i=1;i<M-1;i++)mp[mi[i]=1ll*mi[i-1]*g%m]=i;

for (int i=1,v;i<=siz;i++)(v=read())&&(contain[mp[v]]=1);

for (int i=0;i<M;i++)a[i]=contain[i];

ans[0]=1;

for (;n;n>>=1,NTT.mul(a,a,N))

if (n&1)

NTT.mul(ans,a,N);

printf("%d\n",ans[mp[x]]);

}

0 0 vote

Article Rating

0 评论

Inline Feedbacks

View all comments