Codeforces Round #548 (Div. 2) 1139 D+2021天梯赛l3-3 解题报告（负二项式分布+莫比乌斯容斥+杜教筛（天梯赛））

最新推荐文章于 2022-06-23 22:46:43 发布

nagisa-kun

最新推荐文章于 2022-06-23 22:46:43 发布

阅读量304

点赞数

分类专栏：数论文章标签：莫比乌斯反演 gcd 概率论

本文链接：https://blog.csdn.net/nagisa2019/article/details/116266651

版权

数论专栏收录该内容

25 篇文章 2 订阅

订阅专栏

以下是我这道题的思考过程：
（1）、若选到1，那直接就能结束了。
（2）、若选到其他的数，则我们需要继续选下去，直到选到一个和前面的一个数互质的数。想到互质，可以想到和质数相关。
思考极限情况，我们可以想象这一过程：不停地选数，在这个过程在，我们保持选到的数的都是某一个数x的倍数，直到选到一个数不为x的倍数就结束。
这一过程其实就是概率论中的负二项式分布:

实验包含一系列独立的实验。
每个实验都有成功、失败两种结果。
成功的概率p是恒定的。
实验持续到r次失败，r可以为任意正数。

则X~NB(r,p)，并且 $E(X)={rp\over 1-p}$
(摘自百度百科)
因此，我们可以枚举x，算出p，再对所有的E求和就是答案了。
但是我们发现如果我们令x取遍1～m会存在重复，因为m的倍数和m的因数的倍数之间存在重复，联想之前的质数，我们容易联想到我们可以只枚举质数。但是，仍然有重复，如6同时是2和3的倍数，因此还需要容斥。根据容斥的公式，一个数若含有奇数个次数为1的质因子则+E，若含有偶数个次数为1的因子则-E，这就和莫比乌斯函数有关。
因此，公式如下：
$ANS=1+\sum_{i=2}^m\mu(i)E(i)$
$E(i)={p\over{1-p}}={{{\lfloor{m\over i}\rfloor}\over m}\over{1-{\lfloor{m\over i}\rfloor\over m}}}$
cf的直接套上去就能过了。

#include <bits/stdc++.h>
#define inf 0x3f3f3f3f
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
#define rep(i, a, n) for(int i = a; i <= n; ++ i)
#define per(i, a, n) for(int i = n; i >= a; -- i)
//#define ONLINE_JUDGE
using namespace std;
typedef long long ll;
const int mod=1e9+7;
template<typename T>void write(T x)
{
    if(x<0)
    {
        putchar('-');
        x=-x;
    }
    if(x>9)
    {
        write(x/10);
    }
    putchar(x%10+'0');
}

template<typename T> void read(T &x)
{
    x = 0;char ch = getchar();ll f = 1;
    while(!isdigit(ch)){if(ch == '-')f*=-1;ch=getchar();}
    while(isdigit(ch)){x = x*10+ch-48;ch=getchar();}x*=f;
}

int gcd(int a,int b){return b==0?a:gcd(b,a%b);}
int lcm(int a,int b){return a/gcd(a,b)*b;};
inline int ksm(int a,int n){
	int ans=1;
	while(n){
		if(n&1) ans=(1ll*ans*a)%mod;
		a=1ll*a*a%mod;
		n>>=1;
	}
	return ans%mod;
}
//==============================================================
const int maxn=1e5+10;
int m;
bool vis[maxn];
int pri[maxn],tot;
int mu[maxn];
void init(){
	mu[1]=1;
	for(int i=2;i<maxn;++i){
		if(!vis[i])pri[++tot]=i,mu[i]=mod-1;
		for(int j=1;j<=tot&&i*pri[j]<maxn;++j){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0)break;
			mu[i*pri[j]]=mod-mu[i];
		}
	}
}

int solve(){
	int res=1;
	for(int i=2;i<=m;++i){
		int p=1ll*(m/i)*ksm(m,mod-2)%mod;
		res+=1ll*(mod-mu[i])*p%mod*ksm(1-p,mod-2)%mod;
		res=(res%mod+mod)%mod;
	}
	return res;
}

signed main()
{
	#ifndef ONLINE_JUDGE
	freopen("in.txt","r",stdin);
	freopen("out.txt","w",stdout);
	#endif
	//clock_t c1 = clock();
	//===========================================================
	read(m);
	init();
	cout<<solve()<<endl;
	//===========================================================
	//std::cerr << "Time:" << clock() - c1 << "ms" << std::endl;
	return 0;
}

天体赛的话还需要继续优化。可以看到数据范围高达1e11，八成是要亚线性筛。对于上式，最先想到的应该是数论分块，而还需要一个莫比乌斯函数前缀和，因此，上杜教筛。记得开int128。

#include <bits/stdc++.h>
#define inf 0x3f3f3f3f
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
#define rep(i, a, n) for(int i = a; i <= n; ++ i)
#define per(i, a, n) for(int i = n; i >= a; -- i)
//#define ONLINE_JUDGE
typedef long long ll;
typedef __int128 i128;
#define int ll
ll mod=1e9+7;
using namespace std;
template<typename T>void write(T x)
{
    if(x<0)
    {
        putchar('-');
        x=-x;
    }
    if(x>9)
    {
        write(x/10);
    }
    putchar(x%10+'0');
}

template<typename T> void read(T &x)
{
    x = 0;char ch = getchar();int f = 1;
    while(!isdigit(ch)){if(ch == '-')f*=-1;ch=getchar();}
    while(isdigit(ch)){x = x*10+ch-48;ch=getchar();}x*=f;
}

int gcd(int a,int b){return b==0?a:gcd(b,a%b);}
int lcm(int a,int b){return a/gcd(a,b)*b;};
inline i128 ksm(i128 a,i128 n){
	i128 ans=1;
	while(n){
		if(n&1) ans=((i128)ans*a)%mod;
		a=(i128)a*a%mod;
		n>>=1;
	}
	return ans%mod;
}
//==============================================================
ll n;
inline i128 mul(i128 a,i128 b){
	return (i128)a*b%mod;
}

inline i128 add(i128 a,i128 b){
	return a+b>=mod?a+b-mod:a+b;
}

inline i128 sub(i128 a,i128 b){
	return (a-b+mod)>=mod?a-b:a-b+mod;
}

const int maxn=2e7+10;
bool vis[maxn];
int pri[maxn],tot;
int mu[maxn];
inline void init(){
	mu[1]=1;
	for(int i=2;i<maxn;++i){
		if(!vis[i])pri[++tot]=i,mu[i]=mod-1;
		for(i128 j=1;j<=tot&&i*pri[j]<maxn;++j){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){
				break;
			}
			mu[i*pri[j]]=mod-mu[i];
		}
	}
	for(i128 i=1;i<maxn;++i){
		mu[i]=add(mu[i],mu[i-1]);
	}
}
unordered_map<ll,i128>dpmu;
i128 SumMu(ll n){
    if(n<maxn)return mu[n];
    if(dpmu.count(n))return dpmu[n];
    ll res=1;
    for(ll l=2,r;l<=n;l=r+1){
        ll t=n/l;
        r=n/t;
		res=sub(res,mul(r-l+1,SumMu(n/l)));
        //res=(res-(r-l+1)*SumMu(n/l)%mod+mod)%mod;
    }
    return dpmu[n]=res;
}
inline i128 solve(){
	i128 res=1;
	i128 pre=SumMu(1);
	for(ll l=2,r;l<=n;l=r+1){
		i128 now=n/l;
		r=n/now;
		i128 pos=SumMu(r);
		i128 p=mul((n/l),ksm(n,mod-2));
		res=add(res,mul(mul(((mod-(pos-pre))),p),ksm(1-p,mod-2)));
		pre=pos;
	}
	return (res%mod+mod)%mod;
}

signed main()
{
	#ifndef ONLINE_JUDGE
	freopen("in.txt","r",stdin);
	freopen("out.txt","w",stdout);
	#endif
	//clock_t c1 = clock();
	//===========================================================
	read(n),read(mod);
	init();
	write(solve());
	//===========================================================
	//std::cerr << "Time:" << clock() - c1 << "ms" << std::endl;
	return 0;
}