[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]=\sum\limits_{j=k \oplus p,p为m以内质数} f[i-1][k]\)
就用 \(\oplus\) 表示异或吧…

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

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

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

\[\begin{aligned} FWT(A \oplus B)&=FWT(A) \times FWT(B) \Rightarrow \\ FWT(A \oplus B \oplus C)&=FWT((A \oplus B) \oplus C) \\ &=FWT(A \oplus B) \times FWT(C) \\ &=FWT(A) \times FWT(B) \times 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\n",a[0]);
	}
	
	return 0;
}