bzoj4407: 于神之怒加强版

莫比乌斯反演

---

题目链接

bzoj4407: 于神之怒加强版

题解

求这个东西

$$\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k$$
然后是套路
$$\begin{aligned}\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k &=\sum_{d=1}^{min(n,m)}d^k\sum_{i=1}^n\sum_{j=1}^m\left[(i,j)=d\right]\\ &=\sum_{d=1}^{min(n,m)}d^k\sum_{i=1}^{min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\mu(i)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor\\ &=\sum_{T=1}^{min(n,m)}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{d\mid T}d^k\mu(\frac{T}{d})\end{aligned}$$
显然，$\sum_{d\mid T}d^k\mu(\frac{T}{d})$它是个积性函数
尝试线筛
另$F(T) = \sum_{d\mid T}d^k\mu(\frac{T}{d})$
考虑一个只有一个质因数的情况$F(p_i^{k_i}) = (p_i^{ki})^k - (p_i^{k_i - 1})^k$
所以$F(p_i^{k_i}\times p_i)=F(p_i^{k_i})\times p_i^{k}$
对于不同质因数,因为是积性函数,直接乘起来

代码

#include<cstdio> 
#include<algorithm> 
inline int read() { 
     int x = 0,f = 1 ; 
     char c = getchar(); 
     while(c < '0' || c > '9')c = getchar(); 
     while(c <= '9' && c >= '0') x = x * 10 + c - '0',c = getchar(); 
     return x * f; 
} 
const int mod = 1e9 + 7; 
int t,k; 
const int maxn = 5000007; 
int num,p[maxn],mu[maxn],F[maxn],tmp[maxn]; 
bool np[maxn]; 
int fstpow(int x,int k) {
    int ret = 1; 
    for(;k;k >>= 1,x = 1ll * x * x % mod) 
        if(k & 1) ret = 1ll * ret * x % mod; 
    return ret; 
} 
const int M = 5000000; 
void pre(int k) { 
    F[1] = mu[1] = 1; 
    for(int i = 2;i <= M;++ i) { 
        if(!np[i]) p[++ num] = i, mu[i] = -1, tmp[num] = fstpow(i,k),F[i] = (tmp[num] - 1) % mod; 
        for(int t,j = 1;j <= num && i * p[j] <= M;++ j) { 
            t = i * p[j]; 
            np[t] = 1; 
            if(i % p[j]) mu[t] = -mu[i],F[t] = 1ll * F[i] * F[p[j]] % mod; 
            else {mu[t] = 0,F[t] = 1ll * F[i] * tmp[j] % mod; } 
        } 
    } 
    for(int i = 2;i <= M;++ i) F[i] += F[i - 1],F[i] >= mod && (F[i] -= mod); 
} 
long long solve(int n,int m) { 
    long long ret = 0;  
    for(int i = 1,nxt;i <= std::min(n,m);i = nxt + 1) { 
        nxt = std::min(n / (n / i),m / (m / i)); 
        ret += 1ll * (F[nxt] - F[i - 1] + mod) % mod * (n / i) % mod * (m / i) % mod; 
    } 
    return ret; 
} 
int main() { 
    t = read(),k = read(); 
    pre(k); 
    for(int n,m,i = 1;i <= t;++ i) { 
        n = read(),m = read(); 
        printf("%lld\n",(solve(n,m) % mod + mod) % mod);  
    } 
    return 0; 
}