牛客P19836 裴蜀定理+莫比乌斯反演+杜教筛

题意：

一开始一个人在原点，它拥有$n$个整数$x_i$，并且$x_i\in [1,m]$，他每次可以选择一个$x_i$，来向左走$x_i$或向右走$x_i$步，问有多大概率能够走到点$1$？设答案为$w$，只需要输出$w\cdot m^n$取余$2^{64}$的结果

Solution：

走到$1$即等价于以下方程有解
$$\sum _{i=1}^n{a}_i{x}_i=1,a_i\in Z$$
又裴蜀定理，只需要满足下列等式就是一种可以走到的方案
$$gcd(|a_i|,|a_2|,|a_3|,...,|a_n|)=1$$
令$a_i=|a_i|$，计算符合的方案总数
$$\sum _{a_1=1}^n\sum _{a_2=1}^n...\sum _{a_n=1}^n[gcd(a_1,a_2,...,a_n)=1]$$
莫比乌斯反演即计算
$$\sum _{a_1=1}^n\sum _{a_2=1}^n...\sum _{a_n=1}^n\sum _{d|gcd(a_1,a_2,...,a_n)=1}\mu (d)$$
优先枚举因子$d$，就化为了
$$\sum _{d=1}^n\mu (d)\lfloor \frac{m}{d}{\rfloor }^n$$
右项可以数论分块，而左项在数论分块时需要区间求和，转化为前缀和相减，前缀和只需要用杜教筛即可

最后概率应为
$$\frac{{\sum }_{d=1}^n\mu (d)\lfloor \frac{m}{d}{\rfloor }^n}{m^n}$$
需要输出的恰好把分母乘掉了，取余$2^{64}$次方只需要在unsigned ll下计算就会自动取余了，ull可以直接赋值负数，也会有取余的效果

时间复杂度：

最坏情况下进行$\sqrt{n}$次杜教筛，即$O(n^{\frac{7}{6}})$

// #include
#include
#include
#include
#include
#include
#include
#include
#include
using namespace std;

using ll=long long;
using ull=unsigned long long;
const int N=2e7+5,inf=0x3fffffff;
const long long INF=0x3fffffffffffffff,mod=1e9+9;

ull qpow(ull a,ull b)
{
    ull ret=1,base=a;
    while(b)
    {
        if(b&1) ret=ret*base;
        base=base*base;
        b>>=1;
    }
    return ret;
}

ull n,m,prime[N],cnt,mu[N];
bitset<N>nt;

void make_prime()
{
    mu[1]=1;
    for(ull i=2;i<N;i++)
    {
        if(!nt[i]) prime[++cnt]=i,mu[i]=-1;
        for(ull j=1;j<=cnt&&i*prime[j]<N;j++)
        {
            nt[i*prime[j]]=true;
            if(i%prime[j]==0) break;
            else mu[i*prime[j]]=mu[i]*mu[prime[j]];
        }
    }
    for(int i=1;i<N;i++) mu[i]+=mu[i-1];
}

map<ull,ull>map1;

ull calcmu(ull x)
{
    if(x<N) return mu[x];
    if(map1.count(x)) return map1[x];
    ull ret=1;
    for(ull l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        ret-=(r-l+1)*calcmu(x/l);
    }
    return map1[x]=ret;
}

int main()
{
    #ifdef stdjudge
        freopen("in.txt","r",stdin);
    #endif
    cin>>n>>m;
    make_prime();
    ull ans=0;
    for(ull l=1,r;l<=m;l=r+1)
    {
        r=m/(m/l);
        ans+=(calcmu(r)-calcmu(l-1))*qpow(m/l,n);
    }
    cout<<ans;
    return 0;
}
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