[bzoj4589] Hard Nim

Description

(Claris) 和 (NanoApe) 在玩石子游戏，他们有 (n) 堆石子，规则如下：

1. (Claris) 和 (NanoApe) 两个人轮流拿石子，(Claris) 先拿。
2. 每次只能从一堆中取若干个，可将一堆全取走，但不可不取，拿到最后1颗石子的人获胜。

不同的初始局面，决定了最终的获胜者，有些局面下先拿的 (Claris) 会赢，其余的局面 (Claris) 会负。
(Claris) 很好奇，如果这 (n) 堆石子满足每堆石子的初始数量是不超过 (m) 的质数，而且他们都会按照最优策略玩游戏，那么 (NanoApe) 能获胜的局面有多少种。
由于答案可能很大，你只需要给出答案对 (10^9+7) 取模的值。

Input

输入文件包含多组数据，以 (EOF) 为结尾。
对于每组数据：
共一行两个正整数 (n) 和 (m) 。
每组数据有 (1 leq n leq 10^9, 2 leq m leq 50000) 。
不超过80组数据。

Output

对于每组数据，输出 (NanoApe) 能获胜的局面有多少种。

Sample Input

3 7
4 13

Sample Output

6
120

---

想法

首先，这是一个最基本的 (Nim) 游戏，每个数的 (sg) 值都等于它本身。
根据 (SG) 定理，每堆石子 (sg) 值的异或和若为0，则先手输，不为0则先手赢。
这道题中，我们要求的就是 (n) 堆石子的排列，每堆个数都是 (m) 以内的质数，求它们异或和为0（后手胜）的方案数。

于是我们考虑 (dp)
(f[i][j]) 表示现在有 (i) 堆石子，它们的异或和为 (j) 的方案数。
那么转移就是 (f[i][j]=sumlimits_{j=k oplus p,p为m以内质数} f[i-1][k])
就用 (oplus) 表示异或吧…

看到这里，我想到了矩阵乘法。
但它的复杂度是 (O(logn imes s^3)) ，其中 (s) 表示异或和的最大值，是 (1e5) 级别的，不行
似乎还可以倍增优化转移，但我认真地不明白怎么搞 (QwQ)

但是看到 (j=k oplus p) ，可以想到我昨天刚刚学的 (FWTxor)
把式子改写一下，令 (P[i]) 表示 (i) 是否为质数，是的话为1
则 (f[i][j]=sumlimits_{j=k oplus p} f[i-1][k] imes P[p])
这就是标准的 (xor) 卷积形式了。

但还有一个问题，(n) 是 (1e9) 级别的，不能一个个转移
那就快速幂咯~
先正向 (fwt) ，然后快速幂，最后反向 (fwt) 变回来
大概的原理就是

[egin{aligned} FWT(A oplus B)&=FWT(A) imes FWT(B) Rightarrow \ FWT(A oplus B oplus C)&=FWT((A oplus B) oplus C) \ &=FWT(A oplus B) imes FWT(C) \ &=FWT(A) imes FWT(B) imes FWT(C) end{aligned} ]

最后，写代码时注意边界。
(f[0][0]=1) 或者 (f[1][p]=1,p为m以内质数)

---

代码

#include<cstdio>
#include<iostream>
#include<algorithm>

#define xzy 1000000007
#define ljh 500000004

using namespace std;

const int N = 530005;

int p[N],prime[N],pnum;
void getp(){
	for(int i=2;i<N;i++) p[i]=1;
	for(int i=2;i<N;i++){
		if(p[i]) prime[pnum++]=i;
		for(int j=0;j<pnum && 1ll*prime[j]*i<N;j++){
			p[i*prime[j]]=0;
			if(i%prime[j]==0) break;
		}
	}
}

int w;
void fwt(int *A,int ty){
	for(int i=2;i<=w;i<<=1)
		for(int j=0;j<w;j+=i)
			for(int k=j;k<j+i/2;k++){
				int t=A[k+i/2];
				A[k+i/2]=(A[k]-t+xzy)%xzy;
				A[k]=(A[k]+t)%xzy;
				if(ty==-1) {
					A[k+i/2]=1ll*A[k+i/2]*ljh%xzy;
					A[k]=1ll*A[k]*ljh%xzy;
				}
			}
}

int n,m;
int a[N],b[N];

int main()
{
	getp();
	while(~scanf("%d%d",&n,&m)){
		w=1;
		while(w<=m) w<<=1;
		for(int i=0;i<w;i++) a[i]=0; a[0]=1;
		for(int i=0;i<=m;i++) b[i]=p[i];
		for(int i=m+1;i<w;i++) b[i]=0;
		fwt(a,1); fwt(b,1);
		while(n){
			if(n&1) for(int i=0;i<w;i++) a[i]=1ll*a[i]*b[i]%xzy;
			for(int i=0;i<w;i++) b[i]=1ll*b[i]*b[i]%xzy;
			n>>=1;
		}
		fwt(a,-1);
		printf("%d
",a[0]);
	}
	
	return 0;
}