51nod 1237 最大公约数之和 V3

题目大意：

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

题解：

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

反演1：

∑ni=1∑nj=1gcd(i,j)
=2∗∑ni=1∑ij=1gcd(i,j)−n∗(n+1)/2

∑ni=1∑ij=1gcd(i,j)
=∑nd=1d∗(∑ni=1∑ij=1(gcd(i,j)=d))
=∑nd=1d∗(∑⌊n/d⌋i=1∑ij=1(gcd(i,j)=1))
=∑nd=1d∗φ(⌊n/d⌋)
=∑i=1φ(i)∗(∑⌊n/i⌋j=1j)
=∑i=1φ(i)∗⌊n/i⌋∗(⌊n/i⌋+1)/2

反演2：

必要结论：
∑d|nd∗μ(n/d) = φ(n)
证明：
狄利克雷卷积形式,老套路证明。
只需要证明:∑d|pqd∗μ(pq/d)=φ(pq)
∑d|pqd∗μ(pq/d)
=∑qk=0pk∗μ(pq−k)
因为要使μ(pq−k)不为0，所以缩小循环范围。
=∑qk=q−1pk∗μ(pq−k)
=pq−pq−1
=φ(pq)
那么根据积性函数的性质可得：
所以∑d|nd∗μ(n/d) = φ(n)

∑ni=1∑nj=1gcd(i,j)
=∑nd=1d∗(∑ni=1∑nj=1(gcd(i,j)=d))
中间莫比乌斯反演过程略。
=∑nd=1d∗∑⌊n/d⌋i=1⌊n/(d∗i)⌋2∗μ(i)
设T=d∗i，则
=∑nT=1⌊n/T⌋2∗∑d|Td∗μ(n/d)
=∑nT=1⌊n/T⌋2∗φ(T)

坑爹的发现：

综合反演1和反演2会得到以下式子：
∑ni=1⌊n/i⌋∗φ(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);}