Luogu 1390 - 公約數的和 （尤拉函式/莫比烏斯反演、分塊）
UVa 11426 - GCD - Extreme (II)
題意
給定數字\(n\)，計算下式結果
\[\sum_{i=1}^n\sum_{j=i+1}^n\gcd(i,j) \]
限制
\(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(n/i,j/i)=1\)
即\(g(n,i)\)表示所有小於\(n/i\)的並與\(n/i\)互質的數的個數
根據尤拉函式可得，所有滿足\(gcd(n/i,j/i)=1\)的\(j/i\)個數即為\(\varphi(n/i)\)
所以\(g(n,i)=\varphi(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\)分開來考慮，以“數量*乘積”的方式表示答案
則原式則會變成
\[\sum_{i=1}^n\sum_{j=i+1}^n\gcd(i,j)=\sum_{k=1}^nk\sum_{i=1}^n\sum_{j=i+1}^n[\gcd(i,j)=k] \]
這裡我們需要列舉\(k\)來求和所有\(\gcd(i,j)=k\)的情況
因為\(gcd(a,b)=c\iff gcd(\frac{a}{c},\frac{b}{c})=1\)
我們還能繼續化簡上式，得到
\[\sum_{k=1}^nk\sum_{i=1}^n\sum_{j=i+1}^n[\gcd(i,j)=k]=\sum_{k=1}^nk\sum_{i=1}^{\frac{n}{k}}\sum_{j=i+1}^{\frac{n}{k}}[\gcd(i,j)=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[i*prim[j]]=-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[i*prim[j]]=-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(nlogn)\)
有沒有辦法再繼續降低複雜度呢？當然，我們可以引入分塊的思想
可以發現，雖然列舉的\(d\)在數值小的時候\(\frac n d\)變化很大，以至於我們每次都需要重新計算一次
但是如果\(d\)逐漸變大，由於整除的關係，會有一段連續的\(d\)使得\(\frac n d\)相同
那我們就可以根據這個性質，使得這一整段只計算一次就得出結果
假設當前\(d\)列舉到某個分塊的左邊界，那麼\(d_r=\lfloor \frac n {\lfloor \frac n d \rfloor} \rfloor\)即為這一分塊的右邊界，即滿足\(\frac n d=\frac n {d_r}\)的最大的\(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[i*prim[j]]=-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[i*prim[j]]=-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;
}
思路四（莫比烏斯反演、分塊）（再進階）
啊好……還沒學呢
最後的公式還沒看懂（公式還能再化簡）