跳转至

杜教筛

author: lyn

杜教筛用于在低于线性时间复杂度内求出一个积性函数 \(f\) 的前缀和 \(s(n)\) 。

我们自己搞一个辅助函数 \(g\) ，设 \(h=f\cdot g\) ，则有：

\[g(1)s(n)=\sum_{i=1}^{n}h(i)-\sum_{i=2}^{n}g(i)s(\lfloor\frac{n}{i} \rfloor )\]

杜教筛公式的证明

由狄利克雷卷积的定义有 \(h(i)=\sum_{d|i}g(d)f(\frac{i}{d})\)

对 \(h\) 求和则有

\[ \begin{aligned} \sum_{i=1}^nh(i) &= \sum_{i=1}^n\sum_{d|i}g(d)f(\frac{i}{d}) \\ &= \sum_{d=1}^ng(d)\sum_{d|i}^if(\frac{i}{d}) \\ &= \sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}f(i) \\ &= \sum_{d=1}^ng(d)s(\lfloor \frac{n}{d} \rfloor) \end{aligned} \]

这其中主要变换为：

  1. 将第一步枚举 \(i\) 的因数变为对 \(1\) 到 \(n\) 中每个数分别计算其贡献。
  2. 将第二步枚举 \(i\) 的因数变为在保证不超过范围时枚举 \(i\) 为 \(d\) 的多少倍。
  3. 将对 \(f\) 求和的式子变为 \(s\)。

提出上式中 \(d=1\) 的情况，则有

\[\sum_{i=1}^nh(i)=g(1)s(n)+\sum_{d=2}^ng(d)s(\lfloor \frac{n}{d} \rfloor)\]

我们发现右边的 \(d\) 与左边的 \(i\) 并无直接联系，则可把 \(d\) 变为 \(i\) ，有：

\[g(1)s(n)=\sum_{i=1}^{n}h(i)-\sum_{i=2}^{n}g(i)s(\lfloor\frac{n}{i} \rfloor )\]

这就是杜教筛公式

杜教筛的使用要求

  1. 我们搞的 \(g\) 和 \(h\) 的前缀和要易于计算，否则你就白搞了。
  2. 我们要预处理一部分 \(s\) ，根据数学证明，处理到 \(s(n^{\frac{2}{3}})\) 时跑得最快。

例题

P4213 【模板】 杜教筛

给定一个正整数，求

\[ans_1=\sum_{i=1}^n\varphi(i)\]
\[ans_2=\sum_{i=1}^n \mu(i)\]

对于全部的测试点，保证 \(1 \leq T \leq 10\)， \(1 \leq n \lt 2^{31}\)。

欧拉函数求前缀和

考虑卷积

\[\varphi \cdot 1 = \text{id}\]

取 \(f = \varphi, g = 1, h = \text{id}\)，公式化为

\[s(n) = \frac{n(n+1)}{2} - \sum_{i=2}^ns(\lfloor\frac{n}{i}\rfloor)\]

莫比乌斯函数求前缀和

考虑卷积

\[\mu \cdot 1 = \varepsilon\]

取 \(f = \mu, g = 1, h = \varepsilon\)，公式化为

\[s(n) = 1 - \sum_{i=2}^ns(\lfloor\frac{n}{i}\rfloor)\]
实现 
 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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 5e6 + 5;
ll phi[N];
int mu[N];
int pri[N], cnt;
unordered_map<int, ll> sumphi;
unordered_map<int, int> summu;
ll queryphi(ll x)
{
    if (x < 5e6)
        return phi[x];
    if (sumphi[x])
        return sumphi[x];
    ll ans = x * (x + 1) / 2;
    for (ll l = 2, r; l <= x; l = r + 1)
    {
        r = x / (x / l);
        ans -= (r - l + 1) * queryphi(x / l);
    }
    sumphi[x] = ans;
    return ans;
}
ll querymu(ll x)
{
    if (x < 5e6)
        return mu[x];
    if (summu[x])
        return summu[x];
    ll ans = 1;
    for (ll l = 2, r; l <= x; l = r + 1)
    {
        r = x / (x / l);
        ans -= (r - l + 1) * querymu(x / l);
    }
    summu[x] = ans;
    return ans;
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);
    mu[1] = phi[1] = 1;
    for (int i = 2; i < 5e6; i++)
    {
        if (!phi[i])
        {
            pri[++cnt] = i;
            phi[i] = i - 1;
            mu[i] = -1;
        }
        for (int j = 1; i * pri[j] < 5e6 && j <= cnt; j++)
        {
            int m = i * pri[j];
            if (i % pri[j] == 0)
            {
                phi[m] = phi[i] * pri[j];
                mu[m] = 0;
                break;
            }
            phi[m] = phi[i] * phi[pri[j]];
            mu[m] = -mu[i];
        }
    }
    for (int i = 1; i < 5e6; i++)
    {
        phi[i] = phi[i - 1] + phi[i];
        mu[i] = mu[i - 1] + mu[i];
    }
    int t;
    cin >> t;
    while (t--)
    {
        int x;
        cin >> x;
        cout << queryphi(x) << " " << querymu(x) << '\n';
    }
    return 0;
}
P1587 [NOI2016] 循环之美

牛牛是一个热爱算法设计的高中生。在他设计的算法中，常常会使用带小数的数进行计算。牛牛认为，如果在 \(k\) 进制下，一个数的小数部分是纯循环的，那么它就是美的。现在，牛牛想知道：对于已知的十进制数 \(n\) 和 \(m\)，在 \(k\) 进制下，有多少个数值上互不相等的纯循环小数，可以用分数 \(\frac xy\) 表示,其中 \(1≤x≤n,1≤y≤m\)，且 \(x,y\) 是整数。一个数是纯循环的，当且仅当其可以写成以下形式：

\[a.\dot{c_1} c_2 c_3 \dots c_{p - 1} \dot{c_p}\]

其中，\(a\) 是一个整数，\(p≥1\)；对于 \(1≤i≤p\)， \(c_i\) 是 \(k\) 进制下的一位数字。

例如，在十进制下，\(0.45454545……=0.\dot {4} \dot {5}\) 是纯循环的，它可以用 \(\frac {5}{11}\)、 \(\frac{10}{22}\) 等分数表示；在十进制下，\(0.1666666……=0.1\dot6\) 则不是纯循环的，它可以用 \(\frac 16\) 等分数表示。需要特别注意的是，我们认为一个整数是纯循环的，因为它的小数部分可以表示成 \(0\) 的循环或是 \(k-1\) 的循环；而一个小数部分非 \(0\) 的有限小数不是纯循环的。

对于所有的测试点，保证 \(1\leq n\leq 10^9\)， \(1\leq m \leq 10^9\)， \(2\leq k \leq 2\times 10^3\)。

实现 
 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
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 5e6 + 5;
int n, m, k, mu[N], p[N], v[N], g[2005], f[N], cnt;
int ans;
unordered_map<int, int> ff;
int G(int x) { return x / k * g[k] + g[x % k]; }
void init()
{
    for (int i = 1; i <= k; i++)
        g[i] = g[i - 1] + (__gcd(i, k) == 1);
    mu[1] = f[1] = 1;
    for (int i = 2; i <= N - 5; i++)
    {
        if (!v[i])
        {
            p[++cnt] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && i * p[j] <= N - 5; j++)
        {
            v[i * p[j]] = 1;
            if (i % p[j] == 0)
                break;
            else
                mu[i * p[j]] = -mu[i];
        }
        f[i] = f[i - 1] + mu[i] * (G(i) - G(i - 1));
    }
}
int F(int x)
{
    if (x <= N - 5)
        return f[x];
    if (ff[x])
        return ff[x];
    int ans = 1;
    for (int l = 2, r; l <= x; l = r + 1)
    {
        r = x / (x / l);
        ans -= F(x / l) * (G(r) - G(l - 1));
    }
    ff[x] = ans;
    return ans;
}
signed main()
{
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);
    cin >> n >> m >> k;
    init();
    for (int l = 1, r; l <= min(n, m); l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        ans += (n / l) * G(m / l) * (F(r) - F(l - 1));
    }
    cout << ans;
    return 0;
}