51nod1238 最小公倍数之和V3

题目：求∑ni=1∑nj=1lcm(i,j)(n<=1010)。
先算一个东西。

C(n)=∑j=1ni[i⊥n](⊥表示互质)

证：

C(n)=∑j=1ni[i⊥n]

反演一波

C(n)=∑g|nμ(g)∑i=1ngi⋅g

C(n)=∑g|nμ(g)⋅g⋅ng⋅(ng+1)2

把分数外边的g和分子上的g消掉，把n提出去。

C(n)=n∑g|nμ(g)⋅(ng+1)2

括号里的两项拆开

C(n)=n⋅∑g|nμ(g)⋅ng2+n⋅∑g|nμ(g)2

左边是经典的式子

C(n)=n⋅φ(n)2+n⋅∑g|nμ(g)2

右边只有n=1时有值

C(n)=nφ(n)+[n=1]2

就是说小于等于n的数中于n互质的和是A(n)。

令

A(i)=∑j=1ilcm(i,j)

A(i)=∑d|i∑j=1idij[id⊥j]

A(i)=∑d|ii∑j=1idj[id⊥j]

带入C(i/d)

A(i)=∑d|ii⋅id⋅φ(id)+[id=1]2

把常数提出来，顺便交换d和i/d

A(i)=i2∑d|i(dφ(d)+[d=1])

ans=∑i=1n∑j=1nlcm(i,j)

ans=2∑i=1n∑j=1ilcm(i,j)−∑i=1ni

把A(i)带入

ans=2∑i=1n(i2∑d|i(dφ(d)+[d=1]))−∑i=1ni

2的常数就没了，然后看[d=1]的项，每一个i都只有一个，系数就是i，所以总共是∑ni=1i和后边抵消

ans=∑i=1ni∑d|idφ(d)

设ab=i，把d|i改写一下

ans=∑i=1ni∑ab=iaφ(a)

ans=∑i=1n∑ab=ia2bφ(a)

枚举a，看有多少个b

ans=∑a=1na2φ(a)∑b=1⌊na⌋b

这样左边就可以配一个f(x)=x2杜教筛，右边进行分块。

#include <stdio.h>#include <algorithm>using namespace std;//long longconst 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;}