XJ只有一组数据。
推式子:
[sum_{i=1}^{n} sum_{j=1}^{m} gcd(i,j)^k\
sum_{i=1}^{n}sum_{j=1}^{m} [gcd(i,j) = d] d^k\
sum_{d=1}^{n}d^k sum_{i=1}^{n} sum_{j=1}^{m} [gcd(i,j) == d]\
sum_{d=1}^{n}d^k sum_{i=1}^{frac{n}{d}} sum_{j=1}^{frac{m}{d}} [gcd(i,j) == 1]\
sum_{d=1}^{n}d^k sum_{i=1}^{frac{n}{d}} sum_{j=1}^{frac{m}{d}} sum_{v|gcd(i,j)} mu(v)\
sum_{d=1}^{n}d^k sum_{i=1}^{lfloorfrac{n}{d}
floor}mu(i)lfloorfrac{n}{di}
floorlfloorfrac{m}{di}
floor\
p = di\
sum_{p = 1}^{n} lfloorfrac{n}{p}
floor lfloorfrac{m}{p}
floor sum_{d|p}d^kmu(frac{p}{d}) \
]
前面的部分可以数论分块,后面的我们把其当作积性函数来做。
[h(p) = sum_{d|p}d^kmu(frac{p}{d})\
h(p) = prod_{i=1}^{s} h(p_i^{x_i})\
h(p_i^0) = 1\
h(p_i^1) = p_i^k - 1\
h(p_i^{x_i}) = p_i ^{x_i * k} - p_i^{x_i * k - k}\
h(i * p) = h(i) * p ^ k ( i Mod p = 0)\
]
欧拉筛即可。
#include<bits/stdc++.h>
using namespace std;
const int N = 5e6;
const int mod = 1e9 + 7;
int n,m,k;
long long ans;
int ksm(int x,int y){
int z = 1;
while(y){
if(y & 1) z = 1ll * z * x % mod;
y >>= 1;
x = 1ll * x * x % mod;
}
return z;
}
int prime[N / 10], cnt, p[N + 5];
int h[N + 5];
int main(){
int T; scanf("%d%d",&T,&k);
p[0] = p[1] = 1; h[0] = 0; h[1] = 1;
for(int i = 2; i <= N; ++ i){
if(!p[i]){
prime[++ cnt] = i;
h[i] = (ksm(i,k) - 1 + mod) % mod;
}
for(int j = 1; j <= cnt && 1ll * prime[j] * i <= N; ++ j){
p[prime[j] * i] = 1;
if(i % prime[j] == 0) { h[i * prime[j]] = 1ll * h[i] * (h[prime[j]] + 1) % mod; break; }
else h[i * prime[j]] = 1ll * h[i] * h[prime[j]] % mod;
}
}
for(int i = 1; i <= N; ++ i) h[i] = (h[i] + h[i - 1]) % mod;
while(T --){
scanf("%d%d",&n,&m);
if(n > m) swap(n,m);
ans = 0;
for(int l1 = 1, r1; l1 <= n; l1 = r1 + 1){
r1 = min(n/(n/l1), m/(m/l1));
r1 = min(r1, n);
ans += 1ll * (n/l1) * (m/l1) % mod * ((h[r1] - h[l1 - 1] + mod)%mod) % mod;
ans %= mod;
}
printf("%lld
",ans);
}
return 0;
}