素数的计算-python版本【图文】

目的：

统计指定范围整数内的素数集合

素数的判断：

假设a * b = N,则a 和 b不能同时大于N的算数平方根

如果N不能被≤N的开方的整数，整除，则N为素数或质数

则分母为2...int(N^0.5)

由于range函数为前闭后开，则分母范围为：range(2, int(math.sqrt(N)) + 1)

方法：

1. python内置的推导式写法（for循环）
2. python内置的filter函数-自定义函数
3. python内置的filter函数-lambda函数
4. 素数定义法
5. 埃氏筛法-定义
6. 埃氏筛法-步长优化（整数范围：去除偶数， 从3开始）
7. 埃氏筛法-步长优化（质数的倍数：去除偶数倍的质数）
8. 埃氏筛法-步长优化（合数的起始有质数的3倍优化为质数的平方
9. 埃氏筛法-numba加速

代码实现：

#!/usr/bin/python
# -*- coding:utf-8 -*-

import numpy as np
import time
import math
from numba import jit
import sys
"""
假设a * b = N,则a 和 b不能同时大于N的算数平方根
如果N不能被≤N的开方的整数，整除，则N为素数或质数
则分母为2...int(N^0.5)
由于range函数为前闭后开，则分母范围为：range(2, int(math.sqrt(N)) + 1)
"""
a = 2
b = 100000000
def timestr(cont = None):
    nowtime = time.strftime("%Y-%m-%d %H:%M:%S",time.localtime(time.time()))
    print(nowtime + " " + str(cont))

# 方法1：直接计算
def python_method(a,b):
    # 方法1：直接计算
    t = time.time()
    p = [p for p in range(a, b) if 0 not in [p % d for d in range(2, int(math.sqrt(p)) + 1)]]
    timestr("直接计算分析时间：%f" % (time.time() - t, ))
    return p

def filter_func(num):
    return 1 if 0 not in [num % d for d in range(2, int(math.sqrt(num)) + 1)] else 0
# 方法2：利用python内置函数filter
def python_filter(a,b):
    t = time.time()
    p = filter(filter_func, range(a,b))
    timestr("利用python内置函数filter：%f" % (time.time() - t,))
    return list(p)

# 方法3：利用python内置函数filter和匿名函数lambda
def python_filter_lambda(a,b):
    t = time.time()
    lamb_filter_func = lambda x : 0 not in [x % d for d in range(2, int(math.sqrt(x)) + 1)]
    p = filter(lamb_filter_func, range(a,b))
    timestr("利用python内置函数filter和匿名函数lambda：%f" % (time.time() - t,))
    return list(p)

def python_define(a,b):
    t = time.time()
    plist = []
    for num in range(a,b):
        flag = True
        for p in plist:
            if p > math.sqrt(num):
                break
            if num % p == 0:
                flag = False
                break
        if flag:
            plist.append(num)
    timestr("利用素数/质数的定义：%f" % (time.time() - t,))
    return plist

"""
    要枚举n以内的素数，可以用埃氏筛法。这是一个与辗转相除法一样古老的算法。
    首先，将2到n范围内的所有整数写下来。其中最小的数字2是素数。将表中所有2的倍数都划去。
    表中剩余的最小数字是3，它不能被更小的数整除，所以是素数。
    再将表中所有3的倍数全都划去。依次类推，如果表中剩余的最小数字是m时，m就是素数。
    然后将表中所有m的倍数全部划去。像这样反复操作，就能依次枚举n以内的素数
    埃氏筛的复杂度仅有O（nloglogn）
"""
def eratosthenes_define(a,b):
    t = time.time()
    plist = []                      # 用于存储所有的素数
    is_prime = [1] * (b + 1)        # 用于存储判断某整数num是否为素数（正好对应于索引num的位置）
    for num in range(2, b ):
        if is_prime[num]:
            if num >= a:
                plist.append(num)
            for j in range(num * 2, b, num):  # 则质数的倍数均为合数
                is_prime[j] = 0
    timestr("埃氏筛法定义：%f" % (time.time() - t,))
    return plist

def eratosthenes_opt_rangestep(a,b):
    t = time.time()
    plist = [2]                      # 用于存储所有的素数
    is_prime = [1] * (b + 1)        # 用于存储判断某整数num是否为素数（正好对应于索引num的位置）
    # for num in range(2, b ):
    for num in range(3, b, 2):        # 起点3，步长2，直接把2的倍数都略过
        if is_prime[num]:
            if num >= a:
                plist.append(num)
            for j in range(num * 2, b, num):# 则质数的倍数均为合数
                is_prime[j] = 0
    timestr("埃氏筛法-优化步长:去除范围内的偶数：%f" % (time.time() - t,))
    return plist

def eratosthenes_opt_rangestep_notp(a,b):
    t = time.time()
    plist = [2]                      # 用于存储所有的素数
    is_prime = [1] * (b + 1)        # 用于存储判断某整数num是否为素数（正好对应于索引num的位置）
    # for num in range(2, b ):
    for num in range(3, b, 2):        # 起点3，步长2，直接把2的倍数都略过
        if is_prime[num]:
            if num >= a:
                plist.append(num)
            # for j in range(num * 2, b, num):# 则质数的倍数均为合数
            for j in range(num * 3, b, num * 2):  # 则质数的倍数均为合数 2j，3j， 4j， 5j， 6j...中数字为偶数在上一层循环已经被去除，没有必要在进行遍历和标记
                is_prime[j] = 0
    timestr("埃氏筛法-优化步长:去除范围内的偶数-质数的倍数：%f" % (time.time() - t,))
    return plist

def eratosthenes_opt_rangestep_notp1(a,b):
    """
    任意合数，如果它不是一个整数的平方，那它一定有至少一个比它平方根小的素数因子。所以，当这个较小的素数因子出现时，这个合数就被标记了
    当素数 num 出现时，本应该把 num 的所有倍数都标记，但小于 num 平方的合数都已经被比 num 小的素数标记了。
    num 出现后，只要标记 num 平方及比 num 平方更大的 num 的倍数即可
    合数：num*num, num(num+1), num(num+2), ..., num+1/3/../2n-1均为偶数，不必标记，则步长设为num * 2
    """
    t = time.time()
    plist = [2]                      # 用于存储所有的素数
    is_prime = [1] * (b + 1)        # 用于存储判断某整数num是否为素数（正好对应于索引num的位置）
    # for num in range(2, b ):
    for num in range(3, b, 2):        # 起点3，步长2，直接把2的倍数都略过
        if is_prime[num]:
            if num >= a:
                plist.append(num)
            # for j in range(num * 2, b, num):# 则质数的倍数均为合数
            # for j in range(num * 3, b, num * 2):  # 则质数的倍数均为合数 2j，3j， 4j， 5j， 6j...中数字为偶数在上一层循环已经被去除，没有必要在进行遍历和标记
            for j in range(num ** 2, b, num * 2):
                is_prime[j] = 0
    timestr("埃氏筛法-优化步长:去除范围内的偶数-质数平方及其的倍数：%f" % (time.time() - t,))
    return plist

@jit(nopython=True)
def eratosthenes_opt_rangestep_notp1_numba(a,b):
    """
    任意合数，如果它不是一个整数的平方，那它一定有至少一个比它平方根小的素数因子。所以，当这个较小的素数因子出现时，这个合数就被标记了
    当素数 num 出现时，本应该把 num 的所有倍数都标记，但小于 num 平方的合数都已经被比 num 小的素数标记了。
    num 出现后，只要标记 num 平方及比 num 平方更大的 num 的倍数即可
    合数：num*num, num(num+1), num(num+2), ..., num+1/3/../2n-1均为偶数，不必标记，则步长设为num * 2
    """
    #t = time.time()
    plist = [2]                      # 用于存储所有的素数
    is_prime = [1] * (b + 1)        # 用于存储判断某整数num是否为素数（正好对应于索引num的位置）
    # for num in range(2, b ):
    for num in range(3, b, 2):        # 起点3，步长2，直接把2的倍数都略过
        if is_prime[num]:
            if num >= a:
                plist.append(num)
            # for j in range(num * 2, b, num):# 则质数的倍数均为合数
            # for j in range(num * 3, b, num * 2):  # 则质数的倍数均为合数 2j，3j， 4j， 5j， 6j...中数字为偶数在上一层循环已经被去除，没有必要在进行遍历和标记
            for j in range(num ** 2, b, num * 2):
                is_prime[j] = 0
    #timestr("埃氏筛法-优化步长:去除范围内的偶数-质数平方及其的倍数：%f" % (time.time() - t,))
    return plist

p = python_method(a,b)
print(p[:30])
p = python_filter(a,b)
print(p[:30])
p = python_filter_lambda(a,b)
print(p[:30])
p = python_define(a,b)
print(p[:30])
p = eratosthenes_define(a,b)
print(p[:30])
p = eratosthenes_opt_rangestep(a,b)
print(p[:30])
p = eratosthenes_opt_rangestep_notp(a,b)
print(p[:30])
p = eratosthenes_opt_rangestep_notp1(a,b)
print(p[:30])

t = time.time()
p = eratosthenes_opt_rangestep_notp1_numba(a,b)
timestr("埃氏筛法-优化步长:去除范围内的偶数-质数平方及其的倍数-numba：%f" % (time.time() - t,))
print(p[:30])

更多文章请关注《万象专栏》