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

题意：

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

Solution：

走到 1 1 1即等价于以下方程有解
∑ i = 1 n a i x i = 1 , a i ∈ Z \sum_{i=1}^{n}a_{i}x_{i}=1,a_{i}\in Z i=1∑naixi=1,ai∈Z
又裴蜀定理，只需要满足下列等式就是一种可以走到的方案
g c d ( ∣ a i ∣ , ∣ a 2 ∣ , ∣ a 3 ∣ , . . . , ∣ a n ∣ ) = 1 gcd(|a_{i}|,|a_{2}|,|a_{3}|,...,|a_{n}|)=1 gcd(∣ai∣,∣a2∣,∣a3∣,...,∣an∣)=1
令 a i = ∣ a i ∣ a_{i}=|a_{i}| ai=∣ai∣，计算符合的方案总数
∑ a 1 = 1 n ∑ a 2 = 1 n . . . ∑ a n = 1 n [ g c d ( a 1 , a 2 , . . . , a n ) = 1 ] \sum_{a_{1}=1}^{n}\sum_{a_{2}=1}^{n}...\sum_{a_{n}=1}^{n}[gcd(a_{1},a_{2},...,a_{n})=1] a1=1∑na2=1∑n...an=1∑n[gcd(a1,a2,...,an)=1]
莫比乌斯反演即计算
∑ a 1 = 1 n ∑ a 2 = 1 n . . . ∑ a n = 1 n ∑ d ∣ g c d ( a 1 , a 2 , . . . , a n ) = 1 μ ( d ) \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) a1=1∑na2=1∑n...an=1∑nd∣gcd(a1,a2,...,an)=1∑μ(d)
优先枚举因子 d d d，就化为了
∑ d = 1 n μ ( d ) ⌊ m d ⌋ n \sum_{d=1}^{n}\mu(d)\lfloor \frac{m}{d} \rfloor^n d=1∑nμ(d)⌊dm⌋n
右项可以数论分块，而左项在数论分块时需要区间求和，转化为前缀和相减，前缀和只需要用杜教筛即可

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

时间复杂度：

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

// #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 stdjudgefreopen("in.txt","r",stdin);#endifcin>>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;
}