Project Euler 550

发表于 更新于 分类于 阅读次数： 本文字数： 944 阅读时长 ≈ 1 分钟

Project Euler 550

题目

Divisor game

Two players are playing a game. There are $k$ piles of stones. When it is his turn a player has to choose a pile and replace it by two piles of stones under the following two conditions:

  • Both new piles must have a number of stones more than one and less than the number of stones of the original pile.
  • The number of stones of each of the new piles must be a divisor of the number of stones of the original pile.

The first player unable to make a valid move loses.

Let $f(n,k)$ be the number of winning positions for the first player, assuming perfect play, when the game is played with k piles each having between $2$ and $n$ stones (inclusively).$f(10,5)=40085$.

Find $f(10^7,10^{12})$.Give your answer modulo $987654321$.

解决方案

令$N=10^7,Q=10^{12}$。

可以看出这是一个公平组合游戏，我们将每一堆石头视为独立的一个游戏局面，使用SG定理来解决本题。这里的石头一共有$k$堆，第$i$堆一共有$a_i$个，那么整个游戏的$sg$函数值为$sg(a_1)\oplus sg(a_2)\oplus\dots sg(a_k)$。

考虑目前的某一堆石头，一共有$n$个，每一次操作是用两堆石头$d_1,d_2$来代替它，这两堆石头是两个独立的新局面。因此，可以写出一堆石头的$sg$函数定义：

其中运算$\text{mex}(s)$表示集合$s$中未出现的最小的非负整数。

可以发现，$1\sim N$中的$sg$函数值都很小（不超过$2^6$）。用一个数组$c[i]$表示$1\sim N$中有多少个数的$sg$函数值为$i$。

不过在实现的时候，直接枚举$1\sim N$中的所有因数再进行二重循环的效率是很低的。这里使用了一个优化方法：

对于$n$的分解$n=\prod_{i=1}^k p_i^{e_i}$而言，$sg(n)$只和无序多重集${e_i}$有关，和$p_i$无关。因此，将$1\sim N$中的每个数分解质因数后对集合${e_i}$进行哈希。只要哈希值相同，那么$sg$函数值一定相同。如果发现哈希值是第一次生成，那么才按照上面的二重循环计算$sg$函数值。经过实现，在这个数量级下，只有约$500$个不同的哈希值。

最终问题转化为有多少种方法可以将$sg$函数值填入一个长度为$Q$的序列，并且这$Q$个值的异或和不为$0$。

令状态$g(i,j)(0\le i\le Q,0\le j<2^6)$表示填入长度为$i$的序列，并且异或和为$j$的方案数有多少个。那么不难写出状态转移方程为：

但是，对于$Q$这种数据范围而言，直接进行转移效率非常低。

如果我们已经求出了长度为$a$的序列的各种填法$g(a,\cdot)$和长度为$b$的各种填法$g(b,\cdot)$，不难知道我们可以通过卷积组合出$g(a+b,\cdot)$。

枚举所有的$g(a,i)$和$g(b,j)$，将这两部分值直接合并起来：

$g(a,i)\cdot g(b,j)\rightarrow g(a+b,i\oplus j)$

因此，这给了我们一个方案：依次求出$g(2^0,\cdot),g(2^1,\cdot),g(2^2,\cdot),\dots$。然后针对$Q$，选择这些求出的$g(2^i,\cdot)$进行合并即可。

最终答案为$\sum_{i=1}^{2^6-1} g(Q,i)$。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <bits/stdc++.h>
# include "tools.h"
using namespace std;
typedef long long ll;
const int N=10000000;
const ll Q=1e12;
ll mod=987654321;
int B=log2(N+1)+1;
const int M=1<<7;

int v[N+4],pr[N/10+100],m=0;
int sg[N+4];
bool mex[M];
int p[10],e[10],e2[10],o=0;
unordered_map<ll,int>sg_mem;
void fact(int n){
    o=0;
    for(;n!=1;){
        int x=v[n];
        p[++o]=x;
        e[o]=0;
        for(;n%x==0;n/=x,++e[o]);
    }
}
ll a[M],b[M];
void cal(ll a[M],ll b[M],ll ans[M]){
    ll c[M]={0};
    for(int i=0;i<M;i++)
        for(int j=0;j<M;j++)
            c[i^j]=(c[i^j]+a[i]*b[j])%mod;
    memcpy(ans,c,sizeof(c));
}
int main(){
    for(int i=2;i<=N;i++){
        if(v[i]==0){
            v[i]=i;pr[++m]=i;
        }
        for(int j=1;j<=m;j++){
            if(pr[j]>v[i]||pr[j]>N/i) break;
            v[i*pr[j]]=pr[j];
        }
    }
    for(int i=2;i<=N;i++){
        fact(i);
        memcpy(e2,e,sizeof(e));
        sort(e2+1,e2+o+1);
        ll hs=0;
        for(int j=1;j<=o;j++)
            hs=hs*B+e2[j];
        if(sg_mem.count(hs))
            sg[i]=sg_mem[hs];
        else{
            vector<pl>fact_v;
            for(int j=1;j<=o;j++)
                fact_v.push_back(pl(p[j],e[j]));
            vector<ll>d=divisors(fact_v);
            d.pop_back();
            d.erase(d.begin());
            memset(mex,0,sizeof(mex));
            for(int i=0;i<d.size();i++)
                for(int j=i;j<d.size();j++)
                    mex[sg[d[i]]^sg[d[j]]]=1;
            int j=0;
            for(;mex[j];++j);
            sg[i]=sg_mem[hs]=j;
        }
        ++b[sg[i]];

    }
    a[0]=1;
    for(ll n=Q;n;n>>=1){
        if(n&1) cal(a,b,a);
        cal(b,b,b);
    }
    ll ans=0;
    for(int j=1;j<M;j++)
        ans=(ans+a[j])%mod;
    printf("%lld\n",ans);
}