题目:求
∑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
这样左边就可以配一个 f(x)=x2 杜教筛,右边进行分块。
#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;
}