暴力算法寻找素数的效率是底下的,可以通过素数筛法来在一个自然数表中标记处素数。

Eratosthenes筛法

首先是Eratosthenes筛法,基本方法就是首先排除所有大于2的偶数,然后从3开始在奇数中寻找素数。具体操作就是选取一个素数,然后在数表中删去它的倍数。以3到50为例,寻找所有的素数过程如下

3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51

首先删去3的倍数(3除外)

3 5 7 11 13 17 19 23 25 29 31 35 37 41 43 47 49 51

然后删去5的倍数(5除外)

3 5 7 11 13 17 19 23 29 31 37 41 43 47 49 51

再然后7的倍数(7除外)

3 5 7 11 13 17 19 23 29 31 37 41 43 47 51

直到根号50为止,以上就是Eratosthenes筛法的过程。

程序与优化

首先发现规律,在删除3的倍数的时候,我们删去了:9、15、21、27、35、39、45,是一个9为首项6为公差的等差数列。

在删除5的倍数的时候,我们删除了:25、35,25为首项10为公差。

删除7的倍数的时候,删去了49。

由此可以发现:

每次删去某个素数p的倍数的时候,第一个删去的就是p^2,下一个就是p^2+2p、p^2+4p、...

第一个删去的是p^2而不是2*p是因为小于p^2的合数一定会被某个小于p的素数所删去,例如在删去5的倍数的时候我们从25开始,因为15已经作为3的倍数被删去了。

删去倍数的时候没有考虑p^2+p、p^2+3p、...是因为这些数是偶数,所以每次跳过一个倍数去删除。

为了作为数组用程序去实现,把以上所有数标出数组索引,

0 1 2 3   4   5   6   7   8   9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

3 5 7 11 13 17 19 23 29 31 37 41 43 47 51

以3为例:

删去的倍数有:9、15、21、27、33、39、45

对应的索引是:3、  6、  9、12、15、18、21

由此可知,在这些数中,索引为i的数的值为:value(i) = 2*i+3,反过来已知数value,其索引为 index=(value-3)/2。

我们知道,在删除第i个数的倍数中,删去的第一个元素为i的平方,它的下标是:

index(value(i)) = [(2*i+3)^2-3]/2 = 2i^2+6i+3

删去的相邻两个元素为i^2+2i、i^2+4i、...

需要算出一个数的k倍和(k+2)倍中间相差多少个元素,就需要计算两数索引的差值,

即step=index((k+2)(2*i+3)) - index(k(i*2+3)) = 2*i+3

至此,已经可以写出筛法程序了:

def make_sieve(marker, first, last, factor):
    # 标记为False,表明这个索引的数为合数
    marker[first] = False
    while last - first > factor:
        first += factor
        marker[first] = False


# 求n以内的素数
def prime_table(n):
    if n < 3:
        return [2]
    marker = [True] * n
    last = n
    i = 0  # 第0个素数
    # 最开始从索引3开始,删去3的倍数
    index_square = factor = 3
    while index_square < n:
        if marker[i]:
            # 在数表marker中以索引index_square开始,每隔factor标记一个合数,直到last为止
            # index_square为当前迭代中的第一个合数的索引,即素数p的平方的索引
            # factor即两个合数之间的索引间隔
            make_sieve(marker, index_square, last, factor)
        i += 1
        factor = 2 * i + 3
        index_square = 2 * i * (i + 3) + 3
    res = [2]
    for ind in range((n-2) // 2):
        if marker[ind]:
            res.append((ind+1)*2+1)
    return res




if __name__ == '__main__':
    print(prime_table(101))

打印0-100的素数,结果为:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

最后的优化

    目前仍有可以优化的空间,即factor和index_square的计算。

    首先观察两者在每次迭代中的增量。

    factor(i+1)-factor(i)=2

    index_square(i+1)-index_square(i)=(2i+3)+(2(i+1)+3)=factor(i)+factor(i+1)

    这样的规律就有了优化的空间,不必再在每一次迭代中通过表达式来进行计算了,可以直接通过增量(加法),代替表达式计算中的乘法,通过开销低的运算(加法)等效替代开销高的运算(乘法)。

把:

        i += 1
        factor = 2 * i + 3
        index_square = 2 * i * (i + 3) + 3

修改为:

        i += 1
        index_square += factor
        factor += 2
        index_square += factor

完整代码:

def make_sieve(marker, first, last, factor):
    # 标记为False,表明这个索引的数为合数
    marker[first] = False
    while last - first > factor:
        first += factor
        marker[first] = False


# 求n以内的素数
def prime_table(n):
    if n < 3:
        return [2]
    marker = [True] * n
    last = n
    i = 0  # 第0个素数
    # 最开始从索引3开始,删去3的倍数
    index_square = factor = 3
    while index_square < n:
        if marker[i]:
            # 在数表marker中以索引index_square开始,每隔factor标记一个合数,直到last为止
            # index_square为当前迭代中的第一个合数的索引,即素数p的平方的索引
            # factor即两个合数之间的索引间隔
            make_sieve(marker, index_square, last, factor)
        i += 1
        index_square += factor
        factor += 2
        index_square += factor
    res = [2]
    for ind in range((n-2) // 2):
        if marker[ind]:
            res.append((ind+1)*2+1)
    return res




if __name__ == '__main__':
    print(prime_table(101))

打印0-100的素数,结果为:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

C++程序如下:

#include <iostream>
using namespace std;

void make_sieve(bool* marker, int first, int last, int factor) {
    marker[first] = false;
    while (last - first > factor) {
        first += factor;
        marker[first] = false;
    }
}

int* prime_sieve(int n) {
    if (n < 3)return new int[1] {2};
    bool *marker = new bool[n];
    std::fill(marker, marker + n, true);
    int last = n;
    int i = 0;
    int index_square = 3;
    int factor = 3;
    while (index_square < n) {
        if (marker[i]) {
            make_sieve(marker, index_square, last, factor);
        }
        ++i;
        index_square += factor;
        factor += 2;
        index_square += factor;
    }
    int count = 1;
    for (int j = 0; j < ((n - 2) >> 1); j++) {
        if (marker[j]) {
            count++;
        }
    }
    int ind = 0;
    int *res = new int[count];
    res[ind++] = 2;
    for (int j = 0; j < ((n - 2) >> 1); j++) {
        if (marker[j]) {
            res[ind++] = (j + 1) * 2 + 1;
        }
    }
    return res;
}

int main(){
    int *prime_table = prime_sieve(1000000);
    for (int j = 0; j < 120784; j++) {
        cout << prime_table[j]<<" ";
    }
    if (prime_sieve) {
        delete[]prime_table;
    }
}

Java程序如下:

private static void makeSieve(boolean[] marker, int first, int last, int factor) {
    marker[first] = false;
    while (last - first > factor) {
        first = first + factor;
        marker[first] = false;
    }
}

public static int[] sift(int n) {
    if (n<3) {
        return new int[] {2};
    }
    boolean[] marker = new boolean[n];
    Arrays.fill(marker, true);
    int last = marker.length;
    int i = 0;
    int indexSquare = 3;
    int factor = 3;
    while (indexSquare < n) {
        if (marker[i]) {
            makeSieve(marker, indexSquare, last, factor);
        }
        ++i;
        indexSquare += factor;
        factor += 2;
        indexSquare += factor;
    }
    int count = 1;
    for (int j = 0; j < ((marker.length-2)>>1); j++) {
        if (marker[j]) {
            count++;
        }
    }
    int ind = 0;
    int[] res = new int[count];
    res[ind++] = 2;
    for (int j = 0; j < ((marker.length-2)>>1); j++) {
        if (marker[j]) {
            res[ind++] = (j+1)*2+1;
        }
    }
    return res;
}
public static void main(String[] args) {
    System.out.println(Arrays.toString(sift(100)));
}
11-18 15:03