Luogu3321/BZOJ3992 [SDOI2015]序列统计

题目传送门

题目大意

给定一个集合(S)，(forall xin S,sin [0,m-1])，求从这个集合(S)中选取(n)个数，使得乘积为(x)的方案数。

算法分析

这其实是一类非常经典的问题。

引入例题1

给定一个集合(S)，(forall xin S,sin [1,n])，在这个集合中选取2个数（可以重复），使得加和为(A)的方案数。

不难想到可以从数值角度入手，我们运用桶思想，设在集合中(i)这个数字出现了(h(i))次，则不难得到答案为

[sum_{i=1}^Ah(i) imes h(A-i) ]

显然这个可以(O(n))求解，具体实现不再赘述。

当然这个式子符合多项式卷积的一般形式，所以也可以(O(nlog n))求解。（不过一道普及组难度的题目被搞成这样好吗）

引入例题2

给定一个集合(S)，(forall xin S,sin [1,n])，在这个集合中选取3个数（可以重复），使得加和为(A)的方案数。

类比于上一道例题，不难得到答案为

[sum_{i=1}^Asum_{j=1}^{A-i}h(i) imes h(j) imes h(A-i-j) ]

时间复杂度是(O(n^2))，似乎不是很好。

这个能不能优化呢？

[1+2=3 ]

选取3个数可以看作在选取1个数的集合(S_1)和选取2个数的集合(S_2)中分别选取1个数。

我们可以先求出选取2个数的方案，记从(S)中选取2个数加和为(i)的方案数为(h_2(i))，则不难得到答案为

[sum_{i=1}^Ah(i) imes h_2(A-i) ]

那如何快速求出(h_2)呢？

由上一道例题可得

[h_2(x)=sum_{i=1}^xh(i) imes h(x-j) ]

这是标准的卷积形式，可以使用FFT/NTT快速求出(h_2)，时间复杂度为(O(nlog n))。

推荐题目：

• 【bzoj3513】[MUTC2013]idiots
• 【bzoj3771】Triple

引入例题3

给定一个集合(S)，(forall xin S,sin [1,n])，在这个集合中选取(k)个数（可以重复），使得加和为(A)的方案数。

类比于例题2，不难得到一种方法：设(h_t(i))表示从集合中选取(t)个数加和为(i)的方案数，则有如下的递推式：

[h_t(x)=sum_{i=1}^nh_{t-1}(i) imes h(x-i) ]

这样时间复杂度是(O(kmlog m))，还有一定的优化空间。

考虑倍增，设(h_t(i))表示从集合中选取(2^t)个数加和为(i)的方案数。选取(2^t)个数相当于在选取(2^{t-1})个数的集合中选取2个数，则不难得到如下的递推式：

[h_t(x)=sum_{i=1}^nh_{t-1}(i) imes h_{t-1}(x-i) ]

这样预处理出(h_t)，把(k)二进制拆解，类比于快速幂，可以在(O(mlog mlog k))的时间复杂度下实现。

核心代码如下：

ans[0]=1;//初始化
for(int i=30;~i;i--){
    if((k>>i)&1){
    	mul(h[i],ans,ans);
    }
}

引入例题4

给定一个集合(S)，(forall xin S,sin [0,n-1])，在这个集合中选取(k)个数（可以重复），使得加和在(mod n)意义下为(A)的方案数。

和上一题相比，这一题需要的是模意义下的加法，因此核心算法和上一题一样，只是有一个细节要注意：多项式乘法后，需要把(ge n)的那些项累加到对应模运算以后的位置，因为这些也是可能的方案，对答案有贡献。

可能用代码表示会更加清楚：

void mul(int *f,int *g,int *ans){
	f-->tmp1
    g-->tmp2//复制两个数组，保证f,g本身对应的数组不发生改变
    ntt(tmp1,1);ntt(tmp2,1);
    for(int i=0;i<N;i++)tmp1[i]=tmp1[i]*tmp2[i]%mod;
    ntt(tmp1,-1);
	// 以上为ntt常规操作
    for(int i=0;i<n;i++)ans[i]=tmp1[i];
    for(int i=n;i<N;i++)ans[i%n]=(ans[i%n]+tmp1[i])%mod;//这就是上述的变化
}

回到原题

引入例题4和原题之间只有一步之遥，唯一的变化是加法变成了乘法。（变量名不统一这类的就不考虑了）

这里我们把乘法变成加法！

回想初等函数，我们发现有2种函数可以实现乘法与加法的互化：指数函数和对数函数。

这里略讲一下指数函数和对数函数，防止有些同学没学过。

---

指数函数，形如(f(x)=a^x,age 0,a eq1)，根据幂运算性质，不难得到：(f(x+y)=f(x) imes f(y))。

取对数运算，对于(a^x=N)，则有(log_aN=x)，根据(a^x imes a^y=a^{x+y})，不难得出(log_aa^x+log_aa^y=log_aa^{x+y}=log_a(a^x imes a^y))，推广可得(log_ax+log_ay=log_a(x imes y))。

对数函数，形如(f(x)=log_ax,age 0,a eq 1,x>0)，根据对数运算相关性质，不难得到(f(xy)=f(x)+f(y))。

指数函数可以使加法转化为乘法，对数函数可以使乘法转化为加法。

---

但是指对运算都是在实数域上的，而本题是在模意义下的，怎么处理呢？

模意义下的指数/对数？

求指数还是很容易的，随便以一个数为底（设为(a)），在模意义下求出(a^x,xin[1,m-1])，那么根据对数的定义，则可令(log_a(a^xmod m)=x)，这样建立映射关系。

看上去似乎没有问题？

我们举一个例子，如果(m=7,a=2)，那么(a^1 equiv a^4equiv 2 (mod m))，那么(log_a(2)=?)

我们发现这种映射关系必须是一一对应的，因此并不是随便取一个数为底就可以的。

我们需要找到一个底数(a)，使得(a^1,a^2,cdots,a^{m-1})在(mod m)意义下互不相同。

注意：这里没有包含(a^0)，这是因为在(mod m)意义下不存在一个数(x)使得(a^xequiv 0(mod m))，因此实际上只有(m-1)个可对应位置，在本题中(S)集合中的(0)对答案无贡献，故忽略(a^0)。

这个东西似乎很熟悉啊，这好像就是原根。

对于(p)的原根(g)，满足(g^1,g^2,cdots,g^{p-1})在(mod p)意义下互不相同。

因此我们只需要找到(m)的一个原根，即可把乘法转化为加法。

综上，我们得到这一题的解决方法：

1. 求出(m)的原根，并把集合(S)中的数(x)转化为(log_gx)。
2. 求出(h_t)，含义与上面引入例题的类似，其中(h_0(i))表示集合(S)中满足(log_gx=i)数的数量。
3. 将(n)进行二进制拆分，并做多项式卷积，求出在去过对数的集合中选取(n)个数加和为(i)的方案数为(ans(i))。
4. 对于输入的(X)，答案即为(ans(log_gX))。

时间复杂度为(O(mlog mlog n))

代码实现

#include<bits/stdc++.h>
using namespace std;
#define maxn 100005
#define maxm 2000005
#define inf 0x3f3f3f3f
#define int long long
#define mod 1004535809
#define local
template <typename Tp> void read(Tp &x){
	int fh=1;char c=getchar();x=0;
	while(c>'9'||c<'0'){if(c=='-'){fh=-1;}c=getchar();}
	while(c>='0'&&c<='9'){x=(x<<1)+(x<<3)+(c&15);x%=mod;c=getchar();}x*=fh;
}
int ksm(int B,int P,int Mod){int ret=1;while(P){if(P&1)ret=1ll*ret*B%Mod;B=1ll*B*B%Mod;P>>=1;}return ret;}

namespace Poly{
	const int Gmod=3,invG=334845270;
	int a[maxn],b[maxn];
	int tr[21][maxn];
	int Wn[2][21];
	int LG[maxn];
	int inv[maxn];
	void preprocess(int maxN){
		for(int i=0;(1<<i)<=maxN*2;i++)LG[1<<i]=i;
		for(int i=0;(1<<i)<=maxN*2;i++){
			Wn[0][i]=ksm(invG,(mod-1)/(1<<(i+1)),mod);
			Wn[1][i]=ksm(Gmod,(mod-1)/(1<<(i+1)),mod);
			int N=(1<<i);
			for(int j=0;j<N;j++)tr[i][j]=((tr[i][j>>1]>>1)|((j&1)?(N>>1):0));
		}
		inv[1]=1;
		for(int i=2;i<=maxN*2;i++)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
	}
	void ntt(int *f,int N,int typ){
		for(int i=0;i<N;i++)if(i<tr[LG[N]][i])swap(f[i],f[tr[LG[N]][i]]);
		for(int len=1;len<N;len<<=1){
			int wn=Wn[typ==1][LG[len]];
			for(int i=0;i<N;i+=(len<<1)){
				int buf=1;
				for(int j=0;j<len;j++,buf=1ll*buf*wn%mod){
					int FL=f[i+j],FR=1ll*f[i+j+len]*buf%mod;
					f[i+j]=(FL+FR)%mod;
					f[i+j+len]=(FL-FR)%mod;
				}
			}
		}
		if(typ==-1){
			int invN=inv[N];
			for(int i=0;i<N;i++)f[i]=1ll*f[i]*invN%mod;
		}
		for(int i=0;i<N;i++)f[i]=(f[i]+mod)%mod;
	}
	void mul(int *ff,int *gg,int *ans,int n,int m,int mod_x){
		int N;
		for(N=1;N<n+m-1;N<<=1);
		for(int i=0;i<n;i++)a[i]=ff[i];for(int i=n;i<N;i++)a[i]=0;
		for(int i=0;i<m;i++)b[i]=gg[i];for(int i=m;i<N;i++)b[i]=0;
		ntt(a,N,1);ntt(b,N,1);
		for(int i=0;i<N;i++)a[i]=1ll*a[i]*b[i]%mod;
		ntt(a,N,-1);
		for(int i=0;i<mod_x;i++)ans[i]=a[i];
		for(int i=mod_x;i<N;i++)ans[i%n]=(ans[i%n]+a[i])%mod;//这里即为模意义下的的特殊之处，对应引入例题4
		for(int i=mod_x;i<N;i++)ans[i]=0;
	}//多项式问题要注意及时清零，否则可能会有一些比较难调试出来的错误
}

int n,m;
int P,S;
int h[maxn];
int f[35][maxn],g[maxn];
void Ksm(int p){
	g[0]=1;
	for(int i=32;~i;i--){
		if((p>>i)&1)
		Poly::mul(f[i],g,g,m-1,m-1,m-1);
	}
}//二进制拆分，类比于快速幂
int dtol[maxn],ltod[maxn];
//dtol指真数到对数的映射
//ltod指对数到真数的映射（也是指数到幂的映射）
//本题中ltod没有太大作用
bool check(int gg,int x){
	int tmp[8005]={0};
	for(int i=0,tep=1;i<x-1;i++,tep=1ll*tep*gg%x){
		tmp[tep]++;
		if(tmp[tep]>1)return 0;
	}
	return 1;
}//这个是最暴力的原根判定方法，但因为一个质数最小的原根普遍较小，且本题 m 的范围也比较小，因此暴力判断足矣
void get_G(int x){
	int GGG;
	for(int i=2;i<x;i++){
		if(check(i,x)){
			GGG=i;break;
		}
	}
	for(int i=0,tep=1;i<x-1;i++,tep=1ll*tep*GGG%x){
		dtol[tep]=i;ltod[i]=tep;
	}//构造对数表
}
signed main(){
	int X;
	read(n);read(m);Poly::preprocess(m<<2);
	get_G(m);
	read(X);read(S);
	for(int i=1,a;i<=S;i++){
		read(a);
		if(!a)continue;h[dtol[a]]++;
	}
	for(int i=0;i<m;i++)f[0][i]=h[i];
	for(int i=1;i<=32;i++){
		Poly::mul(f[i-1],f[i-1],f[i],m-1,m-1,m-1);
	}//倍增求出f（即为分析中的ht）
	Ksm(n);
	printf("%lld
",g[dtol[X]]);
	return 0;
}