Luogu P3768 简单的数学题

题目

Luogu P3768

给定 \(p,n\) ，求：

\[(\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j))\bmod p \]

分析

\[\begin{aligned} &\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\\ =&\sum_{i=1}^n\sum_{j=1}^n\sum_{d=1}^n ijd[\gcd(i,j)=d]\\ =&\sum_{d=1}^n d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ijd^2[\gcd(i,j)=1]\\ =&\sum_{d=1}^n d^3 \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor} ij\sum_{k\mid i,k\mid j} \mu(k)\\ =&\sum_{d=1}^n d^3\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i[k\mid i]\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}j[k\mid j] \end{aligned} \]

可以发现后面两个和式 \(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i[k\mid i]=k(1+2+\cdots+\lfloor\frac{n}{dk}\rfloor)\) ，用 \(sum(x)\) 表示 \((1+2+\cdots +x)\) 那么式子可以写成：

\[\sum_{d=1}^n d^3\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}k^2\mu(k)sum^2(\lfloor\frac{n}{dk}\rfloor) \]

设 \(T=dk\) 则有：

\[\begin{aligned} &\sum_{T=1}^n sum^2(\lfloor\frac{n}{T}\rfloor)\sum_{d\mid T} d^3(\frac{T}{d})^2\mu(\frac{T}{d})\\ =&\sum_{T=1}^n sum^2(\lfloor\frac{n}{T}\rfloor)T^2 (\mu \ast\operatorname{id})(T)\\ =&\sum_{T=1}^n sum^2(\lfloor\frac{n}{T}\rfloor)T^2 \varphi(T) \end{aligned} \]

由于 \(sum(\lfloor\frac{n}{T}\rfloor)\) 可以分块求和，所以考虑如何求 \(T^2\varphi(T)\) 的前缀和，设 \(f(x)=\operatorname{id_2}(x)\varphi(x),g(x)=\operatorname{id_2}(x)\) 那么有：

\[\begin{aligned} (f\ast g )(x)&=\sum_{d\mid x}\varphi(d)d^2(\frac{x}{d})^2\\ &=x^2 (\varphi\ast1)(x)\\ &=x^3 \end{aligned} \]

那么 \(f\ast g\) 的前缀和也容易得出为 \(sum^2(x)\) ，\(g(x)\) 的前缀和为 \(\frac {n(n+1)(2n+1)}{6}\) 都可以 \(O(1)\) 计算，所以用杜教筛：

\[g(1)S(n)=\sum_{i=1}^n (f\ast g)(i)-\sum_{i=2}^n g(i)S(\lfloor\frac{n}{i}\rfloor) \]

这样就解决了本题

代码

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

const int MAX_23 = 8000000;
bool is_prime[MAX_23 + 5];
int prime[MAX_23 + 5], s1[MAX_23 + 5];
unordered_map<ll, ll> s2;
int cnt;
ll inv6, inv2;
ll p, n;

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

ll sum(ll x)
{
    return x % p * (x + 1) % p * inv2 % p;
}

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

void sieve(int n)
{
    cnt = 0;
    memset(is_prime, 1, sizeof(is_prime));
    is_prime[1] = false;
    s1[1] = 1;
    for(int i = 2; i <= n; i++) {
        if(is_prime[i]) {
            prime[++cnt] = i;
            s1[i] = i - 1;
        }
        for(int j = 1; j <= cnt && i * prime[j] <= n; j++) {
            int t = i * prime[j];
            is_prime[t] = false;
            if(i % prime[j]) {
                s1[t] = s1[i] * s1[prime[j]];
            } else {
                s1[t] = s1[i] * prime[j];
                break;
            }
        }
    }
    for(int i = 1; i <= n; i++) {
        s1[i] = 1ll * i * i % p * s1[i] % p;
        s1[i] = (s1[i - 1] + s1[i]) % p;
    }
}


ll get_sum(ll n)
{
    if(n <= MAX_23)
        return s1[n];
    if(s2[n])
        return s2[n];
    ll res = sum(n);
    res = res * res % p;
    for(ll l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        res = (res - (sum2(r) - sum2(l - 1)) * get_sum(n / l) % p + p) % p;
    }
    s2[n] = res;
    return res;
}

int main()
{
    cin >> p >> n;
    sieve(MAX_23);
    inv6 = qpow(6, p - 2);
    inv2 = qpow(2, p - 2);
    ll ans = 0;
    for(ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ll t = sum(n / l);
        t = t * t % p;
        ll g = (get_sum(r) - get_sum(l - 1)) % p;
        ans = (ans + g * t % p) % p;
    }
    cout << (ans + p) % p << endl;
    return 0;
}