Chaos Slover 补档计划 - 数论函数变换(1)

2015.11.22 16:10 Sun| 2 visits oi_2016| ChaosSlover补档计划| Text

又开了这么一个说大不大说小不小补得慢忘得快看到题就发虚两个半个小时推公式一分钟写暴力然后忘记开 long long 挂得死去活来的坑。

我相信我写的这玩意看起来最不需要动脑了。别人(尤其是大爷,大爷就是大爷)写的 PPT 看起来那个费劲啊,还不如自己推。

莫比乌斯反演

大概就是一个函数 mu:

这个函数的两个性质

另外这个函数是一个积性函数,所以可以线性筛。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
void getmu(int n)
{
    mu[1] = 1;
    for (int i = 2; i <= n; ++i)
    {
        if (!v[i]) p[++tot] = i, mu[i] = -1;
        for (int j = 1; i * p[j] <= n; ++j)
        {
            v[i * p[j]] = true;
            if (i % p[j]) mu[i * p[j]] = -mu[i];
            else { mu[i * p[j]] = 0; break; }
        }
    }
}

还有一个定理

对于函数 $F(n)=\sum_{d|n}f(d)$,

这个定理的另一种形式:

对于函数 $F(n)=\sum_{n|d}f(d)$,

妈呀,看 PPT 的时候没看清楚那一个前提条件(这一定是在逗我)!

然后这类题目在做题的时候要先弄明白函数 f() 和 F() 的含义,然后通过求出 F() 函数的取值来反演出 f() 的取值。常常需要用到分块优化求解。

BZOJ2301 [HAOI2011]Problem b

简单地利用容斥原理转化为求满足 gcd(x, y)=k,1<=x<=n,1<=y<=m 的数对 (x, y) 的个数。

设 f(k) 表示满足 gcd(x, y)=k 的数对个数,F(k) 表示满足 k|gcd(x, y) 的数对个数。由上面的公式,可以得出(嗯,我真的比较喜欢右面的那个形式,因为好像更好记好理解一点):

很显然:

由此可见,仅当 $i*k\le\min(n,m)$ 时,函数 F() 的取值不为零。这样我们就可以简单地写出一个时间复杂度为 O(n/k) 的代码:

1
2
3
4
5
6
7
long long calc(int n, int m, int k)
{
    long long re = 0;
    for (int i = 1; i * k <= n && i * k <= m; ++i)
        re += 1ll * mu[i] * (n / (i * k)) * (m / (i * k));
    return re;
}

观察到很多时候,$\lfloor\frac{m}{x}\rfloor*\lfloor\frac{n}{x}\rfloor$ 的取值都是相同的。简单计算就可以发现,它最多有 O(sqrt(n)) 的取值。因此我们可以枚举它的每一个取值,计算相应的 i 的取值范围,之后乘上 mu 函数在这一段区间里面的和。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
long long calc(int n, int m, int k)
{
    long long re = 0;
    if (n > m) swap(n, m);
    for (int i = 1, last = 0; i * k <= n; i = last + 1)
    {
        last = min(n / (n / (i * k) * k), m / (m / (i * k) * k));
        re += 1ll * (smu[last] - smu[i - 1]) * (n / (i * k)) * (m / (i * k));
    }
    return re;
}

在这一段代码中,last 表示的是使得那一坨式子的乘号左半边和右半边取值都不改变的最大的 i。

// 万峰源那个家伙没扒明白大爷的代码,还吐槽代码慢。明明不 RE 就烧高香了啊!

// 妈呀,BZOJ 的 cin cout 变 RE 是有多该死啊!

// 我不想再吐槽我的这个博客 markdown 编辑器了。在行文中内嵌的 mathjax 代码里面如果有小星星哪怕离那么远都会被识别成强调然后变红高亮斜体显示……

 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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include <bits/stdc++.h>
using namespace std;

#define N 50005

int n, a, b, c, d, k;
int mu[N], v[N], p[N], tot, smu[N];

inline int read()
{
    int x = 0, f = 1; char ch = getchar();
    while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
    return x * f;
}

void getmu(int n)
{
    mu[1] = 1;
    for (int i = 2; i <= n; ++i)
    {
        if (!v[i]) p[++tot] = i, mu[i] = -1;
        for (int j = 1; i * p[j] <= n; ++j)
        {
            v[i * p[j]] = true;
            if (i % p[j]) mu[i * p[j]] = -mu[i];
            else { mu[i * p[j]] = 0; break; }
        }
    }
    for (int i = 1; i <= n; ++i)
        smu[i] = smu[i - 1] + mu[i];
}

long long calc(int n, int m, int k)
{
    long long re = 0;
    if (n > m) swap(n, m);
    for (int i = 1, last = 0; i * k <= n; i = last + 1)
    {
        last = min(n / (n / (i * k) * k), m / (m / (i * k) * k));
        re += 1ll * (smu[last] - smu[i - 1]) * (n / (i * k)) * (m / (i * k));
    }
    return re;
}

int main()
{
    getmu(50000);
    n = read();
    for (int i = 1; i <= n; ++i)
    {
        a = read(); b = read(); c = read(); d = read(); k = read();
        printf("%lld\n", calc(b, d, k) - calc(a - 1, d, k)
            - calc(b, c - 1, k) + calc(a - 1, c - 1, k));
    }
    return 0;
}

BZOJ2820 YY的GCD

和上一题很类似嘛。枚举那个素数,然后公式什么的和上一题一模一样。当然,i 和 p 的乘积不会大于 min(n,m) 啦,这里就直接写在上面了。

上式中左面的求和符号下面的 p 就代表素数啦!求和符号下面写一大堆什么的实在是太麻烦了呢。

然后把求和符号调换一下顺序。p 的限制条件就不需要了,所以公式就变成了下面的样子。

然而在左面的求和符号下面写两个变量看起来太山寨了(意会就好了……),还是把 i*p 替换成 k 吧。

好啦,发现右面的 F(k) 的取值已经和 p 没有关系了,直接提出来。

到此为止,我们发现这个公式的形式和刚刚的那一道题非常相近,只需要把上一道题的 mu 函数替换为右面的和式, k 的取值设成 1,就可以一样地求解了。设 $g(i)=\sum_{p|i}\mu(\frac{i}{p})$,最后一步就是如何快速求得 g(i) 的取值。

其实,函数 g(i) 是可以线性筛的。连我这种数学这么差的人都能自己推出来……

听说大爷证明了这个枚举每一个素数函数之后暴力求也是 O(n) 的。

然后这道题就结束了。写这篇题解同时理清思路用了整整五十分钟啊!到此为止,也就算是结束了莫比乌斯反演的入门吧!

 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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include <bits/stdc++.h>
using namespace std;

#define N 10000005

int t, n, m;
int tot, p[N], mu[N], sum[N], v[N];

inline int read()
{
    int x = 0, f = 1; char ch = getchar();
    while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
    return x * f;
}

void getmu(int n)
{
    mu[1] = 1;
    for (int i = 2; i <= n; ++i)
    {
        if (!v[i]) p[++tot] = i, mu[i] = -1;
        for (int j = 1; i * p[j] <= n; ++j)
        {
            v[i * p[j]] = true;
            if (i % p[j]) mu[i * p[j]] = -mu[i];
            else { mu[i * p[j]] = 0; break; }
        }
    }
    for (int i = 2; i <= n; ++i)
    {
        if (!v[i]) sum[i] = mu[1];
        for (int j = 1; i * p[j] <= n; ++j)
        {
            if (i % p[j]) sum[i * p[j]] = -sum[i] + mu[i];
            else { sum[i * p[j]] = mu[i]; break; }
        }
    }
    for (int i = 2; i <= n; ++i)
        sum[i] += sum[i - 1];
}

long long calc(int n, int m)
{
    long long re = 0;
    if (n > m) swap(n, m);
    for (int i = 1, last = 0; i <= n; i = last + 1)
    {
        last = min(n / (n / i), m / (m / i));
        re += 1ll * (sum[last] - sum[i - 1]) * (n / i) * (m / i);
    }
    return re;
}

int main()
{
    getmu(10000000);
    t = read();
    while (t--)
    {
        n = read(); m = read();
        printf("%lld\n", calc(n, m));
    }
    return 0;
}