—— Summon, Light!
0%

hdu 5901 Count primes

发表于 更新于 Valine：

题意

输出[1…n]的质数个数 (1 <= n <= 1e11) 。时间限制6s，空间限制64M 。

题解

很显然，要用线性筛的话时间和空间都不够。

有一种Meissel–Lehmer算法可以解决该问题，不过本文要介绍的是一种dp解法。

定义$SR(n,p)$为，$2…n$被小于等于$p$的质数筛后剩下的数的个数；也就是说，在n的范围内，是质数，或者是只由大于$p$的质数相乘得到的数的个数。
记小于等于$n$的质数个数为$\pi(n)$，那么显然有$SR(n,n)=\pi(n)$ 。

对$SR(n,p)$分两类讨论：

  • 当$p$不是质数或者$n<p^2$时，有$SR(n,p)=SR(n,p-1)$；
  • 当$p$是质数且$n\ge p^2$时，$SR(n,p)$ 也可由 $SR(n,p-1)$ 推得：$$SR(n,p)= SR(n,p−1)− SR(\frac{n}{p}, p−1)+ SR(p−1,p−1)$$ 表示要从 $2…p-1$ 筛后剩下的数中去掉那些质因子均大于等于 $p$ 且含 $p$ 的数，因为若有小于 $p$ 的质因子则该数已经被筛去了，同理，该数除以p也一定不会被 $2…p-1$ 筛去，所以减去 $SR(\frac{n}{p}, p−1)$，与此同时 $2…p-1$ 的质数也被减掉了，所以加上 $SR(p-1,p-1)$ 。

当然因为空间限制肯定不能直接这样dp，考虑到整个过程中只用到了 $p$ 和 $\frac{n}{p}$，我们直接分为两段dp， 用H[k]表示 $k\le \sqrt{n}$ 时 $SR(\frac{n}{k},p)$ 的值，用L[k]表示 $k \le \sqrt{n}$ 时 $SR(k,p)$ 的值。

先放出代码，原作者*adkroxx*：

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
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
const int N=320000;
ll H[N],L[N];

ll SR(ll n,ll p)
{
    ll m;
    for (m=1;m*m<=n;++m) {
        H[m]=n/m-1;
    }
    for (ll i=1;i<=m;++i) {
        L[i]=i-1;
    }
    for (ll i=2;i<=m;++i){
        if (L[i]==L[i-1]) continue;
        for (ll j=1;j<=min(m-1,n/i/i);++j) {
            if(i*j<m) H[j]-=H[i*j]-L[i-1];
            else H[j]-=L[n/i/j]-L[i-1];
        }
        for (ll j=m;j>=i*i;--j) {
            L[j]-=L[j/i]-L[i-1];
        }
    }
    return H[1];
}

int main()
{
    ll n;
    while (scanf("%lld",&n)!=EOF) {
        printf("%lld\n",SR(n,n));
    }
}

L[],H[]初始化为筛内总的个数(不含1)。代码中的i用来枚举质数(L[i]==L[i-1]表示i不是质数)，然后对于 $j>i^2$ 的 $j$ 都更新一遍dp值，只是放在了两个不同的数组而已：对于 $j\le \sqrt{n}$ 的$j$更新L[]，对于 $j>\sqrt{n}$ 的 $j$ 更新H[]。只要理解了递推式代码还是很好理解的。

这个算法的时间复杂度为$O(n^{\frac{3}{4}})$（证明可以看参考资料内的下一条评论），虽然不如Meissel–Lehmer算法，但是dp的思想非常巧妙，值得学习。

后记

看到csdn和cnblogs上很多人都是把这份代码当作模板直接贴出来不讲解，不仅不尊重原作者，也有悖大家学习数据结构与算法的初衷，是以我不仅要翻译，还要加上一些自己的理解~~

参考资料

Editorial of Educational Codeforces Round 12

咖啡，亦我所需也