51nod 1237 最大公约数之和 V3

题目大意：

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

题解：

杜教筛。
这里提供两种反演姿势，之后我们会发现一个神奇的东西。
在此之前你得会杜教筛φ。

反演1：

$\sum_{i = 1}^n \sum_{j = 1}^n gcd(i, j)$
$= 2 *\sum_{i = 1}^n \sum_{j = 1}^i gcd(i, j) - n * (n + 1) / 2$

$\sum_{i = 1}^n \sum_{j = 1}^i gcd(i, j)$
$= \sum_{d = 1}^n d * (\sum_{i = 1}^n \sum_{j = 1}^ i (gcd(i,j) = d))$
$= \sum_{d = 1}^n d * (\sum_{i = 1}^{\lfloor n / d \rfloor} \sum_{j = 1}^i (gcd(i,j) = 1))$
$= \sum_{d = 1}^n d * φ(\lfloor n / d \rfloor)$
$= \sum_{i = 1} φ(i) * (\sum_{j = 1}^{\lfloor n / i \rfloor} j)$
$= \sum_{i = 1} φ(i) * \lfloor n / i \rfloor * (\lfloor n / i \rfloor + 1) / 2$

反演2：

必要结论：
$\sum_{d | n} d * μ(n / d)$ = $φ(n)$
证明：
狄利克雷卷积形式,老套路证明。
只需要证明:$\sum_{d | p^q} d * μ(p^q / d) = φ(p ^ q)$
$\sum_{d | p^q} d * μ(p^q / d)$
$= \sum_{k = 0}^q p^k * μ(p^{q - k})$
因为要使$μ(p^{q - k})$不为0，所以缩小循环范围。
$= \sum_{k = q - 1}^q p ^ k *μ(p^{q - k})$
$= p^q - p^{q - 1}$
$= φ(p ^ q)$
那么根据积性函数的性质可得：
所以$\sum_{d | n} d * μ(n / d)$ = $φ(n)$

$\sum_{i = 1}^n \sum_{j = 1}^n gcd(i, j)$
$= \sum_{d = 1}^n d * (\sum_{i = 1}^n \sum_{j = 1}^n (gcd(i, j)=d))$
中间莫比乌斯反演过程略。
$= \sum_{d = 1}^n d * \sum_{i = 1}^{\lfloor n / d \rfloor} \lfloor n / (d * i) \rfloor^2 * μ(i)$
设$T = d * i$，则
$= \sum_{T = 1}^n \lfloor n / T \rfloor ^2 * \sum_{d | T} d * μ(n / d)$
$= \sum_{T = 1}^n \lfloor n / T \rfloor ^2 * φ(T)$

坑爹的发现：

综合反演1和反演2会得到以下式子：
$\sum_{i = 1}^n \lfloor n / i \rfloor * φ(i) = n * (n +1)/ 2$
显然我不知道它有什么用，当做一个比较吊的结论吧。

反演2Code：

#include<cmath>
#include<cstdio>
#define ll long long
#define fo(i, x, y) for(ll i = x; i <= y; i ++)
using namespace std;

const ll Maxn = 5000000, M = 12461987, ni_2 = 5e8 + 4, mo = 1e9 + 7;
ll p[700000], phi[Maxn + 5], h[M][2];
bool bz[Maxn + 5];

ll hash(ll x) {
    ll y = x % M;
    while(h[y][0] != 0 && h[y][0] != x)
        y = (y + 1) % M;
    return y;
}

ll dg(ll n) {
    if(n <= Maxn) return phi[n];
    ll y = hash(n);
    if(h[y][0] == n) return h[y][1];
    ll t = n % mo;
    h[y][1] = t * (t + 1) % mo * ni_2 % mo;
    fo(i, 2, n) {
        ll j = n / (n / i);
        h[y][1] -= dg(n / i) * (j - i + 1) % mo;
        i = j;
    }
    h[y][0] = n; h[y][1] = (h[y][1] % mo + mo) % mo;
    return h[y][1];
}

ll n;

int main() {
    phi[1] = 1;
    fo(i, 2, Maxn) {
        if(!bz[i]) p[++ p[0]] = i, phi[i] = i - 1;
        fo(j, 1, p[0]) {
            ll k = i * p[j];
            if(k > Maxn) break;
            bz[k] = 1;
            if(i % p[j] == 0) {
                phi[k] = phi[i] * p[j];
                break;
            }
            phi[k] = phi[i] * (p[j] - 1);
        }
    }
    fo(i, 1, Maxn) phi[i] = (phi[i - 1] + phi[i]) % mo;
    ll ans = 0;
    scanf("%lld", &n);
    fo(i, 1, n)  {
        ll j = n / (n / i), k = n / i % mo;
        k = (k * k) % mo;
        ans += k * (dg(j) - dg(i - 1) + mo) % mo;
        ans %= mo;
        i = j;
    }
    printf("%lld", ans);
}