【51nod 1237】最大公约数之和 V3-CSDN博客

本文介绍了一种利用杜教筛算法解决数学问题的方法，具体为计算小于等于给定数值N的所有整数对的最大公约数之和。通过巧妙的数学转换，将原始问题转化为求解欧拉函数的过程，并给出了高效的实现方案。

description

给出一个数N（<= 10^10），输出小于等于N的所有数，两两之间的最大公约数之和。

相当于计算这段程序（程序中的gcd(i,j)表示i与j的最大公约数）：
由于结果很大，输出Mod 1000000007的结果。

G=0;
for(i=1;i< N;i++)
for(j=1;j<=N;j++)
{
G = (G + gcd(i,j)) % 1000000007;
}

Solution

这道题想出来后被人告知是杜教筛……
前面的都很顺利，非常快的推出式子。

a n s = \sum i = 1 n \sum j = 1 n g c d (i, j) = \sum d = 1 n d \sum i = 1 n / d ⌊ n i d ⌋ * ⌊ n i d ⌋ * μ (i) = \sum T = 1 n ⌊ n T ⌋ * ⌊ n T ⌋ * \sum d | T d * μ (T / d)

$ans=\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)=\sum_{d=1}^{n}d\sum_{i=1}^{n/d}{{\lfloor{n\over id }\rfloor}*{\lfloor{n\over id }\rfloor}*\mu(i)}=\sum_{T=1}^{n}{\lfloor{n\over T }\rfloor}*{\lfloor{n\over T }\rfloor}*\sum_{d|T}d*\mu(T/d)$
推到这里，我们发现一个很喜欢的式子

∑d|Td∗μ(T/d) $\sum_{d|T}d*\mu(T/d)$ ，因为我们知道

ϕ(n)=∑d|nd∗μ(n/d) $\phi(n)=\sum_{d|n}d*\mu(n/d)$ 。

a n s = \sum T = 1 n ⌊ n T ⌋ * ⌊ n T ⌋ * ϕ (T) 。

$ans=\sum_{T=1}^{n}{\lfloor{n\over T }\rfloor}*{\lfloor{n\over T }\rfloor}*\phi(T)。$
然后考虑分块，所以我们要求一段数的欧拉函数，这和之前做的一道 51nod1239 欧拉函数之和简直一模一样，这里就不再推导了，直接上结论。

设f(x)=∑xi=1ϕ(i) $设f(x)=\sum_{i=1}^{x}\phi(i)$
那么

f (x) = (1 + x) * x / 2 - \sum T = 2 x f (⌊ x / T ⌋) 。

$f(x)=(1+x)*x/2-\sum_{T=2}^{x}f(\lfloor x/T \rfloor)。$ 所以ans中欧拉函数再分块中用

f(j)−f(i−1) $f(j)-f(i-1)$ 算一下即可。复杂度O(

n2/3 $n^{2/3}$ )。

Code

#include<iostream>
#include<cmath>
#include<cstring>
#include<cstdio>
#include<algorithm>
#define ll long long
using namespace std;
const ll maxn=1e6+5,mo=1e9+7,mo2=5e8+4;
ll d[maxn],bz[maxn],p[maxn],h[maxn*10],f[maxn*10];
ll n,i,t,j,k,l,x,y,z,ans;
int hash(ll x){
    ll t=x%maxn;
    while (h[t] && h[t]!=x) t=(t+1)%maxn;
    return t;
}
ll dg(ll n){
    if (n<=maxn) return p[n];
    ll t=n%mo,k=t*(t+1)%mo*mo2%mo,i=2,x=hash(n),y;
    if (h[x]) return f[x];
    while (i<=n){
        t=n/(n/i);
        k-=dg(n/i)*(t-i+1)%mo;i=t+1;
    }
    k=(k%mo+mo)%mo;h[x]=n;f[x]=k;return k;
}
int main(){
    //freopen("data.in","r",stdin);
    scanf("%lld",&n);p[1]=1;
    for (i=2;i<=maxn;i++){
        if (!bz[i]) d[++d[0]]=i,p[i]=i-1;
        for (j=1;j<=d[0];j++){
            if (i*d[j]>maxn) break;
            bz[i*d[j]]=1;
            if (i%d[j]==0){
                p[i*d[j]]=p[i]*d[j];break;
            }else p[i*d[j]]=p[i]*p[d[j]];
        }
    }
    for (i=1;i<=maxn;i++)
        p[i]+=p[i-1];
    i=1;
    while (i<=n){
        t=n/(n/i);
        k=(dg(t)-dg(i-1)+mo)%mo;l=n/i%mo;
        ans=(ans+l*l%mo*k%mo)%mo;
        i=t+1;
    }
    printf("%lld\n",ans);
}