[杜教筛][莫比乌斯反演] LOJ #6229. 这是一道简单的数学题

$Solution$

求

$$\sum_{i=1}^n\sum_{j=1}^i{ij\over(i,j)^2}$$ 可能算是 裸题?
相当于求 $${1\over2}\sum_{d=1}^n\big(d[d=1]+d^2\varphi(d)\big)\lfloor{n\over d}\rfloor$$
关于 $id^n\cdot\varphi$可以直接杜教筛 $(id^n\cdot\varphi)*id^n=id^{n+1}$。

#include <bits/stdc++.h>
#define show(x) cerr << #x << " = " << x << endl
using namespace std;
typedef long long ll;
typedef pair<int, int> Pairs;

const int N = 1010101;
const int MOD = 1000000007;
const int INV2 = (MOD + 1) / 2;
const int INV6 = (MOD + 1) / 6;

int n, lim, ans, pnt;
int phi[N], phid[N];
int prime[N], vis[N];
map<int, int> mp;

inline void Add(int &x, int a) {
    x += a; while (x >= MOD) x -= MOD;
}
inline int sum2(int n) {
    return (ll)n * (n + 1) % MOD * (2 * n + 1) % MOD * INV6 % MOD;
}
inline int sum2(int l, int r) {
    return (sum2(r) - sum2(l - 1) + MOD) % MOD;
}
inline int sum3(int n) {
    int s = (ll)n * (n + 1) / 2 % MOD;
    return (ll)s * s % MOD;
}
void pre(int n) {
    int x; phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) {
            prime[++pnt] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= pnt; j++) {
            if ((x = prime[j] * i) > n) break;
            if (i % prime[j]) phi[x] = phi[i] * phi[prime[j]];
            else { phi[x] = phi[i] * prime[j]; break; }
        }
    }
    for (int i = 1; i <= n; i++)
        Add(phid[i], phid[i - 1] + (ll)phi[i] * i * i % MOD);
}
int S(int n) {
    if (n <= lim) return phid[n];
    if (mp.count(n)) return mp[n];
    int res = sum3(n), pos;
    for (int i = 2; i <= n; i = pos + 1) {
        pos = n / (n / i);
        Add(res, MOD - (ll)S(n / i) * sum2(i, pos) % MOD);
    }
    return mp[n] = res;
}

int main(void) {
    freopen("1.in", "r", stdin);
    freopen("1.out", "w", stdout);
    cin >> n; pre(lim = pow(n, 0.666));
    ans = n;
    int lst = 0, cur, pos;
    for (int i = 1; i <= n; i = pos + 1) {
        pos = n / (n / i);
        cur = S(pos);
        Add(ans, (ll)(cur - lst + MOD) * (n / i) % MOD);
        lst = cur;
    }
    ans = (ll)ans * INV2 % MOD;
    cout << ans << endl;
    return 0;
}