bzoj4589 Hard Nim

http://www.lydsy.com/JudgeOnline/problem.php?id=4589

题意：石子堆数为N且每堆石子的数量都是不大于M的质数的Nim游戏，求先手必败的局面数量模10^9+7。N<=10^18，M<=50,000，数据组数<=80。

也即，有依次N个不大于M的质数，求异或和为0的方案数。

考虑暴力dp，令f[i][j]为前i堆石子异或和为j的方案数，则f[0][0]=1，f[0][j]=0(j>0)，f[i][j]=sum{f[i-1][k]+g[j xor k]}(i>0)。时间复杂度O(NM)，TLE。

上述式子中，当且仅当k为不大于M的质数时g[k]=1否则g[k]=0。

考虑一种特殊的基于异或的多项式乘法，一般的多项式乘法a*b=c是a[i]*b[j]对c[i+j]有贡献，也即c[i]=sum{a[j]*b[i-j]}。

这里比较特殊，a[i]*b[j]对c[i xor j]有贡献，也即c[i]=sum{a[j]*b[i xor j]}。

容易发现，答案就是(g^n)[0]。注意到对于任意x>0，f[x][0]就等于(g^x)[]，用归纳法易证。

这种特殊的多项式乘法显然是满足结合律的，于是我们计算乘方的时候可以用快速幂来做，而乘法是暴力乘。时间复杂度O(M^2logN)，还是TLE。

然后我们来讲讲FWT，它是一个可以在O(NlogN)的时间复杂度内计算规模为N的基于异或的多项式乘法。

考虑N=2的情况，此时下标只有0和1，我们用(x,y)来表示a[0]=x且a[1]=y的多项式a。

我们考虑FFT的思想，构造一种特殊的变换trans，使得(trans(a*b))[i]=(trans(a))[i]*(trans(b))[i]。

这样就可以很容易的计算出trans(ans)，再用逆变换得到ans即可。

这里当N=2时，变换即为trans((x,y))=(x-y,x+y)，同时trans((z,w))=(z-w,z+w)。另外，(x,y)*(z,w)=(xz+yw,xw+yz)。

验证一下，trans((xz+yw,xw+yz))=(xz+yw-xw-yz,xz+yw+xw+yz)=((x-y)*(z-w),(x+y)*(z+w))，无误。

逆变换应该是trans'((x,y))=((x+y)/2,(y-x)/2)，注意这里做除法的时候应该乘以逆元。

考虑更大的N=8，下标为0~7。变换的时候，先对01,23,45,67做，再对02,13,46,57做，最后对04,15,26,37做。逆变换把顺序反过来就好了。

事实上写成代码的话FWT的结构和FFT很像。

于是我们照样快速幂，乘法的时候用FWT，时间复杂度O(MlogMlogN)，似乎还是会TLE。

注意到我们每次乘法完的时候都将其逆变换，然后进行下一次乘法的时候又将其变换，容易发现中间的这些变换都是多余的。

也就是说，我们只需要对g[]变换一下，然后对每一项求快速幂，再逆变换就行了。时间复杂度只有O(MlogM+MlogN)。

代码：

#include<cstdio>
#define ll long long
#define N 66666
const int mo=1000000007;
const int iv=mo+1>>1;
void fwt(int*a,int n){
int i,j,k,x;
for(k=2;k<=n;k<<=1){
for(i=0;i<n;i+=k){
for(j=i;j<i+(k>>1);j++){
x=a[j]+a[j+(k>>1)];if(x>=mo) x-=mo;
a[j]=a[j]-a[j+(k>>1)];if(a[j]<0) a[j]+=mo;
a[j+(k>>1)]=x;
}
}
}
}
void twf(int*a,int n){
int i,j,k,x;
for(k=n;k>=2;k>>=1){
for(i=0;i<n;i+=k){
for(j=i;j<i+(k>>1);j++){
x=a[j]+a[j+(k>>1)];
a[j+(k>>1)]=(int)(1LL*(a[j+(k>>1)]-a[j]+mo)*iv%mo);
a[j]=(int)(1LL*x*iv%mo);
}
}
}
}
ll pw(ll x,ll y){
if(!y) return 1;
ll temp=pw(x,y>>1);
(temp*=temp)%=mo;
if(y&1) (temp*=x)%=mo;
return temp;
}
int pr[N],i,j,m,len,a[N];ll n;
int main(){
pr[1]=0;for(i=2;i<=50000;i++) pr[i]=1;
for(i=2;i<=50000;i++) if(pr[i]){
for(j=2;i*j<=50000;j++) pr[i*j]=0;
}
while(scanf("%lld%d",&n,&m)==2){
len=1;while(len<=m) len<<=1;
for(i=0;i<len;i++) a[i]=pr[i]&(i<=m);
fwt(a,len);
for(i=0;i<len;i++) a[i]=(int)pw(a[i],n);
twf(a,len);
printf("%d\n",a[0]);
}
}