杜教筛——省选前的学习1

BZOJ 3944 ——Sum 

题目要求给定一个数$N$,$N leq 2^{32} - 1$,求

$$ans1 = sum_{i = 1}^N phi (i), ans2 = sum_{i = 1}^N mu (i) $$

——蛤?这个怎么做?我只知道线性筛。。。

先介绍$O(N)$的算法——线性筛。

由于Euler函数和Mobius函数都是积性函数(即对于$(a, b) = 1$, $f(ab) = f(a)f(b)$),可以利用类似贾志鹏线性筛的方法$O(n)$筛出1到$N$所有数的函数值。

 1 #include<cstdio>
 2 #include<cstring>
 3 typedef long long LL;
 4 const int N = 1000050;
 5 LL phi[N], mobius[N];
 6 int pr[300000];
 7 void calc(int n) {                    //计算[1, n]的phi值及mobius值 
 8   memset(phi, -1, sizeof phi);
 9   memset(mobius, -1, sizeof mobius);
10   int prnum = 0, i, j;
11   phi[1] = mobius[1] = 1;
12   for (i = 2; i < n; ++i) {
13     if (phi[i] < 0) {
14       pr[prnum++] = i;
15       phi[i] = i - 1;
16       mobius[i] = -1;
17     }
18     for (j = 0; j < prnum && (LL)i * pr[j] <= n; ++j) {
19       if (i % pr[j]) {
20         phi[i * pr[j]] = phi[i] * (pr[j] - 1);
21         mobius[i * pr[j]] = -mobius[i];
22       } else {
23         phi[i * pr[j]] = phi[i] * pr[j];
24         mobius[i * pr[j]] = 0;
25         break;
26       }
27     }
28   }
29 }
30 int main() {                          //O(n) - O(1)
31   calc(1000000);
32   for (int i = 2; i <= 1000000; ++i) {
33     phi[i] += phi[i - 1];
34     mobius[i] += mobius[i - 1];
35   }
36   long long ans1 = 0, ans2 = 0;
37   int T, n;
38   scanf("%d", &T);
39   while (T--) {
40     scanf("%d", &n);
41     printf("%lld %lld", phi[n], mobius[n]);
42   }
43   return 0;
44 } 
线性筛

更一般的,如果我们有任意积性函数$f(i)$,那么我们可以利用贾志鹏线性筛求出每个数的最小质因数及其幂次,从而(在p为素数时$f(p^k)$可以快速计算的情况下)线性递推出所有[1,N]的函数值。

那么,下面我们来介绍重点——杜教筛。

首先介绍狄利克雷卷积:若f,g均为$mathbb Z  ightarrow mathbb C$的函数(称为数论函数),定义其狄利克雷卷积为

$$(f * g)(n) = sum_{d mid n} f(d)gleft(frac n d ight)$$

数论函数集与狄利克雷卷积形成群,这意味着狄利克雷卷积满足(设数论函数集为$mathbf I$):

1. 结合律: $forall f, g, k in mathbf I, (f * g) * k = f * (g * k)$

2. 单位元: $exists epsilon in mathbf I, forall f in mathbf I, f * epsilon = epsilon * f = f$

3. 逆元:    $forall f in mathbf I, exists f in mathbf I, f * g = epsilon$

4. 封闭性: $forall f, g in mathbf I, f * g in mathbf I$

另外,狄利克雷卷积还满足交换律:

5. 交换律: $forall f, g in mathbf I, f * g = g * f$

性质2的单位函数$epsilon (n) = left[n = 1 ight] = egin{cases} 1 & n = 1\ ;0 & n > 1\ end{cases}$

有了狄利克雷卷积,我们来看一看杜教筛:

假设我们要计算$S(n) = sum_{i = 1}^N f(i)$,其中$f in mathbf I$。

那么,假设$g in mathbf I, g(1) ot= 0$, 那么

$$egin{aligned} sum_{i = 1}^N (f * g)(i) & = sum_{i = 1}^N sum_{d mid i} g(d) fleft(frac i d ight) \
& = sum_{d = 1}^N g(d) sum_{1 leq i leq N, d | i} fleft(frac i d ight) \
& = sum_{d = 1}^N g(d) sum_{1 leq i leq lfloor frac n d floor} f(i) \
& = sum_{d = 1}^N g(d) Sleft(leftlfloor frac n d ight floor ight) end{aligned}$$

$$sum_{i = 1}^N (f * g)(i)  = sum_{d = 1}^N g(d) Sleft(leftlfloor frac n d ight floor ight)$$

$$ herefore egin{equation} label{recursion} g(1)S(n) = sum_{i = 1}^N (f * g)(i)  - sum_{i = 2}^N g(i) Sleft(leftlfloor frac n i ight floor ight)end{equation}$$

如果我们能够快速地对$g$和$f * g$求和,那么由于所有的$leftlfloor frac n i ight floor$只有$O(sqrt n)$种,则对于所有$i leq sqrt n$ 和$i > sqrt n$, 计算$Sleft(lfloor frac n i floor ight)$的时间复杂度为

$$Oleft(sum_{i = 1}^{leftlfloor sqrt n ight floor} sqrt i ight) + Oleft(sum_{i = 1}^{leftlfloor sqrt n ight floor} sqrt{leftlfloor frac n i ight floor}   ight)$$

而第二项明显大一些, 它等于

$$Oleft(sum_{i = 1}^{lfloor sqrt n floor} lfloor sqrt{frac n i} floor  ight) approx Oleft(int_0^{sqrt n}sqrt{frac n x} dx ight) = O(n^{frac34})$$

如果$f$为积性函数,还可以通过线性筛预处理$S(i)$的前$n^{frac23}$项就能将复杂度降到$Oleft(n^{frac23} ight)$。

终于写完了,累死我了。蛤?还有?

考虑原题,定义$mathbf 1: mathbb Z ightarrow mathbb Z, mathbf 1(n) = 1$,易证$mathbf 1 * mu = epsilon$(貌似这就是Mobius函数的定义之一)。

代入式$( ef{recursion})$, 得

$$ans2(n) = 1  - sum_{i = 2}^n ans2left(leftlfloor frac n i ight floor ight)$$

所以只要预处理出$S$的前$n^{frac23}$项就可以$Oleft(n^{frac23} ight)$计算了。

同理,将$mathbf 1$和$ans1$代入$( ef{recursion})$,得

$$ans1(n) = frac {n(n+1)}2  - sum_{i = 2}^n ans1left(leftlfloor frac n i ight floor ight)$$

 AC代码:

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<map>
 4 typedef long long LL;
 5 typedef std::map<int, LL> Map;
 6 Map _phi, _mu;
 7 int n;
 8 const int N = 5000000;
 9 LL phi[N], mu[N];
10 int pr[1000000];
11 void init(int n) {                    //计算[1, n]的phi值及mu值 
12   memset(phi, -1, sizeof phi);
13   memset(mu, -1, sizeof mu);
14   int prnum = 0, i, j;
15   phi[1] = mu[1] = 1;
16   phi[0] = mu[0] = 0;
17   for (i = 2; i <= n; ++i) {
18     if (phi[i] < 0) {
19       pr[prnum++] = i;
20       phi[i] = i - 1;
21       mu[i] = -1;
22     }
23     for (j = 0; j < prnum && (LL)i * pr[j] <= n; ++j) {
24       if (i % pr[j]) {
25         phi[i * pr[j]] = phi[i] * (pr[j] - 1);
26         mu[i * pr[j]] = -mu[i];
27       } else {
28         phi[i * pr[j]] = phi[i] * pr[j];
29         mu[i * pr[j]] = 0;
30         break;
31       }
32     }
33   }
34   //for (i = 1; i <= 30; ++i) {
35   //  printf("%d %d
", i, phi[i]);
36   //}
37   for (i = 2; i <= n; ++i) {
38     phi[i] += phi[i - 1];
39     mu[i] += mu[i - 1];
40   }
41 }
42 LL CalcPhi(LL n) {
43   Map::iterator it; 
44   if (n < N)
45     return phi[n];
46   if ((it = _phi.find(n)) != _phi.end())
47     return it->second;
48   LL i, last, ans = (LL)n * (n + 1) >> 1;
49   for (i = 2; i <= n; i = last + 1) {
50     last = n / (n / i);
51     ans -= (last - i + 1) * CalcPhi(n / i);
52   }
53   return _phi[n] = ans;
54 }
55 LL CalcMu(LL n) {
56   Map::iterator it; 
57   if (n < N)
58     return mu[n];
59   if ((it = _mu.find(n)) != _mu.end())
60     return it->second;
61   LL i, last, ans = 1;
62   for (i = 2; i <= n; i = last + 1) {
63     last = n / (n / i);
64     ans -= (last - i + 1) * CalcMu(n / i);
65   }
66   return _mu[n] = ans;
67 }
68 int main() {
69   //freopen("sum.in", "r", stdin);
70   //freopen("sum.out", "w", stdout);
71   init(N - 1);
72   int T;
73   scanf("%d", &T);
74   while (T--) {
75     scanf("%d", &n);
76     printf("%lld %lld
", CalcPhi(n), CalcMu(n));
77   }
78   return 0;
79 } 
杜教筛
原文地址:https://www.cnblogs.com/y-clever/p/6514901.html