Luogu3768 简单数学题

#include <cmath>
#include <cstdio>

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

using namespace std;

typedef long long ll;

const int N = 5000010;

struct getter {
    int& operator[](ll ind) const;
} s;

int p, pc, r6, bd;
ll n;

bool vis[N];
int phi[N], sum[N], sum2[N], ps[1000010];

void euler();
void du();
int rev(int x);
void ex_gcd(int a, int b, int& x, int& y);
int sq(ll x);
int s1(ll x);
int s2(ll x);
int s3(ll x);

int main() {
    int ans = 0;
    scanf("%d%lld", &p, &n);
    r6 = rev(6);
    bd = pow(n, 2. / 3);
    euler();
    du();
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans = (ans + (ll)s3(n / l) * ((s[r] - s[l - 1]) % p) % p) % p;
        if (ans < 0)
            ans += p;
    }
    printf("%d\n", ans);
    return 0;
}

inline int sq(ll x) {
    x %= p;
    return x * x % p;
}

inline int s1(ll x) {
    x %= p;
    return x * (x + 1) / 2 % p;
}

inline int s3(ll x) {
    return sq(s1(x));
}

inline int s2(ll x) {
    x %= p;
    return x * (x + 1) % p * (x * 2 + 1) % p * r6 % p;
}

int& getter::operator[](ll ind) const {
    if (ind <= bd)
        return sum[ind];
    return sum2[n / ind];
}

int rev(int a) {
    int x, y;
    ex_gcd(a, p, x, y);
    if (x < 0)
        x += p;
    return x;
}

void ex_gcd(int a, int b, int& x, int& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return;
    }
    ex_gcd(b, a % b, y, x);
    y -= a / b * x;
}

void euler() {
    phi[1] = 1;
    for (int i = 2; i <= bd; ++i) {
        if (!vis[i]) {
            ps[++pc] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; i * ps[j] <= bd; ++j) {
            vis[i * ps[j]] = true;
            if (i % ps[j] == 0) {
                phi[i * ps[j]] = phi[i] * ps[j];
                break;
            } else
                phi[i * ps[j]] = phi[i] * (ps[j] - 1);
        }
    }
    for (int i = 1; i <= bd; ++i)
        sum[i] = (sum[i - 1] + (ll)sq(i) * phi[i]) % p;
}

void du() {
    for (ll l0 = bd, r0; l0 <= n; l0 = r0 + 1) {
        r0 = n / (n / l0);
        s[r0] = s3(r0);
        for (ll l = 2, r; l <= r0; l = r + 1) {
            r = r0 / (r0 / l);
            s[r0] = (s[r0] - (ll)s[r0 / l] * ((s2(r) - s2(l - 1)) % p) % p) % p;
            if (s[r0] < 0)
                s[r0] += p;
        }
    }
}