【BZOJ4407】于神之怒加强版

【题目链接】

• BZOJ4407

【前置技能】

• 线性筛
• 反演

【题解】

• $\sum _{i=1}^n\sum _{j=1}^{m}gcd(i,j{)}^k$
  $=\sum _x{x}^k\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=x]$
  $=\sum _x{x}^k\sum _{i=1}^{\lfloor \frac{n}{x}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{x}\rfloor }ϵ(gcd(i,j))$
  $=\sum _x{x}^k\sum _{i=1}^{\lfloor \frac{n}{x}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{x}\rfloor }\sum _{d|gcd(i,j)}\mu (d)$
  $=\sum _x{x}^k\sum _d\mu (d)\sum _{i=1}^{\lfloor \frac{n}{x}\rfloor }[d|i]\sum _{j=1}^{\lfloor \frac{m}{x}\rfloor }[d|j]$
  $=\sum _x{x}^k\sum _d\mu (d)\lfloor \frac{n}{xd}\rfloor \lfloor \frac{m}{xd}\rfloor$
  $=\displaystyle \sum i \lfloor\frac n i \rfloor \lfloor\frac m i \rfloor \sum{d|i} {d^k * \mu (\frac{i}{d})} $
• 令$F(i)=\sum _{d|i}{d}^k\ast \mu (\frac{i}{d})$，原式$=\sum _i\lfloor \frac{n}{i}\rfloor \lfloor \frac{m}{i}\rfloor F(i)$
• 发现$F$是两个积性函数的卷积，也是一个积性函数，可以线性筛预处理出$F$的值。然后按照$\lfloor \frac{n}{i}\rfloor \lfloor \frac{m}{i}\rfloor$的值分成$O(\sqrt{n})$组计算。
• 时间复杂度$O(N+Q\sqrt{N})$(NM同阶)

【代码】

#include<bits/stdc++.h>
#define INF 0x3f3f3f3f
#define LL  long long
#define MAXN    5000010 
#define mod 1000000007
using namespace std;
int Q, n, m, k, p[MAXN], pri[350000], cnt;
LL f[MAXN], q[MAXN];
 
template <typename T> void chkmin(T &x, T y){x = min(x, y);}
template <typename T> void chkmax(T &x, T y){x = max(x, y);}
template <typename T> void read(T &x){
    x = 0; int f = 1; char ch = getchar();
    while (!isdigit(ch)) {if (ch == '-') f = -1; ch = getchar();}
    while (isdigit(ch)) {x = x * 10 + ch - '0'; ch = getchar();}
    x *= f;
}
 
LL qpow(LL a, int b){
    LL ret = 1;
    while (b){
        if (b & 1) ret = ret * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return ret;
}
 
void init(){
    for (int i = 2; i <= 5000000; ++i){
        if (!p[i]) p[i] = i, pri[++cnt] = i;
        for (int j = 1; j <= cnt && pri[j] <= p[i] && pri[j] * i <= 5000000; ++j)
            p[pri[j] * i] = pri[j];
    }
    f[1] = 1;
    for (int i = 2; i <= 5000000; ++i){
        if (p[i / p[i]] == p[i]) q[i] = q[i / p[i]] * p[i];
        else q[i] = p[i];
        if (p[i] == i) f[i] = (mod + qpow(i, k) - 1) % mod;
        else if (q[i] == i) f[i] = f[i / p[i]] * qpow(p[i], k) % mod;
        else f[i] = f[q[i]] * f[i / q[i]] % mod;
    }
    for (int i = 2; i <= 5000000; ++i)
        f[i] = (f[i] + f[i - 1]) % mod;
}
 
LL work(){
    LL ret = 0;
    for (int i = 1; i <= n && i <= m; ++i){
        int nxt = n / (n / i), mxt = m / (m / i);
        if (nxt <= mxt){
            ret += 1ll * (n / i) * (m / i) % mod * (f[nxt] - f[i - 1]) % mod;
            ret %= mod;
            i = nxt;
        } else {
            ret += 1ll * (n / i) * (m / i) % mod * (f[mxt] - f[i - 1]) % mod;
            ret %= mod;
            i = mxt;
        }
    }
    if (ret < 0) ret += mod;
    return ret;
}
 
int main(){
    read(Q), read(k);
    init();
    while (Q--){
        read(n), read(m);
        LL ans = work();
        printf("%lld\n", ans);
    }
    return 0;
}