BZOJ3309 DZY Loves Math <莫比乌斯反演>

Problem

DZY Loves Math

Time  Limit:  20  Sec\mathrm{Time\;Limit:\;20\;Sec}
Memory  Limit:  512  MB\mathrm{Memory\;Limit:\;512\;MB}

Description

对于正整数nn,定义f(n)f(n)nn所含质因子的最大幂指数。例如f(1960)=f(23×51×7×2)=3f(1960)=f(2^3\times5^1\times7\times2)=3, f(10007)=1f(10007)=1, f(1)=0f(1)=0
给定正整数a,ba,b,求i=1aj=1bf(gcd(i,j))\sum_{i=1}^{a}\sum_{j=1}^{b}{f(\gcd(i,j))}

Input

第一行一个数TT,表示询问数。
接下来TT行,每行两个数a,ba,b,表示一个询问。

Output

对于每一个询问,输出一行一个非负整数作为回答。

Sample Input

1
2
3
4
5
4
7558588 9653114
6514903 4451211
7425644 1189442
6335198 4957

Sample Output

1
2
3
4
35793453939901
14225956593420
4332838845846
15400094813

HINT

T104T\le 10^4
1a,b1071\le a,b\le 10^7

标签:莫比乌斯反演

Solution

好题,getget线筛新姿势。

首先套路转化出莫比乌斯函数:

Ans=i=1aj=1bf(gcd(i,j))=d=1min(a,b)i=1aj=1bf(d)[gcd(i,j)=d]=d=1min(a,b)f(d)i=1adj=1bdkgcd(i,j)μ(k)=d=1min(a,b)f(d)k=1min(a,b)dμ(k)×ad×k×bd×k=t=1min(a,b)atbtdtμ(td)×f(d)\begin{aligned} Ans&=\sum_{i=1}^{a}\sum_{j=1}^{b}f(\gcd(i,j))\\ &=\sum_{d=1}^{\min(a,b)}\sum_{i=1}^{a}\sum_{j=1}^{b}f(d)[\gcd(i,j)=d]\\ &=\sum_{d=1}^{\min(a,b)}f(d)\sum_{i=1}^{\lfloor\frac{a}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{b}{d}\rfloor}\sum_{k|\gcd(i,j)}\mu(k)\\ &=\sum_{d=1}^{\min(a,b)}f(d)\sum_{k=1}^{\lfloor\frac{\min(a,b)}{d}\rfloor}\mu(k)\times\lfloor\frac{a}{d\times k}\rfloor\times\lfloor\frac{b}{d\times k}\rfloor\\ &=\sum_{t=1}^{\min(a,b)}\lfloor\frac{a}{t}\rfloor\lfloor\frac{b}{t}\rfloor\sum_{d|t}\mu(\frac{t}{d})\times f(d)\\ \end{aligned}

这时会发现前面用根号分块很好处理,而后面的部分需要O(1)O(1)计算,所以需要线性筛预处理。

g(t)=dtμ(td)×f(d)g(t)=\sum_{d|t}\mu(\frac{t}{d})\times f(d),考虑通过μ\muff的性质找到其积性关系。

t=p1a1×p2a2××pkakt=p_1^{a_1}\times p_2^{a_2}\times\cdots\times p_k^{a_k}td=p1a1×p2a2××pkak\frac{t}{d}=p_1^{a_1'}\times p_2^{a_2'}\times\cdots\times p_k^{a_k'}d=p1a1a1×p2a2a2××pkakakd=p_1^{a_1-a_1'}\times p_2^{a_2-a_2'}\times\cdots\times p_k^{a_k-a_k'}
那么一定有0a1,a2,,ak10\le a_1',a_2',\cdots,a_k'\le 1,否则μ(td)=0\mu(\lfloor\frac{t}{d}\rfloor)=0,不计入总贡献。

  1. a1=a2==ak=max{a}a_1=a_2=\cdots=a_k=\max\{a\}
    • 对于f(d)=max{a}1f(d)=\max\{a\}-1的情况,只有一种,即a1=a2==ak=0a_1'=a_2'=\cdots=a_k'=0。而μ(td)=(1)k\mu(\frac{t}{d})=(-1)^k。故贡献为(a1)×(1)k=a×(1)k(1)k(a-1)\times(-1)^k=a\times(-1)^k-(-1)^k
    • 对于f(d)=max{a}f(d)=\max\{a\}的情况,根据组合原理,有i=0k1(1)i(ik)\sum_{i=0}^{k-1}(-1)^i\binom{i}{k},而又由二项式基本定理知i=0k(1)k(ik)=0\sum_{i=0}^{k}(-1)^k\binom{i}{k}=0,因而贡献为(1)k(kk)=(1)k(-1)^k\binom{k}{k}=(-1)^k
    • 故此情况g(d)=a×(1)k(1)k+(1)k=a×(1)kg(d)=a\times(-1)^k-(-1)^k+(-1)^k=a\times(-1)^k
  2. i,j\exists i,j使得ij,  aiaji\ne j,\;a_i\ne a_j
    • 不论f(d)=max{a}f(d)=\max\{a\}还是f(d)=max{a}1f(d)=\max\{a\}-1,都存在至少一个质因数prp_r使得ara_r'不论取00还是11f(d)f(d)的取值都没有影响。然而ara_r'0011会使得μ(td)\mu(\frac{t}{d})取到1-111,此处的贡献为f(d)+(f(d))=0f(d)+(-f(d))=0,一定全部被抵消。
    • 故此情况g(d)=0g(d)=0

综上,线性筛预处理g(t)g(t)需要知道每个数最小的质因数的次数numnum和最小质因数的幂指数次幂spsp,这样看是否有num[i×pri[j]]=num[i/sp[i]]num[i\times pri[j]]=num[i/sp[i]]即可知i×pri[j]i\times pri[j]的最小与次小质因数的次数是否相等,由此可判断是情况11还是情况22。这样先线性筛预处理后,对每次询问O(min(a,b))O(\sqrt{\min(a,b)})进行根号分块,即可达到O(Tmin(a,b))O(T\sqrt{\min(a,b)})的复杂度。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#include <bits/stdc++.h>
#define MAX_N 10000000
using namespace std;
typedef long long lnt;
template <class T> inline void read(T &x) {
x = 0; int c = getchar(), f = 1;
for (; !isdigit(c); c = getchar()) if (c == 45) f = -1;
for (; isdigit(c); c = getchar()) (x *= 10) += f*(c-'0');
}
int g[MAX_N+5], sp[MAX_N+5], num[MAX_N+5];
int pri[MAX_N+5], cnt; lnt ans;
bool NotPri[MAX_N+5];
void init(int n) {
NotPri[0] = NotPri[1] = true;
for (int i = 2; i <= n; i++) {
if (!NotPri[i]) pri[cnt++] = sp[i] = i, g[i] = num[i] = 1;
for (int j = 0; j < cnt; j++) {
if (i*pri[j] > n) break; NotPri[i*pri[j]] = true;
if (i%pri[j]) sp[i*pri[j]] = pri[j], num[i*pri[j]] = 1, g[i*pri[j]] = num[i] == 1 ? -g[i] : 0;
else sp[i*pri[j]] = sp[i]*pri[j], num[i*pri[j]] = num[i]+1,
g[i*pri[j]] = sp[i] == i ? 1 : (num[i/sp[i]] == num[i*pri[j]] ? -g[i/sp[i]] : 0);
if (i%pri[j] == 0) break;
}
}
for (int i = 1; i <= n; i++) g[i] += g[i-1];
}
int main() {
int T; read(T), init(MAX_N);
while (T--) {
int a, b; read(a), read(b), ans = 0;
for (int l = 1, r; l <= min(a,b); l = r+1)
r = min(a/(a/l), b/(b/l)), ans += 1LL*(a/l)*(b/l)*(g[r]-g[l-1]);
printf("%lld\n", ans);
}
return 0;
}