51nod 1238 最小公倍数之和 V3

题目描述：

求：

∑ni=1∑nj=1lcm(i,j)

题解：

∑ni=1∑nj=1lcm(i,j)

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

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

=(∑nd=1d∗∑n/di=1i∗∑ij=1j(gcd(i,j)=1))∗2−n∗(n+1)/2

=∑nd=1d∗∑n/di=1i∗i∗φ(i)

=∑ni=1φ(i)∗i2∗⌊n/i⌋∗(⌊n/i⌋+1)/2

我们需要杜教筛∑ni=1φ(i)∗i2

设s(n)=∑ni=1φ(i)∗i2

∑ni=1i∗i∑j|iφ(j)

=∑i=1i3

=∑ni=1i2∑⌊n/i⌋j=1φ(j)∗j2

=∑ni=1i2∗s(⌊n/i⌋)

s(n)

=∑ni=1i3−∑ni=2i2∗s(⌊n/i⌋)

=n2∗(n+1)2/4−∑ni=2i2∗s(⌊n/i⌋)

像往常一样筛即可。

Code：

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

const ll mo = 1e9 + 7, ni_2 = 5e8 + 4, ni_4 = 250000002, ni_6 = 166666668;
const int Maxn = 5000000, M = 12131891;

bool bz[Maxn + 1]; int p[Maxn / 10];
ll phi[Maxn + 1];
ll h[M][2], n;

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

ll sa(ll x) {
x %= mo; return x * (x + 1) % mo * ni_2 % mo;
}

ll ss(ll x) {
x %= mo; return x * (x + 1) % mo * (2 * x + 1) % mo * ni_6 % mo;
}

ll dg(ll n) {
if(n <= Maxn) return phi
;
int y = hash(n);
if(h[y][0] != 0) return h[y][1];
h[y][0] = n;
ll k = sa(n) * sa(n) % mo;
for(ll i = 2; i <= n; i ++) {
ll j = n / (n / i);
k -= dg(n / i) * (ss(j) - ss(i - 1) + mo) % mo;
i = j;
}
h[y][1] = (k % mo + mo) % mo;
return h[y][1];
}

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]) {
int 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] * phi[p[j]];
}
}
fo(i, 1, Maxn)
phi[i] = (phi[i - 1] + (ll) i * i % mo * phi[i] % mo) % mo;
scanf("%lld", &n);
ll ans = 0;
for(ll i = 1; i <= n; i ++) {
ll j = n / (n / i);
ll x = (n / i) % mo;
ans += x * (x + 1) % mo * (dg(j) - dg(i - 1)) % mo;
i = j;
}
ans = (ans % mo + mo) * ni_2 % mo;
printf("%lld", ans);
}