NJU 1017 [JSCPC2016]Heresy 莫比乌斯反演

题目链接

题意

求

$$\sum_{i=1}^{n}\sum_{j=1}^{m}i^2j^2gcd(i,j)$$
（题面上写 $n\leq 106$，事实上是 $n\leq 1e6$.
十分感谢 WuBaizhe，不然我就一直RE死不瞑目了…
莫比乌斯反演这块也是一篇一篇看着 WuBaizhe的blog学的，对初学者十分友好，每篇的推导都很详细，非常感谢原Po）

初级版(280ms)

推导

令 $gcd(i,j)=k$，则 $gcd(\frac{i}{k},\frac{j}{k})=1$，则

$$\begin{align} 原式&=\sum_{i=1}^{n}\sum_{j=1}^{m}i^2j^2k[gcd(\frac{i}{k},\frac{j}{k})=1]\\ &=\sum_{k=1}^{min(n,m)}k\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}(ik)^2(jk)^2[gcd(i,j)=1]\\ &=\sum_{k=1}^{min(n,m)}k^5\sum_{d=1}^{min(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor)}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{kd}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{kd}\rfloor}(id)^2(jd)^2\\ &=\sum_{k=1}^{min(n,m)}k^5\sum_{d=1}^{min(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor)}d^4\mu(d)\sum_{i=1}^{\lfloor\frac{n}{kd}\rfloor}i^2\sum_{j=1}^{\lfloor\frac{m}{kd}\rfloor}j^2 \end{align}$$
$O(n)$ 线性筛出 $d^4\mu(d)$，算出 $k^5$ 与 $i^2$ 的前缀和，再两次分块 $\sqrt n*\sqrt n=n$就能搞定。
预处理 $O(n)$，每组数据 $O(n)$。

Code

#include <cstdio>
#include <iostream>
#define maxn 1000010
#define maxm maxn + 10
using namespace std;
typedef long long LL;
const LL mod = 1e9 + 7;
bool check[maxm];
int prime[maxm], kas;
LL h[maxm], pre[maxm], sum5[maxm], sum2[maxm];
LL poww(LL a, LL b) {
    LL ret = 1;
    while (b) {
        if (b & 1) ret = ret * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return ret;
}
void init() {
    h[1] = 1; int tot = 0;
    for (int i = 2; i <= maxn; ++i) {
        if (!check[i]) {
            prime[tot++] = i;
            h[i] = (-poww(i, 4) + mod) % mod;
        }
        for (int j = 0; j < tot; ++j) {
            if (i * prime[j] > maxn) break;
            check[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                h[i * prime[j]] = 0;
                break;
            }
            h[i * prime[j]] = (-h[i] * poww(prime[j], 4) % mod + mod) % mod;
        }
    }
    for (int i = 1; i <= maxn; ++i) {
        pre[i] = (pre[i - 1] + h[i]) % mod;
        sum5[i] = (sum5[i - 1] + poww(i, 5) % mod) % mod;
        sum2[i] = (sum2[i - 1] + poww(i, 2) % mod) % mod;
    }
}
LL calc(int a, int b) {
    int lim = min(a, b), le, ri;
    LL ret = 0;
    for (int i = 1; i <= lim; i = ri + 1) {
        le = i, ri = min(a / (a / i), b / (b / i));
        ret = (ret + (pre[ri] - pre[le - 1] + mod) * sum2[a / i] % mod * sum2[b / i] % mod + mod) % mod;
    }
    return ret;
}
void work() {
    int n, m;
    scanf("%d%d", &n, &m);
    int lim = min(n, m), le, ri;
    LL ans = 0;
    for (int i = 1; i <= lim; i = ri + 1) {
        le = i, ri = min(n / (n / i), m / (m / i));
        ans = (ans + (sum5[ri] - sum5[le - 1] + mod) % mod * calc(n / i, m / i) % mod + mod) % mod;
    }
    printf("Case #%d: %lld\n", ++kas, ans);
}
int main() {
    init();
    int T;
    scanf("%d", &T);
    while (T--) work();
    return 0;
}

升级版(96ms)

推导

对上面的式子继续推导，令 $T=kd$，

$$\begin{align} 上式&=\sum_{k=1}^{min(n,m)}k^5\sum_{d=1}^{min(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor)}d^4\mu(d)\sum_{i=1}^{\lfloor\frac{n}{kd}\rfloor}i^2\sum_{j=1}^{\lfloor\frac{m}{kd}\rfloor}j^2\\ &=\sum_{T=1}^{min(n,m)}\sum_{i=1}^{\lfloor\frac{n}{T}\rfloor}i^2\sum_{j=1}^{\lfloor\frac{m}{T}\rfloor}j^2\sum_{d|T}(\frac{T}{d})^5d^4\mu(d)\\ &=\sum_{T=1}^{min(n,m)}\sum_{d|T}\frac{T^5}{d}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{T}\rfloor}i^2\sum_{j=1}^{\lfloor\frac{m}{T}\rfloor}j^2 \end{align}$$
令 $f(T)=\sum_{d|T}\frac{T^5}{d}\mu(d)$，显然也是可以 $O(n)$ 线性筛出的，接下来再一次分块就行。
预处理 $O(n)$，每组数据 $O(\sqrt n)$。

Code

#include <cstdio>
#include <iostream>
#define maxn 1000010
#define maxm maxn + 10
using namespace std;
typedef long long LL;
const LL mod = 1e9 + 7;
bool check[maxm];
int prime[maxm], kas;
LL h[maxm], pre[maxm], sum5[maxm], sum2[maxm];
LL poww(LL a, LL b) {
    LL ret = 1;
    while (b) {
        if (b & 1) ret = ret * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return ret;
}
void init() {
    h[1] = 1; int tot = 0;
    for (int i = 2; i <= maxn; ++i) {
        if (!check[i]) {
            prime[tot++] = i;
            h[i] = (poww(i, 5) - poww(i, 4) + mod) % mod;
        }
        for (int j = 0; j < tot; ++j) {
            if (i * prime[j] > maxn) break;
            check[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                h[i * prime[j]] = h[i] * poww(prime[j], 5) % mod;
                break;
            }
            h[i * prime[j]] = h[i] * h[prime[j]] % mod;
        }
    }
    for (int i = 1; i <= maxn; ++i) {
        pre[i] = (pre[i - 1] + h[i]) % mod;
        sum2[i] = (sum2[i - 1] + poww(i, 2) % mod) % mod;
    }
}
void work() {
    int n, m;
    scanf("%d%d", &n, &m);
    int lim = min(n, m), le, ri;
    LL ans = 0;
    for (int i = 1; i <= lim; i = ri + 1) {
        le = i, ri = min(n / (n / i), m / (m / i));
        ans = (ans + (pre[ri] - pre[le - 1] + mod) % mod * sum2[n / i] % mod * sum2[m / i] % mod) % mod;
    }
    printf("Case #%d: %lld\n", ++kas, ans);
}
int main() {
    init();
    int T;
    scanf("%d", &T);
    while (T--) work();
    return 0;
}