• 杜教筛——省选前的学习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 } 
    杜教筛
  • 相关阅读:
    git 知道这些就够了
    接私活可用的 Springboot + Vue 快速开发框架
    Vue 组件传值
    Vue实现点击按钮复制功能
    vue 获取组件高度
    git commit 提交的时候报错husky > pre-commit hook failed (add --no-verify to bypass)(解决办法)
    vue中异步函数async和await的用法
    JS设置浏览器缩放比例
    CSS修改滚动条的样式
    JS代码查看浏览器页面放大比例
  • 原文地址:https://www.cnblogs.com/y-clever/p/6514901.html
Copyright © 2020-2023  润新知