Luogu 3768 简单的数学题

莫比乌斯反演 + 杜教筛

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

$$= sum_{d = 1}^{n}dsum_{i = 1}^{n}sum_{j = 1}^{n}ij[gcd(i, j) == d]$$

$$= sum_{d = 1}^{n}d^3sum_{i = 1}^{left lfloor frac{n}{d} ight floor}sum_{j = 1}^{left lfloor frac{n}{d} ight floor}ij[gcd(i, j) == 1]$$

$$= sum_{d = 1}^{n}d^3sum_{i = 1}^{left lfloor frac{n}{d} ight floor}sum_{j = 1}^{left lfloor frac{n}{d} ight floor}ijsum_{t | i, t | j}mu(t)$$

$$= sum_{d = 1}^{n}d^3sum_{t = 1}^{left lfloor frac{n}{d} ight floor}mu(t)sum_{i = 1}^{left lfloor frac{n}{d} ight floor}i[t | i]sum_{j = 1}^{left lfloor frac{n}{d} ight floor}j[t | j]$$

$$= sum_{d = 1}^{n}d^3sum_{t = 1}^{left lfloor frac{n}{d} ight floor}mu(t)(tsum_{i = 1}^{left lfloor frac{n}{td} ight floor}i)^2$$

$$= sum_{d = 1}^{n}d^3sum_{t = 1}^{left lfloor frac{n}{d} ight floor}t^2mu(t)(frac{left lfloor frac{n}{td} ight floor(left lfloor frac{n}{td} ight floor + 1)}{2})^2$$

记$T = td$，

$$sum_{T = 1}^{n}(frac{left lfloor frac{n}{T} ight floor(left lfloor frac{n}{T} ight floor + 1)}{2})^2sum_{d | T}d^3mu(frac{T}{d})frac{T^2}{d^2}$$

$$= sum_{T = 1}^{n}(frac{left lfloor frac{n}{T} ight floor(left lfloor frac{n}{T} ight floor + 1)}{2})^2T^2sum_{d | T}dmu(frac{T}{d})$$

注意到$sum_{d | n}frac{n}{d}mu(d) = varphi(n)$，

$$= sum_{T = 1}^{n}(frac{left lfloor frac{n}{T} ight floor(left lfloor frac{n}{T} ight floor + 1)}{2})^2(T^2varphi(T))$$

记$h(i) = i^2varphi(i)$，

$$= sum_{T = 1}^{n}(frac{left lfloor frac{n}{T} ight floor(left lfloor frac{n}{T} ight floor + 1)}{2})^2h(T)$$

前面那部分已经可以分块计算了，现在的主要矛盾就是$h(T)$怎么算。

考虑杜教筛的套路，我们让$h$函数卷上$g(n) = n^2$，得到

$$(h * g)(n) = sum_{d | n}d^2varphi(d)(frac{n}{d})^2 = n^2sum_{d | n}varphi(d) = n^3$$

而

$$S(n) = sum_{i = 1}^{n}i^3 - sum_{i = 2}^{n}i^2S(left lfloor frac{n}{i} ight floor)$$

根据小学奥数的知识

$$sum_{i = 1}^{n}i^3 = (sum_{i = 1}^{n}i)^2 = frac{n^2(n + 1)^2}{4}$$

$$sum_{i = 1}^{n}i^2 = frac{n(n + 1)(2n + 1)}{6}$$

然后就可以愉快地杜教筛了。

并不会算时间复杂度，反正能过。

另外注意这题的$n$有可能大于等于模数$P$，所以我疯狂取模之后效率很低。

Code：

// luogu-judger-enable-o2
#include <cstdio>
#include <cstring>
#include <unordered_map>
using namespace std;
typedef long long ll;

const int N = 4e6 + 5;
const int Maxn = 4e6;

int pCnt = 0, pri[N], phi[N];
ll P, inv2, inv4, inv6, sum[N], h[N];
bool np[N];
unordered_map <ll, ll> sH;

template <typename T>
inline void inc(T &x, T y) {
    x += y;
    if (x >= P) x -= P;
}

template <typename T>
inline void sub(T &x, T y) {
    x -= y;
    if (x < 0) x += P; 
}

inline ll fpow(ll x, ll y) {
    ll res = 1LL;
    for (; y > 0; y >>= 1) {
        if (y & 1) res = res * x % P;
        x = x * x % P;
    }
    return res;
}

inline void sieve() {
    phi[1] = 1;
    for (int i = 2; i <= Maxn; i++) {
        if (!np[i]) pri[++pCnt] = i, phi[i] = i - 1;
        for (int j = 1; j <= pCnt && i * pri[j] <= Maxn; j++) {
            np[i * pri[j]] = 1;
            if (i % pri[j] == 0) {
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            }
            phi[i * pri[j]] = phi[i] * phi[pri[j]];
        }
    }
    
    for (int i = 1; i <= Maxn; i++) {
        h[i] = 1LL * i * i % P * phi[i] % P;
        inc(sum[i], sum[i - 1]), inc(sum[i], h[i]);
    }
}

/*inline ll getPhi(ll n) {
    if (n <= 0) return 0;
    if (n <= Maxn) return sum[n];
    if (sPhi[n]) return sPhi[n];
    
    ll res = n * (n + 1) % P * inv2 % P;
    for (ll l = 2, r; l <= n; l = r + 1) {
        r = (n / (n / l));
        sub(res, (r - l + 1) * getPhi(n / l) % P);
    }
    return sPhi[n] = res;
}   */

inline ll squSum(ll n) {
    if (n <= 0) return 0;
    return n % P * ((n + 1) % P) % P * ((2LL * n % P + 1) % P) % P * inv6 % P;
}

inline ll squ(ll n) {
    n %= P;
    return n * n % P;
}

ll getH(ll n) {
    if (n <= 0) return 0LL;
    if (n <= Maxn) return sum[n];
    if (sH[n]) return sH[n];
    
    ll res = squ(n + 1) * squ(n) % P * inv4 % P;
    for (ll l = 2, r; l <= n; l = r + 1) {
        r = (n / (n / l));
        sub(res, (squSum(r) - squSum(l - 1) + P) % P * getH(n / l) % P);
//        res = (res - (squSum(r) - squSum(l - 1) + P) % P * getH(n / l) % P + P) % P;
    }
    return sH[n] = res;
}

int main() {
//    freopen("testdata.in", "r", stdin);
    
    ll n;
    scanf("%lld%lld", &P, &n);
    
    inv2 = fpow(2, P - 2), inv4 = fpow(4, P - 2), inv6 = fpow(6, P - 2);
    sieve();
    
    ll ans = 0;
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = (n / (n / l));
        inc(ans, squ(n / l) * squ(n / l + 1) % P * ((getH(r) - getH(l - 1) + P) % P) % P * inv4 % P);
    }
    
    printf("%lld
", ans);
    return 0;
}