[BZOJ 2820] YY的gcd(莫比乌斯反演)
[BZOJ 2820] YY的gcd(莫比乌斯反演+数论分块)
题面
给定N, M,求1≤x≤N,1≤y≤M且gcd(x, y)为质数的(x, y)有多少对。q组询问
分析
我们要求的是
∑p∈Pn∑i=1m∑j=1[gcd(i,j)=p]
(大写P表示质数集合)
根据kgcd(i,j)=gcd(ki,kj),
原式=∑p∈P⌊n/p⌋∑i=1⌊m/p⌋∑j=1[gcd(i,j)=1]
又根据莫比乌斯反演里的一个常用结论(证明见 BZOJ 2301)
n∑i=1m∑j=1[gcd(i,j)=1]=min(n,m)∑d=1μ(d)⌊nd⌋⌊md⌋
原式=∑p∈Pmin(⌊n/p⌋,⌊m/p⌋)∑d=1μ(d)⌊npd⌋⌊mpd⌋
令T=pd,则d=Tp
改变求和顺序,
原式=min(n,m)∑T=1∑p|t ∩ p∈P⌊nT⌋⌊mT⌋μ(Tp)
=min(n,m)∑T=1)⌊nT⌋⌊mT⌋∑p|t ∩ p∈Pμ(Tp)
令g(n)=∑p|n ∩ p∈Pμ(np)
原式=min(n,m)∑T=1)⌊nT⌋⌊mT⌋g(T)
前面的部分可以数论分块求解,考虑如何快速求出g(T)
对于每个质数p,我们从1开始枚举j,并保证jp≤n,然后用μ(j)更新g(jp)的值。
由于1/1+1/2+1/3+...+1/n=O(logn),每次更新的复杂度是均摊O(logn)的,而1~n的质数约有nlnn个,所以预处理g函数的总时间复杂度为O(n)
总时间复杂度O(n+q√n)
代码
#include<iostream> #include<cstdio> #include<algorithm> #include<cstring> #define maxn 10000000 using namespace std; typedef long long ll; int t,n,m; int cnt; int prime[maxn+5]; bool vis[maxn+5]; int mu[maxn+5]; ll g[maxn+5]; ll sumg[maxn+5]; void sieve(int n){ mu[1]=1; for(int i=2;i<=n;i++){ if(!vis[i]){ mu[i]=-1; prime[++cnt]=i; } for(int j=1;j<=cnt&&i*prime[j]<=n;j++){ vis[i*prime[j]]=1; if(i%prime[j]==0){ mu[i*prime[j]]=0; break; }else mu[i*prime[j]]=-mu[i]; } } for(int i=1;i<=cnt;i++){ for(int j=1;j*prime[i]<=n;j++){ g[prime[i]*j]+=mu[j]; } } for(int i=1;i<=n;i++){ sumg[i]=sumg[i-1]+g[i]; } } int cas; ll calc(int n,int m){ int nn=min(n,m); ll ans=0; for(int l=1,r;l<=nn;l=r+1){ r=min(n/(n/l),m/(m/l)); ans+=(sumg[r]-sumg[l-1])*(n/l)*(m/l); } return ans; } int main(){ sieve(maxn); scanf("%d",&cas); while(cas--){ scanf("%d %d",&n,&m); printf("%lld\n",calc(n,m)); } }
版权声明:因为我是蒟蒻,所以请大佬和神犇们不要转载(有坑)的文章,并指出问题,谢谢
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 你所不知道的 C/C++ 宏知识
· 聊一聊 操作系统蓝屏 c0000102 的故障分析
· SQL Server 内存占用高分析
· .NET Core GC计划阶段(plan_phase)底层原理浅谈
· .NET开发智能桌面机器人:用.NET IoT库编写驱动控制两个屏幕
· 我干了两个月的大项目,开源了!
· 推荐一款非常好用的在线 SSH 管理工具
· 千万级的大表,如何做性能调优?
· 聊一聊 操作系统蓝屏 c0000102 的故障分析
· .NET周刊【1月第1期 2025-01-05】