题解P3768 简单的数学题

P3768 简单的数学题

非常典型的一道杜教筛求前缀和题目

首先我们先顺着题目给出的式子往下面推

$$\sum\limits_{i=1}^n\sum\limits_{j=1}^nijgcd(i,j)=\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^nijd[gcd(i,j)=d]\\=\sum\limits_{d=1}^nd^3\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{n}{d}}ij[gcd(i,j)=1]$$

用莫比乌斯反演,令 $g(n)=\sum\limits_{i=1}^n\sum\limits_{j=1}^nij$

$$\Rightarrow ans=\sum\limits_{d=1}^nd^3\sum\limits_{i=1}^{\frac{n}{d}}\mu(i)i^2g(\frac{n}{id})$$

令 $T=id$,尝试交换循环顺序

$$ans=\sum\limits_{T=1}^ng(\frac{n}{T})T^2\sum\limits_{d|T}d\mu(\frac{T}{d})$$

将 $\sum\limits_{d|T}d\mu(\frac{T}{d})$ 写成迪利克雷卷积的形式 $\mu*id=\phi$

$$\Rightarrow ans=\sum\limits_{T=1}^ng(\frac{n}{T})T^2\phi(T)$$

令 $f(n)=n^2\phi(n)$

$$ans=\sum\limits_{T=1}^ng(\frac{n}{T})f(T)$$

对于 $g$ 可以使用数论分块,考虑如何快速求 $S(n)=\sum\limits_{i=1}^nf(i)$

首先搬出杜教筛式子 $g(1)S(n)=\sum\limits_{i=1}^n(f*g)(i)-\sum\limits_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)$

$$(f*g)(n)=\sum\limits_{d|n}g(d)f(\frac{n}{d})=\sum\limits_{d|n}\frac{n^2}{d^2}\phi(\frac{n}{d})g(d)$$

不妨令 $g(i)=i^2$

$$(f*g)(n)=\sum\limits_{d|n}n^2\phi(\frac{n}{d})=n^2\sum\limits_{d|n}\phi(\frac{n}{d})=n^3$$

$$\Rightarrow S(n)=\sum\limits_{i=1}^ni^3-\sum\limits_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)$$

$$ans=\sum\limits_{T=1}^ng(\frac{n}{T})T^2\phi(T)$$

前面那坨可以数论分块,后面那一坨可以杜教筛 $O(n^{\frac{2}{3}})$ 处理前缀和,注意取模

#include<iostream>
#include<cstdio>
#include<unordered_map>
using namespace std;
typedef long long ll;
const int N = 8000001;

bool st[N];
int prime[N], tot;
ll sum[N], p, n, phi[N], inv2, inv6;
unordered_map<ll, ll> S;

void init()
{
    phi[1] = 1;
    for (int i = 2; i < N; i++)
    {
        if (!st[i])
            prime[++tot] = i, phi[i] = i - 1;
        for (int j = 1; i * prime[j] < N; j++)
        {
            st[i * prime[j]] = true;
            if (i % prime[j] == 0)
            {
                phi[i * prime[j]] = (phi[i] * prime[j]) % p;
                break;
            }
            phi[i * prime[j]] = (phi[i] * (prime[j] - 1)) % p;
        }
    }
    for (int i = 1; i < N; i++)
        sum[i] = (sum[i - 1] + (ll)phi[i] * i % p * i % p) % p;
}

ll Sum1(ll x)
{
    ll res = (((x + 1) % p * x % p) % p) * inv2 % p;
    return res * res % p;
}

ll Sum2(ll x)
{
    ll r1 = (x % p * (x + 1) % p) % p;
    return ((r1 * ((x + x + 1) % p)) % p * inv6) % p;
}

ll Djf(ll x)
{
    if (x < N)
        return sum[x];
    if (S[x])
        return S[x];
    ll res = Sum1(x);
    for (ll l = 2, r; l <= x; l = r + 1)
    {
        r = min(x, (x / (x / l)));
        res -= (Sum2(r) - Sum2(l - 1))* (Djf(x / l)) % p;
        res %= p;
    }
    return S[x] = (res + p) % p;
}

ll qmi(ll a, ll b)
{
    ll res = 1;
    while (b)
    {
        if (b & 1)
            res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res;
}

int main()
{
    scanf("%lld%lld", &p, &n);
    init();
    inv2 = qmi(2, p - 2), inv6 = qmi(6, p - 2);
    ll ans = 0;
    for (ll l = 1, r; l <= n; l = r + 1)
    {
        r = min(n, n / (n / l));
        ans = (ans + (Sum1(n / l) % p * ((Djf(r) - Djf(l - 1)) % p)) % p + p) % p;
    }
    printf("%lld", ans);
    return 0;
}