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

最新推荐文章于 2025-11-25 00:06:26 发布

原创最新推荐文章于 2025-11-25 00:06:26 发布 · 321 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#算法

该博客探讨了一个数学问题：一个人从原点出发，拥有若干整数，每次可以按这些整数的值向左或向右走，求到达点1的概率。博主利用裴蜀定理和数论方法，通过计算满足条件的方案数，得出概率的表达式，并通过数论分块和杜教筛优化算法，实现了高效计算。最后，博主提供了C++代码实现整个计算过程。

题意：

一开始一个人在原点，它拥有 $n$ 个整数 $x_{i}$ ，并且 $xi∈[1,m]x_{i}\in[1,m]$ ，他每次可以选择一个 $x_{i}$ ，来向左走 $x_{i}$ 或向右走 $x_{i}$ 步，问有多大概率能够走到点 $1$ ？设答案为 $w$ ，只需要输出 $w⋅mnw\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可以直接赋值负数，也会有取余的效果

时间复杂度：

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

// #include<bits/stdc++.h>
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<vector>
#include<bitset>
#include<map>
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;
}