Loj #572. 「LibreOJ Round #11」Misaka Network 与求和（莫比乌斯反演 + 杜教筛 + min_25筛（递推版））

---

直接反演一下：$$\sum _{i=1}^n\sum _{i=1}^{n}f(gcd(i,j){)}^k$$$$=\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }f(d{)}^k\ast [gcd(i,j)=1]$$$$=\sum _{d=1}^{n}f(d{)}^k\sum _{p=1}^{\lfloor \frac{n}{d}\rfloor }\mu (p)\ast {\lfloor \frac{n}{p\ast d}\rfloor }^2$$$$=\sum _{T=1}^n{\lfloor \frac{n}{T}\rfloor }^2\sum _{d|T}f(d{)}^k\mu (\frac{T}{d})$$这个式子可以分块，需要预处理右半边前缀和。
令$g(T)={\sum }_{d|T}f(d{)}^k\mu (\frac{T}{d})$，则$g(T)=f^k\ast \mu$，$\ast$代表卷积
按杜教筛的套路，将$g$函数卷上 $I$ 函数进行求和，令 $h=g\ast I=f^k$
$$\sum _{i=1}^{n}h(i)$$$$=\sum _{i=1}^n\sum _{d|i}I(d)\ast g(\frac{i}{d})$$$$=\sum _{d=1}^{n}I(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }g(i)$$$$=\sum _{d=1}^{n}I(d)S(\lfloor \frac{n}{d}\rfloor )=\sum _{d=1}^{n}S(\lfloor \frac{n}{d}\rfloor )$$提出第一项 $I(1)S(n)=S(n)$ 并移项：$$S(n)=\sum _{i=1}^{n}h(i)-\sum _{i=2}^{n}S(\lfloor \frac{n}{i}\rfloor )$$$$S(n)=\sum _{i=1}^{n}f(i{)}^k-\sum _{i=2}^{n}S(\lfloor \frac{n}{i}\rfloor )$$$f(i{)}^k$可以用min_25筛筛出来，筛法和Uoj 188相同，但是这题$f(i)$质数部分贡献不为0，先筛出合数部分再加上质数部分即可。由于递归版没有记忆化，杜教筛里每次都要计算会T，因此需要用递推版，并且由于模数是 $2^{32}$，可以用 unsigned int 类型自然溢出来取代取模运算减少常数。

---

代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 2e5 + 10;
const ll mod = 1ll << 32;
bool ispri[maxn];
ll pri[maxn], num, tot;
ll pw[maxn];
ll g[maxn], id1[maxn], id2[maxn], w[maxn], sqr;
ll sum[maxn];
ll n, k;
ll mp[maxn];
ll fpow(ll a, ll b) {
    ll r = 1;
    while (b) {
        if (b & 1)
            r = r * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return r;
}
void sieve(int n) {
    ispri[0] = ispri[1] = true;
    num = 0;
    for (int i = 2; i <= n; i++) {
        if (!ispri[i])
            pri[++num] = i;
        for (int j = 1; j <= num && i * pri[j] <= n; j++) {
            ispri[i * pri[j]] = true;
            if (i % pri[j] == 0)
                break;
        }
    }
}
/*ll S(ll x, ll y) {			//递归版，T飞
    if (pri[y] > x)
        return 0;
    ll z = x <= sqr ? id1[x] : id2[n / x];
    ll res = 0;
    for (ll i = y + 1; i <= num && pri[i] * pri[i] <= x; i++) {
        for (ll e = 1, pe = pri[i]; pe * pri[i] <= x; e++, pe *= pri[i]) {
            ll t = x / pe;
            ll k = t <= sqr ? id1[t] : id2[n / t];
            res += (S(t, i) + (g[k] - i + 1) % mod * pw[i] % mod) % mod;
            res %= mod;
        }
    }
    return res;
}*/
ll getsum(ll x) {						//杜教筛
    ll t = x <= sqr ? id1[x] : id2[n / x];
    if (mp[t] != -1)
        return mp[t];
    ll res = (sum[t] + g[t]) % mod;
    for (ll i = 2, j; i <= x; i = j + 1) {
        j = x / (x / i);
        res -= getsum(x / i) * (j - i + 1) % mod;
        if (res < 0)
            res += mod;
    }
    return mp[t] = res;
}
int main() {
    sieve(maxn - 10);
    scanf("%lld%lld", &n, &k);
    memset(mp, -1, sizeof mp);
    for (int i = 1; i <= num; i++) {
        pw[i] = fpow(pri[i], k);
    }
    sqr = sqrt(n);
    for (ll i = 1, j; i <= n; i = j + 1) {		//min_25筛离散化
        j = n / (n / i);
        w[++tot] = n / i;
        g[tot] = (w[tot] - 1) % mod;
        if (n / i <= sqr)
            id1[n / i] = tot;
        else
            id2[j] = tot;
    }
    for (int i = 1; i <= num; i++) {			//min_25筛预处理 g
        for (int j = 1; j <= tot && pri[i] * pri[i] <= w[j]; j++) {
            ll t = w[j] / pri[i];
            ll k = t <= sqr ? id1[t] : id2[n / t];
            g[j] -= (g[k] - i + 1) % mod;
            if (g[j] < 0)
                g[j] += mod;
        }
    }
    for (int i = num; i >= 1; i--) {			//min_25筛 递推版
        for (int j = 1; j <= tot && pri[i] * pri[i] <= w[j]; j++) {
            for (ll e = 1, pe = pri[i]; pe * pri[i] <= w[j]; pe = pe * pri[i]) {
                ll t = w[j] / pe;
                ll k = t <= sqr ? id1[t] : id2[n / t];
                sum[j] += (sum[k] + (g[k] - i + 1) % mod * pw[i] % mod) % mod;
                sum[j] %= mod;
            }
        }
    }
    ll ans = 0;
    for (ll i = 1, j; i <= n; i = j + 1) {
        j = n / (n / i);
        ll p = n / i % mod;
        ans += (getsum(j) - getsum(i - 1) + mod) % mod * p * p % mod;
        ans %= mod;
    }
    printf("%lld\n", ans);
    return 0;
}