重学了一下莫反,感觉收获了好多。(
原来你还记得自己学过莫反)
原题链接:https://www.acwing.com/problem/content/217/
题目描述:
达达正在破解一段密码,他需要回答很多类似的问题:
对于给定的整数 a,b 和 d,有多少正整数对 x,y,满足 x≤a,y≤b,并且 gcd(x,y)=d。作为达达的同学,达达希望得到你的帮助。
我们分析一下,题目很明了,要求的就是$gcd(x , y) = d$的个数,转化一下:
$$
x’ = x / d , y’ = y / d
$$
$$
[gcd(x , y) = d] = [gcd(x’ , y’) = 1]
$$
这两个是等价的,但是看的时候我还是愣了一下(还是比较笨的),但是仔细思考下,感觉还是很容易证明的。
证明:
首先,x , y 的最大公约数是d,那么我们除去最大公约数,两个数就互质了,这个毋庸置疑,对于每个$gcd(x , y) = d$都如此操作,那么可以证明,这两个的数量是一一对应的。
(这么显然,我还要思考一会,真滴笨!OVQ)
那么我们的问题转化为:$x \le a/d$ , $y \le b/d$ , $gcd(x , y)= 1$的个数,我们考虑容斥原理:互质的数 = 总数 - 不是互质的数,
总数很容易,$a’ = a / d , b’ = b / d$ , 那么总数就是$a’ * b’$,下面分析下我们要减去的不互质的数:
$$
∑^{min(a,b)}_{i=2} [\frac{a’}{i}]∗[\frac{b’}{i}]∗mobius[i]
$$
式子上午看的,人是上午没的,我直接喵喵喵?? 愣了一会,发现显然(我是笨比).
证明:
我们要减去的是什么呢?是不互质的数,不互质数是什么呢?是最大公因数不为1的数,不为1,那可以为几呢?不是1就行(废话),我们枚举最大公因数不为1的个数即可。枚举每一个可能的最大公因数,并且计算个数。
当最大公因数是只有一个质因子的时候,我们可以列出:
$$
[\frac{a’}{2}] * [\frac{b’}{2}] + [\frac{a’}{3}] * [\frac{b’}{3}] + .....
$$
最大公因数是2的两对数,实际上就等2的倍数个数相乘, 这两个个数是等同的,同理3也是,但是我们要考虑重复的情况,如6,我们2,3,都枚举了一次,那么会重复,根据容斥原理,我们要减去,最后我们发现,其实符号就是$Mobius$函数(这个证明是显然的)。以此类推,枚举到$a’$和$b’$中较小的就可以了,因为大于较小的,则式子为0,我们就能推出上面的式子了。(被自己蠢哭了)
那么答案呼之欲出:
$$
a’b’ + ∑^{min(a,b)}_{i=2} [\frac{a’}{i}]∗[\frac{b’}{i}]∗mobius[i]\\
$$
我们把这玩意合体一下
$$
∑^{min(a,b)}_{i=1} [\frac{a’}{i}]∗[\frac{b’}{i}]∗mobius[i]
$$
(当当当当!)
然后我们发现,高兴的还是有点早,每次都是 $O(n)$,总体是$O(n^2)$,时间复杂度很大,我们忍受不了。
然后神奇的东西就出现了,我们称之为 “优雅的暴力——分块”。
下面来分析下分块,这个还好!没怎么卡。
首先,我们发现式子里面有很多取整,那么实际上数据应该是一段一段的,我们就可以考虑通过这点来加速,既然这一段都是一样的,我们直接一段一段的算,发现效率可观,可以达到$O(n\sqrt{n}) )$,下面写下分块的原理,我们定义一个 get(x) 函数,这个函数和lowbit函数一样,虽然短小,但是精悍,他主要做的事是:传回当前一样的一段中,下标最大的那个,也就是 $\frac{a}{get(x) + 1}$就会变小。
int get(int n , int x) {
return n / (n / x);
}
我们说这玩意想着就是一段一段的,那么到底多少呢?答案是:$2\sqrt{a}$
证明: 我们把1到n分成:1到$\sqrt{a}$和$\sqrt{a} + 1$到$n$,我们发现$\frac{a}{1 \sim \sqrt{a}}$,只有 $\sqrt{a}$项,并且$\frac{a}{\sqrt{a} + 1 \sim n}$取值都小于 $\sqrt{a}$,综上所述,只有$2\sqrt{a}$个。
下面证明一下:$[\frac{a}{x}] =[\frac{a}{get(x)}]$
$$
get(x) = [\frac{a}{[\frac{a}{x}]}] \le [\frac{a}{\frac{a}{x}}] = x
$$
$$
[\frac{a}{x}] \ge [\frac{a}{get(x)}]
$$
$$
[\frac{a}{get(x)}] = [\frac{a}{ [\frac{a}{[\frac{a}{x}]}] }] \ge [\frac{a}{ \frac{a}{[\frac{a}{x}]}}] = [\frac{a}{x}]
$$
因为 $[\frac{a}{x}] \le [\frac{a}{get(x)}]$且$[\frac{a}{x}] \ge [\frac{a}{get(x)}]$,那么$[\frac{a}{x}] = [\frac{a}{get(x)}]$得证。’
还有一个需要证明,这个不证明,时间复杂度没法保证:
$$
[\frac{a}{x}] \ge [\frac{a}{get(x) + 1}]
$$
设:$a = kx + r, 1 \le r < x$
即证明:
$$
k \ge [\frac{a}{[\frac{a}{k}]+ 1}]
$$
$$
k*[ [\frac{a}{k}]+ 1] \ge a
$$
设:$a = pk + q, 1 \le q < k$
得:
$$
k(p + 1) > pk + q
$$
$$
kp + k > pk + q
$$
$$
k > q
$$
得证.(呼)
那么我们就可以愉快的写代码了~
如下
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long LL;
const int N = 50010;
int tt , n , a , b , d;
int primes[N] , mobius[N] , sum[N] , cnt;
bool st[N];
void init(int n) {
mobius[1] = 1;
for (int i = 2;i <= n;i ++) {
if (!st[i]) {
primes[cnt ++ ] = i;
mobius[i] = -1;
}
for (int j = 0; primes[j] <= n / i;j ++) {
st[primes[j] * i] = true;
if(i % primes[j] == 0) {
mobius[i * primes[j]] = 0;
break;
}
mobius[i * primes[j]] = -mobius[i];
}
}
for (int i = 1;i <= n;i ++ ) sum[i] = sum[i - 1] + mobius[i];
}
int get(int a , int x) {
return a / (a / x);
}
int main() {
init(N - 1);
ios::sync_with_stdio(false);
cin.tie(0); cin >> tt;
for (int i = 0;i < tt;i ++) {
LL ans = 0;
cin >> a >> b >> d;
a = a / d , b = b / d;
n = min(a , b);
for(int l = 1 , r;l <= n;l = r + 1) {
r = min(n , min(get(a , l) , get(b , l)));
ans += (sum[r] - sum[l - 1]) * (LL)(a / l) * (LL)(b / l);
}
printf("%lld\n" , ans);
}
return 0;
}
当然,我们也可以用莫比乌斯反演直接写出来公式
$$ 若F(n) = \sum\limits_{n|d}f(d),则f(n)=\sum\limits_{n|d}\mu(\cfrac{d}{n})*F(d) $$
我们定义
$$ F(n) = \sum_{x = 1} ^ a \sum_{y = 1} ^ b [n | (x , y)] $$
$$ f(n)= \sum_{x = 1} ^ a\sum_{y = 1}[(x , y) = n] $$
可以直接反演出来:
$$
∑^{min(a,b)}_{i=1} [\frac{a’}{i}]∗[\frac{b’}{i}]∗mobius[i]
$$
数学,真是奇妙呢
orz
有个疑问,分块如果有两项相乘的话,那还能保证整个分块的时间复杂度是O(n√n)的呢?
orz
学废了
细节解释得太棒了Orz
这题最好的题解,没有之一且不接受反驳
代码上面倒数第五个数学式子应该是证:
$$ k > \lfloor \frac{a}{\lfloor \frac{a}{k} \rfloor + 1} \rfloor $$
才能得出:
$$ k * \lfloor \lfloor \frac{a}{k} \rfloor + 1 \rfloor > a $$
原文中的上式是推不出下式的。
五个月前写的,已经忘记具体细节了,不过感谢指出。
$gcd$
$get$
$\to$ $\gcd$\gcd
学到了
这篇才是最强的qwq
orz
orz