题目
题目链接：https://gmoj.net/senior/#main/show/4496
求
[sum^{n}_{i=1}sum_{j|i}gcd(j,frac{i}{j})
]
(nleq 10^{11})。
思路
[sum^{n}_{i=1}sum_{j|i}gcd(j,frac{i}{j})
]
[=sum^{n}_{i=1}sum^{lfloorfrac{n}{i}
floor}_{j=1}gcd(i,j)
]
枚举 (gcd=d)
[sum^{n}_{d=1}dsum^{lfloorfrac{n}{d^2}
floor}_{i=1}sum^{lfloorfrac{n}{d^2i}
floor}_{j=1}[gcd(i,j)=1]
]
令 (g(n)=sum^{n}_{i=1}sum^{lfloorfrac{n}{i} floor}_{j=1}[gcd(i,j)=1])，则原式 (=sum^{n}_{d=1}g(lfloorfrac{n}{d^2} floor))。
[g(n)=sum^{sqrt{n}}_{d=1}mu(d)sum^{lfloorfrac{n}{d^2}
floor}_{i=1}lfloorfrac{n}{d^2i}
floor
]
整除分块即可。记忆化一下可以卡进 (1s)。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=350010,M=10000000;
int m,prm[N],mu[N];
bool v[N];
ll n,ans,f[M];
void findprm(int n)
{
mu[1]=1;
for (int i=2;i<=n;i++)
{
if (!v[i]) prm[++m]=i,mu[i]=-1;
for (int j=1;j<=m;j++)
{
if (i>n/prm[j]) break;
v[i*prm[j]]=1; mu[i*prm[j]]=-mu[i];
if (!(i%prm[j])) { mu[i*prm[j]]=0; break; }
}
}
}
ll calc2(ll n)
{
if (n<M && f[n]) return f[n];
ll res=0;
for (ll l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
res+=(r-l+1)*(n/l);
}
if (n<M) f[n]=res;
return res;
}
ll calc1(ll n)
{
ll res=0;
for (ll i=1;i*i<=n;i++)
res+=mu[i]*calc2(n/i/i);
return res;
}
int main()
{
freopen("gcd.in","r",stdin);
freopen("gcd.out","w",stdout);
findprm(N-1);
scanf("%lld",&n);
for (ll i=1;i*i<=n;i++)
ans+=i*calc1(n/i/i);
cout<<ans;
return 0;
}