[BZOJ4407] 于神之怒加强版

Portal

我们先来化式子:
\[ Ans = \sum_{i = 1}^{n} \sum_{j = 1}^{m} (i, j) ^{k} \pmod {1e9 + 7}\\]
首先是套路枚举GCD + 反演
\[ Ans = \sum_{i}\sum_{j}\sum_{l = 1} ^ {min(n, m)} l ^ {k} ~\cdot~ [(i, j) == l] \= \sum_{l = 1} ^ {min(n, m)} l ^ {k} \sum_{i}\sum_{j} [(\frac{i}{l}, \frac{j}{l})= 1] \= \sum_{l = 1} ^ {min(n, m)} l ^ {k} \sum_{i}^{\frac{n}{l}}\sum_{j}^{\frac{m}{l}} [(i,j)= 1] \= \sum_{l = 1} ^ {min(n, m)} l ^ {k} \sum_{i}\sum_{j} \sum_{d} [d | i][d | j] \mu(d) \= \sum_{l}l^k \sum_{d} \mu(d)\frac{n}{ld}\frac{m}{ld} \\]
然后就是合并\(\sum\)的套路: 观察到式子中有两个sum结合起来的部分.我们可以枚举结合部,然后再枚举另一部.
\[ Ans = \sum_{l}\sum_{d} l^k\mu(d)\frac{n}{ld}\frac{m}{ld}\=\sum_{c}\sum_{p | c} p ^ k \frac{}{} \mu(\frac{c}{p}) \frac{}{} \frac{n}{c} \frac{m}{c} \]
定义\[f(n) = \sum_{d | n} d ^ k \frac{}{} \mu(\frac{n}{d}) \],发现这就是\(Id * \mu\),那么f为积性函数, 可以采用Min25线性筛筛一下前缀和.然后整除分块就做完了.

以下是线筛筛积性函数的通用办法.

明显的, 每一个数都只会被它自己的最小素因子筛掉, 那么我们记录一个mp[i]表示i质因数分解中的\(p_1 ^ {a_1}\)

那么就很好做了:

对于某一个\(prime_j\):

\((1)\) 如果 i % prime[j] != 0, F[i * prime[j]] = F[i] * F[prime[j]]

否则\((2)\) 当mp[i] == i, 表示i是一个质数的幂,我们可以通过递推或其他操作得到.否则 F[i * prime[j]] = F[i / mp[i]] * F[mp[i] * prime[j]]

这样的话,把最小质因子除掉了,他们就是互质的.

Code

#include<bits/stdc++.h>
using namespace std;
#define rep(i, a, b) for(int i = (a), i##_end_ = (b); i <= i##_end_; ++i)
#define drep(i, a, b) for(int i = (a), i##_end_ = (b); i >= i##_end_; --i)
#define clar(a, b) memset((a), (b), sizeof(a))
#define debug(...) fprintf(stderr, __VA_ARGS__)
typedef long long LL;
typedef long double LD;
int read() {
    char ch = getchar();
    int x = 0, flag = 1;
    for (;!isdigit(ch); ch = getchar()) if (ch == ‘-‘) flag *= -1;
    for (;isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
    return x * flag;
}
void write(int x) {
    if (x < 0) putchar(‘-‘), x = -x;
    if (x >= 10) write(x / 10);
    putchar(x % 10 + 48);
}

const int Mod = 1e9 + 7;
const int Maxn = 5e6 + 9;
int fpm(int base, int tims) {
    int r = 1; base %= Mod;
    while (tims) {
        if (tims & 1) r = 1ll * base * r % Mod;
        base = 1ll * base * base % Mod;
        tims >>= 1;
    }
    return r;
}

int T, k;

int prime[Maxn], isnprime[Maxn], mu[Maxn], tot, F[Maxn], low[Maxn], pk[Maxn];
int sumfixF[Maxn];
inline int getPk(int u) { return pk[u] ? pk[u] : (pk[u] = fpm(u, k)); }
void linearSieve() {
    mu[1] = 1; F[1] = 1;
    rep (i, 2, Maxn - 1) {
        if (!isnprime[i]) prime[++tot] = i, mu[i] = -1, F[i] = (getPk(i) - 1ll + Mod) % Mod, low[i] = i;
        for (int k, j = 1; j <= tot && (k = prime[j] * i) < Maxn; ++j) {
            isnprime[k] = 1;
            if (i % prime[j] == 0) {
                low[k] = low[i] * prime[j], mu[k] = 0;
                if (low[k] == k) F[k] = 1ll * F[i] * getPk(prime[j]) % Mod; 
/**/            else F[k] = 1ll * F[k / low[k]] * F[low[k]] % Mod;
                break;
            } else {
                mu[k] = mu[i] * mu[prime[j]];
                F[k] = 1ll * F[i] * F[prime[j]] % Mod;
                low[k] = prime[j];
            }
        }
    }

    rep (i, 1, Maxn - 1) sumfixF[i] = (1ll * sumfixF[i - 1] + F[i]) % Mod; 
    
}

vector <pair<int, int> > s;
void init() {
    T = read(), k = read();
    rep (i, 1, T) {
        int u = read(), v = read();
        s.push_back(make_pair(u, v));
    }
    
    linearSieve();
}

void solve() {
    rep (i, 1, T) {
        int n = s[i - 1].first, m = s[i - 1].second;
        int Limit = min(n, m); LL ans = 0;
        for (int r, l = 1; l <= Limit; l = r + 1) {
            r = min(min(n / (n / l), m / (m / l)), Limit);
            (ans += 1ll * (1ll * sumfixF[r] - sumfixF[l - 1] + Mod) % Mod * (n / l) % Mod * (m / l) % Mod) %= Mod;
        }
        printf("%d\n", ans);
    }
}

int main() {
    freopen("BZOJ4407.in", "r", stdin);
    freopen("BZOJ4407.out", "w", stdout);

    init();
    solve();

#ifdef Qrsikno
    debug("\nRunning time: %.3lf(s)\n", clock() * 1.0 / CLOCKS_PER_SEC);
#endif
    return 0;
}

[BZOJ4407] 于神之怒加强版

原文：https://www.cnblogs.com/qrsikno/p/10197221.html