大素数

Posted on Edited on In Word count in article: 605 Reading time ≈ 2 mins.
十分高效的素数个数的求解

记号声明

首先，声明几个记号：

  • \(p_i\) 表示第 \(i\) 个素数
  • \(\pi(x)\) 表示 \([1,x]\) 范围内的素数个数
  • \(\Phi(m,n)\) 表示不大于 \(m\) 且不能整除任何一个 \(p_i\quad i\in[1,n]\) 的自然数个数

运算方式

那么显然有：

  • \[ \Phi(m,n)=\Phi(m-n-1)-\Phi(\frac m{p_i}, n-1) \]
  • \[ \Phi(m,0)=\left\lfloor m\right\rfloor \]

给定一个自然数 \(m\)，如果 \(n=\pi(\sqrt[3]{m})\quad\mu=\pi(\sqrt m)-n\)，那么有：

\[ \pi(m)=\Phi(m,n)+n(\mu+1)+\frac{\mu^2-\mu}2-1-\sum_{k=1}^\mu\pi\left(\frac m{p_{n+k}}\right) \]

这个方法大概可以计算 \(1e8\) 数量级的素数个数 但是有一种更快的处理方式，对于实数 \(m\)，和自然数 \(n\) 和 \(k\), 定义 \(P_k(m,n)\) 表示不大于 \(m\)，且恰好有 \(k\) 个大于 \(p_n\) 的素因子的整数个数，更进一步，设定： \(P_0(m,n)=1\)，那么有：

\[ \Phi(m,n)=\sum_{k\ge0}P_k(m,n) \]

这个和 实际上只有有限个非零的项。设 \(y\) 为一个整数，使得 \(\sqrt[3]m\le y\le\sqrt m\)，并设 \(n=\pi(y)\)。那么当 \(k \ge 3\) 时， \(P_1(m,n)=\pi(m)-n\) 且 \(P_k(m,n)=0\)。因此有：

\[ \pi(m)=\Phi(m,n)+n-1-P_2(m,n) \]

\(P_2(m,n)\) 的计算方式可以用这种方法来获得：

\[ P_2(m,n)=\sum_{y<p\le \sqrt m}\left(\pi\left(\frac mp\right)-\pi(p)+1\right) \]

这样，我们就可以在 \(O(n^{\frac 23})\) 的时间复杂度内完成 \(1\sim n\) 的素数个数的计算

例题LOJ6235

板子题。。。

Code

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
#include <cstdio>
#include <cmath>

typedef long long ll;
const int maxn = 1e5 + 10;
const int maxp = maxn * 10;
const int maxm = 105;

namespace prime
{
    int phi[maxn][maxm];
    int p[maxn], pi[maxn * 10], cnt;
    bool vis[maxn * 10];

    inline void init()
    {
        for (int i = 2; i < maxn * 10; ++i) {
            if (!vis[i]) p[++cnt] = i;
            pi[i] = cnt;
            for (int j = 1; j <= cnt && i * p[j] < maxn * 10; ++j) {
                vis[i * p[j]] = 1;
                if (i % p[j] == 0) break;
            }
        }

        for (int n = 0; n < maxm; ++n)
            for (int m = 0; m < maxn; ++m)
                if (!n) phi[m][n] = m;
                else phi[m][n] = phi[m][n - 1] - phi[m / p[n]][n - 1];
    }

    inline ll Phi(ll m, int n)
    {
        if (n == 0) return m;
        else if (p[n] >= m) return 1;
        else if (m < maxn && n < maxm) return phi[m][n];
        return Phi(m, n - 1) - Phi(m / p[n], n - 1);
    }

    inline ll Pi(ll m)
    {
        if (m < maxn) return pi[m];
        int s = sqrt(m + 0.9), y = cbrt(m + 0.9), n = pi[y];
        ll res = Phi(m, n) + n - 1;
        for (int i = n + 1; p[i] <= s; ++i) res -= Pi(m / p[i]) - Pi(p[i]) + 1;
        return res;
    }
}

ll n;
int main()
{
    prime::init();
    while (~scanf("%lld", &n))
        printf ("%lld\n", prime::Pi(n));
}