题目大意
给定 n , m , k n,m,k n,m,k,计算 ∑ i = 1 n ∑ j = 1 m g c d ( i , j ) k \sum_{i=1}^n \sum_{j=1}^mgcd(i,j)^k ∑i=1n∑j=1mgcd(i,j)k
反演部分
方便起见设 n ≤ m , g c d ( i , j ) ⇔ ( i , j ) n\leq m,gcd(i,j) \Leftrightarrow (i,j) n≤m,gcd(i,j)⇔(i,j)
常用的转化思路,枚举
g
c
d
gcd
gcd:
∑
i
=
1
n
∑
j
=
1
m
g
c
d
(
i
,
j
)
k
=
∑
i
=
1
n
∑
j
=
1
m
∑
d
=
1
n
d
k
[
g
c
d
(
i
,
j
)
=
d
]
\sum_{i=1}^n \sum_{j=1}^mgcd(i,j)^k=\sum_{i=1}^n \sum_{j=1}^m\sum_{d=1}^nd^k[gcd(i,j)=d]
∑i=1n∑j=1mgcd(i,j)k=∑i=1n∑j=1m∑d=1ndk[gcd(i,j)=d]
然后将后面部分提到前面:
∑ d = 1 n d k ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = d ] \sum_{d=1}^nd^k\sum_{i=1}^n \sum_{j=1}^m[gcd(i,j)=d] ∑d=1ndk∑i=1n∑j=1m[gcd(i,j)=d]
根据 g c d ( i , j ) = d ⇔ g c d ( ⌊ i d ⌋ , ⌊ i d ⌋ ) = 1 gcd(i,j)=d \Leftrightarrow gcd(\lfloor \frac{i}{d}\rfloor,\lfloor \frac{i}{d}\rfloor)=1 gcd(i,j)=d⇔gcd(⌊di⌋,⌊di⌋)=1,并转化到两个 ∑ \sum ∑上,得:
∑ d = 1 n d k ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ [ g c d ( i , j ) = 1 ] \sum_{d=1}^nd^k\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor} \sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}[gcd(i,j)=1] ∑d=1ndk∑i=1⌊dn⌋∑j=1⌊dm⌋[gcd(i,j)=1]
然后根据 [ g c d ( i , j ) = 1 ] ⇔ [ 1 n ] ⇔ ∑ d ∣ n μ ( d ) [gcd(i,j)=1] \Leftrightarrow [\frac{1}{n}] \Leftrightarrow \sum_{d|n} \mu(d) [gcd(i,j)=1]⇔[n1]⇔∑d∣nμ(d),得:
∑ d = 1 n d k ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ ∑ t ∣ ( i , j ) μ ( t ) \sum_{d=1}^nd^k\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor} \sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}\sum_{t|(i,j)} \mu(t) ∑d=1ndk∑i=1⌊dn⌋∑j=1⌊dm⌋∑t∣(i,j)μ(t)
然后是关键的一步,改为枚举 t t t,然后联系整除分块得:
∑ d = 1 n d k ∑ t = 1 ⌊ n d ⌋ μ ( t ) ∑ i = 1 , t ∣ i ⌊ n d ⌋ ∑ j = 1 , t ∣ j ⌊ m d ⌋ 1 = ∑ d = 1 n d k ∑ t = 1 ⌊ n d ⌋ μ ( t ) ⌊ n d t ⌋ ⌊ m d t ⌋ \sum_{d=1}^nd^k \sum_{t=1}^{\lfloor \frac{n}{d}\rfloor}\mu(t)\sum_{i=1,t|i}^{\lfloor \frac{n}{d}\rfloor} \sum_{j=1,t|j}^{\lfloor \frac{m}{d}\rfloor}1 = \sum_{d=1}^nd^k \sum_{t=1}^{\lfloor \frac{n}{d}\rfloor}\mu(t)\lfloor \frac{n}{dt}\rfloor\lfloor \frac{m}{dt}\rfloor ∑d=1ndk∑t=1⌊dn⌋μ(t)∑i=1,t∣i⌊dn⌋∑j=1,t∣j⌊dm⌋1=∑d=1ndk∑t=1⌊dn⌋μ(t)⌊dtn⌋⌊dtm⌋
令 T = d t T=dt T=dt,然后改为枚举 T T T,最后化简为 ∑ T = 1 n ⌊ n T ⌋ ⌊ m T ⌋ ∑ d ∣ T d k μ ( T d ) \sum_{T=1}^n \lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor \sum_{d|T}d^k\mu(\frac{T}{d}) ∑T=1n⌊Tn⌋⌊Tm⌋∑d∣Tdkμ(dT)
线性筛部分
对于上述式子的最后一部分: f ( n ) = ∑ d ∣ T d k μ ( T d ) f(n)=\sum_{d|T}d^k\mu(\frac{T}{d}) f(n)=∑d∣Tdkμ(dT),因为幂函数 i d k ( n ) id_k(n) idk(n)和莫比乌斯函数 μ ( n ) \mu(n) μ(n)都是积性函数,根据狄利克雷卷积,那么二者的乘积 f = i d k ∗ μ f=id_k*\mu f=idk∗μ也是一个积性函数,积性函数可以考虑欧拉筛埃氏筛,对于 5 e 6 5e6 5e6的数据范围, n l o g n nlogn nlogn的筛法会 T L E TLE TLE,只能想办法线性筛
蒟蒻第一次筛奇怪神奇的积性函数,参考了这篇博客
- 当 n = p , p ∈ p r i m e s n=p,p\in primes n=p,p∈primes, f ( p ) = 1 k ∗ μ ( p ) + p k ∗ μ ( 1 ) = p k − 1 f(p)=1^k*\mu(p)+p^k*\mu(1)=p^k-1 f(p)=1k∗μ(p)+pk∗μ(1)=pk−1
- 当
n
=
p
x
,
p
∈
p
r
i
m
e
s
n=p^x,p \in primes
n=px,p∈primes,注意到莫比乌斯函数只要含有平方项以上的素数就为0,那么
f
(
p
)
=
(
p
k
)
x
−
1
∗
μ
(
p
)
+
(
p
k
)
x
∗
μ
(
1
)
=
(
p
k
)
x
−
1
(
p
k
−
1
)
f(p)=(p^k)^{x-1}*\mu(p)+(p^k)^x*\mu(1)=(p^k)^{x-1}(p^k-1)
f(p)=(pk)x−1∗μ(p)+(pk)x∗μ(1)=(pk)x−1(pk−1)
再考虑 f ( p x + 1 ) f(p^{x+1}) f(px+1),显然: f ( p x + 1 ) = ( p k ) x ( p k − 1 ) = f ( p x ) p k f(p^{x+1})=(p^k)^{x}(p^k-1)=f(p^x)p^k f(px+1)=(pk)x(pk−1)=f(px)pk
有了上面的结论,设 x , y ; x ∈ Z + , y ∈ p r i m e s x,y; x\in Z_+,y\in primes x,y;x∈Z+,y∈primes,考虑如何线性筛:
- 当 y ∣ x y|x y∣x时,取 x x x中 y y y的幂次 c c c,由 f f f的积性有 f ( x y ) = f ( x y y c + 1 ) f ( y c + 1 ) f(xy)=f(\frac{xy}{y^{c+1}})f(y^{c+1}) f(xy)=f(yc+1xy)f(yc+1),又因为 f ( p c + 1 ) = f ( p c ) p k f(p^{c+1})=f(p^c)p^k f(pc+1)=f(pc)pk,得 f ( x y ) = f ( x y c ) f ( y c ) y k = f ( x ) y k f(xy)=f(\frac{x}{y^c})f(y^c)y^k=f(x)y^k f(xy)=f(ycx)f(yc)yk=f(x)yk
- 当 y ∤ x y \nmid x y∤x时,由积性函数的性质得 f ( x y ) = f ( x ) f ( y ) f(xy)=f(x)f(y) f(xy)=f(x)f(y)
代码
//
// Created by Happig on 2020/9/25
//
#include <bits/stdc++.h>
#include <unordered_map>
#include <unordered_set>
using namespace std;
#define fi first
#define se second
#define pb push_back
#define ins insert
#define Vector Point
#define ENDL "\n"
#define lowbit(x) (x&(-x))
#define mkp(x, y) make_pair(x,y)
#define mem(a, x) memset(a,x,sizeof a);
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
typedef pair<double, double> pdd;
const double eps = 1e-8;
const double pi = acos(-1.0);
const int inf = 0x3f3f3f3f;
const double dinf = 1e300;
const ll INF = 1e18;
const int Mod = 1e9 + 7;
const int maxn = 5e6 + 10;
vector<int> prime;
bitset<maxn> vis;
ll f[maxn], sum[maxn];
int k;
ll qkp(ll x, ll n, ll p) {
x %= p;
ll ans = 1;
while (n) {
if (n & 1) ans = ans * x % p;
x = x * x % p;
n >>= 1;
}
return ans;
}
void init() {
memset(sum, 0, sizeof sum);
f[1] = 1;
vis.reset();
for (int i = 2; i < maxn; i++) {
if (!vis[i]) {
prime.push_back(i);
f[i] = (qkp(i, k, Mod) - 1 + Mod) % Mod;
}
for (int j = 0; j < prime.size() && i * prime[j] < maxn; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j]) {
f[i * prime[j]] = f[i] * f[prime[j]] % Mod;
} else {
f[i * prime[j]] = f[i] * qkp(prime[j], k, Mod) % Mod;
break;
}
}
}
for (int i = 1; i < maxn; i++) {
sum[i] = (sum[i - 1] + f[i]) % Mod;
}
}
int main() {
//freopen("in.txt","r",stdin);
//freopen("out.txt", "w", stdout);
ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);
int t, n, m;
cin >> t >> k;
init();
while (t--) {
cin >> n >> m;
ll up = min(n, m), ans = 0;
for (ll l = 1, r; l <= up; l = r + 1) {
r = min(n / (n / l), m / (m / l));
if (r > up) r = up;
ans += (sum[r] - sum[l - 1] + Mod) % Mod * (n / l) % Mod * (m / l) % Mod;
ans %= Mod;
}
cout << ans << ENDL;
}
return 0;
}