LOJ6229 这是一道简单的数学题

#include <cmath>
#include <cstdio>

#define LOG(FMT...) fprintf(stderr, FMT)

using namespace std;

const int N = 1e9 + 10, N23 = 1e6 + 10;

typedef long long ll;

struct gter {
    int a[N23], b[N23];

    const int& operator[](int k) const {
        return k < N23 ? a[k] : b[N / k];
    }

    int& operator[](int k) {
        return k < N23 ? a[k] : b[N / k];
    }
};

const int P = 1e9 + 7, REV2 = 500000004, REV3 = 333333336,
          REV6 = (ll)REV2 * REV3 % P;

int n, n23, pc;
bool vis[N23];
int f[N23], g[N23], mu[N23], minp[N23];
int p[N23];
gter sf, sg, smu;

void euler();
void dyh();

int sq(int x) { return (ll)x * x % P; }
int sum1(int x) { return (ll)x * (x + 1) % P * REV2 % P; }
int sum2(int x) { return (ll)x * (x + 1) % P * (2 * x % P + 1) % P * REV6 % P; }

int main() {
    int ans = 0;
    scanf("%d", &n);
    n23 = pow(n, 2. / 3);
    euler();
    dyh();

    for (int l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans = (ans + (ll)P + (ll)(sf[r] - sf[l - 1]) % P * sq(sum1(n / l)) % P) % P;
    }
    ans = (ans + n) % P * (ll)REV2 % P;
    printf("%d\n", ans);
    return 0;
}

void euler() {
    mu[1] = 1;
    f[1] = 1;
    g[1] = 1;
    for (int x = 2; x <= n23; ++x) {
        if (!vis[x]) {
            p[++pc] = x;
            mu[x] = P - 1;
            g[x] = (ll)(P - 1) * sq(x) % P;
            f[x] = (1 + g[x]) % P;
            minp[x] = x;
        }
        for (int i = 1; x * p[i] <= n23; ++i) {
            int ind = x * p[i];
            vis[ind] = true;
            if (x % p[i] == 0) {
                // g[ind] = 0;
                // mu[ind] = 0;
                minp[ind] = minp[x] * p[i];
                if (minp[ind] != ind)
                    f[ind] = (ll)f[ind / minp[ind]] * f[minp[ind]] % P;
                else
                    f[ind] = f[x];
                break;
            } else {
                minp[ind] = p[i];
                mu[ind] = (P - mu[x]) % P;
                g[ind] = (ll)g[x] * g[p[i]] % P;
                f[ind] = (ll)f[x] * f[p[i]] % P;
            }
        }
    }
    for (int i = 1; i <= n23; ++i) {
        sf[i] = (sf[i - 1] + f[i]) % P;
        sg[i] = (sg[i - 1] + g[i]) % P;
        smu[i] = (smu[i - 1] + mu[i]) % P;
    }
}

void dyh() {
    for (int l0 = n23 + 1, r0; l0 <= n; l0 = r0 + 1) {
        r0 = n / (n / l0);
        smu[r0] = 1;
        sg[r0] = 1;
        for (int l = 2, r; l <= r0; l = r + 1) {
            r = r0 / (r0 / l);
            smu[r0] = (smu[r0] - (r - l + 1) * (ll)smu[r0 / l] % P) % P;
            if (smu[r0] < 0)
                smu[r0] += P;
            sg[r0] = (sg[r0] - (ll)(P + sum2(r) - sum2(l - 1)) % P * sg[r0 / l] % P);
            if (sg[r0] < 0)
                sg[r0] += P;
        }
        sf[r0] = sg[r0];
        for (int l = 2, r; l <= r0; l = r + 1) {
            r = r0 / (r0 / l);
            sf[r0] = (P + sf[r0] - (ll)(smu[r] - smu[l - 1]) % P * sf[r0 / l] % P) % P;
        }
    }
}