题目要求给定一个数$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 }