51nod1238 最小公倍数之和V3

题目：求$\sum^n_{i=1}\sum^n_{j=1}lcm(i,j)(n<=10^{10})$。
先算一个东西。

$$C(n)=\sum^n_{j=1}i[i\bot n](\bot 表示互质)$$
证：
$$C(n)=\sum^n_{j=1}i[i\bot n]$$
反演一波
$$C(n)=\sum_{g|n}\mu(g) \sum^{n\over g}_{i=1} i\cdot g$$
$$C(n)=\sum_{g|n}\mu(g)\cdot g \cdot {{{\frac n g}\cdot (\frac n g + 1)} \over 2}$$
把分数外边的g和分子上的g消掉，把n提出去。
$$C(n)=n\sum_{g|n}{\mu(g)\cdot ({n \over g}+1) \over 2}$$
括号里的两项拆开
$$C(n)={n\cdot \sum_{g|n}\mu(g)\cdot {n \over g} \over 2 }+{n\cdot \sum_{g|n}\mu(g)\over2}$$
左边是经典的式子
$$C(n)={n\cdot\varphi(n)\over 2}+{n\cdot\sum_{g|n}\mu(g)\over2}$$
右边只有n=1时有值
$$C(n)={n\varphi(n)+[n=1]\over2}$$
就是说小于等于n的数中于n互质的和是A(n)。

令

$$A(i)=\sum^i_{j=1}lcm(i,j)$$
$$A(i)=\sum_{d|i}\sum^{i\over d}_{j=1}ij[{i\over d}\bot j]$$
$$A(i)=\sum_{d|i}i\sum^{i\over d}_{j=1}j[{i\over d} \bot j]$$
带入C(i/d)
$$A(i)=\sum_{d|i}i\cdot {{i\over d}\cdot {\varphi({i\over d})}+[{i\over d}=1]\over2}$$
把常数提出来，顺便交换d和i/d
$$A(i)={i\over2 }\sum_{d|i}(d\varphi(d)+[d=1])$$
$$ans=\sum^n_{i=1}\sum^n_{j=1}lcm(i,j)$$
$$ans=2\sum^n_{i=1}\sum^i_{j=1}lcm(i,j)-\sum^n_{i=1}i$$
把A(i)带入
$$ans=2\sum^n_{i=1}({i\over2 }\sum_{d|i}(d\varphi(d)+[d=1]))-\sum^n_{i=1}i$$
2的常数就没了，然后看$[d=1]$的项，每一个i都只有一个，系数就是i，所以总共是$\sum^n_{i=1}i$和后边抵消
$$ans=\sum^n_{i=1}i\sum_{d|i}d\varphi(d)$$
设$ab=i$，把$d|i$改写一下
$$ans=\sum^n_{i=1}i\sum_{ab=i}a\varphi(a)$$
$$ans=\sum^n_{i=1}\sum_{ab=i}a^2b\varphi(a)$$
枚举$a$，看有多少个b
$$ans=\sum^n_{a=1}a^2\varphi(a)\sum^{\lfloor {n\over a}\rfloor}_{b=1}b$$
这样左边就可以配一个$f(x)=x^2$杜教筛，右边进行分块。

#include <stdio.h>
#include <algorithm>

using namespace std;

//long long
const long long maxn = 5100000;
const long long mod = 1000000007;
const long long mod1 = 1000007;
const long long rev2 = 500000004;
const long long rev6 = 166666668;
struct link {
    long long x , y;
    link *next;
} pool[maxn] , *g[mod1 + 13];
long long top;
long long n , N;
long long p[maxn] , prime[maxn] , tot;
long long f[maxn];
long long ans;
void predo () {
    long long i , j;
    N = 4700000;
    f[1] = 1;
    for ( i = 2 ; i <= N ; i++ ) {
        if ( p[i] == 0 ) {
            prime[++tot] = i;
            f[i] = i - 1;
        }
        for ( j = 1 ; j <= tot ; j++ ) {
            if ( i * prime[j] > N ) break;
            p[i*prime[j]] = 1;
            f[i*prime[j]] = f[i] * (prime[j] - 1);
            if ( i % prime[j] == 0 ) {
                f[i*prime[j]] = f[i] * prime[j];
                break;
            }
        }
    }
    for ( i = 1 ; i <= N ; i++ ) {
        f[i] = (((f[i]*i)%mod)*i)%mod;
    }
    for ( i = 1 ; i <= N ; i++ ) f[i] = (f[i] + f[i-1]) % mod;
}
long long linsum ( long long x ) {
    long long ret;
    ret = ((x%mod)*((x+1)%mod))%mod;
    ret = (ret*rev2) % mod;
    return ret;
}
long long sqrsum ( long long x ) {
    long long ret;
    ret = ( (x%mod) * ((x+1)%mod) ) % mod;
    ret = (ret * ((2*x+1)%mod) ) % mod;
    ret = (ret * rev6) % mod;
    return ret;
}
long long get ( long long x ) {
    if ( x <= N ) return f[x];
    for ( link *j = g[x%mod1] ; j ; j = j -> next ) if ( j -> x == x ) return j -> y;
    long long l , r , ret = ((((x%mod)*((x+1)%mod))%mod)*rev2)%mod;
    ret = (ret*ret) % mod;
    for ( l = 2 ; l <= x ; l = r + 1 ) {
        r = x/(x/l);
        ret -= ( ((sqrsum(r)-sqrsum(l-1)+mod)%mod) * get ( x / l ) ) % mod;
        ret += mod;
        ret %= mod;
    }
    link *tmp = &pool[++top];
    tmp -> x = x; tmp -> y = ret; tmp -> next = g[x%mod1]; g[x%mod1] = tmp;
    return ret;
}
void work () {
    long long l , r;
    scanf ( "%lld" , &n );
    for ( l = 1 ; l <= n ; l = r + 1 ) {
        r = n/(n/l);
        ans += ( ((get(r)-get(l-1)+mod)%mod) * linsum (n/l) ) % mod;
        ans %= mod;
    }
    printf ( "%lld\n" , ans );
}
int main () {
    predo ();
    work ();
    return 0;
}