Luogu 1390 - 公约数的和 （欧拉函数/莫比乌斯反演、分块）
UVa 11426 - GCD - Extreme (II)
题意
给定数字\(n\)，计算下式结果
限制
\(Case=1,\ n\leq2\times10^6,\ 1000ms\) （洛谷）
\(Case\leq100,\ n\leq4\times10^6,\ 10000ms\) （UVa）
思路一（欧拉函数）
刘汝佳蓝书 P124-P125
设\(f(n)=\gcd(1,n)+\gcd(2,n)+...+\gcd(n-1,n)=\sum_{i=1}^{n-1}\gcd(i,n)\)
则所求答案即为\(\sum_{i=1}^n\sum_{j=i+1}^n\gcd(i,j)=f(2)+f(3)+...+f(n)=\sum_{i=2}^nf(i)\)
洛谷仅单测试点，故最后直接\(O(n)\)累加所有\(f(i)\)输出即可
UVa有多测试点，需要预处理前缀和后再输出
设\(g(n,i)\)为满足\(\gcd(n,j)=i,\ j<n\)的所有满足条件的\(j\)的个数
即\(g(n,i)=\sum_{j=1}^{n-1}[\gcd(n,j)=i]\)
又因为\(\gcd(n,j)=i\iff \gcd(\frac n i,\frac j i)=1\)
即\(g(n,i)\)表示所有小于\(n/i\)的并与\(n/i\)互质的数的个数
根据欧拉函数可得，所有满足\(\gcd(\frac n i,\frac j i)=1\)的\(\frac j i\)个数即为\(\varphi(\frac n i)\)
所以\(g(n,i)=\varphi(\frac n i)\)，可以根据欧拉筛法预处理求出欧拉函数先
然后再考虑\(f(n)=\sum i\times g(n,i),\ i|n\)
要想得到\(f(n)\)，首先就得找到所有小于\(n\)的因子（不包括\(n\)）
所以可以根据埃氏筛法，从小到大枚举因子，将答案加给因子的倍数即可
代码一 - Luogu 1390
Case Max (291ms/1000ms)
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=2000000;
ll eular[N+50],prim[N+50];
ll f[N+50];
bool vis[N+50];
void init(int n)
{
memset(vis,false,sizeof vis);
int p=0;
for(int i=2;i<=n;i++)
{
if(!vis[i])
{
prim[p++]=i;
eular[i]=i-1;
}
for(int j=0;j<p;j++)
{
int k=i*prim[j];
if(k>n)
break;
vis[k]=true;
if(i%prim[j]==0)
{
eular[k]=eular[i]*prim[j];
break;
}
eular[k]=eular[i]*eular[prim[j]];
}
}
}
int main()
{
int n;
scanf("%d",&n);
init(n);
for(int i=1;i<=n;i++)
{
for(int j=i*2;j<=n;j+=i)
f[j]+=eular[j/i]*i;
}
ll ans=0;
for(int i=2;i<=n;i++)
ans+=f[i];
printf("%lld\n",ans);
return 0;
}
代码一 - UVa 11426
(630ms/10000ms)
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=4000000;
ll eular[N+50],prim[N+50];
ll f[N+50],ans[N+50];
bool vis[N+50];
void init(int n)
{
memset(vis,false,sizeof vis);
int p=0;
for(int i=2;i<=n;i++)
{
if(!vis[i])
{
prim[p++]=i;
eular[i]=i-1;
}
for(int j=0;j<p;j++)
{
int k=i*prim[j];
if(k>n)
break;
vis[k]=true;
if(i%prim[j]==0)
{
eular[k]=eular[i]*prim[j];
break;
}
eular[k]=eular[i]*eular[prim[j]];
}
}
}
int main()
{
init(N);
for(int i=1;i<=N;i++)
{
for(int j=i*2;j<=N;j+=i)
f[j]+=eular[j/i]*i;
}
for(int i=2;i<=N;i++)
ans[i]=ans[i-1]+f[i];
int n;
while(scanf("%d",&n)!=EOF&&n)
{
printf("%lld\n",ans[n]);
}
return 0;
}
思路二（莫比乌斯反演）（简单）
类似于思路一，我们将所有不同的\(\gcd(i,j)=k\)分开来考虑，以“数量*乘积”的方式表示答案
则原式则会变成
这里我们需要枚举\(k\)来求和所有\(\gcd(i,j)=k\)的情况
因为\(gcd(a,b)=c\iff gcd(\frac{a}{c},\frac{b}{c})=1\)
我们还能继续化简上式，得到
明显发现，该式子后半部分即为莫比乌斯反演的模板题——求两区间内互质数的对数
直接套板子即可
代码二 - Luogu 1390
Case Max(139ms/1000ms)
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=2000000;
ll mu[N+50],prim[N+50];
bool vis[N+50];
void init(int n)
{
memset(vis,false,sizeof vis);
mu[1]=1;
int p=0;
for(int i=2;i<=n;i++)
{
if(!vis[i])
{
prim[p++]=i;
mu[i]=-1;
}
for(int j=0;j<p;j++)
{
int k=i*prim[j];
if(k>n)
break;
vis[k]=true;
if(i%prim[j]==0)
{
mu[k]=0;
break;
}
else
mu[k]=-mu[i];
}
}
}
int main()
{
int n;
scanf("%d",&n);
init(n);
ll ans=0;
for(int d=1;d<=n;d++)
{
ll ansd=0;
int m=n/d;
for(int i=1;i<=m;i++)
ansd+=mu[i]*(m/i)*(m/i);
ans+=(ansd>>1)*d;
}
printf("%lld\n",ans);
return 0;
}
代码二 - UVa 11426
(2310ms/10000ms)
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=4000000;
ll mu[N+50],prim[N+50];
bool vis[N+50];
void init(int n)
{
memset(vis,false,sizeof vis);
mu[1]=1;
int p=0;
for(int i=2;i<=n;i++)
{
if(!vis[i])
{
prim[p++]=i;
mu[i]=-1;
}
for(int j=0;j<p;j++)
{
int k=i*prim[j];
if(k>n)
break;
vis[k]=true;
if(i%prim[j]==0)
{
mu[k]=0;
break;
}
else
mu[k]=-mu[i];
}
}
}
int main()
{
init(N);
int n;
while(scanf("%d",&n)&&n)
{
ll ans=0;
for(int d=1;d<=n;d++)
{
ll ansd=0;
int m=n/d;
for(int i=1;i<=m;i++)
ansd+=mu[i]*(m/i)*(m/i);
ans+=(ansd>>1)*d;
}
printf("%lld\n",ans);
}
return 0;
}
思路三（莫比乌斯反演、分块）（进阶）
在思路二中，我们将式子化简到了直接套板子的情况
但是会发现，我们每次枚举\(d\)时，对于\(\gcd=d\)的种类数就必须计算\(\frac n d\)次才能出结果
总时间复杂度严格为\(O(n\sqrt n)\)
有没有办法再继续降低复杂度呢？当然，我们可以引入分块的思想
可以发现，虽然枚举的\(d\)在数值小的时候\(\frac n d\)变化很大，以至于我们每次都需要重新计算一次
但是如果\(d\)逐渐变大，由于整除的关系，会有一段连续的\(d\)使得\(\frac n d\)相同
那我们就可以根据这个性质，使得这一整段只计算一次就得出结果
假设当前\(d\)枚举到某个分块的左边界，那么\(d_r=\lfloor \frac n {\lfloor \frac n d \rfloor} \rfloor\)即为这一分块的右边界，即满足\(\lfloor\frac n d\rfloor=\lfloor\frac n {d_r}\rfloor\)的最大的\(d_r\)
代码三 - Luogu 1390
Case Max(79ms/1000ms)
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=2000000;
ll mu[N+50],prim[N+50];
bool vis[N+50];
void init(int n)
{
memset(vis,false,sizeof vis);
mu[1]=1;
int p=0;
for(int i=2;i<=n;i++)
{
if(!vis[i])
{
prim[p++]=i;
mu[i]=-1;
}
for(int j=0;j<p;j++)
{
int k=i*prim[j];
if(k>n)
break;
vis[k]=true;
if(i%prim[j]==0)
{
mu[k]=0;
break;
}
else
mu[k]=-mu[i];
}
}
}
int main()
{
int n;
scanf("%d",&n);
init(N);
ll ans=0;
for(int d=1;d<=n;)
{
int m=n/d;
ll ansd=0;
for(int i=1;i<=m;i++)
ansd+=mu[i]*(m/i)*(m/i);
ansd>>=1;
int nxt=n/(n/d); //计算这一块的右边界
ans+=ansd*d;
for(int i=d+1;i<=nxt;i++)
ans+=ansd*i; //需要遍历乘上权值
d=nxt+1; //d跳到下一分块的左边界
}
printf("%lld\n",ans);
return 0;
}
代码三 - UVa 11426
(1650ms/10000ms)
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=4000000;
ll mu[N+50],prim[N+50];
bool vis[N+50];
void init(int n)
{
memset(vis,false,sizeof vis);
mu[1]=1;
int p=0;
for(int i=2;i<=n;i++)
{
if(!vis[i])
{
prim[p++]=i;
mu[i]=-1;
}
for(int j=0;j<p;j++)
{
int k=i*prim[j];
if(k>n)
break;
vis[k]=true;
if(i%prim[j]==0)
{
mu[k]=0;
break;
}
else
mu[k]=-mu[i];
}
}
}
int main()
{
init(N);
int n;
while(scanf("%d",&n)&&n)
{
ll ans=0;
for(int d=1;d<=n;)
{
int m=n/d;
ll ansd=0;
for(int i=1;i<=m;i++)
ansd+=mu[i]*(m/i)*(m/i);
ansd>>=1;
int nxt=n/(n/d); //计算这一块的右边界
ans+=ansd*d;
for(int i=d+1;i<=nxt;i++)
ans+=ansd*i; //需要遍历乘上权值
d=nxt+1; //d跳到下一分块的左边界
}
printf("%lld\n",ans);
}
return 0;
}
