洛谷3768 简单的数学题

题目大意：求$\sum{i = 1}^{n} \sum{j = 1}^{n} ij\gcd(i, j) \mod p$。

题解

遇到gcd的题目不仅仅可以反演，也可以采用一些其他方法。
这道题目就是，直接反演非常麻烦，所以采用下面的变换：

$$\begin{aligned} Ans &= \sum_{d=1}^{n} d \sum_{i = 1}^{n} \sum_{j = 1}^{n} ij[\gcd(i, j)=d] \\ &= \sum_{d=1}^{n} d^3 \sum_{i = 1}^{n} \sum_{j = 1}^{n} \frac{i}{d}\frac{j}{d}[\gcd(\frac{i}{d}, \frac{j}{d})=1] \\ &= \sum_{d=1}^{n} d^3 \sum_{i = 1}^{\lfloor \frac{n}{d}\rfloor} \sum_{j = 1}^{\lfloor \frac{n}{d}\rfloor} ij[\gcd(i, j)=1] \\ &= \sum_{d=1}^{n} d^3 \sum_{i = 1}^{\lfloor \frac{n}{d}\rfloor}i \sum_{j = 1}^{\lfloor \frac{n}{d}\rfloor} j[\gcd(i, j)=1] \\ &= \sum_{d=1}^{n} d^3 \sum_{i = 1}^{\lfloor \frac{n}{d}\rfloor}i\cdot i\phi(i) \\ &= \sum_{d=1}^{n} d^3 \sum_{i = 1}^{\lfloor \frac{n}{d}\rfloor}i^2 \phi(i) \end{aligned}$$

设$S(n)=\sum_{i = 1}^{n}i^2 \phi(i)$，那么有

$$\begin{aligned} S(n) &= \sum_{i = 1}^{n}i^2 (i - \sum_{d\mid i, d < i}\phi(d)) \\ &= \sum_{i = 1}^{n}i^3 - \sum_{i = 1}^{n}\sum_{d\mid i, d < i}\phi(d) \\ &= \sum_{i = 1}^{n}i^3 - \sum_{i = 2}^{n}\sum_{d=1}^{\lfloor \frac{n}{i}\rfloor}(id)^2 \phi(d) \\ &= \sum_{i = 1}^{n}i^3 - \sum_{i = 2}^{n} i^2 \sum_{d=1}^{\lfloor \frac{n}{i}\rfloor}d^2 \phi(d) \\ &= \sum_{i = 1}^{n}i^3 - \sum_{i = 2}^{n} i^2 S(\lfloor \frac{n}{i}\rfloor) \end{aligned}$$

用杜教筛，时间复杂度$O(n^{\frac{2}{3}})$。

貌似这种做法的底层原理还是莫比乌斯反演加上狄利克雷卷积，这个其实我没有太多的考证（其实是还没学到）。。。
不过运气比较好，这道题的两个维度都是$n$，如果是$m$和$n$就很难办了。。。

#include <bits/stdc++.h>
#define INF 2000000000
#define MAXN 4641600
using namespace std;
typedef long long ll;
ll n;
int phi[MAXN + 5], tot = 0, prime[MAXN >> 1], p, inv2, inv4, inv6;
int sum[MAXN + 5] = {0}, sum2[MAXN + 5];
bool vis[MAXN + 5] = {0};
map<ll, int> mp;
int poww(int a, int b, int c){
    int res = 1;
    while(b){
        if(b & 1) res = 1ll * res * a % c;
        a = 1ll * a * a % c, b >>= 1;
    }
    return res;
}
inline int poww3(int r){
    int res = 1ll * r * r % p;
    int res2 = 1ll * (r + 1) * (r + 1) % p;
    res = 1ll * res * res2 % p;
    res = 1ll * res * inv4 % p;
    return res;
}
inline int poww2(int r){
    int res = 1ll * r * (r + 1) % p;
    int res2 = 1ll * (r + r + 1) * inv6 % p;
    res = 1ll * res * res2 % p;
    return res;
}
inline int q3(int l, int r){
    int res = poww3(r) - poww3(l - 1);
    if(res < 0) res += p;
    return res;
}
inline int q2(int l, int r){
    int res = poww2(r) - poww2(l - 1);
    if(res < 0) res += p;
    return res;
}
void getPhi(int N){
    phi[1] = sum[1] = sum2[1] = 1;
    for(int i = 2; i <= N; ++i){
        if(!vis[i])
            prime[++tot] = i, phi[i] = i - 1;
        for(int j = 1; j <= tot; ++j){
            ll t = 1ll * i * prime[j];
            if(t > 1ll * N) break;
            vis[t] = true, phi[t] = phi[i] * (prime[j] - 1);
            if(i % prime[j] == 0){
                phi[t] += phi[i];
                break;
            }
        }
        sum[i] = sum[i - 1] + phi[i];
        if(sum[i] >= p) sum[i] -= p;
        int term = 1ll * i * i % p;
        term = 1ll * term * phi[i] % p;
        sum2[i] = sum2[i - 1] + term;
        if(sum2[i] >= p) sum2[i] -= p;
    }
}
int solve_s(ll N){
    if(N <= MAXN) return sum2[N];
    if(mp.count(N)) return mp[N];
    int ans = poww3(N % p);
    for(ll i = 2, lst = 0; i <= N; i = lst + 1){
        lst = N / (N / i);
        int distract = 1ll * q2(i % p, lst % p) * solve_s(N / i) % p;
        ans -= distract;
        if(ans < 0) ans += p;
    }
    mp[N] = ans;
    return ans;
}
void init(){
    scanf("%d%lld", &p, &n);
    getPhi(MAXN);
    inv2 = poww(2, p - 2, p);
    inv4 = poww(4, p - 2, p);
    inv6 = poww(6, p - 2, p);
}
void solve(){
    int ans = 0;
    for(ll i = 1, lst = 0; i <= n; i = lst + 1){
        lst = n / (n / i);
        int term = 1ll * q3(i % p, lst % p) * solve_s(n / i) % p;
        ans += term;
        if(ans >= p) ans -= p; 
    }
    printf("%d\n", ans);
}
int main(){
    init();
    solve();
    return 0;
}