51nod 1238 最小公倍数之和 V3

题目描述：

求：
$\sum_{i = 1}^n \sum_{j = 1}^n lcm(i, j)$

题解：

$\sum_{i = 1}^n \sum_{j = 1}^n lcm(i, j)$
$=(\sum_{i = 1}^n \sum_{j = 1}^i lcm(i,j)) * 2 - n * (n + 1) / 2$
$=(\sum_{i = 1}^n \sum_{j = 1}^i i * j / gcd(i, j)) * 2 - n * (n + 1) / 2$
$=(\sum_{d = 1}^n d * \sum_{i = 1}^{n /d} i * \sum_{j = 1}^i j( gcd(i, j) = 1))* 2 - n * (n + 1) / 2$
$=\sum_{d = 1}^n d * \sum_{i = 1}^{n /d} i * i * φ(i)$
$= \sum_{i = 1}^n φ(i) * i ^ 2 * \lfloor n / i \rfloor * (\lfloor n / i \rfloor + 1) / 2$

我们需要杜教筛$\sum_{i = 1}^n φ(i) * i ^ 2$

设$s(n) = \sum_{i = 1}^n φ(i) * i ^ 2$
$\sum_{i = 1}^n i * i \sum_{j | i} φ(j)$
$= \sum_{i = 1}i ^ 3$
$= \sum_{i = 1}^n i ^ 2 \sum_{j = 1}^{\lfloor n / i \rfloor} φ(j) * j^2$
$= \sum_{i = 1}^n i ^ 2 *s(\lfloor n / i \rfloor)$
$s(n)$
$= \sum_{i = 1}^n i ^ 3 - \sum_{i = 2}^n i ^ 2*s(\lfloor n / i \rfloor)$
$= n^2 *(n + 1)^2 / 4 - \sum_{i = 2}^n i ^ 2*s(\lfloor n / i \rfloor)$

像往常一样筛即可。

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[n];
    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);
}