RError.com

RError.com Logo RError.com Logo

RError.com Navigation

  • 主页

Mobile menu

Close
  • 主页
  • 系统&网络
    • 热门问题
    • 最新问题
    • 标签
  • Ubuntu
    • 热门问题
    • 最新问题
    • 标签
  • 帮助
主页 / 问题 / 962194
Accepted
MaxU - stop genocide of UA
MaxU - stop genocide of UA
Asked:2020-03-28 00:10:04 +0000 UTC2020-03-28 00:10:04 +0000 UTC 2020-03-28 00:10:04 +0000 UTC

找到所有小于 N 的素数的快速方法

  • 772

除了众所周知的之外,找到所有小于 N 的素数的最快算法是什么:

  • 埃拉托色尼筛
  • 阿特金筛
  • 筛子孙达拉玛

什么是一些快速的 Python 实现(Vanilla、Numpy 等)?

Eratosthenes 筛子的幼稚实现:

def eratosthenes2(n):
    multiples = set()
    for i in range(2, n+1):
        if i not in multiples:
            yield i
            multiples.update(range(i*i, n+1, i))

SO英文版中的类似问题

PS 创建这个问题的想法与不断涌现的新问题有关,即如何有效地实现查找所有小于 N 的素数。

python
  • 4 4 个回答
  • 10 Views

4 个回答

  • Voted
  1. Best Answer
    MaxU - stop genocide of UA
    2020-03-28T00:10:04Z2020-03-28T00:10:04Z

    作为答案,我决定比较这个问题的答案中的最佳实现并进行时间测量(为此,我使用了@jfs答案中的模块reporttime:

    素数.py

    import numpy
    from math import sqrt, ceil
    
    
    def rwh_primes(n):
        """ Returns  a list of primes < n """
        # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
        sieve = [True] * n
        for i in range(3,int(n**0.5)+1,2):
            if sieve[i]:
                sieve[i*i::2*i]=[False]*((n-i*i-1)//(2*i)+1)
        return [2] + [i for i in range(3,n,2) if sieve[i]]
    
    
    def rwh_primes1(n):
        """ Returns  a list of primes < n """
        # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
        sieve = [True] * (n//2)
        for i in range(3,int(n**0.5)+1,2):
            if sieve[i//2]:
                sieve[i*i//2::i] = [False] * ((n-i*i-1)//(2*i)+1)
        return [2] + [2*i+1 for i in range(1,n//2) if sieve[i]]
    
    
    def rwh_primes2(n):
        """ Input n>=6, Returns a list of primes, 2 <= p < n """
        # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
        n, correction = n-n%6+6, 2-(n%6>1)
        sieve = [True] * (n//3)
        for i in range(1,int(n**0.5)//3+1):
          if sieve[i]:
            k=3*i+1|1
            sieve[      k*k//3      ::2*k] = [False] * ((n//6-k*k//6-1)//k+1)
            sieve[k*(k-2*(i&1)+4)//3::2*k] = [False] * ((n//6-k*(k-2*(i&1)+4)//6-1)//k+1)
        return [2,3] + [3*i+1|1 for i in range(1,n//3-correction) if sieve[i]]
    
    
    def primesfrom3to(n):
        """ Returns a array of primes, 3 <= p < n """
        # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
        sieve = numpy.ones(n//2, dtype=numpy.bool)
        for i in range(3,int(n**0.5)+1,2):
            if sieve[i//2]:
                sieve[i*i//2::i] = False
        return 2*numpy.nonzero(sieve)[0][1::]+1
    
    
    def primesfrom2to(n):
        """ Input n>=6, Returns a array of primes, 2 <= p < n """
        # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
        sieve = numpy.ones(n//3 + (n%6==2), dtype=numpy.bool)
        for i in range(1,int(n**0.5)//3+1):
            if sieve[i]:
                k=3*i+1|1
                sieve[       k*k//3     ::2*k] = False
                sieve[k*(k-2*(i&1)+4)//3::2*k] = False
        return numpy.r_[2,3,((3*numpy.nonzero(sieve)[0][1:]+1)|1)]
    
    
    def sieveOfEratosthenes(n):
        """sieveOfEratosthenes(n): return the list of the primes < n."""
        # Code from: <dickinsm@gmail.com>, Nov 30 2006
        # http://groups.google.com/group/comp.lang.python/msg/f1f10ced88c68c2d
        if n <= 2:
            return []
        sieve = list(range(3, n, 2))
        top = len(sieve)
        for si in sieve:
            if si:
                bottom = (si*si - 3) // 2
                if bottom >= top:
                    break
                sieve[bottom::si] = [0] * -((bottom - top) // si)
        return [2] + [el for el in sieve if el]
    
    
    def sieveOfAtkin(end):
        """sieveOfAtkin(end): return a list of all the prime numbers <end
        using the Sieve of Atkin."""
        # Code by Steve Krenzel, <Sgk284@gmail.com>, improved
        # Code: https://web.archive.org/web/20080324064651/http://krenzel.info/?p=83
        # Info: http://en.wikipedia.org/wiki/Sieve_of_Atkin
        assert end > 0
        lng = ((end-1) // 2)
        sieve = [False] * (lng + 1)
    
        x_max, x2, xd = int(sqrt((end-1)/4.0)), 0, 4
        for xd in range(4, 8*x_max + 2, 8):
            x2 += xd
            y_max = int(sqrt(end-x2))
            n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
            if not (n & 1):
                n -= n_diff
                n_diff -= 2
            for d in range((n_diff - 1) << 1, -1, -8):
                m = n % 12
                if m == 1 or m == 5:
                    m = n >> 1
                    sieve[m] = not sieve[m]
                n -= d
    
        x_max, x2, xd = int(sqrt((end-1) / 3.0)), 0, 3
        for xd in range(3, 6 * x_max + 2, 6):
            x2 += xd
            y_max = int(sqrt(end-x2))
            n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
            if not(n & 1):
                n -= n_diff
                n_diff -= 2
            for d in range((n_diff - 1) << 1, -1, -8):
                if n % 12 == 7:
                    m = n >> 1
                    sieve[m] = not sieve[m]
                n -= d
    
        x_max, y_min, x2, xd = int((2 + sqrt(4-8*(1-end)))/4), -1, 0, 3
        for x in range(1, x_max + 1):
            x2 += xd
            xd += 6
            if x2 >= end: y_min = (((int(ceil(sqrt(x2 - end))) - 1) << 1) - 2) << 1
            n, n_diff = ((x*x + x) << 1) - 1, (((x-1) << 1) - 2) << 1
            for d in range(n_diff, y_min, -8):
                if n % 12 == 11:
                    m = n >> 1
                    sieve[m] = not sieve[m]
                n += d
    
        primes = [2, 3]
        if end <= 3:
            return primes[:max(0,end-2)]
    
        for n in range(5 >> 1, (int(sqrt(end))+1) >> 1):
            if sieve[n]:
                primes.append((n << 1) + 1)
                aux = (n << 1) + 1
                aux *= aux
                for k in range(aux, end, 2 * aux):
                    sieve[k >> 1] = False
    
        s  = int(sqrt(end)) + 1
        if s  % 2 == 0:
            s += 1
        primes.extend([i for i in range(s, end, 2) if sieve[i >> 1]])
    
        return primes
    
    
    def sundaram3(max_n):
        # https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/2073279#2073279
        numbers = list(range(3, max_n+1, 2))
        half = (max_n)//2
        initial = 4
    
        for step in range(3, max_n+1, 2):
            for i in range(initial, half, step):
                numbers[i-1] = 0
            initial += 2*(step+1)
    
            if initial > half:
                return [2] + list(filter(None, numbers))

    时间测量:

    import os
    import pandas as pd
    
    os.chdir(r'C:\work\project\primes_timing')
    
    import primes
    from reporttime import get_functions_with_prefix, measure
    
    funcs = get_functions_with_prefix('', module=primes)
        
    d = {}
    for n in [10**3, 10**4, 10**5, 10**6]:
        m = measure(funcs, args=[n])
        d[n] = pd.Series({b:a for a,b in m}) 
    r = pd.DataFrame(d)
    

    打印输出:

    name                     time ratio comment
    primesfrom3to         21 usec  1.00 [1000]
    sieveOfEratosthenes 40.1 usec  1.91 [1000]
    rwh_primes          50.2 usec  2.39 [1000]
    primesfrom2to       63.6 usec  3.02 [1000]
    sundaram3           91.8 usec  4.37 [1000]
    rwh_primes1          110 usec  5.25 [1000]
    rwh_primes2          129 usec  6.12 [1000]
    sieveOfAtkin         286 usec 13.61 [1000]
    
    name                     time ratio comment
    primesfrom3to       69.3 usec  1.00 [10000]
    primesfrom2to        172 usec  2.48 [10000]
    sieveOfEratosthenes  375 usec  5.41 [10000]
    rwh_primes           447 usec  6.45 [10000]
    rwh_primes2          503 usec  7.26 [10000]
    rwh_primes1          577 usec  8.33 [10000]
    sundaram3           1.22 msec 17.66 [10000]
    sieveOfAtkin         2.2 msec 31.72 [10000]
    
    name                     time ratio comment
    primesfrom3to        306 usec  1.00 [100000]
    primesfrom2to        371 usec  1.21 [100000]
    sieveOfEratosthenes 3.82 msec 12.47 [100000]
    rwh_primes          3.94 msec 12.84 [100000]
    rwh_primes2         4.19 msec 13.68 [100000]
    rwh_primes1         5.37 msec 17.53 [100000]
    sundaram3           15.2 msec 49.57 [100000]
    sieveOfAtkin        21.3 msec 69.47 [100000]
    
    name                     time  ratio comment
    primesfrom2to       2.29 msec   1.00 [1000000]
    primesfrom3to       3.05 msec   1.33 [1000000]
    rwh_primes2         38.7 msec  16.91 [1000000]
    rwh_primes1         48.3 msec  21.12 [1000000]
    rwh_primes          52.5 msec  22.96 [1000000]
    sieveOfEratosthenes   61 msec  26.64 [1000000]
    sieveOfAtkin         205 msec  89.64 [1000000]
    sundaram3            245 msec 107.27 [1000000]
    

    带有结果的数据框:

    In [67]: r
    Out[67]:
                          1000      10000     100000    1000000
    primesfrom2to        0.000064  0.000172  0.000371  0.002288
    primesfrom3to        0.000021  0.000069  0.000306  0.003046
    rwh_primes           0.000050  0.000447  0.003936  0.052516
    rwh_primes1          0.000110  0.000577  0.005373  0.048307
    rwh_primes2          0.000129  0.000503  0.004191  0.038678
    sieveOfAtkin         0.000286  0.002197  0.021293  0.205074
    sieveOfEratosthenes  0.000040  0.000375  0.003823  0.060955
    sundaram3            0.000092  0.001223  0.015194  0.245407
    

    日程:

    import matplotlib.pyplot as plt
    import matplotlib
    
    matplotlib.style.use('ggplot')
    
    r.T.plot(loglog=True, grid=True, figsize=(8, 6))
    plt.xlabel('N - generate all primes below N')
    plt.ylabel('Time in seconds')
    plt.show()
    

    在此处输入图像描述

    • 10
  2. MaxU - stop genocide of UA
    2020-03-29T22:08:23Z2020-03-29T22:08:23Z

    对我所知的Vanilla Python中 Eratosthenes 筛的最快实现工作的分析,即 无需使用额外的模块。

    (c) Bruno Astrolino E Silva - 我刚刚添加了调试信息:

    def primes(n, verbose=0):
        """ Returns  a list of primes < n """
        # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
        def pr(*args):
            if verbose > 0:
                print(*args)
        sieve = [True] * n
        pr("все чётные числа игнорируются и будут пропущены при возврате...\n")
        for i in range(3,int(n**0.5)+1,2):
            if sieve[i]:
                pr('содержимое решета:\t{}'.format([x for x in range(3,n,2) if sieve[x]]))
                pr(f'i:{i} вычёркиваем все числа кратные "{i}", начиная с "{i}^2": {list(range(i*i, n, 2*i))}')
                sieve[i*i::2*i]=[False]*((n-i*i-1)//(2*i)+1)
                pr(f'sieve[{i}*{i}::2*{i}]=[False]*(({n-i}*{i-1})//(2*{i})+1)')
                pr('содержимое решета:\t{}'.format([x for x in range(3,n,2) if sieve[x]]))
                pr('*' * 60)
        return [2] + [i for i in range(3,n,2) if sieve[i]]
    

    50PS这个实现非常有效地使用切片,并且由于这一点,最小化了循环迭代的次数 - 在下面的例子中,只需要三个循环迭代就可以找到所有小于 的素数。

    输出n=50:

    In [165]: primes(50, verbose=1)
    все чётные числа игнорируются и будут пропущены при возврате...
    
    содержимое решета:      [3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 43, 45, 47, 49]
    i:3 вычёркиваем все числа кратные "3", начиная с "3^2": [9, 15, 21, 27, 33, 39, 45]
    sieve[3*3::2*3]=[False]*((47*2)//(2*3)+1)
    содержимое решета:      [3, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, 37, 41, 43, 47, 49]
    ************************************************************
    содержимое решета:      [3, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, 37, 41, 43, 47, 49]
    i:5 вычёркиваем все числа кратные "5", начиная с "5^2": [25, 35, 45]
    sieve[5*5::2*5]=[False]*((45*4)//(2*5)+1)
    содержимое решета:      [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49]
    ************************************************************
    содержимое решета:      [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49]
    i:7 вычёркиваем все числа кратные "7", начиная с "7^2": [49]
    sieve[7*7::2*7]=[False]*((43*6)//(2*7)+1)
    содержимое решета:      [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
    ************************************************************
    Out[165]: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47
    
    • 3
  3. Vikk317
    2020-04-03T20:09:09Z2020-04-03T20:09:09Z

    实际上,例如,直到 2000000 的所有素数之和的输出不到 0.01 秒,我编写了这段代码,实现了特定任务的 Eratosthenes 筛,但是在给定范围内查找素数是一个相当普遍的过程:

    #include <stdio.h>
    #include <malloc.h>
    #include <math.h>
    
    main (){
    
    int * pointer;
    int i, ii;
    int size=2000000;
    long long int sum=2;
    pointer = (int*) malloc (size*sizeof(int));
    for (i=0;i<=size;i++)
        *(pointer+i)=i;
    *(pointer+1)=0;
    for (i=2;i<=size;i=i+2)
        *(pointer+i)=0;
    for (i=3;i*i<size;i++){
        if (*(pointer+i)!=0) {
            for (ii=i*i;ii<size;ii=ii+i){
            *(pointer+ii)=0;
            } 
        }
    }
    for (i=0;i<size;i++){
        sum=sum+*(pointer+i);
    }
    printf ("Sum for %d elem is %lld", size, sum);
    free (pointer);
    return 0;
    }
    
    • 2
  4. Evgeny
    2020-02-16T21:53:01Z2020-02-16T21:53:01Z

    如果你需要一个list素数列表进行进一步处理,那么你可以这样做,计算速度令人满意。该算法基于 2 个原则: 1. 如果一个数不能被先前的素数整除,则它是素数。2.一个合数在k / 2内有一个除数。这是我的解决方案:

    my_list = [2]
    for k in range(3, 100000, 2):
        i = 0
        while k > my_list[i]:
            if k % my_list[i] == 0:
                break
            if my_list[i] * 2 > k:
               my_list.append(k)
               break
            i +=1
    print(my_list)
    
    • 1

相关问题

Sidebar

Stats

  • 问题 10021
  • Answers 30001
  • 最佳答案 8000
  • 用户 6900
  • 常问
  • 回答
  • Marko Smith

    根据浏览器窗口的大小调整背景图案的大小

    • 2 个回答
  • Marko Smith

    理解for循环的执行逻辑

    • 1 个回答
  • Marko Smith

    复制动态数组时出错(C++)

    • 1 个回答
  • Marko Smith

    Or and If,elif,else 构造[重复]

    • 1 个回答
  • Marko Smith

    如何构建支持 x64 的 APK

    • 1 个回答
  • Marko Smith

    如何使按钮的输入宽度?

    • 2 个回答
  • Marko Smith

    如何显示对象变量的名称?

    • 3 个回答
  • Marko Smith

    如何循环一个函数?

    • 1 个回答
  • Marko Smith

    LOWORD 宏有什么作用?

    • 2 个回答
  • Marko Smith

    从字符串的开头删除直到并包括一个字符

    • 2 个回答
  • Martin Hope
    Alexandr_TT 2020年新年大赛! 2020-12-20 18:20:21 +0000 UTC
  • Martin Hope
    Alexandr_TT 圣诞树动画 2020-12-23 00:38:08 +0000 UTC
  • Martin Hope
    Air 究竟是什么标识了网站访问者? 2020-11-03 15:49:20 +0000 UTC
  • Martin Hope
    Qwertiy 号码显示 9223372036854775807 2020-07-11 18:16:49 +0000 UTC
  • Martin Hope
    user216109 如何为黑客设下陷阱,或充分击退攻击? 2020-05-10 02:22:52 +0000 UTC
  • Martin Hope
    Qwertiy 并变成3个无穷大 2020-11-06 07:15:57 +0000 UTC
  • Martin Hope
    koks_rs 什么是样板代码? 2020-10-27 15:43:19 +0000 UTC
  • Martin Hope
    Sirop4ik 向 git 提交发布的正确方法是什么? 2020-10-05 00:02:00 +0000 UTC
  • Martin Hope
    faoxis 为什么在这么多示例中函数都称为 foo? 2020-08-15 04:42:49 +0000 UTC
  • Martin Hope
    Pavel Mayorov 如何从事件或回调函数中返回值?或者至少等他们完成。 2020-08-11 16:49:28 +0000 UTC

热门标签

javascript python java php c# c++ html android jquery mysql

Explore

  • 主页
  • 问题
    • 热门问题
    • 最新问题
  • 标签
  • 帮助

Footer

RError.com

关于我们

  • 关于我们
  • 联系我们

Legal Stuff

  • Privacy Policy

帮助

© 2023 RError.com All Rights Reserve   沪ICP备12040472号-5