Description
(Claris) 和 (NanoApe) 在玩石子游戏，他们有 (n) 堆石子，规则如下：
- (Claris) 和 (NanoApe) 两个人轮流拿石子，(Claris) 先拿。
- 每次只能从一堆中取若干个，可将一堆全取走，但不可不取，拿到最后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) 变回来
大概的原理就是
最后，写代码时注意边界。
(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;
}