# 题目链接

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}

$F(T) = \sum_{d\mid T}d^k\mu(\frac{T}{d})$

# 代码

#include<cstdio>
#include<algorithm>
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() {