【题解】最大公约数之和 V3 51nod 1237 杜教筛

题目传送门 http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1237

数学题真是做的又爽又痛苦，爽在于只要推出来公式基本上就是AC，痛苦就在于推公式。。。

题意很简单，求

$\Large\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}gcd(i,j)$

其中$n\le 10^{10}$

这个题有很多做法，除了普及组的$O(n^2\log n)$做法，还有用莫比乌斯反演+分块优化的$O(n)$做法

然而因为这个题两个维度的限制是相等的，都是$n$，所以可以用杜教筛做到$O(n^{\frac{2}{3}})$

然后推一波公式

$\Large\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}gcd(i,j)$

按照套路，可以枚举最大公因数$g$，于是有

$\Large\sum\limits_{g=1}^{n}\sum\limits_{g=gcd(i,j)}g$

继续变形

$\Large\sum\limits_{g=1}^{n}g\cdot(\sum\limits_{g=gcd(i,j)}1)$

$\Large\sum\limits_{g=1}^{n}g\cdot(\sum\limits_{1=gcd(i,j)}^{i\le\lfloor\frac{n}{g}\rfloor,j\le\lfloor\frac{n}{g}\rfloor}1)$

到这里，我们就可以用莫比乌斯反演的套路做到$O(n)$的复杂度了，但是我们换种形式继续推式子

我们先关注这一部分的变形，暂时不管前面的那一部分

$\Large\sum\limits_{1=gcd(i,j)}^{i\le\lfloor\frac{n}{g}\rfloor,j\le\lfloor\frac{n}{g}\rfloor}1$

展开，得到

$\Large\sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{g}\rfloor}[gcd(i,j)=1]\cdot 1$

其中$[gcd(i,j)=1]$的意思是，当$gcd(i,j)=1$的时候这个东西的值为$1$，否则为$0$

然后用$\varphi()$函数进行化简，于是式子就变成了

$\Large 2\cdot(\sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor}\varphi(i))-1$

注意到中间的一部分是$\varphi()$函数的前缀和，于是我们可以用杜教筛算

为了方便，设函数

$\Large S(n)=\sum\limits_{i=1}^{n}\varphi(i)$

然后再看一眼总的式子

$\Large\sum\limits_{g=1}^{n}g\cdot(2\cdot S(\lfloor\frac{n}{g}\rfloor)-1)$

于是$\lfloor\frac{n}{g}\rfloor$是整除，可以用分块优化，这样的话，除去杜教筛求$S()$的部分，复杂度是$O(n^{\frac{1}{2}})$

因为杜教筛有记忆化，复杂度和分块优化的部分是分离的，所以总复杂度是

$\Large O(n^{\frac{1}{2}}+n^{\frac{2}{3}})=O(n^{\frac{2}{3}})$

单点时限$5s$，我杜教筛记忆化用的$map$，没有自己手写$hash$，于是复杂度就多了$O(\log n)$，但是仍然每个点跑到了$1s$以内，还是挺不错的

上代码：

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <map>

using namespace std;
typedef long long ll;
const int MOD = 1e9+7;
const int LIMIT = 4641589;

bool vis[LIMIT] = {0};
ll prime[LIMIT], pidx = 0, phi[LIMIT] = {0}, pfphi[LIMIT] = {0};
void prelude_phi() { // 线性筛预处理一部分前缀和
    phi[1] = 1;
    for( ll i = 2; i < LIMIT; ++i ) {
        if( !vis[i] ) {
            prime[pidx++] = i;
            phi[i] = i-1;
        }
        for( int j = 0; j < pidx; ++j ) {
            ll k = i * prime[j];
            if( k >= LIMIT ) break;
            vis[k] = true;
            if( i % prime[j] ) phi[k] = phi[i] * phi[prime[j]] % MOD;
            else {
                phi[k] = phi[i] * prime[j] % MOD;
                break;
            }
        }
    }
    for( int i = 1; i < LIMIT; ++i )
        pfphi[i] = (pfphi[i-1] + phi[i]) % MOD;
}

ll inv2;
ll pow_mod( ll a, ll b ) {
    if( !b ) return 1;
    ll rtn = pow_mod(a,b>>1);
    rtn = rtn * rtn % MOD;
    if( b&1 ) rtn = rtn * a % MOD;
    return rtn;
}
ll inv( ll x ) {
    return pow_mod(x,MOD-2);
}

map<ll,ll> mp; // 记忆化用
ll S( ll n ) { // 杜教筛求前缀和
    if( n < LIMIT ) return pfphi[n];
    if( mp.count(n) ) return mp[n];
    ll ans = (n % MOD) * ((n+1) % MOD) % MOD * inv2 % MOD;
    for( ll i = 2, j; i <= n; i = j+1 ) {
        j = n/(n/i);
        ans = (ans - S(n/i) * ((j-i+1) % MOD) % MOD + MOD) % MOD;
    }
    mp[n] = ans;
    return ans;
}

ll n;

void solve() {
    ll ans = 0;
    for( ll i = 1, j; i <= n; i = j+1 ) {
        j = n/(n/i); // 分块优化
        ans = (ans + (((i+j) % MOD) * ((j-i+1) % MOD) % MOD * inv2 % MOD) * ((2 * S(n/i) - 1) % MOD) % MOD ) % MOD;
    }
    printf( "%lld\n", ans );
}

int main() {
    prelude_phi();
    inv2 = inv(2);
    scanf( "%lld", &n );
    solve();
    return 0;
}