「Bzoj 4407」于神之怒加强版(莫比乌斯反演)

BZOJ 4407
题意：
给定$n,m,k$，求
$$
\sum_{i=1}^n\sum_{j=1}^m \gcd(i,j)^k \mod (10^9+7)
$$

按照套路反演以后得到

$$
\sum_{d=1}^n d^k \sum_{i=1}^{\lfloor \frac nd \rfloor} \mu(i) \lfloor \frac n{di} \rfloor \lfloor \frac m{di}\rfloor
$$

此时我们发现可以分块套分块达到复杂度$O(n)$单次询问。
但是还是不够，那么继续变形

设$T=di$，则

$$
\sum_{T=1}^n \lfloor \frac nT \rfloor \lfloor \frac mT \rfloor \sum_{d|T} d^k \mu(\frac Td)
$$

设$f(T)=\sum_{d|T}d^k \mu(\frac Td)$，发现$d^k$和$\mu(k)$都是积性函数，所以乘积也是积性函数，所以我们可以线性筛$O(n)$算这个积性函数值。

如何线性筛？
线性筛积性函数考虑两个
第一个是$f(p), p$为质数时的取值怎么$O(1)$算
第二个是$f(p^x), p$为质数时的取值怎么$O(1)$算，也就是说我们需要从$f(p^x)$转到$f(p^{x+1})$怎么求

考虑将$T$质因数分解，那么对于每个$p_i$

$$
f(p_i^{x_i})=\sum_{d|p_i^{x_i}}d^k\mu(\frac{p_i^{x_i}}d)
$$

因为有$\mu$的存在，只需要考虑$d=p_i^{x_i}$和$d=p_i^{x_i-1}$的情况，则原式化为

$$
f(p_i^{x_i})=\mu(p_i) \cdot p_i^{k(x_i-1)} + \mu(1) \cdot p_i^{x_ik}
$$

即

$$
f(p_i^{x_i}) = p_i^{x_ik} - p_i^{k(x_i-1)}
$$

$$
f(p_i^{x_i}) = p_i^{k(x_i-1)}(p_i^k-1)
$$

这个式子可以方便解决$f(p^x), p$为质数时的取值

对于第一个，$f(p)=\mu(1) \cdot p^k + \mu(p) \cdot 1^k=p^k-1$

那么这样就可以线筛了。

#include<cstdio> 
#include<cstring>
#include<algorithm>
#include<iostream>
#include<vector>
#include<set>
#include<cmath>
#define ms(i, j) memset(i, j, sizeof i)
#define LL long long
#define db double
#define fir first
#define sec second
#define mp make_pair
using namespace std;

namespace flyinthesky {

    const int MAXN = 5000000 + 5;
    const LL MO = 1000000007;

    int T, k, cnt, pri[MAXN], vis[MAXN];
    LL f[MAXN];

    LL ksm(LL a, LL b) {
        LL bs = a % MO, ans = 1;
        while (b) {
            if (b & 1) ans = (ans * bs) % MO;
            bs = (bs * bs) % MO;
            b >>= 1;
        }
        return ans;
    }

    void clean() {
    }
    int solve() {
        scanf("%d%d", &T, &k);
        clean();
        cnt = 0, vis[1] = 1, f[1] = 1;
        for (int i = 2; i <= 5000000; ++i) {
            if (!vis[i]) pri[++cnt] = i, f[i] = ksm(i, k) - 1;
            for (int j = 1; j <= cnt && (LL)i * pri[j] <= 5000000ll; ++j) {
                vis[i * pri[j]] = 1;
                if (i % pri[j] == 0) {
                    f[i * pri[j]] = (f[i] * ksm(pri[j], k)) % MO;
                    break ;
                }
                f[i * pri[j]] = (f[i] * f[pri[j]]) % MO;
            }
        }
        for (int i = 2; i <= 5000000; ++i) f[i] += f[i - 1];
        while (T--) {
            int n, m; scanf("%d%d", &n, &m);
            if (n > m) swap(n, m);
            int l = 1;
            LL ans = 0ll;
            while (l <= n) {
                int r = min(n / (n / l), m / (m / l));
                ans = (ans + (LL)(n / l) * (m / l) % MO * (f[r] - f[l - 1]) % MO) % MO;
                l = r + 1;
            }
            printf("%lld\n", ans);
        }
        return 0; 
    }
}
int main() {
    flyinthesky::solve();
    return 0;
}