题目链接http://www.51nod.com/Challenge/Problem.html#!#problemId=1237
题意:求$\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)$ ,$1\leq{n}\leq10^{10}$.
知识提要:$n=\sum_{d|n}\varphi(d)$ 该公式证明https://blog.csdn.net/chen1352/article/details/50695930
根据莫比乌斯反演可以得到 $\varphi(n)=\sum_{d|n}\mu(\frac{n}{d})d$
解析:试着化简这个柿子,枚举最大公因数d
$$ \quad\sum_{d=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[(i,j)=1]d\\$$
$$=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[(i,j)=1]$$
右边根据莫比乌斯反演化简
$$\sum_{d=1}^{n}d\sum_{d'=1}^{\lfloor\frac{n}{d}\rfloor}\mu(d')\lfloor{\frac{n}{dd'}}\rfloor\lfloor{\frac{n}{dd'}}\rfloor$$
令T=dd'
$$\sum_{T=1}^{n}\lfloor{\frac{n}{T}}\rfloor\lfloor{\frac{n}{T}}\rfloor\sum_{d|T}\mu(\frac{T}{d})d$$
把右边替换成$\varphi(T)$得到
$$\sum_{T=1}^{n}\lfloor{\frac{n}{T}}\rfloor\lfloor{\frac{n}{T}}\rfloor\varphi(T)$$
所以求出来欧拉函数的前缀和 这个式子也就求出来了
AC代码
#include <bits/stdc++.h> #define pb push_back #define mp make_pair #define fi first #define se second #define all(a) (a).begin(), (a).end() #define fillchar(a, x) memset(a, x, sizeof(a)) #define huan printf("\n"); #define debug(a,b) cout<<a<<" "<<b<<" "; using namespace std; const int maxn=1e6+10,inf=0x3f3f3f3f; typedef long long ll; const ll mod = 1000000007; typedef pair<int,int> pii; int check[maxn],prime[maxn],phi[maxn],sum[maxn]; void Phi(int N)//莫比乌斯函数线性筛 { int pos=0;sum[0]=0; sum[1]=phi[1]=1; for(int i = 2 ; i <= N ; i++) { if (!check[i]) prime[pos++] = i,phi[i]=i-1; for (int j = 0 ; j < pos && i*prime[j] <= N ; j++) { check[i*prime[j]] = 1; if (i % prime[j] == 0) { phi[i*prime[j]]=phi[i]*prime[j]; break; } else phi[i*prime[j]]=phi[i]*(prime[j]-1); } sum[i]=(sum[i-1]+phi[i])%mod; } } unordered_map<ll,ll> ma; ll inv2=500000004; ll solve(ll n) { if(n<=1e6) return sum[n]; else if(ma.count(n)) return ma[n]; ll temp = ((n%mod)*((n+1)%mod)%mod)*inv2%mod; for(ll i=2,j;i<=n;i=j+1) { j=n/(n/i); temp = (temp-solve(n/i)*(j-i+1)%mod+mod)%mod; } return ma[n]=temp; } int main() { ll n; Phi(1e6); scanf("%lld",&n); ll ans=0; for(ll i=1,j;i<=n;i=j+1) { j=n/(n/i); ll k=(n/i)%mod; k = k*k%mod; ll temp=((solve(j)-solve(i-1)+mod)%mod)*k%mod; ans=(ans+temp)%mod; } printf("%lld\n",ans); }