[Luogu3768]简单的数学题

Decription

给定$n$，求

$$\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)$$

Solution

$$\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)$$
$$=\sum d^3\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor}ij[gcd(i,j=1)]$$

设$S(n)=\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j=1)]$

则原式

$$=\sum d^3G(\lfloor \frac{n}{d}\rfloor)$$

$$S(n)=\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j=1)]$$
$$=\sum_{i=1}^ni\sum_{j=1}^nj[gcd(i,j=1)]$$
$$=2\sum_{i=1}^ni\sum_{j=1}^ij[gcd(i,j=1)]-1$$
$$=2\sum_{i=1}^ni\frac{i\phi(i)+[i=1]}{2}-1$$
$$=\sum_{i=1}^ni^2i\phi(i)$$

设$f(i)=i^2\phi(i),g(i)=i^2$,则$S(n)=\sum_{i=1}^nf(i)$

最后用杜教筛求出$S(n)$

$$S(n)=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)S(i)=\sum_{i=1}^ni^3-\sum_{i=2}^ng(i)S(i)=(\frac{n*(n+1)}{2})^2-\sum_{i=2}^ng(i)S(i)$$

#include <bits/stdc++.h>
using namespace std;

const int maxn = 5000005, Max = 5000000;
typedef long long lint;

int p, inv1, inv2;
lint n;

inline int pow(int x, int k)
{
    int ret = 1;
    while (k) {
        if (k & 1) ret = (lint)ret * x % p;
        x = (lint)x * x % p; k >>= 1;
    }
    return ret;
}

int s[maxn], phi[maxn], pr[maxn], pcnt;
bool isp[maxn];
void prepare(int n)
{
    memset(isp, -1, sizeof(isp));
    isp[1] = 0; phi[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (isp[i]) {phi[i] = i - 1; pr[++pcnt] = i;}
        for (int j = 1; i * pr[j] <= n && j <= pcnt; ++j) {
            isp[i * pr[j]] = 0;
            if (i % pr[j]) phi[i * pr[j]] = (lint)phi[i] * phi[pr[j]] % p;
            else phi[i * pr[j]] = (lint)phi[i] * pr[j] % p;
        }
    }
    for (int i = 1; i <= n; ++i) {
        s[i] = (lint)i * i % p * phi[i] % p;
        s[i] += s[i - 1];
        if (s[i] >= p) s[i] -= p;
    }
}

bool vis[maxn];
lint S[maxn];

inline lint calc1(lint x)
{
    x %= p;
    return (lint)x * (x + 1) % p * (2 * x + 1) % p * inv1 % p;
}

inline lint calc2(lint x)
{
    x %= p;
    x = x * (x + 1) % p * inv2 % p;
    return x * x % p;
}

lint dfs(lint x)
{
    if (x < Max) return s[x];
    int m = n / x;
    if (vis[m]) return S[m];
    vis[m] = true;
    lint sum = calc2(x);
    for (lint i = 2, j; i <= x; i = j + 1) {
        j = x / (x / i);
        (sum -= (calc1(j) - calc1(i - 1)) * dfs(x / j) % p) %= p;
    }
    if (sum < 0) sum += p;
    return S[m] = sum;
}

int main()
{
    scanf("%d%lld", &p, &n); inv1 = pow(6, p - 2); inv2 = pow(2, p - 2);

    prepare(Max);
    dfs(n);

    lint ans = 0;
    for (lint i = 1, j; i <= n; i = j + 1) {
        j = n / (n / i);
        if (n / i <= Max) ans += (lint)(calc2(j) - calc2(i - 1) + p) % p * s[n / i] % p;
        else ans += (lint)(calc2(j) - calc2(i - 1) + p) % p * S[i] % p;
        if (ans >= p) ans -= p;
    }

    printf("%lld\n", ans);

    return 0;
}