题目大意
给定一个集合\(S\),\(\forall x\in S,s\in [0,m-1]\),求从这个集合\(S\)中选取\(n\)个数,使得乘积为\(x\)的方案数。
算法分析
这其实是一类非常经典的问题。
引入例题1
给定一个集合\(S\),\(\forall x\in S,s\in [1,n]\),在这个集合中选取2个数(可以重复),使得加和为\(A\)的方案数。
不难想到可以从数值角度入手,我们运用桶思想,设在集合中\(i\)这个数字出现了\(h(i)\)次,则不难得到答案为
显然这个可以\(O(n)\)求解,具体实现不再赘述。
当然这个式子符合多项式卷积的一般形式,所以也可以\(O(n\log n)\)求解。(不过一道普及组难度的题目被搞成这样好吗)
引入例题2
给定一个集合\(S\),\(\forall x\in S,s\in [1,n]\),在这个集合中选取3个数(可以重复),使得加和为\(A\)的方案数。
类比于上一道例题,不难得到答案为
时间复杂度是\(O(n^2)\),似乎不是很好。
这个能不能优化呢?
选取3个数可以看作在选取1个数的集合\(S_1\)和选取2个数的集合\(S_2\)中分别选取1个数。
我们可以先求出选取2个数的方案,记从\(S\)中选取2个数加和为\(i\)的方案数为\(h_2(i)\),则不难得到答案为
那如何快速求出\(h_2\)呢?
由上一道例题可得
这是标准的卷积形式,可以使用FFT/NTT
快速求出\(h_2\),时间复杂度为\(O(n\log n)\)。
推荐题目:
【bzoj3513】[MUTC2013]idiots
【bzoj3771】Triple
引入例题3
给定一个集合\(S\),\(\forall x\in S,s\in [1,n]\),在这个集合中选取\(k\)个数(可以重复),使得加和为\(A\)的方案数。
类比于例题2,不难得到一种方法:设\(h_t(i)\)表示从集合中选取\(t\)个数加和为\(i\)的方案数,则有如下的递推式:
这样时间复杂度是\(O(km\log m)\),还有一定的优化空间。
考虑倍增,设\(h_t(i)\)表示从集合中选取\(2^t\)个数加和为\(i\)的方案数。选取\(2^t\)个数相当于在选取\(2^{t-1}\)个数的集合中选取2个数,则不难得到如下的递推式:
这样预处理出\(h_t\),把\(k\)二进制拆解,类比于快速幂,可以在\(O(m\log m\log k)\)的时间复杂度下实现。
核心代码如下:
ans[0]=1;//初始化
for(int i=30;~i;i--){
if((k>>i)&1){
mul(h[i],ans,ans);
}
}
引入例题4
给定一个集合\(S\),\(\forall x\in S,s\in [0,n-1]\),在这个集合中选取\(k\)个数(可以重复),使得加和在\(\bmod \ 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,a\ge 0,a\neq1\),根据幂运算性质,不难得到:\(f(x+y)=f(x)\times f(y)\)。
取对数运算,对于\(a^x=N\),则有\(\log_aN=x\),根据\(a^x\times a^y=a^{x+y}\),不难得出\(\log_aa^x+\log_aa^y=\log_aa^{x+y}=\log_a(a^x\times a^y)\),推广可得\(\log_ax+\log_ay=\log_a(x\times y)\)。
对数函数,形如\(f(x)=\log_ax,a\ge 0,a\neq 1,x>0\),根据对数运算相关性质,不难得到\(f(xy)=f(x)+f(y)\)。
指数函数可以使加法转化为乘法,对数函数可以使乘法转化为加法。
但是指对运算都是在实数域上的,而本题是在模意义下的,怎么处理呢?
模意义下的指数/对数?
求指数还是很容易的,随便以一个数为底(设为\(a\)),在模意义下求出\(a^x,x\in[1,m-1]\),那么根据对数的定义,则可令\(\log_a(a^x\bmod m)=x\),这样建立映射关系。
看上去似乎没有问题?
我们举一个例子,如果\(m=7,a=2\),那么\(a^1 \equiv a^4\equiv 2 \ (\bmod\ m)\),那么\(\log_a(2)=?\)
我们发现这种映射关系必须是一一对应的,因此并不是随便取一个数为底就可以的。
我们需要找到一个底数\(a\),使得\(a^1,a^2,\cdots,a^{m-1}\)在\(\bmod \ m\)意义下互不相同。
注意:这里没有包含\(a^0\),这是因为在\(\bmod \ m\)意义下不存在一个数\(x\)使得\(a^x\equiv 0(\bmod \ m)\),因此实际上只有\(m-1\)个可对应位置,在本题中\(S\)集合中的\(0\)对答案无贡献,故忽略\(a^0\)。
这个东西似乎很熟悉啊,这好像就是原根。
对于\(p\)的原根\(g\),满足\(g^1,g^2,\cdots,g^{p-1}\)在\(\bmod \ p\)意义下互不相同。
因此我们只需要找到\(m\)的一个原根,即可把乘法转化为加法。
综上,我们得到这一题的解决方法:
- 求出\(m\)的原根,并把集合\(S\)中的数\(x\)转化为\(\log_gx\)。
- 求出\(h_t\),含义与上面引入例题的类似,其中\(h_0(i)\)表示集合\(S\)中满足\(\log_gx=i\)数的数量。
- 将\(n\)进行二进制拆分,并做多项式卷积,求出在去过对数的集合中选取\(n\)个数加和为\(i\)的方案数为\(ans(i)\)。
- 对于输入的\(X\),答案即为\(ans(\log_gX)\)。
时间复杂度为\(O(m\log m\log 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\n",g[dtol[X]]);
return 0;
}