Moyen le plus rapide de répertorier tous les nombres premiers sous N

358

C'est le meilleur algorithme que j'ai pu trouver.

def get_primes(n):
    numbers = set(range(n, 1, -1))
    primes = []
    while numbers:
        p = numbers.pop()
        primes.append(p)
        numbers.difference_update(set(range(p*2, n+1, p)))
    return primes

>>> timeit.Timer(stmt='get_primes.get_primes(1000000)', setup='import   get_primes').timeit(1)
1.1499958793645562

Peut-il être rendu encore plus rapide?

Ce code a un défaut: puisqu'il numberss'agit d'un ensemble non ordonné, il n'y a aucune garantie qui numbers.pop()supprimera le numéro le plus bas de l'ensemble. Néanmoins, cela fonctionne (au moins pour moi) pour certains numéros d'entrée:

>>> sum(get_primes(2000000))
142913828922L
#That's the correct sum of all numbers below 2 million
>>> 529 in get_primes(1000)
False
>>> 529 in get_primes(530)
True
jbochi
la source
Le code sniplet en question est beaucoup plus rapide si les nombres déclarés comme les nombres = set (range (n, 2, -2)). Mais ne peut pas battre sundaram3. Merci pour la question.
Shekhar
3
Ce serait bien s'il pouvait y avoir des versions Python 3 des fonctions dans les réponses.
Michael Foukarakis du
Il y a sûrement une bibliothèque pour le faire, nous n'avons donc pas à rouler notre propre> xkcd promis que Python est aussi simple que import antigravity. N'y a-t-il rien comme require 'prime'; Prime.take(10)(Ruby)?
Colonel Panic
2
@ColonelPanic En l'occurrence, j'ai mis à jour github.com/jaredks/pyprimesieve pour Py3 et ajouté à PyPi. C'est certainement plus rapide que ces derniers mais pas des ordres de grandeur - plus comme ~ 5 fois plus rapide que les meilleures versions numpy.
Jared
3
@ColonelPanic: Je pense que la modification des anciennes réponses pour noter qu'elles ont vieilli est appropriée, car cela en fait une ressource plus utile. Si la réponse «acceptée» n'est plus la meilleure, modifiez peut-être une note dans la question avec une mise à jour de 2015 pour orienter les gens vers la meilleure méthode actuelle.
Peter Cordes

Réponses:

366

Attention: timeit résultats peuvent varier en raison de différences de matériel ou de version de Python.

Voici un script qui compare un certain nombre d'implémentations:

Un grand merci à stephan d' avoir porté sieve_wheel_30 à mon attention. Le mérite revient à Robert William Hanks pour primesfrom2to, primesfrom3to, rwh_primes, rwh_primes1 et rwh_primes2.

Parmi les méthodes Python simples testées, avec psyco , pour n = 1000000, rwh_primes1 a été la plus rapide testée.

+---------------------+-------+
| Method              | ms    |
+---------------------+-------+
| rwh_primes1         | 43.0  |
| sieveOfAtkin        | 46.4  |
| rwh_primes          | 57.4  |
| sieve_wheel_30      | 63.0  |
| rwh_primes2         | 67.8  |    
| sieveOfEratosthenes | 147.0 |
| ambi_sieve_plain    | 152.0 |
| sundaram3           | 194.0 |
+---------------------+-------+

Parmi les méthodes Python simples testées, sans psyco , pour n = 1000000, rwh_primes2 était la plus rapide.

+---------------------+-------+
| Method              | ms    |
+---------------------+-------+
| rwh_primes2         | 68.1  |
| rwh_primes1         | 93.7  |
| rwh_primes          | 94.6  |
| sieve_wheel_30      | 97.4  |
| sieveOfEratosthenes | 178.0 |
| ambi_sieve_plain    | 286.0 |
| sieveOfAtkin        | 314.0 |
| sundaram3           | 416.0 |
+---------------------+-------+

De toutes les méthodes testées, permettant numpy , pour n = 1000000, primesfrom2to a été la plus rapide testée.

+---------------------+-------+
| Method              | ms    |
+---------------------+-------+
| primesfrom2to       | 15.9  |
| primesfrom3to       | 18.4  |
| ambi_sieve          | 29.3  |
+---------------------+-------+

Les délais ont été mesurés à l'aide de la commande:

python -mtimeit -s"import primes" "primes.{method}(1000000)"

avec {method}remplacé par chacun des noms de méthode.

primes.py:

#!/usr/bin/env python
import psyco; psyco.full()
from math import sqrt, ceil
import numpy as np

def rwh_primes(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    sieve = [True] * n
    for i in xrange(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 xrange(3,n,2) if sieve[i]]

def rwh_primes1(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    sieve = [True] * (n/2)
    for i in xrange(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 xrange(1,n/2) if sieve[i]]

def rwh_primes2(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns a list of primes, 2 <= p < n """
    correction = (n%6>1)
    n = {0:n,1:n-1,2:n+4,3:n+3,4:n+2,5:n+1}[n%6]
    sieve = [True] * (n/3)
    sieve[0] = False
    for i in xrange(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+4*k-2*k*(i&1))/3::2*k]=[False]*((n/6-(k*k+4*k-2*k*(i&1))/6-1)/k+1)
    return [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) if sieve[i]]

def sieve_wheel_30(N):
    # http://zerovolt.com/?p=88
    ''' Returns a list of primes <= N using wheel criterion 2*3*5 = 30

Copyright 2009 by zerovolt.com
This code is free for non-commercial purposes, in which case you can just leave this comment as a credit for my work.
If you need this code for commercial purposes, please contact me by sending an email to: info [at] zerovolt [dot] com.'''
    __smallp = ( 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, 101, 103, 107, 109, 113, 127, 131, 137, 139,
    149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
    229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311,
    313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401,
    409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491,
    499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
    601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683,
    691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
    809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887,
    907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997)

    wheel = (2, 3, 5)
    const = 30
    if N < 2:
        return []
    if N <= const:
        pos = 0
        while __smallp[pos] <= N:
            pos += 1
        return list(__smallp[:pos])
    # make the offsets list
    offsets = (7, 11, 13, 17, 19, 23, 29, 1)
    # prepare the list
    p = [2, 3, 5]
    dim = 2 + N // const
    tk1  = [True] * dim
    tk7  = [True] * dim
    tk11 = [True] * dim
    tk13 = [True] * dim
    tk17 = [True] * dim
    tk19 = [True] * dim
    tk23 = [True] * dim
    tk29 = [True] * dim
    tk1[0] = False
    # help dictionary d
    # d[a , b] = c  ==> if I want to find the smallest useful multiple of (30*pos)+a
    # on tkc, then I need the index given by the product of [(30*pos)+a][(30*pos)+b]
    # in general. If b < a, I need [(30*pos)+a][(30*(pos+1))+b]
    d = {}
    for x in offsets:
        for y in offsets:
            res = (x*y) % const
            if res in offsets:
                d[(x, res)] = y
    # another help dictionary: gives tkx calling tmptk[x]
    tmptk = {1:tk1, 7:tk7, 11:tk11, 13:tk13, 17:tk17, 19:tk19, 23:tk23, 29:tk29}
    pos, prime, lastadded, stop = 0, 0, 0, int(ceil(sqrt(N)))
    # inner functions definition
    def del_mult(tk, start, step):
        for k in xrange(start, len(tk), step):
            tk[k] = False
    # end of inner functions definition
    cpos = const * pos
    while prime < stop:
        # 30k + 7
        if tk7[pos]:
            prime = cpos + 7
            p.append(prime)
            lastadded = 7
            for off in offsets:
                tmp = d[(7, off)]
                start = (pos + prime) if off == 7 else (prime * (const * (pos + 1 if tmp < 7 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 11
        if tk11[pos]:
            prime = cpos + 11
            p.append(prime)
            lastadded = 11
            for off in offsets:
                tmp = d[(11, off)]
                start = (pos + prime) if off == 11 else (prime * (const * (pos + 1 if tmp < 11 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 13
        if tk13[pos]:
            prime = cpos + 13
            p.append(prime)
            lastadded = 13
            for off in offsets:
                tmp = d[(13, off)]
                start = (pos + prime) if off == 13 else (prime * (const * (pos + 1 if tmp < 13 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 17
        if tk17[pos]:
            prime = cpos + 17
            p.append(prime)
            lastadded = 17
            for off in offsets:
                tmp = d[(17, off)]
                start = (pos + prime) if off == 17 else (prime * (const * (pos + 1 if tmp < 17 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 19
        if tk19[pos]:
            prime = cpos + 19
            p.append(prime)
            lastadded = 19
            for off in offsets:
                tmp = d[(19, off)]
                start = (pos + prime) if off == 19 else (prime * (const * (pos + 1 if tmp < 19 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 23
        if tk23[pos]:
            prime = cpos + 23
            p.append(prime)
            lastadded = 23
            for off in offsets:
                tmp = d[(23, off)]
                start = (pos + prime) if off == 23 else (prime * (const * (pos + 1 if tmp < 23 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 29
        if tk29[pos]:
            prime = cpos + 29
            p.append(prime)
            lastadded = 29
            for off in offsets:
                tmp = d[(29, off)]
                start = (pos + prime) if off == 29 else (prime * (const * (pos + 1 if tmp < 29 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # now we go back to top tk1, so we need to increase pos by 1
        pos += 1
        cpos = const * pos
        # 30k + 1
        if tk1[pos]:
            prime = cpos + 1
            p.append(prime)
            lastadded = 1
            for off in offsets:
                tmp = d[(1, off)]
                start = (pos + prime) if off == 1 else (prime * (const * pos + tmp) )//const
                del_mult(tmptk[off], start, prime)
    # time to add remaining primes
    # if lastadded == 1, remove last element and start adding them from tk1
    # this way we don't need an "if" within the last while
    if lastadded == 1:
        p.pop()
    # now complete for every other possible prime
    while pos < len(tk1):
        cpos = const * pos
        if tk1[pos]: p.append(cpos + 1)
        if tk7[pos]: p.append(cpos + 7)
        if tk11[pos]: p.append(cpos + 11)
        if tk13[pos]: p.append(cpos + 13)
        if tk17[pos]: p.append(cpos + 17)
        if tk19[pos]: p.append(cpos + 19)
        if tk23[pos]: p.append(cpos + 23)
        if tk29[pos]: p.append(cpos + 29)
        pos += 1
    # remove exceeding if present
    pos = len(p) - 1
    while p[pos] > N:
        pos -= 1
    if pos < len(p) - 1:
        del p[pos+1:]
    # return p list
    return p

def sieveOfEratosthenes(n):
    """sieveOfEratosthenes(n): return the list of the primes < n."""
    # Code from: <[email protected]>, Nov 30 2006
    # http://groups.google.com/group/comp.lang.python/msg/f1f10ced88c68c2d
    if n <= 2:
        return []
    sieve = 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, <[email protected]>, 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 xrange(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 xrange((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 xrange(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 xrange((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 xrange(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 xrange(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 xrange(5 >> 1, (int(sqrt(end))+1) >> 1):
        if sieve[n]:
            primes.append((n << 1) + 1)
            aux = (n << 1) + 1
            aux *= aux
            for k in xrange(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 xrange(s, end, 2) if sieve[i >> 1]])

    return primes

def ambi_sieve_plain(n):
    s = range(3, n, 2)
    for m in xrange(3, int(n**0.5)+1, 2): 
        if s[(m-3)/2]: 
            for t in xrange((m*m-3)/2,(n>>1)-1,m):
                s[t]=0
    return [2]+[t for t in s if t>0]

def sundaram3(max_n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/2073279#2073279
    numbers = range(3, max_n+1, 2)
    half = (max_n)//2
    initial = 4

    for step in xrange(3, max_n+1, 2):
        for i in xrange(initial, half, step):
            numbers[i-1] = 0
        initial += 2*(step+1)

        if initial > half:
            return [2] + filter(None, numbers)

################################################################################
# Using Numpy:
def ambi_sieve(n):
    # http://tommih.blogspot.com/2009/04/fast-prime-number-generator.html
    s = np.arange(3, n, 2)
    for m in xrange(3, int(n ** 0.5)+1, 2): 
        if s[(m-3)/2]: 
            s[(m*m-3)/2::m]=0
    return np.r_[2, s[s>0]]

def primesfrom3to(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns a array of primes, p < n """
    assert n>=2
    sieve = np.ones(n/2, dtype=np.bool)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i/2]:
            sieve[i*i/2::i] = False
    return np.r_[2, 2*np.nonzero(sieve)[0][1::]+1]    

def primesfrom2to(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
    sieve[0] = False
    for i in xrange(int(n**0.5)/3+1):
        if sieve[i]:
            k=3*i+1|1
            sieve[      ((k*k)/3)      ::2*k] = False
            sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
    return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]

if __name__=='__main__':
    import itertools
    import sys

    def test(f1,f2,num):
        print('Testing {f1} and {f2} return same results'.format(
            f1=f1.func_name,
            f2=f2.func_name))
        if not all([a==b for a,b in itertools.izip_longest(f1(num),f2(num))]):
            sys.exit("Error: %s(%s) != %s(%s)"%(f1.func_name,num,f2.func_name,num))

    n=1000000
    test(sieveOfAtkin,sieveOfEratosthenes,n)
    test(sieveOfAtkin,ambi_sieve,n)
    test(sieveOfAtkin,ambi_sieve_plain,n) 
    test(sieveOfAtkin,sundaram3,n)
    test(sieveOfAtkin,sieve_wheel_30,n)
    test(sieveOfAtkin,primesfrom3to,n)
    test(sieveOfAtkin,primesfrom2to,n)
    test(sieveOfAtkin,rwh_primes,n)
    test(sieveOfAtkin,rwh_primes1,n)         
    test(sieveOfAtkin,rwh_primes2,n)

L'exécution du script teste que toutes les implémentations donnent le même résultat.

unutbu
la source
4
Si vous êtes intéressé par du code Python non pur, alors vous devriez vérifier gmpy- il a un assez bon support pour les nombres premiers, via la next_primeméthode de son mpztype.
Alex Martelli
1
Si vous utilisez pypy, ces benchmarks (les psyco) semblent assez éteints. Étonnamment, j'ai trouvé que sieveOfEratosthenes et ambi_sieve_plain étaient les plus rapides avec Pypy. C'est ce que j'ai trouvé pour les non-engourdis gist.github.com/5bf466bb1ee9e5726a52
Ehsan Kia
1
Si quelqu'un se demande comment les fonctions ici se comparent à PG7.8 de Wikibooks pour python pur sans psyco ni pypy: pour n = 1000000: PG7.8: 4,93 s par boucle; rwh_primes1: 69 ms par boucle; rwh_primes2: 57,1 ms par boucle
génial
8
Pouvez-vous mettre à jour cela avec PyPy, maintenant que psyco est mort et que PyPy l'a remplacé?
noɥʇʎԀʎzɐɹƆ
3
Ce serait génial si ces fonctions et temporisations pouvaient être mises à jour pour python3.
cs95
135

Code Python pur plus rapide et plus mémoire:

def primes(n):
    """ Returns  a list of primes < n """
    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]]

ou en commençant par un demi-tamis

def primes1(n):
    """ Returns  a list of primes < n """
    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]]

Code numpy plus rapide et plus efficace en mémoire:

import numpy
def primesfrom3to(n):
    """ Returns a array of primes, 3 <= p < n """
    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

une variation plus rapide à partir d'un tiers de tamis:

import numpy
def primesfrom2to(n):
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    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)]

Une version (difficile à coder) en python pur du code ci-dessus serait:

def primes2(n):
    """ Input n>=6, Returns a list of primes, 2 <= p < n """
    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]]

Malheureusement, pure-python n'adopte pas la méthode numpy plus simple et plus rapide pour effectuer l'affectation, et appeler len()à l'intérieur de la boucle comme dans [False]*len(sieve[((k*k)//3)::2*k])est trop lent. J'ai donc dû improviser pour corriger les entrées (et éviter plus de mathématiques) et faire de la magie mathématique extrême (et douloureuse).

Personnellement, je pense que c'est dommage que numpy (qui est si largement utilisé) ne fasse pas partie de la bibliothèque standard Python, et que les améliorations de la syntaxe et de la vitesse semblent être complètement ignorées par les développeurs Python.

Robert William Hanks
la source
2
Numpy est désormais compatible avec Python 3. Le fait qu'il ne soit pas dans la bibliothèque standard est bon, de cette façon, ils peuvent avoir leur propre cycle de publication.
Adam
pour simplement stocker les valeurs binaires dans un tableau, je suggère bitarray- tel qu'utilisé ici (pour le tamis principal le plus simple; pas un concurrent dans la course ici!) stackoverflow.com/questions/31120986/…
hiro protagonist
Lors de la coulée dans la primesfrom2to()méthode, la division doit-elle être à l'intérieur des supports?
355durch113
3
Pour une version pure python compatible avec python 3, suivez ce lien: stackoverflow.com/a/33356284/2482582
Moebius
1
Holy buttsnacks cette ventouse est rapide.
Scott
42

Il y a un échantillon assez soigné du livre de recettes Python ici - la version la plus rapide proposée sur cette URL est:

import itertools
def erat2( ):
    D = {  }
    yield 2
    for q in itertools.islice(itertools.count(3), 0, None, 2):
        p = D.pop(q, None)
        if p is None:
            D[q*q] = q
            yield q
        else:
            x = p + q
            while x in D or not (x&1):
                x += p
            D[x] = p

ce qui donnerait

def get_primes_erat(n):
  return list(itertools.takewhile(lambda p: p<n, erat2()))

Mesurant à l'invite du shell (comme je préfère le faire) avec ce code dans pri.py, j'observe:

$ python2.5 -mtimeit -s'import pri' 'pri.get_primes(1000000)'
10 loops, best of 3: 1.69 sec per loop
$ python2.5 -mtimeit -s'import pri' 'pri.get_primes_erat(1000000)'
10 loops, best of 3: 673 msec per loop

il semble donc que la solution Cookbook soit deux fois plus rapide.

Alex Martelli
la source
1
@jbochi, vous êtes les bienvenus - mais regardez cette URL, y compris les crédits: il nous a fallu dix d'entre nous pour affiner collectivement le code jusqu'à ce point, y compris des luminaires de performance Python tels que Tim Peters et Raymond Hettinger (j'ai écrit le texte final de la recette depuis que j'ai édité le livre de cuisine imprimé, mais en termes de codage, ma contribution était comparable à celle des autres) - au final, c'est un code vraiment subtil et finement réglé, et ce n'est pas surprenant! -)
Alex Martelli
@Alex: Sachant que votre code est "seulement" deux fois plus rapide que le mien, je suis alors assez fier. :) L'URL était également très intéressante à lire. Merci encore.
jbochi
Et cela peut être rendu encore plus rapide avec un changement mineur: voir stackoverflow.com/questions/2211990/…
tzot
1
... Et cela peut être rendu encore plus rapide avec une accélération supplémentaire de ~ 1,2x-1,3x, une réduction drastique de l'empreinte mémoire de O (n) à O (sqrt (n)) et une amélioration de la complexité temporelle empirique, en différant l'ajout de amorce le dict jusqu'à ce que leur carré soit vu dans l'entrée. Testez-le ici .
Will Ness
28

En utilisant Sieve de Sundaram , je pense avoir battu le record de Python pur:

def sundaram3(max_n):
    numbers = range(3, max_n+1, 2)
    half = (max_n)//2
    initial = 4

    for step in xrange(3, max_n+1, 2):
        for i in xrange(initial, half, step):
            numbers[i-1] = 0
        initial += 2*(step+1)

        if initial > half:
            return [2] + filter(None, numbers)

Comparaison:

C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.get_primes_erat(1000000)"
10 loops, best of 3: 710 msec per loop

C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.daniel_sieve_2(1000000)"
10 loops, best of 3: 435 msec per loop

C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.sundaram3(1000000)"
10 loops, best of 3: 327 msec per loop
jbochi
la source
1
J'ai réussi à accélérer votre fonction d'environ 20% en ajoutant "zéro = 0" en haut de la fonction, puis en remplaçant le lambda dans votre filtre par "zéro .__ sous__". Pas le plus joli code du monde, mais un peu plus rapide :)
truppo
1
@truppo: Merci pour votre commentaire! Je viens de réaliser que passer à la Noneplace de la fonction d'origine fonctionne et c'est encore plus rapide quezero.__sub__
jbochi
7
Saviez-vous que si vous réussissez, sundaram3(9)il reviendra [2, 3, 5, 7, 9]? Il semble faire cela avec de nombreux - peut - être tous - nombres impairs (même quand ils ne sont pas premiers)
wrhall
1
il a un problème: sundaram3 (7071) inclut 7071 alors qu'il n'est pas le premier
bigOther
18

L'algorithme est rapide, mais il a un sérieux défaut:

>>> sorted(get_primes(530))
[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, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163,
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 527, 529]
>>> 17*31
527
>>> 23*23
529

Vous supposez que numbers.pop()cela retournerait le plus petit nombre de l'ensemble, mais ce n'est pas du tout garanti. Les ensembles ne sont pas ordonnés et pop()supprime et renvoie un élément arbitraire , il ne peut donc pas être utilisé pour sélectionner le premier nombre parmi les nombres restants.

qch
la source
17

Pour une solution vraiment plus rapide avec un N suffisamment grand serait de télécharger une liste pré-calculée de nombres premiers , de la stocker en tant que tuple et de faire quelque chose comme:

for pos,i in enumerate(primes):
    if i > N:
        print primes[:pos]

Si N > primes[-1] seulement alors calculez plus de nombres premiers et enregistrez la nouvelle liste dans votre code, la prochaine fois, c'est aussi rapide.

Sortez toujours des sentiers battus.

Kimvais
la source
9
Pour être juste, cependant, vous devez compter le temps de téléchargement, de décompression et de formatage des nombres premiers et les comparer avec le temps de générer des nombres premiers à l'aide d'un algorithme - n'importe lequel de ces algorithmes pourrait facilement écrire les résultats dans un fichier pour plus tard utilisation. Je pense que dans ce cas, avec suffisamment de mémoire pour calculer réellement tous les nombres premiers inférieurs à 982 451 653, la solution numpy serait encore plus rapide.
Daniel G
3
@Daniel correct. Cependant, le magasin ce que vous avez et continuez chaque fois que nécessaire est toujours debout ...
Kimvais
@Daniel GI pense que le temps de téléchargement n'est pas pertinent. Il ne s'agit pas vraiment de générer des chiffres, vous devriez donc prendre en compte l'algorithme utilisé pour créer la liste que vous téléchargez. Et toute complexité temporelle ignorerait une fois le transfert de fichiers lui étant donné O (n).
Ross
La FAQ de la page Prime UTM suggère que le calcul de petits nombres premiers est plus rapide que de les lire sur un disque (la question est de savoir ce que signifie petits).
Batman
12

Si vous ne voulez pas réinventer la roue, vous pouvez installer la bibliothèque symbolique de mathématiques sympy (oui c'est compatible Python 3)

pip install sympy

Et utilisez la fonction primerange

from sympy import sieve
primes = list(sieve.primerange(1, 10**6))
Colonel Panic
la source
9

Si vous acceptez itertools mais pas numpy, voici une adaptation de rwh_primes2 pour Python 3 qui tourne environ deux fois plus vite sur ma machine. Le seul changement substantiel consiste à utiliser un bytearray au lieu d'une liste pour le booléen, et à utiliser compress au lieu d'une compréhension de liste pour construire la liste finale. (J'ajouterais ceci comme un commentaire comme moarningsun si je pouvais.)

import itertools
izip = itertools.zip_longest
chain = itertools.chain.from_iterable
compress = itertools.compress
def rwh_primes2_python3(n):
    """ Input n>=6, Returns a list of primes, 2 <= p < n """
    zero = bytearray([False])
    size = n//3 + (n % 6 == 2)
    sieve = bytearray([True]) * size
    sieve[0] = False
    for i in range(int(n**0.5)//3+1):
      if sieve[i]:
        k=3*i+1|1
        start = (k*k+4*k-2*k*(i&1))//3
        sieve[(k*k)//3::2*k]=zero*((size - (k*k)//3 - 1) // (2 * k) + 1)
        sieve[  start ::2*k]=zero*((size -   start  - 1) // (2 * k) + 1)
    ans = [2,3]
    poss = chain(izip(*[range(i, n, 6) for i in (1,5)]))
    ans.extend(compress(poss, sieve))
    return ans

Comparaisons:

>>> timeit.timeit('primes.rwh_primes2(10**6)', setup='import primes', number=1)
0.0652179726976101
>>> timeit.timeit('primes.rwh_primes2_python3(10**6)', setup='import primes', number=1)
0.03267321276325674

et

>>> timeit.timeit('primes.rwh_primes2(10**8)', setup='import primes', number=1)
6.394284538007014
>>> timeit.timeit('primes.rwh_primes2_python3(10**8)', setup='import primes', number=1)
3.833829450302801
Jason
la source
Mise en œuvre très cool. :)
Krish
7

Il est instructif d'écrire votre propre code de recherche principal, mais il est également utile d'avoir une bibliothèque rapide et fiable à portée de main. J'ai écrit un wrapper autour de la bibliothèque C ++ primesieve , je l' ai nommé primesieve-python

Essayez-le pip install primesieve

import primesieve
primes = primesieve.generate_primes(10**8)

Je serais curieux de voir la vitesse comparée.

Colonel Panic
la source
Ce n'est pas exactement ce que l'OP a ordonné mais je ne vois pas pourquoi le downvote. C'est une solution 2.8sec contrairement à certains autres modules externes. J'ai remarqué dans la source qu'il est fileté, avez-vous des tests sur la façon dont il évolue?
ljetibo
@ljetibo cheers. Le goulot d'étranglement semble copier le vecteur C ++ dans la liste Python, donc la count_primesfonction est beaucoup plus rapide quegenerate_primes
Colonel Panic
Sur mon ordinateur, il peut facilement générer des nombres premiers jusqu'à 1e8 (il donne MemoryError pour 1e9) et compter des nombres premiers jusqu'à 1e10. @HappyLeapSecond ci-dessus compare les algorithmes pour 1e6
Colonel Panic
7

Voici deux versions mises à jour (Python 3.6 pur) d'une des fonctions les plus rapides,

from itertools import compress

def rwh_primes1v1(n):
    """ Returns  a list of primes < n for n > 2 """
    sieve = bytearray([True]) * (n//2)
    for i in range(3,int(n**0.5)+1,2):
        if sieve[i//2]:
            sieve[i*i//2::i] = bytearray((n-i*i-1)//(2*i)+1)
    return [2,*compress(range(3,n,2), sieve[1:])]

def rwh_primes1v2(n):
    """ Returns a list of primes < n for n > 2 """
    sieve = bytearray([True]) * (n//2+1)
    for i in range(1,int(n**0.5)//2+1):
        if sieve[i]:
            sieve[2*i*(i+1)::2*i+1] = bytearray((n//2-2*i*(i+1))//(2*i+1)+1)
    return [2,*compress(range(3,n,2), sieve[1:])]
Bruno Astrolino
la source
1
Dans Python 3, j'ai utilisé cette fonction stackoverflow.com/a/3035188/7799269 mais remplacé / par // et xrange avec plage et ils semblaient beaucoup plus rapides que cela.
samerivertwice
4

Une mise en œuvre déterministe du test de primalité de Miller-Rabin sur l'hypothèse que N <9 080 191

import sys
import random

def miller_rabin_pass(a, n):
    d = n - 1
    s = 0
    while d % 2 == 0:
        d >>= 1
        s += 1

    a_to_power = pow(a, d, n)
    if a_to_power == 1:
        return True
    for i in xrange(s-1):
        if a_to_power == n - 1:
            return True
        a_to_power = (a_to_power * a_to_power) % n
    return a_to_power == n - 1


def miller_rabin(n):
    for a in [2, 3, 37, 73]:
      if not miller_rabin_pass(a, n):
        return False
    return True


n = int(sys.argv[1])
primes = [2]
for p in range(3,n,2):
  if miller_rabin(p):
    primes.append(p)
print len(primes)

Selon l'article sur Wikipedia ( http://en.wikipedia.org/wiki/Miller–Rabin_primality_test ), tester N <9 080 191 pour a = 2,3,37, et 73 est suffisant pour décider si N est composite ou non.

Et j'ai adapté le code source de l'implémentation probabiliste du test de Miller-Rabin original trouvé ici: http://en.literateprograms.org/Miller-Rabin_primality_test_(Python)

Ruggiero Spearman
la source
1
Merci pour le test de primalité de Miller-Rabin, mais ce code est en fait plus lent et ne fournit pas les résultats corrects. 37 est premier et ne réussit pas le test.
jbochi
Je suppose que 37 est l'un des cas spéciaux, mon mauvais. J'avais cependant bon espoir pour la version déterministe :)
Ruggiero Spearman
Il n'y a pas de cas particulier pour rabin miller.
Misguided
2
Vous avez mal lu l'article. C'est 31, pas 37. C'est pourquoi votre implémentation échoue.
Logan
4

Si vous contrôlez N, le moyen le plus rapide de répertorier tous les nombres premiers est de les précalculer. Sérieusement. Le précalcul est un moyen d'optimisation négligé.

Dave W. Smith
la source
3
Ou téléchargez-les ici primes.utm.edu/lists/small/millions , mais l'idée est de tester la limite de python et de voir si un beau code émerge de l'optimisation.
jbochi
4

Voici le code que j'utilise normalement pour générer des nombres premiers en Python:

$ python -mtimeit -s'import sieve' 'sieve.sieve(1000000)' 
10 loops, best of 3: 445 msec per loop
$ cat sieve.py
from math import sqrt

def sieve(size):
 prime=[True]*size
 rng=xrange
 limit=int(sqrt(size))

 for i in rng(3,limit+1,+2):
  if prime[i]:
   prime[i*i::+i]=[False]*len(prime[i*i::+i])

 return [2]+[i for i in rng(3,size,+2) if prime[i]]

if __name__=='__main__':
 print sieve(100)

Il ne peut pas rivaliser avec les solutions plus rapides publiées ici, mais au moins c'est du pur python.

Merci d'avoir posté cette question. J'ai vraiment beaucoup appris aujourd'hui.

MAK
la source
3

Pour le code le plus rapide, la solution numpy est la meilleure. Pour des raisons purement académiques, cependant, je publie ma version en python pur, qui est un peu moins de 50% plus rapide que la version du livre de cuisine publiée ci-dessus. Puisque je fais la liste entière en mémoire, vous avez besoin de suffisamment d'espace pour tout contenir, mais cela semble assez bien évoluer.

def daniel_sieve_2(maxNumber):
    """
    Given a number, returns all numbers less than or equal to
    that number which are prime.
    """
    allNumbers = range(3, maxNumber+1, 2)
    for mIndex, number in enumerate(xrange(3, maxNumber+1, 2)):
        if allNumbers[mIndex] == 0:
            continue
        # now set all multiples to 0
        for index in xrange(mIndex+number, (maxNumber-3)/2+1, number):
            allNumbers[index] = 0
    return [2] + filter(lambda n: n!=0, allNumbers)

Et les résultats:

>>>mine = timeit.Timer("daniel_sieve_2(1000000)",
...                    "from sieves import daniel_sieve_2")
>>>prev = timeit.Timer("get_primes_erat(1000000)",
...                    "from sieves import get_primes_erat")
>>>print "Mine: {0:0.4f} ms".format(min(mine.repeat(3, 1))*1000)
Mine: 428.9446 ms
>>>print "Previous Best {0:0.4f} ms".format(min(prev.repeat(3, 1))*1000)
Previous Best 621.3581 ms
Daniel G
la source
3

Une implémentation légèrement différente d'un demi-tamis utilisant Numpy:

http://rebrained.com/?p=458

importer des mathématiques
importer numpy
def prime6 (jusqu'à):
    nombres premiers = numpy.arange (3, jusqu'à + 1,2)
    isprime = numpy.ones ((jusqu'à-1) / 2, dtype = bool)
    pour le facteur en nombres premiers [: int (math.sqrt (upto))]:
        si isprime [(facteur 2) / 2]: isprime [(facteur * 3-2) / 2: (jusqu'à 1) / 2: facteur] = 0
    return numpy.insert (nombres premiers [isprime], 0,2)

Quelqu'un peut-il comparer cela avec les autres horaires? Sur ma machine, il semble assez comparable aux autres demi-tamis Numpy.

nolfonzo
la source
upto=10**6: primesfrom2to()- 7 ms; prime6()- 12 ms ideone.com/oDg2Y
jfs
3

Tout est écrit et testé. Il n'est donc pas nécessaire de réinventer la roue.

python -m timeit -r10 -s"from sympy import sieve" "primes = list(sieve.primerange(1, 10**6))"

nous donne un record de 12,2 ms !

10 loops, best of 10: 12.2 msec per loop

Si ce n'est pas assez rapide, vous pouvez essayer PyPy:

pypy -m timeit -r10 -s"from sympy import sieve" "primes = list(sieve.primerange(1, 10**6))"

ce qui se traduit par:

10 loops, best of 10: 2.03 msec per loop

La réponse avec 247 votes positifs répertorie 15,9 ms pour la meilleure solution. Comparez cela !!!

lifolofi
la source
3

J'ai testé certaines fonctions d'unutbu , je les ai calculées avec des millions de millions

Les gagnants sont les fonctions qui utilisent la bibliothèque numpy,

Remarque : Il serait également intéressant de faire un test d'utilisation de la mémoire :)

Résultat du temps de calcul

Exemple de code

Code complet sur mon dépôt github

#!/usr/bin/env python

import lib
import timeit
import sys
import math
import datetime

import prettyplotlib as ppl
import numpy as np

import matplotlib.pyplot as plt
from prettyplotlib import brewer2mpl

primenumbers_gen = [
    'sieveOfEratosthenes',
    'ambi_sieve',
    'ambi_sieve_plain',
    'sundaram3',
    'sieve_wheel_30',
    'primesfrom3to',
    'primesfrom2to',
    'rwh_primes',
    'rwh_primes1',
    'rwh_primes2',
]

def human_format(num):
    # /programming/579310/formatting-long-numbers-as-strings-in-python?answertab=active#tab-top
    magnitude = 0
    while abs(num) >= 1000:
        magnitude += 1
        num /= 1000.0
    # add more suffixes if you need them
    return '%.2f%s' % (num, ['', 'K', 'M', 'G', 'T', 'P'][magnitude])


if __name__=='__main__':

    # Vars
    n = 10000000 # number itereration generator
    nbcol = 5 # For decompose prime number generator
    nb_benchloop = 3 # Eliminate false positive value during the test (bench average time)
    datetimeformat = '%Y-%m-%d %H:%M:%S.%f'
    config = 'from __main__ import n; import lib'
    primenumbers_gen = {
        'sieveOfEratosthenes': {'color': 'b'},
        'ambi_sieve': {'color': 'b'},
        'ambi_sieve_plain': {'color': 'b'},
         'sundaram3': {'color': 'b'},
        'sieve_wheel_30': {'color': 'b'},
# # #        'primesfrom2to': {'color': 'b'},
        'primesfrom3to': {'color': 'b'},
        # 'rwh_primes': {'color': 'b'},
        # 'rwh_primes1': {'color': 'b'},
        'rwh_primes2': {'color': 'b'},
    }


    # Get n in command line
    if len(sys.argv)>1:
        n = int(sys.argv[1])

    step = int(math.ceil(n / float(nbcol)))
    nbs = np.array([i * step for i in range(1, int(nbcol) + 1)])
    set2 = brewer2mpl.get_map('Paired', 'qualitative', 12).mpl_colors

    print datetime.datetime.now().strftime(datetimeformat)
    print("Compute prime number to %(n)s" % locals())
    print("")

    results = dict()
    for pgen in primenumbers_gen:
        results[pgen] = dict()
        benchtimes = list()
        for n in nbs:
            t = timeit.Timer("lib.%(pgen)s(n)" % locals(), setup=config)
            execute_times = t.repeat(repeat=nb_benchloop,number=1)
            benchtime = np.mean(execute_times)
            benchtimes.append(benchtime)
        results[pgen] = {'benchtimes':np.array(benchtimes)}

fig, ax = plt.subplots(1)
plt.ylabel('Computation time (in second)')
plt.xlabel('Numbers computed')
i = 0
for pgen in primenumbers_gen:

    bench = results[pgen]['benchtimes']
    avgs = np.divide(bench,nbs)
    avg = np.average(bench, weights=nbs)

    # Compute linear regression
    A = np.vstack([nbs, np.ones(len(nbs))]).T
    a, b = np.linalg.lstsq(A, nbs*avgs)[0]

    # Plot
    i += 1
    #label="%(pgen)s" % locals()
    #ppl.plot(nbs, nbs*avgs, label=label, lw=1, linestyle='--', color=set2[i % 12])
    label="%(pgen)s avg" % locals()
    ppl.plot(nbs, a * nbs + b, label=label, lw=2, color=set2[i % 12])
print datetime.datetime.now().strftime(datetimeformat)

ppl.legend(ax, loc='upper left', ncol=4)

# Change x axis label
ax.get_xaxis().get_major_formatter().set_scientific(False)
fig.canvas.draw()
labels = [human_format(int(item.get_text())) for item in ax.get_xticklabels()]

ax.set_xticklabels(labels)
ax = plt.gca()

plt.show()
Bruno Adelé
la source
2
pour comparer les performances algorithmiques , il est préférable de tracer à une échelle log-log .
Will Ness
3

Pour Python 3

def rwh_primes2(n):
    correction = (n%6>1)
    n = {0:n,1:n-1,2:n+4,3:n+3,4:n+2,5:n+1}[n%6]
    sieve = [True] * (n//3)
    sieve[0] = False
    for i in range(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+4*k-2*k*(i&1))//3::2*k]=[False]*((n//6-(k*k+4*k-2*k*(i&1))//6-1)//k+1)
    return [2,3] + [3*i+1|1 for i in range(1,n//3-correction) if sieve[i]]
SmartManoj
la source
3

Tamis primaire le plus rapide en Pure Python :

from itertools import compress

def half_sieve(n):
    """
    Returns a list of prime numbers less than `n`.
    """
    if n <= 2:
        return []
    sieve = bytearray([True]) * (n // 2)
    for i in range(3, int(n ** 0.5) + 1, 2):
        if sieve[i // 2]:
            sieve[i * i // 2::i] = bytearray((n - i * i - 1) // (2 * i) + 1)
    primes = list(compress(range(1, n, 2), sieve))
    primes[0] = 2
    return primes

J'ai optimisé Sieve of Eratosthenes pour la vitesse et la mémoire.

Référence

from time import clock
import platform

def benchmark(iterations, limit):
    start = clock()
    for x in range(iterations):
        half_sieve(limit)
    end = clock() - start
    print(f'{end/iterations:.4f} seconds for primes < {limit}')

if __name__ == '__main__':
    print(platform.python_version())
    print(platform.platform())
    print(platform.processor())
    it = 10
    for pw in range(4, 9):
        benchmark(it, 10**pw)

Production

>>> 3.6.7
>>> Windows-10-10.0.17763-SP0
>>> Intel64 Family 6 Model 78 Stepping 3, GenuineIntel
>>> 0.0003 seconds for primes < 10000
>>> 0.0021 seconds for primes < 100000
>>> 0.0204 seconds for primes < 1000000
>>> 0.2389 seconds for primes < 10000000
>>> 2.6702 seconds for primes < 100000000
MrSeeker
la source
2

La première fois que vous utilisez python, certaines des méthodes que j'utilise peuvent sembler un peu lourdes. Je viens de convertir mon code c ++ en python et c'est ce que j'ai (bien qu'un peu lent en python)

#!/usr/bin/env python
import time

def GetPrimes(n):

    Sieve = [1 for x in xrange(n)]

    Done = False
    w = 3

    while not Done:

        for q in xrange (3, n, 2):
            Prod = w*q
            if Prod < n:
                Sieve[Prod] = 0
            else:
                break

        if w > (n/2):
            Done = True
        w += 2

    return Sieve



start = time.clock()

d = 10000000
Primes = GetPrimes(d)

count = 1 #This is for 2

for x in xrange (3, d, 2):
    if Primes[x]:
        count+=1

elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"

pythonw Primes.py

664579 nombres premiers trouvés en 12.799119 secondes!

#!/usr/bin/env python
import time

def GetPrimes2(n):

    Sieve = [1 for x in xrange(n)]

    for q in xrange (3, n, 2):
        k = q
        for y in xrange(k*3, n, k*2):
            Sieve[y] = 0

    return Sieve



start = time.clock()

d = 10000000
Primes = GetPrimes2(d)

count = 1 #This is for 2

for x in xrange (3, d, 2):
    if Primes[x]:
        count+=1

elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"

pythonw Primes2.py

664579 nombres premiers trouvés en 10.230172 secondes!

#!/usr/bin/env python
import time

def GetPrimes3(n):

    Sieve = [1 for x in xrange(n)]

    for q in xrange (3, n, 2):
        k = q
        for y in xrange(k*k, n, k << 1):
            Sieve[y] = 0

    return Sieve



start = time.clock()

d = 10000000
Primes = GetPrimes3(d)

count = 1 #This is for 2

for x in xrange (3, d, 2):
    if Primes[x]:
        count+=1

elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"

python Primes2.py

664579 nombres premiers trouvés en 7.113776 secondes!

smac89
la source
2

Je sais que le concours est fermé depuis quelques années. …

Néanmoins, c'est ma suggestion pour un tamis primaire en python pur, basé sur l'omission des multiples de 2, 3 et 5 en utilisant les étapes appropriées lors du traitement du tamis vers l'avant. Néanmoins, il est en réalité plus lent pour N <10 ^ 9 que @Robert William Hanks solutions supérieures rwh_primes2 et rwh_primes1. En utilisant un tableau de tamis ctypes.c_ushort supérieur à 1,5 * 10 ^ 8, il s'adapte en quelque sorte aux limites de la mémoire.

10 ^ 6

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (1000000)" 10 boucles, au mieux de 3: 46,7 msec par boucle

à comparer: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (1000000)" 10 boucles, meilleur de 3: 43,2 msec par boucle à comparer: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2 (1000000) "10 boucles, meilleur de 3: 34,5 msec par boucle

10 ^ 7

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (10000000)" 10 boucles, au mieux de 3: 530 msec par boucle

à comparer: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (10000000)" 10 boucles, au mieux de 3: 494 msec par boucle à comparer: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2 (10000000) "10 boucles, meilleur de 3: 375 ms par boucle

10 ^ 8

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (100000000)" 10 boucles, au mieux de 3: 5,55 s par boucle

à comparer: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (100000000)" 10 boucles, meilleur de 3: 5,33 sec par boucle à comparer: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2 (100000000) "10 boucles, meilleur de 3: 3,95 s par boucle

10 ^ 9

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (1000000000)" 10 boucles, au mieux de 3: 61,2 s par boucle

à comparer: $ python -mtimeit -n 3 -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (1000000000)" 3 boucles, le meilleur de 3: 97,8 s par boucle

à comparer: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2 (1000000000)" 10 boucles, au mieux de 3: 41,9 s par boucle

Vous pouvez copier le code ci-dessous dans ubuntus primeSieveSpeedComp pour revoir ces tests.

def primeSieveSeq(MAX_Int):
    if MAX_Int > 5*10**8:
        import ctypes
        int16Array = ctypes.c_ushort * (MAX_Int >> 1)
        sieve = int16Array()
        #print 'uses ctypes "unsigned short int Array"'
    else:
        sieve = (MAX_Int >> 1) * [False]
        #print 'uses python list() of long long int'
    if MAX_Int < 10**8:
        sieve[4::3] = [True]*((MAX_Int - 8)/6+1)
        sieve[12::5] = [True]*((MAX_Int - 24)/10+1)
    r = [2, 3, 5]
    n = 0
    for i in xrange(int(MAX_Int**0.5)/30+1):
        n += 3
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 2
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 1
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 2
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 1
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 2
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 3
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 1
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
    if MAX_Int < 10**8:
        return [2, 3, 5]+[(p << 1) + 1 for p in [n for n in xrange(3, MAX_Int >> 1) if not sieve[n]]]
    n = n >> 1
    try:
        for i in xrange((MAX_Int-2*n)/30 + 1):
            n += 3
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 2
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 1
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 2
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 1
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 2
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 3
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 1
            if not sieve[n]:
                r.append((n << 1) + 1)
    except:
        pass
    return r
ABri
la source
pour visualiser vos résultats de test, tracez-les sur une échelle log-log, pour voir et comparer les ordres empiriques de croissance .
Will Ness
@ Will merci pour l'entrée,
j'aurai
1

Voici une version numpy de Sieve of Eratosthenes ayant à la fois une bonne complexité (inférieure au tri d'un tableau de longueur n) et une vectorisation. Comparé à @unutbu, c'est aussi rapide que les packages avec 46 microsecondes pour trouver tous les nombres premiers inférieurs à un million.

import numpy as np 
def generate_primes(n):
    is_prime = np.ones(n+1,dtype=bool)
    is_prime[0:2] = False
    for i in range(int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i**2::i]=False
    return np.where(is_prime)[0]

Calendrier:

import time    
for i in range(2,10):
    timer =time.time()
    generate_primes(10**i)
    print('n = 10^',i,' time =', round(time.time()-timer,6))

>> n = 10^ 2  time = 5.6e-05
>> n = 10^ 3  time = 6.4e-05
>> n = 10^ 4  time = 0.000114
>> n = 10^ 5  time = 0.000593
>> n = 10^ 6  time = 0.00467
>> n = 10^ 7  time = 0.177758
>> n = 10^ 8  time = 1.701312
>> n = 10^ 9  time = 19.322478
Peter Mølgaard Pallesen
la source
1

J'ai mis à jour une grande partie du code pour Python 3 et je l'ai lancé sur perfplot (un de mes projets) pour voir lequel est réellement le plus rapide. Il s'avère que, pour les gros n, primesfrom{2,3}toprenez le gâteau:

entrez la description de l'image ici


Code pour reproduire l'intrigue:

import perfplot
from math import sqrt, ceil
import numpy as np
import sympy


def rwh_primes(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    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):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    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):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """Input n>=6, Returns a list of primes, 2 <= p < n"""
    assert n >= 6
    correction = n % 6 > 1
    n = {0: n, 1: n - 1, 2: n + 4, 3: n + 3, 4: n + 2, 5: n + 1}[n % 6]
    sieve = [True] * (n // 3)
    sieve[0] = False
    for i in range(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 + 4 * k - 2 * k * (i & 1)) // 3::2 * k] = [False] * (
                (n // 6 - (k * k + 4 * k - 2 * k * (i & 1)) // 6 - 1) // k + 1
            )
    return [2, 3] + [3 * i + 1 | 1 for i in range(1, n // 3 - correction) if sieve[i]]


def sieve_wheel_30(N):
    # http://zerovolt.com/?p=88
    """ Returns a list of primes <= N using wheel criterion 2*3*5 = 30

Copyright 2009 by zerovolt.com
This code is free for non-commercial purposes, in which case you can just leave this comment as a credit for my work.
If you need this code for commercial purposes, please contact me by sending an email to: info [at] zerovolt [dot] com."""
    __smallp = (
        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,
        101,
        103,
        107,
        109,
        113,
        127,
        131,
        137,
        139,
        149,
        151,
        157,
        163,
        167,
        173,
        179,
        181,
        191,
        193,
        197,
        199,
        211,
        223,
        227,
        229,
        233,
        239,
        241,
        251,
        257,
        263,
        269,
        271,
        277,
        281,
        283,
        293,
        307,
        311,
        313,
        317,
        331,
        337,
        347,
        349,
        353,
        359,
        367,
        373,
        379,
        383,
        389,
        397,
        401,
        409,
        419,
        421,
        431,
        433,
        439,
        443,
        449,
        457,
        461,
        463,
        467,
        479,
        487,
        491,
        499,
        503,
        509,
        521,
        523,
        541,
        547,
        557,
        563,
        569,
        571,
        577,
        587,
        593,
        599,
        601,
        607,
        613,
        617,
        619,
        631,
        641,
        643,
        647,
        653,
        659,
        661,
        673,
        677,
        683,
        691,
        701,
        709,
        719,
        727,
        733,
        739,
        743,
        751,
        757,
        761,
        769,
        773,
        787,
        797,
        809,
        811,
        821,
        823,
        827,
        829,
        839,
        853,
        857,
        859,
        863,
        877,
        881,
        883,
        887,
        907,
        911,
        919,
        929,
        937,
        941,
        947,
        953,
        967,
        971,
        977,
        983,
        991,
        997,
    )
    # wheel = (2, 3, 5)
    const = 30
    if N < 2:
        return []
    if N <= const:
        pos = 0
        while __smallp[pos] <= N:
            pos += 1
        return list(__smallp[:pos])
    # make the offsets list
    offsets = (7, 11, 13, 17, 19, 23, 29, 1)
    # prepare the list
    p = [2, 3, 5]
    dim = 2 + N // const
    tk1 = [True] * dim
    tk7 = [True] * dim
    tk11 = [True] * dim
    tk13 = [True] * dim
    tk17 = [True] * dim
    tk19 = [True] * dim
    tk23 = [True] * dim
    tk29 = [True] * dim
    tk1[0] = False
    # help dictionary d
    # d[a , b] = c  ==> if I want to find the smallest useful multiple of (30*pos)+a
    # on tkc, then I need the index given by the product of [(30*pos)+a][(30*pos)+b]
    # in general. If b < a, I need [(30*pos)+a][(30*(pos+1))+b]
    d = {}
    for x in offsets:
        for y in offsets:
            res = (x * y) % const
            if res in offsets:
                d[(x, res)] = y
    # another help dictionary: gives tkx calling tmptk[x]
    tmptk = {1: tk1, 7: tk7, 11: tk11, 13: tk13, 17: tk17, 19: tk19, 23: tk23, 29: tk29}
    pos, prime, lastadded, stop = 0, 0, 0, int(ceil(sqrt(N)))

    # inner functions definition
    def del_mult(tk, start, step):
        for k in range(start, len(tk), step):
            tk[k] = False

    # end of inner functions definition
    cpos = const * pos
    while prime < stop:
        # 30k + 7
        if tk7[pos]:
            prime = cpos + 7
            p.append(prime)
            lastadded = 7
            for off in offsets:
                tmp = d[(7, off)]
                start = (
                    (pos + prime)
                    if off == 7
                    else (prime * (const * (pos + 1 if tmp < 7 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 11
        if tk11[pos]:
            prime = cpos + 11
            p.append(prime)
            lastadded = 11
            for off in offsets:
                tmp = d[(11, off)]
                start = (
                    (pos + prime)
                    if off == 11
                    else (prime * (const * (pos + 1 if tmp < 11 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 13
        if tk13[pos]:
            prime = cpos + 13
            p.append(prime)
            lastadded = 13
            for off in offsets:
                tmp = d[(13, off)]
                start = (
                    (pos + prime)
                    if off == 13
                    else (prime * (const * (pos + 1 if tmp < 13 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 17
        if tk17[pos]:
            prime = cpos + 17
            p.append(prime)
            lastadded = 17
            for off in offsets:
                tmp = d[(17, off)]
                start = (
                    (pos + prime)
                    if off == 17
                    else (prime * (const * (pos + 1 if tmp < 17 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 19
        if tk19[pos]:
            prime = cpos + 19
            p.append(prime)
            lastadded = 19
            for off in offsets:
                tmp = d[(19, off)]
                start = (
                    (pos + prime)
                    if off == 19
                    else (prime * (const * (pos + 1 if tmp < 19 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 23
        if tk23[pos]:
            prime = cpos + 23
            p.append(prime)
            lastadded = 23
            for off in offsets:
                tmp = d[(23, off)]
                start = (
                    (pos + prime)
                    if off == 23
                    else (prime * (const * (pos + 1 if tmp < 23 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 29
        if tk29[pos]:
            prime = cpos + 29
            p.append(prime)
            lastadded = 29
            for off in offsets:
                tmp = d[(29, off)]
                start = (
                    (pos + prime)
                    if off == 29
                    else (prime * (const * (pos + 1 if tmp < 29 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # now we go back to top tk1, so we need to increase pos by 1
        pos += 1
        cpos = const * pos
        # 30k + 1
        if tk1[pos]:
            prime = cpos + 1
            p.append(prime)
            lastadded = 1
            for off in offsets:
                tmp = d[(1, off)]
                start = (
                    (pos + prime)
                    if off == 1
                    else (prime * (const * pos + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
    # time to add remaining primes
    # if lastadded == 1, remove last element and start adding them from tk1
    # this way we don't need an "if" within the last while
    if lastadded == 1:
        p.pop()
    # now complete for every other possible prime
    while pos < len(tk1):
        cpos = const * pos
        if tk1[pos]:
            p.append(cpos + 1)
        if tk7[pos]:
            p.append(cpos + 7)
        if tk11[pos]:
            p.append(cpos + 11)
        if tk13[pos]:
            p.append(cpos + 13)
        if tk17[pos]:
            p.append(cpos + 17)
        if tk19[pos]:
            p.append(cpos + 19)
        if tk23[pos]:
            p.append(cpos + 23)
        if tk29[pos]:
            p.append(cpos + 29)
        pos += 1
    # remove exceeding if present
    pos = len(p) - 1
    while p[pos] > N:
        pos -= 1
    if pos < len(p) - 1:
        del p[pos + 1 :]
    # return p list
    return p


def sieve_of_eratosthenes(n):
    """sieveOfEratosthenes(n): return the list of the primes < n."""
    # Code from: <[email protected]>, 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 sieve_of_atkin(end):
    """return a list of all the prime numbers <end using the Sieve of Atkin."""
    # Code by Steve Krenzel, <[email protected]>, 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 ambi_sieve_plain(n):
    s = list(range(3, n, 2))
    for m in range(3, int(n ** 0.5) + 1, 2):
        if s[(m - 3) // 2]:
            for t in range((m * m - 3) // 2, (n >> 1) - 1, m):
                s[t] = 0
    return [2] + [t for t in s if t > 0]


def sundaram3(max_n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/2073279#2073279
    numbers = 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] + filter(None, numbers)


# Using Numpy:
def ambi_sieve(n):
    # http://tommih.blogspot.com/2009/04/fast-prime-number-generator.html
    s = np.arange(3, n, 2)
    for m in range(3, int(n ** 0.5) + 1, 2):
        if s[(m - 3) // 2]:
            s[(m * m - 3) // 2::m] = 0
    return np.r_[2, s[s > 0]]


def primesfrom3to(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns an array of primes, p < n """
    assert n >= 2
    sieve = np.ones(n // 2, dtype=np.bool)
    for i in range(3, int(n ** 0.5) + 1, 2):
        if sieve[i // 2]:
            sieve[i * i // 2::i] = False
    return np.r_[2, 2 * np.nonzero(sieve)[0][1::] + 1]


def primesfrom2to(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns an array of primes, 2 <= p < n """
    assert n >= 6
    sieve = np.ones(n // 3 + (n % 6 == 2), dtype=np.bool)
    sieve[0] = False
    for i in range(int(n ** 0.5) // 3 + 1):
        if sieve[i]:
            k = 3 * i + 1 | 1
            sieve[((k * k) // 3)::2 * k] = False
            sieve[(k * k + 4 * k - 2 * k * (i & 1)) // 3::2 * k] = False
    return np.r_[2, 3, ((3 * np.nonzero(sieve)[0] + 1) | 1)]


def sympy_sieve(n):
    return list(sympy.sieve.primerange(1, n))


perfplot.save(
    "prime.png",
    setup=lambda n: n,
    kernels=[
        rwh_primes,
        rwh_primes1,
        rwh_primes2,
        sieve_wheel_30,
        sieve_of_eratosthenes,
        sieve_of_atkin,
        # ambi_sieve_plain,
        # sundaram3,
        ambi_sieve,
        primesfrom3to,
        primesfrom2to,
        sympy_sieve,
    ],
    n_range=[2 ** k for k in range(3, 25)],
    logx=True,
    logy=True,
    xlabel="n",
)
Nico Schlömer
la source
0

Je suppose que le plus rapide de toutes les façons est de coder en dur les nombres premiers dans votre code.

Alors pourquoi ne pas simplement écrire un script lent qui génère un autre fichier source contenant tous les nombres, puis importer ce fichier source lorsque vous exécutez votre programme réel.

Bien sûr, cela ne fonctionne que si vous connaissez la limite supérieure de N au moment de la compilation, mais c'est donc le cas pour (presque) tous les problèmes d'Euler du projet.

 

PS: Je peux me tromper, même si l'analyse de la source avec des nombres premiers câblés est plus lente que leur calcul en premier lieu, mais pour autant que je sache, Python s'exécute à partir de .pycfichiers compilés , donc la lecture d'un tableau binaire avec tous les nombres premiers jusqu'à N devrait être sanglante rapide dans ce cas.

akuhn
la source
0

Désolé de déranger mais erat2 () a un sérieux défaut dans l'algorithme.

Lors de la recherche du prochain composite, nous devons uniquement tester les nombres impairs. q, p les deux sont impairs; alors q + p est pair et n'a pas besoin d'être testé, mais q + 2 * p est toujours impair. Cela élimine le test «si même» dans la condition de boucle while et économise environ 30% du temps d'exécution.

Pendant que nous y sommes: au lieu de l'élégant «D.pop (q, None)», la méthode get et delete utilise «if q in D: p = D [q], del D [q]» qui est deux fois plus rapide ! Au moins sur ma machine (P3-1Ghz). Je suggère donc cette implémentation de cet algorithme intelligent:

def erat3( ):
    from itertools import islice, count

    # q is the running integer that's checked for primeness.
    # yield 2 and no other even number thereafter
    yield 2
    D = {}
    # no need to mark D[4] as we will test odd numbers only
    for q in islice(count(3),0,None,2):
        if q in D:                  #  is composite
            p = D[q]
            del D[q]
            # q is composite. p=D[q] is the first prime that
            # divides it. Since we've reached q, we no longer
            # need it in the map, but we'll mark the next
            # multiple of its witnesses to prepare for larger
            # numbers.
            x = q + p+p        # next odd(!) multiple
            while x in D:      # skip composites
                x += p+p
            D[x] = p
        else:                  # is prime
            # q is a new prime.
            # Yield it and mark its first multiple that isn't
            # already marked in previous iterations.
            D[q*q] = q
            yield q
user1016274
la source
pour une addition différée de nombres premiers dans le dict (jusqu'à ce que le carré d'un nombre premier soit vu dans l'entrée), voir stackoverflow.com/a/10733621/849891 .
Will Ness
0

La méthode la plus rapide que j'ai essayée jusqu'à présent est basée sur la fonction livre de recettes Pythonerat2 :

import itertools as it
def erat2a( ):
    D = {  }
    yield 2
    for q in it.islice(it.count(3), 0, None, 2):
        p = D.pop(q, None)
        if p is None:
            D[q*q] = q
            yield q
        else:
            x = q + 2*p
            while x in D:
                x += 2*p
            D[x] = p

Voir cette réponse pour une explication de l'accélération.

tzot
la source
0

Je suis peut-être en retard à la fête, mais je devrai ajouter mon propre code pour cela. Il utilise environ n / 2 dans l'espace car nous n'avons pas besoin de stocker des nombres pairs et j'utilise également le module bitarray python, réduisant encore plus la consommation de mémoire et permettant de calculer tous les nombres premiers jusqu'à 1 000 000 000

from bitarray import bitarray
def primes_to(n):
    size = n//2
    sieve = bitarray(size)
    sieve.setall(1)
    limit = int(n**0.5)
    for i in range(1,limit):
        if sieve[i]:
            val = 2*i+1
            sieve[(i+i*val)::val] = 0
    return [2] + [2*i+1 for i, v in enumerate(sieve) if v and i > 0]

python -m timeit -n10 -s "import euler" "euler.primes_to(1000000000)"
10 loops, best of 3: 46.5 sec per loop

Il a été exécuté sur un MAC OSX 10.8.3 64 bits 2,4 GHz.

cobie
la source
1
poster un chronométrage pour une machine inconnue ne dit rien. La réponse acceptée ici dit "sans psyco, pour n = 1000000, rwh_primes2 était le plus rapide". Donc , si vous fournissez vos horaires pour ce code ainsi que le vôtre, sur la même machine, et à 2, 4, 10 millions et, alors ce serait beaucoup plus instructif.
Will Ness
-1, Ce code dépend des caractéristiques spéciales du bitarray implémenté en C, c'est pourquoi le code est rapide car la plupart du travail se fait en code natif dans l'affectation de tranche. Le paquet BitArray CASSE la définition standard pour les tranches convenables (indexé sur une plage) pour des séquences mutables en ce qu'il permet d' attribuer une valeur booléenne unique 0/1 ou Vrai / Faux à tous les éléments de la tranche, alors que le comportement standard pour python pur semble ne pas autoriser cela et autoriser uniquement la valeur d'affectation de 0, auquel cas elle est traitée comme une suppression de tous les éléments de tranche de la séquence / du tableau.
GordonBGood
suite: Si l'on compare le code natif non standard, nous pouvons aussi bien écrire un package de générateur de séquence "fastprimes" basé sur du code C tel que celui de la prime-tam de Kim Walisch et générer tous les nombres premiers dans les quatre milliards plus 32 -intervalle de numéros de bits en quelques secondes avec un seul appel au générateur de séquence. Cela n'utiliserait également presque pas de mémoire car le code lié est basé sur un tamis segmenté d'Ératosthène et n'utilise donc que quelques dizaines de kilo-octets de RAM, et si une séquence était générée, aucun stockage de liste ne serait nécessaire.
GordonBGood
0

J'ai collecté plusieurs tamis de nombre premier au fil du temps. Le plus rapide sur mon ordinateur est le suivant:

from time import time
# 175 ms for all the primes up to the value 10**6
def primes_sieve(limit):
    a = [True] * limit
    a[0] = a[1] = False
    #a[2] = True
    for n in xrange(4, limit, 2):
        a[n] = False
    root_limit = int(limit**.5)+1
    for i in xrange(3,root_limit):
        if a[i]:
            for n in xrange(i*i, limit, 2*i):
                a[n] = False
    return a

LIMIT = 10**6
s=time()
primes = primes_sieve(LIMIT)
print time()-s
Stefan Gruenwald
la source
0

Je suis lent à répondre à cette question, mais cela semblait être un exercice amusant. J'utilise numpy qui pourrait tricher et je doute que cette méthode soit la plus rapide mais elle devrait être claire. Il filtre un tableau booléen se référant uniquement à ses indices et obtient des nombres premiers à partir des indices de toutes les valeurs True. Aucun modulo nécessaire.

import numpy as np
def ajs_primes3a(upto):
    mat = np.ones((upto), dtype=bool)
    mat[0] = False
    mat[1] = False
    mat[4::2] = False
    for idx in range(3, int(upto ** 0.5)+1, 2):
        mat[idx*2::idx] = False
    return np.where(mat == True)[0]
Alan James Salmoni
la source
c'est incorrect par exemple, ajs_primes3a(10)-> array([2, 3, 5, 7, 9]). 9n'est pas un premier
jfs
Vous avez repéré un cas de bord que je n'avais pas - bien fait! Le problème était dans 'pour idx dans la plage (3, int (jusqu'à ** 0,5), 2):' qui devrait être 'pour idx dans la plage (3, int (jusqu'à ** 0,5) + 1, 2):'. Merci mais ça marche maintenant.
Alan James Salmoni
La raison en était que la boucle idx est passée à «jusqu'à ** 05», ce qui, pour les cas allant jusqu'à 15. inclusivement. Il s'agissait d'un ensemble de cas de bord pour lesquels je n'avais pas testé. Ajouter 1 signifie que cela devrait fonctionner pour tous les nombres.
Alan James Salmoni
Ça semble fonctionner maintenant. Il s'agit des numpysolutions parmi les plus lentes qui retournent un tableau. Remarque: aucune véritable implémentation de Sieve of Eratosthenes n'utilise modulo - pas besoin de le mentionner. Vous pouvez utiliser à la mat[idx*idx::idx]place de mat[idx*2::idx]. Et np.nonzero(mat)[0]au lieu de np.where(mat == True)[0].
jfs
Merci JF. J'ai testé contre prime6 () et obtenu un résultat plus rapide jusqu'à (IIRC) environ 250k lorsque prime6 () a pris le relais. primesfrom2to () était plus rapide. Jusqu'à 20 m, ajs_primes3a () a pris 0,034744977951ms, prime6 () a pris 0,0222899913788ms et primesfrom2to () a pris 0,0104751586914ms (même machine, même charge, au meilleur des 10 temps). C'est honnêtement mieux que ce que je pensais!
Alan James Salmoni
0

Voici une technique intéressante pour générer des nombres premiers (mais pas les plus efficaces) en utilisant les compréhensions de liste de python:

noprimes = [j for i in range(2, 8) for j in range(i*2, 50, i)]
primes = [x for x in range(2, 50) if x not in noprimes]

Vous pouvez trouver l'exemple et quelques explications ici

Alexandre
la source