51Nod 1237 最大公约数之和 V3 杜教筛

题目链接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);
}