Partage

Crible d'Ératosthène

1 2 >>

6 novembre 2011 à 11:11:10

Suite à un autre topic et sur conseil de candide, j'en ouvre un sur le crible d'Ératosthène.

Le crible d'Ératosthène est un procédé qui permet de trouver tous les nombres premiers inférieurs à un certain entier naturel donné N.

Donc, voici mon code en Python3:

def sieveOfEratosthenes(x):
    """Return a list of prime numbers from 2 to x
    build with the sieve of Eratosthenes"""
    prime_list = [2]

    # build a list with odd numbers after 2
    for i in range(3, x, 2):
        prime_list.append(i)

    # remove non prime numbers
    for i in prime_list:
        for n in prime_list:
            if n % i == 0 and i != n:
                prime_list.remove(n)

    return prime_list

ps : si je code en anglais, ce n'est pas pour faire hype mais pour pouvoir poster le code sur project euler et parce-que l'anglais est un peu une langue universelle, surtout en informatique.

Francky__laH

6 novembre 2011 à 11:41:33

Ça n'a pas l'air très efficace. Désolé.
Par tâtonnement, essaye de voir jusqu'où tu peux aller en une minute maxi.
Je vais te donner des pistes pour l'améliorer.
On verra une version simple, efficace et réutilisable (pour d'autres problèmes).

Ensuite, je te montrerai un crible qui déchire tout.
Objectif : aller le plus loin en une minute, on va dire pour obtenir la somme des nombres premiers inférieurs à limit.

from time import time
def primes(limit):
  # à compléter

st=time()
print(sum(p for p in primes(10**6)), time()-st)
# 10**6 c'est pour commencer, on devrait atteindre 10**9

Édit : j'ai la somme jusk 10⁹ en 56s et 3Go de RAM (avec mon crible qui déchire)
La somme jusk 10⁶ en 50ms. (et 135ms avec un crible simple)

-----

Édit : pour commencer, les données prendront beaucoup de place, il faut un truc simple qui prend moins de place.
Conseil : prime = [True]*(limit+1)
te permettra de stocker la primalité de chaque entier jusk limit.
Tu commences ensuite à changer prime[:2] = [False, False]

Ton idée de sauter un nombre sur deux après 3 est bonne, on s'en occupera après, et on fera mieux même.

candide

6 novembre 2011 à 13:13:24

Citation : carrion crow

def sieveOfEratosthenes(x):
    """Return a list of prime numbers from 2 to x
    build with the sieve of Eratosthenes"""
    prime_list = [2]

    # build a list with odd numbers after 2
    for i in range(3, x, 2):
        prime_list.append(i)

    # remove non prime numbers
    for i in prime_list:
        for n in prime_list:
            if n % i == 0 and i != n:
                prime_list.remove(n)

    return prime_list

Je n'ai pas testé la rapidité de ton code mais je veux bien croire qu'il est peu efficace comme le suggère la Hache.

Il est très lent de retirer un élément d'une liste Python (une liste Python est beaucoup moins une liste chaînée qu'un tableau C). L'idée de base et qui m'a toujours suffi sur Euler Project est de partir d'une liste P de longueur N remplie de 0 et de mettre 1 à la place de 0 lorsque dans ton crible tu COCHES (ou tu rayes, comme tu préfères). À la fin de ton crible, ta liste P est formée de 0 et de 1 avec la propriété que P[i]=0 si et seulement si i est premier.

Francky__laH

6 novembre 2011 à 13:22:10

Prime[p]=1 si et seulement si p est premier est conseillé. (1 ou True).
Édit : ou comme le suggère Candide ci dessous Compo[p]=False # ou 0, ssi p est premier.

J'attends ton nouveau code avant de donner une solution.
On va y aller progressivement. À la fin ça risque de piquer les yeux.

candide

6 novembre 2011 à 15:08:18

Citation : la Hache

Prime[p]=1 si et seulement si p est premier est conseillé. (1 ou True).

Je pense que c'est une affaire de convenance personnelle. Pour moi, j'initialise la liste à 0 car c'est naturel d'initialiser ainsi et je passe le 0 en 1 quand le nombre est «marqué» ie criblé. Mais je comprends qu'on puisse faire l'inverse.

Francky__laH

6 novembre 2011 à 23:07:50

Si ce fil devait inactif, sachez que vous trouverez des algos de cribles qui déchirent tout, en Python, en C++ (crible segmenté), ou autre langage, en résolvant le problème 10 du projet Euler.

Mais il vaut mieux commencer tranquillement.
Si quelqu'un demande à apprendre progressivement il faut nourrir ce fil.

carrion crow

6 novembre 2011 à 23:21:00

Bon, je tourne un peu en rond. Je ne sais pas ce qu'il me manque, des connaissances des librairies de Python ou de l'intelligence. J'ai essayé de faire un algo, résultat peu probant, il ne donne rien à 10^6 (j'ai bien attendu 5 min), puis un autre qui me donne 103,58s à 10^6, c'est mieux mais je n'atteint pas ton objectif "la Hache".
J'ai essayé de comprendre ce code, enfin plutôt les explications du site, car le code en C# je ne comprends pas tout.
Je suis ensuite tombé sur cette page et ai découvert que si on fait liste[2::2] ça sélectionne les multiples de 2 à partir de 2, j'ai ensuite refait le même genre d'algo mais en tâtonnant un peu et en essayant de comprendre. Du coup, j'ai fait sieveOfEratosthene3 qui met 1,54s à 10^6 et 18,88s à 10^7, par contre rien à 10^8 et memory error à 10^9

(j'ai 4Go de ram, et un core i540m, d'ailleurs dommage que python n'utilise qu'un seul thread au lieu des 4)

Sans me donner de réponse, sur quoi dois-je continuer ? Je pense qu'il y a surement une librairie Python codée en C qui permettrait de gagner en vitesse, ou un fonctionnement de Python qui m'échappe, genre ce type de truc est mieux par rapport au fonctionnement de python que cet autre type de truc.

def sieveOfEratosthenes(x):
    """Return a list of prime numbers from 2 to x
    build with the sieve of Eratosthenes"""
    prime_list = [2]

    # build a list with odd numbers after 2
    for i in range(3, x, 2):
        prime_list.append(i)

    # remove non prime numbers
    for i in prime_list:
        for n in prime_list:
            if n % i == 0 and i != n:
                prime_list.remove(n)

    return prime_list


def sieveOfEratosthene2(limit):
    """Return a list of numbers from 0 to limit
    build with the sieve of Eratosthenes. Second method
    True if index is prime"""
    prime_list = [True]*(limit+1)

    prime_list[:2] = [False, False]    # for 0 and 1

    # even number after 2 are not prime
    for i in range(4, limit+1, 2):
        prime_list[i] = False

    for i in range(3, limit+1, 2):
        if i and i**2 <= limit:    # only numbers that are still True
            for n in range(2, limit+1):
                if n*i <= limit:
                    prime_list[n*i] = False

    return prime_list


def sieveOfEratosthene3(limit):
    """Return a list of numbers from 0 to limit
    build with the sieve of Eratosthenes. Third method
    True if index is prime"""
    prime_list = [True]*(limit+1)

    prime_list[:2] = [False, False]    # for 0 and 1

    for i in range(2, limit+1):
        if i**2 <= limit:
            prime_list[i*2::i] = [False] * len(prime_list[i*2::i])

    return prime_list

Francky__laH

6 novembre 2011 à 23:31:55

Le trois est très bien, on va l'améliorer, c'est une bonne base.

Le point noir, c'est que tu fais len(bidule) et Python va créer la liste juste pour le plaisir d'avoir la longueur !!!
Essaye de calculer à la main cette longueur, avec i, limit.

Édit : autre point noir. Dans ta version 3, change

if i**2 <= limit:
  §§§§

if i*i <= limit:
  §§§§
else:
  break

i**2 est lent à côté de i*i (je sais c'est nul).
le break est salutaire, tu verras !!!

**********
Tu as réussi l'étape 1.

Solution étape1:

def primes(lim):
  """Générateur de nombres premier jusqu'à 'lim'."""
  # Initialisation
  compo=[True, True]+[False]*(lim-1) # indices de 0 à lim

  # On lance le crible
  for p in range(2, int(lim**0.5)+1):
    if not compo[p]: # p est premier
      compo[2*p::p]=[True]*(lim//p - 1)
      # mais les multiples suivants de p seront composés

  # Le crible est fini, on génère
  for p in range(2, lim+1):
    if not compo[p]: yield p
  return # fini

st=time() ; print(sum(primes(10**6)), time()-st)

Je l'ai écrit avec un générateur, il faut que tu comprennes le yield.
Ça permet d'envoyer un seul nombre premier à la fois, il est ici inutile de tous les stocker pour en faire la somme ; on a pas beaucoup de place !!!
Toute critique est la bienvenue. C'est l'étape 1. Objectif réutilisable à 100% pour d'autres situations.

On va améliorer ça, sans aide de librairies, rien de rien que de l'algo !!!
Et tu auras ton 10⁹ dans la minute.

Avec ça, tu devrais faire un bon temps pour 10⁶.

carrion crow

6 novembre 2011 à 23:43:14

Merci, je vais regarder ça.

Générateur et yield, voilà deux trucs que j'ai déjà croisé dans les livres/docs/codes, mais qu'il va falloir que je comprenne parfaitement.

Je ne savais pas que Python faisait une liste quand on fait len().

edit (je n'ai pas encore regardé ta correction):
J'ai remplacer la ligne :

prime_list[i*2::i] = [False] * len(prime_list[i*2::i])

par :

prime_list[i*2::i] = [False] * ((limit//i)-1)

Je fais 14,07s à 10**7 maintenant mais toujours pas de résultats à 10**8 (ou je n'attends pas assez longtemps)

Francky__laH

6 novembre 2011 à 23:52:38

Pour avoir la longueur d'un truc, il faut bien savoir ce qu'est truc.

---
Avec l'étape 1, tu as déjà un crible redoutable en performance à côté de nombre propositions.

*****
Étape 2: tu as remarqué dans ton premier post que la "moitié" était inutile, c'est de la place gâchée !!! On va gérer ça, au prix d'indices moins évidents.

À l'étape 1 : compo[p] = True ssi p est composé, ie non premier.
Maintenant, on va prendre que les nombres impair >=3.
Ils sont de la forme p=2*i+1.
On va créer un tableau compo, avec compo[i] = True ssi p=2*i+1 est composé.

voici le canevas, à toi de chercher à le compléter.

def primes(lim):
  if lim<2: return
  yield 2  # lui il faut pas l'oublier
  n = (lim-1)//2 # on a besoin de moins d'indices
  compo=[True]+[False]*n # initialisation
  for i in range(  ???? ):
    if not compo[i]:
      p = 2*i+1
      compo[ ??? ]=[True]*( ??? )
  for i in range(n+1):
    if not compo[i]: yield 2*i+1
  return


print(  list(prime(20))  ) # pour t'aider à débugger
st = time() ; print(sum(primes(10**6)), time()-st) # perf !!!

Si tu réussis, tu observeras un gain de temps 50% par rapport à l'étape 1.
Pour passer à l'étape 3, ce sera bien plus dur, et on ne gagnera qu'à peine 30%.
Si ensuite on poursuit l'aventure (étape 4), ce sera une autre direction : le crible segmenté pour diminuer l'empreinte mémoire, car il faut pas gâcher !!!

carrion crow

7 novembre 2011 à 11:52:04

Alors, pour le yield, quand tu fais sum(primes(x)), ça appel primes() en reprenant chaque fois à yield, et la boucle for fait le .next() ? J'ai bon ?

En remplaçant i**2 par i*i et en mettant le break, je passe à 4,52s pour 10**7 et 55,17s pour 10**8 (que je n'arrivais pas à calculer avant)

Voici mon quatrième code, basé sur ton canevas. J'ai modifié la condition, prime_list[i] = True ssi i*2+1 est premier (ça me semble plus logique de mettre True pour un premier)
Je mets 24,97s à 10**8 :-°

Mais ça me colle un MemoryError à la création de prime_list.

J'ai eu un peu de mal avec la ligne suivante, le coup du p=2*i+1 me perturbait au niveau des index de liste.

prime_list[i+p::p] = [False] * ((n-i)//p)

def sieveOfEratosthene4(limit):
    """Return a list of numbers from 0 to limit
    build with the sieve of Eratosthenes. Fourth method
    True if index is prime"""

    if limit < 2: return    # no primes numbers below 2
    yield 2    # first and only even number

    n = (limit-1) // 2

    prime_list = [False] + [True]*n    # create list entries for odd numbers
    # first entry is False because slice step cannot be 0

    for i in range(1, int((n**0.5)+1)):    # index[1] is number 3
        if prime_list[i]:
            p = 2*i+1    # equation for the odd numbers
            prime_list[i+p::p] = [False] * ((n-i)//p)
            # i+p to avoid [i] and ::p to have a p step

    for i in range(n+1):
        if prime_list[i]: yield 2*i+1

    return

edit :
Changement de :

for i in range(1, n+1):

par :

for i in range(1, int((n**0.5)+1)):

je passe à 13,74s à 10**8

edit2:
D'un point de vue mémoire True, False ou 1, 0 ça prend la même place, c'est ça ?

edit3: orthographe

cerium50

7 novembre 2011 à 13:47:43

A priori oui:

>>> import sys
>>> sys.getsizeof(0)
24
>>> sys.getsizeof(False)
24

Edit: pourquoi tu ne yield pas les nombres premiers quand tu les trouve (ligne 15/16), du coup tu pourrais diminuer la taille de ta boucle ligne 20 (j'ai pas vérifié mais il n'y a pas de raison que ça ne marche pas!) ?

Francky__laH

7 novembre 2011 à 13:51:26

@ carrion crow : first and not thirst ! (je corrige le reste ce soir)
C'est déjà très bien. Bravo.
Si tu veux, tu peux préparer une étape suivante (optionnelle):
faire un double crible : pour les 6k±1, avec deux tableaux donc.
Je mettrais le canevas ce soir.

****************

Citation : cerium50

A priori oui:
>>> import sys
>>> sys.getsizeof(0)
24
>>> sys.getsizeof(False)
24
Edit: pourquoi tu ne yield pas les nombres premiers quand tu les trouve (ligne 15/16), du coup tu pourrais diminuer la taille de ta boucle ligne 20 (j'ai pas vérifié mais il n'y a pas de raison que ça ne marche pas!) ?

Ce sont des bits, j'imagine, donc 24 bits = 3 octets, et ça corrobore le fait qu'un tableau de 10⁹ valeurs prenne 3Go.
Avec le crible segmenté, on pourra diminuer considérablement la place mémoire.

Sinon, l'utilisation de bitarray (un bit par valeur True ou False) permet de diviser la taille mémoire par 24, c'est un petit plus.
À suivre...

Édit : c'est par ce que là, on n'a que le début du tableau, on s'arrête à sqrt n,
pour tout le tableau, il vaut mieux le reprendre en entier.

*******

Édit : bon je suis plutôt d'accord (dès le départ) pour prime initialisé à True, plutôt que compo inialisé à False. C'est plus pratique à la fin : if prime[i]: yield .... Candide nous donneras son avis, il y a peut être une autre raison objective.

Sinon, je soupçonne que tu as vu mon code chez Euler :-°

, les indices sont simples, mais pas si évident à trouver !!!

Nombre : p=2i+1 | 3*p=6i+3=2*(3i+1)+1 | 5*p=10i+5=2*(5i+2)+1
Indice :    i   |    3i+1 = i+p       |     5i+2 = i+2p

Ce qui explique honnêtement le [i+p::p]
Le calcul de la longueur de ce tableau est donc (n-i)//p
On en déduit le prime[i+p::p] = [False]*( (n-i)//p ) qui est le killer de l'algo.

Voici donc ma correction de l'étape 2

def primes(lim):
  """Générateur de nombres premier jusqu'à 'lim'."""
  # cas particulier
  if lim<2: return
  yield 2
  # Initialisation
  n = (lim-1)//2
  prime = [True]*(n+1)
  prime[0] = False # 2*0+1 = 1 n'est pas premier
  # Crible
  for i in range((int(lim**0.5)+1)//2):
    if prime[i]:
      p = 2*i+1
      prime[p+i::p]=[False]*( (n-i)//p )
  # Génération
  for i in range(n+1):
    if prime[i]: yield 2*i+1
  return

for i in range(30):  # tests
  print(i, list(primes(i)))

st = time() ; print(sum(primes(10**6)), time()-st)

On constate un gain de -50% du temps d'exécution.

***************

Canevas pour l'étape 3 :
Préparer deux tableaux un pour les 6k-1, un autre pour les 6k+1.
On commencera donc de -1 et 1 (non premiers) à ...
Il faudra étudier les produits des nombres (6i±1)(6j±1) ; (4 cas) afin de remplir les deux tableaux.
Pour la fin il faut jouer fin, pour ne pas en générer un ou deux de trop !!!

def primes(lim):
  # cas particuliers
  if lim<5:
    if lim<2: return
    yield 2
    if lim<3: return
    yield 3
    return

  n = ???
  primeM=[True]*(n+1) # p=2*i-1
  primeP=[True]*(n+1) # p=2*i+1
  primeM[0]=False # 2*0-1 pas premier
  primeP[0]=False # 2*0+1 pas premier
  for i in range( ???? ):
    if primeM[i]: # 6*i-1 est premier
      p = 6*i-1
      primeM[ ???? ]=[False]*( ???? )
      primeP[ ???? ]=[False]*( ???? )
    if primeP[i]: # 6*i+1 est premier
      p = 6*i+1
      primeM[ ???? ]=[False]*( ???? )
      primeP[ ???? ]=[False]*( ???? )

  yield 2
  yield 3
  for i in range(n): # un rang avant la fin en n
    if primeM[i]: yield 6*i-1
    if primeP[i]: yield 6*i+1
  # pour mieux gérer la fin
  if primeM[n] and 6*n-1<=lim: yield 6*n-1
  if primeP[n] and 6*n+1<=lim: yield 6*n+1
  return

for i in range(30): # tests
  print(i, list(prime(i)))

st = time() ; print(sum(primes(10**6)), time()-st)

Ce code qui commence à devenir un peu lourd permet de gagner encore 30% en temps.
Avec ça, on gratte un mauvais programme en C++, à l'aise !
La barre des 10⁹ est en tout cas accessible dans la minute (mais avec 4Go de RAM au total), on poursuivra donc dans une autre direction pour gagner en performance, l'empreinte mémoire !

carrion crow

7 novembre 2011 à 18:11:18

Citation : la Hache

@ carrion crow : first and not thirst ! (je corrige le reste ce soir)
C'est déjà très bien. Bravo.

ok, corrigé, je sais pas à quoi je pensais (mais ça me fait apprendre un mot, thirst c'est la soif)
Merci

Citation : la Hache

Si tu veux, tu peux préparer une étape suivante (optionnelle):
faire un double crible : pour les 6k±1, avec deux tableaux donc.
Je mettrais le canevas ce soir.

...
Édit : bon je suis plutôt d'accord (dès le départ) pour prime initialisé à True, plutôt que compo inialisé à False. C'est plus pratique à la fin : if prime[i]: yield .... Candide nous donneras son avis, il y a peut être une autre raison objective.

Sinon, je soupçonne que tu as vu mon code chez Euler , les indices sont simples, mais pas si évident à trouver !!!

Même pas. J'ai juste essayé dans tous les sens, en réfléchissant quand même de temps en temps

J'en suis au problème 4 à PE, mais je ne l'ai pas commencé [HS]Il y a une méthode pour trouver ou c'est en commençant à 999*999 et en décrémentant ?[/HS]

Citation : la Hache

Canevas pour l'étape 3 :
Préparer deux tableaux un pour les 6k-1, un autre pour les 6k+1.
On commencera donc de -1 et 1 (non premiers) à ...
Il faudra étudier les produits des nombres (6i±1)(6j±1) ; (4 cas) afin de remplir les deux tableaux.
Pour la fin il faut jouer fin, pour ne pas en générer un ou deux de trop !!!
def primes(lim):
  if lim<2: return
  yield 2
  if lim<3: return
  yield 3
  if lim<5: return
  n = ???
  primeM=[True]*(n+1) # p=2*i-1
  primeP=[True]*(n+1) # p=2*i+1
  primeM[0]=False # 2*0-1 pas premier
  primeP[0]=False # 2*0+1 pas premier
  for i in range( ???? ):
    if primeM[i]: # 6*i-1 est premier
      p = 6*i-1
      primeM[ ???? ]=[False]*( ???? )
      primeP[ ???? ]=[False]*( ???? )
    if primeP[i]: # 6*i+1 est premier
      p = 6*i+1
      primeM[ ???? ]=[False]*( ???? )
      primeP[ ???? ]=[False]*( ???? )
  for i in range(n): # un rang avant la fin en n
    if primeM[i]: yield 6*i-1
    if primeP[i]: yield 6*i+1
  # pour mieux gérer la fin
  if primeM[n] and 6*n-1<=lim: yield 6*n-1
  if primeP[n] and 6*n+1<=lim: yield 6*n+1
  return

for i in range(30): # tests
  print(i, list(prime(i)))

st = time() ; print(sum(primes(10**6)), time()-st)
Ce code qui commence à devenir un peu lourd permet de gagner encore 30% en temps.
Avec ça, on gratte un mauvais programme en C++, à l'aise !
La barre des 10⁹ est en tout cas accessible dans la minute (mais avec 4Go de RAM au total), on poursuivra donc dans une autre direction pour gagner en performance, l'empreinte mémoire !

Oui, ça commence à être ardu. Je ne savais pas que les nombres premiers sont de la forme 6n±1.
Je vais m'y mettre.

Francky__laH

7 novembre 2011 à 18:39:18

[HS]PE004, oui tu peux commencer à 999*999 et décrementer, mais il te faut une condition d'arrêt réfléchie, sinon brute force sur tout, ça passe. [/HS]

Les nombres premiers hors 2 et 3 sont de la forme 6k±1. La preuve :
6k+0 : multiple de 2 et 3
6k+1 : candidat à la primalité
6k+2 : multiple de 2
6k+3 : multiple de 3
6k+4 : multiple de 2
6k+5 : candidat, de la forme 6k'-1

carrion crow

8 novembre 2011 à 21:55:33

Salut,
Je n'ai pas abandonné, c'est juste que je suis occupé.
J'ai essayé un peu quand même mais j'ai un peu de mal avec les slices et leur gestion sur deux listes.

carrion crow

15 novembre 2011 à 11:56:54

Bonjour,

J'ai eu le temps de regarder un peu ça. Et c'est pas joyeux, j'ai essayé de résoudre le problème en traçant les deux tableaux, en faisant pas par pas à la main, et en faisant moult essais. Pas vraiment propre comme méthode, ça ne me satisfait pas :colere2:

.
Surtout que j'ai fait ça avec 60 comme limite et que ça ne fonctionne pas avec le test de vitesse.

def sieveOfEratosthene5(limit):
    """Return a list of numbers from 0 to limit
    build with the sieve of Eratosthenes. Fourth method
    True if index is prime"""

    if limit < 5:
        if limit < 2: return
        yield 2
        if limit < 3: return
        yield 3
        return

    n = limit // 4    # only the half of the half

    primeM = [False] + [True]*n    # p=2*i-1
    primeP = [False] + [True]*n    # p=2*i+1

    for i in range(1, int((n**0.5)+1)):    # avoid index[0]
        if primeM[i]: # 6*i-1 is prime
            p = 6*i-1
            primeM[ p+i::p ] = [False] * ( (n-i)//p )
            primeP[ p-i::p ] = [False] * ( (n+i)//p )
        if primeP[i]: # 6*i+1 is prime
            p = 6*i+1
            primeM[ p-i+7::p ] = [False] * ( (n+i)//(p+7) )
            primeP[ p+i::p ] = [False] * ( (n-i)//p )

    yield 2
    yield 3

    for i in range(n): # un rang avant la fin en n
        if primeM[i]: yield 6*i-1
        if primeP[i]: yield 6*i+1

    # pour mieux gérer la fin
    if primeM[n] and 6*n-1 <= limit: yield 6*n-1
    if primeP[n] and 6*n+1 <= limit: yield 6*n+1

    return

Francky__laH

15 novembre 2011 à 13:36:01

Tu n'étais pas loin, tu n'aurais juste pas dû galérer avec l'indice zéro.

Voilà le corrigé, et là tu vas pouvoir faire péter le chrono. ( À toi les 10⁹ !!!)
C'est le crible le plus plus rapide à ma connaissance en python, et c'est le mien.

def prime(lim):
  # cas particuliers
  if lim<5:
    if lim<2: return
    yield 2
    if lim<3: return
    yield 3
    return

  yield 2
  yield 3
  n = (lim+1)//6
  primeP = [True]*(n+1)
  primeM = [True]*(n+1)
  primeP[0] = False
  primeM[0] = False
  for i in range((int(lim**0.5)+1)//6+1):
    if primeM[i]: # 6*i-1 est premier
      p = 6*i-1
      primeM[p+i::p] = [False]*((n-i)//p)
      primeP[p-i::p] = [False]*((n+i)//p)
    if primeP[i]: # 6*i+1 est premier
      p = 6*i+1
      primeP[p+i::p] = [False]*((n-i)//p)
      primeM[p-i::p] = [False]*((n+i)//p)
  for i in range(n):
    if primeM[i]: yield 6*i-1
    if primeP[i]: yield 6*i+1
  if primeM[n] and 6*n-1<=lim: yield 6*n-1
  if primeP[n] and 6*n+1<=lim: yield 6*n+1

Pour l'étape suivante, on va faire du crible segmenté.
On reprendra l'étape 1 (plus simple) et on découpera le tableau différemment.
Ça nous permettra d'aller plus loin, mais un peu moins vite en pratique.
À suivre...

carrion crow

15 novembre 2011 à 14:24:47

Citation : la Hache

Tu n'étais pas loin, tu n'aurais juste pas dû galérer avec l'indice zéro.

Voilà le corrigé, et là tu vas pouvoir faire péter le chrono. ( À toi les 10⁹ !!!)
C'est le crible le plus plus rapide à ma connaissance en python, et c'est le mien.
def prime(lim):
  # cas particuliers
  if lim<5:
    if lim<2: return
    yield 2
    if lim<3: return
    yield 3
    return

  yield 2
  yield 3
  n = (lim+1)//6
  primeP = [True]*(n+1)
  primeM = [True]*(n+1)
  primeP[0] = False
  primeM[0] = False
  for i in range((int(lim**0.5)+1)//6+1):
    if primeM[i]: # 6*i-1 est premier
      p = 6*i-1
      primeM[p+i::p] = [False]*((n-i)//p)
      primeP[p-i::p] = [False]*((n+i)//p)
    if primeP[i]: # 6*i+1 est premier
      p = 6*i+1
      primeP[p+i::p] = [False]*((n-i)//p)
      primeM[p-i::p] = [False]*((n+i)//p)
  for i in range(n):
    if primeM[i]: yield 6*i-1
    if primeP[i]: yield 6*i+1
  if primeM[n] and 6*n-1<=lim: yield 6*n-1
  if primeP[n] and 6*n+1<=lim: yield 6*n+1
Pour l'étape suivante, on va faire du crible segmenté.
On reprendra l'étape 1 (plus simple) et on découpera le tableau différemment.
Ça nous permettra d'aller plus loin, mais un peu moins vite en pratique.
À suivre...

Merci. J'étais quand même loin

Pourquoi tu divise ta limite par 6 ? On ne doit pas diviser par 2 pour avoir deux listes puis par 2 pour enlever les nombres pairs ?

Je fais ~9,8s à 10**8, sacrée diminution, mais toujours "memory error" à 10**9.

Francky__laH

15 novembre 2011 à 15:21:47

Citation : carrion crow

Je fais ~9,8s à 10**8, sacrée diminution, mais toujours "memory error" à 10**9.

Je fais 5s, avec un Athlon64 à 3GHz.

Je pense que tu as 4Go de mémoire (comme moi), mais tu dois avoir un OS privateur venant de Richmond, qui est possiblement goinfre en mémoire. :-°

J'utilise un système GNU/Linux (modèle grand luxe) qui prend 500Mo de mémoire (certains plus spartiates prennent moins de 100Mo), ceci explique peut-être cela.

Essaye à 8×10⁸, on verra déjà.
Le crible segmenté nous permettra d'aller bien plus loin.

carrion crow

15 novembre 2011 à 16:00:00

Citation : la Hache

Citation : carrion crow

Je fais ~9,8s à 10**8, sacrée diminution, mais toujours "memory error" à 10**9.

Je fais 5s, avec un Athlon64 à 3GHz.

Je pense que tu as 4Go de mémoire (comme moi), mais tu dois avoir un OS privateur venant de Richmond, qui est possiblement goinfre en mémoire.
J'utilise un système GNU/Linux (modèle grand luxe) qui prend 500Mo de mémoire (certains plus spartiates prennent moins de 100Mo), ceci explique peut-être cela.

Essaye à 8×10⁸, on verra déjà.
Le crible segmenté nous permettra d'aller bien plus loin.

Oui, j'ai leur dernier OS privateur en 64bits (je n'ai que 2,1Go de dispo avec ce que j'ai qui tourne)
8*10**8 ne passe pas. 4*10**8, oui, en 36,41

J'ai un petit soucis, si je veux la somme des premiers en dessous de 2 millions, je fais sum(sieveOfEratosthene5(2000000)), mais ça ne me donne pas la bonne réponse o_O

Francky__laH

15 novembre 2011 à 16:22:57

Avec ton algo peut-être, moi j'ai la bonne réponse en 90ms :

142913828922

(merci de ne pas indiquer le numéro, ni le nom de tu sais quoi)

carrion crow

15 novembre 2011 à 17:23:15

Citation : la Hache

Avec ton algo peut-être, moi j'ai la bonne réponse en 90ms :
142913828922

(merci de ne pas indiquer le numéro, ni le nom de tu sais quoi)

Oui, c'est de ma faute, j'ai inversé deux signes. Je ne comprenais pas pourquoi tes deux if primeM/primeP n'étaient pas symétriques, et en fait tu as inversé primeM/primeP par rapport au canevas et je ne l'avais pas vu.

Francky__laH

2 janvier 2012 à 22:00:25

J'ai réécrit mes cribles, avec une roue à 1, 2, 6 ou 30 dents.
(même 210 aussi, mais ça devient monstrueux), tests effectués pour la somme des premiers inférieurs à 10⁸.

erat1 : 9.25s
erat2 : 4.5s
erat6 : 2.9s
erat30 : 2.06s

from time import time

def erat1(lim=10**8):
  fin=lim
  prime = [True]*(fin+1)
  prime[0] = False
  prime[1] = False
  p = 2
  l = 4
  while l<=fin:
    if prime[p]: prime[l::p]=[False]*(1+(fin-l)//p)
    p+=1
    l+=2*p-1
  for i in range(fin+1):
    if prime[i]: yield i

def erat2(lim=10**8):
  fin = (lim-1)//2
  prime = [True]*(fin+1)
  prime[0] = False
  i = 1
  p=3
  l = 4
  while l<=fin:
    if prime[i]: prime[l::p]=[False]*(1+(fin-l)//p)
    i+=1
    p+=2
    l += 4*i
  yield 2
  for i in range(fin+1):
    if prime[i]: yield 2*i+1

def erat6(lim = 10**8):
  fin = (lim-1)//6
  prime1 = [True]*(fin+1)
  prime5 = [True]*(fin+1)
  prime1[0] = False
  i = 0
  l = 0
  while l<=fin:
    if prime1[i]:
      p = 6*i+1
      l = i*(p+1)
      prime1[l::p] = [False]*(1+(fin-l)//p)
      l += 4*i
      prime5[l::p] = [False]*(1+(fin-l)//p)

    if prime5[i]:
      p = 6*i+5
      l = i*(p+5)+4
      prime1[l::p] = [False]*(1+(fin-l)//p)
      l += 2*i+1
      prime5[l::p] = [False]*(1+(fin-l)//p)
    i += 1

  yield 2
  yield 3
  for i in range(fin):
    if prime1[i]: yield 6*i+1
    if prime5[i]: yield 6*i+5
  if prime1[fin] and 6*fin+1<=lim: yield 6*fin+1
  if prime5[fin] and 6*fin+5<=lim: yield 6*fin+5



def erat30(lim=10**8):
  # cas particuliers
  if lim<7:
    if lim<2: return
    yield 2
    if lim<3: return
    yield 3
    if lim<5: return
    yield 5
    return

  n = (lim-1)//30
  prime1 = [True]*(n+1)
  prime7 = [True]*(n+1)
  prime11 = [True]*(n+1)
  prime13 = [True]*(n+1)
  prime17 = [True]*(n+1)
  prime19 = [True]*(n+1)
  prime23 = [True]*(n+1)
  prime29 = [True]*(n+1)
  prime1[0] = False
  i = l = 0
  while l<=n:
    if prime1[i]:
      p = 30*i+1
      l = i*(p+1)
      prime1[l::p] = [False]*(1+(n-l)//p)
      l += i*6
      prime7[l::p] = [False]*(1+(n-l)//p)
      l += i*4
      prime11[l::p] = [False]*(1+(n-l)//p)
      l += i*2
      prime13[l::p] = [False]*(1+(n-l)//p)
      l += i*4
      prime17[l::p] = [False]*(1+(n-l)//p)
      l += i*2
      prime19[l::p] = [False]*(1+(n-l)//p)
      l += i*4
      prime23[l::p] = [False]*(1+(n-l)//p)
      l += i*6
      prime29[l::p] = [False]*(1+(n-l)//p)
    if prime7[i]:
      p = 30*i+7
      l = i*(p+7)+1
      prime19[l::p] = [False]*(1+(n-l)//p)
      l += i*4+1
      prime17[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime1[l::p] = [False]*(1+(n-l)//p)
      l += i*4
      prime29[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime13[l::p] = [False]*(1+(n-l)//p)
      l += i*4+1
      prime11[l::p] = [False]*(1+(n-l)//p)
      l += i*6+1
      prime23[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime7[l::p] = [False]*(1+(n-l)//p)
    if prime11[i]:
      p = 30*i+11
      l = i*(p+11)+4
      prime1[l::p] = [False]*(1+(n-l)//p)
      l += i*2
      prime23[l::p] = [False]*(1+(n-l)//p)
      l += i*4+2
      prime7[l::p] = [False]*(1+(n-l)//p)
      l += i*2
      prime29[l::p] = [False]*(1+(n-l)//p)
      l += i*4+2
      prime13[l::p] = [False]*(1+(n-l)//p)
      l += i*6+2
      prime19[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime11[l::p] = [False]*(1+(n-l)//p)
      l += i*6+2
      prime17[l::p] = [False]*(1+(n-l)//p)
    if prime13[i]:
      p = 30*i+13
      l = i*(p+13)+5
      prime19[l::p] = [False]*(1+(n-l)//p)
      l += i*4+2
      prime11[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime7[l::p] = [False]*(1+(n-l)//p)
      l += i*4+1
      prime29[l::p] = [False]*(1+(n-l)//p)
      l += i*6+3
      prime17[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime13[l::p] = [False]*(1+(n-l)//p)
      l += i*6+3
      prime1[l::p] = [False]*(1+(n-l)//p)
      l += i*4+1
      prime23[l::p] = [False]*(1+(n-l)//p)
    if prime17[i]:
      p = 30*i+17
      l = i*(p+17)+9
      prime19[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime23[l::p] = [False]*(1+(n-l)//p)
      l += i*4+3
      prime1[l::p] = [False]*(1+(n-l)//p)
      l += i*6+3
      prime13[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime17[l::p] = [False]*(1+(n-l)//p)
      l += i*6+3
      prime29[l::p] = [False]*(1+(n-l)//p)
      l += i*4+3
      prime7[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime11[l::p] = [False]*(1+(n-l)//p)
    if prime19[i]:
      p = 30*i+19
      l = i*(p+19)+12
      prime1[l::p] = [False]*(1+(n-l)//p)
      l += i*4+2
      prime17[l::p] = [False]*(1+(n-l)//p)
      l += i*6+4
      prime11[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime19[l::p] = [False]*(1+(n-l)//p)
      l += i*6+4
      prime13[l::p] = [False]*(1+(n-l)//p)
      l += i*4+2
      prime29[l::p] = [False]*(1+(n-l)//p)
      l += i*2+2
      prime7[l::p] = [False]*(1+(n-l)//p)
      l += i*4+2
      prime23[l::p] = [False]*(1+(n-l)//p)
    if prime23[i]:
      p = 30*i+23
      l = i*(p+23)+17
      prime19[l::p] = [False]*(1+(n-l)//p)
      l += i*6+5
      prime7[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime23[l::p] = [False]*(1+(n-l)//p)
      l += i*6+5
      prime11[l::p] = [False]*(1+(n-l)//p)
      l += i*4+3
      prime13[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime29[l::p] = [False]*(1+(n-l)//p)
      l += i*4+4
      prime1[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime17[l::p] = [False]*(1+(n-l)//p)
    if prime29[i]:
      p = 30*i+29
      l = i*(p+29)+28
      prime1[l::p] = [False]*(1+(n-l)//p)
      l += i*2+1
      prime29[l::p] = [False]*(1+(n-l)//p)
      l += i*6+6
      prime23[l::p] = [False]*(1+(n-l)//p)
      l += i*4+4
      prime19[l::p] = [False]*(1+(n-l)//p)
      l += i*2+2
      prime17[l::p] = [False]*(1+(n-l)//p)
      l += i*4+4
      prime13[l::p] = [False]*(1+(n-l)//p)
      l += i*2+2
      prime11[l::p] = [False]*(1+(n-l)//p)
      l += i*4+4
      prime7[l::p] = [False]*(1+(n-l)//p)
    i+=1
  yield 2
  yield 3
  yield 5
  for i in range(n):
    if prime1[i]: yield 30*i+1
    if prime7[i]: yield 30*i+7
    if prime11[i]: yield 30*i+11
    if prime13[i]: yield 30*i+13
    if prime17[i]: yield 30*i+17
    if prime19[i]: yield 30*i+19
    if prime23[i]: yield 30*i+23
    if prime29[i]: yield 30*i+29
  if prime1[n] and (30*n+1)<=lim: yield 30*n+1
  if prime7[n] and (30*n+7)<=lim: yield 30*n+7
  if prime11[n] and (30*n+11)<=lim: yield 30*n+11
  if prime13[n] and (30*n+13)<=lim: yield 30*n+13
  if prime17[n] and (30*n+17)<=lim: yield 30*n+17
  if prime19[n] and (30*n+19)<=lim: yield 30*n+19
  if prime23[n] and (30*n+23)<=lim: yield 30*n+23
  if prime29[n] and (30*n+29)<=lim: yield 30*n+29



t1 = time()
print(sum(erat1()))
print(time()-t1)

print()
t1 = time()
print(sum(erat2()))
print(time()-t1)

print()
t1 = time()
print(sum(erat6()))
print(time()-t1)

print()
t1 = time()
print(sum(erat30()))
print(time()-t1)

Pour écrire ces cribles, j'ai opté pour une méthode radicale pour éviter les erreurs, un programme qui écrit mon code, en effet si on choisit une roue à 210 dents, le code fait plus de 7000 lignes, je ne l'aurai pas écrit sans ça.

q = 3 # the size of the wheel ; 3-> 2×3×5 = 30

prem = [2, 3, 5, 7]#, 11, 13, 17, ... the first primes
prod = 1
for p in prem[:q]: prod*=p
mod = list(range(prod*2))
for p in prem[:q]:
  for i in range(prod*2):
    if i in mod and i%p==0: mod.remove(i)
lng = len(mod)//2

print("def erat", prod, "(lim=10**8):", sep="")
print("  n = (lim-1)//", prod, sep="")
for i in range(lng):
  print("  prime", mod[i], " = [1]*(n+1)", sep="")
print("  prime1[0] = 0")
print("  i = l = 0")
print("  while l<=n:")
for i in range(lng):
  p1 = mod[i]
  print("    if prime", p1, "[i]:", sep="")
  print("      p = ", prod, "*i+", p1, sep="")


  p2 = mod[i]
  a0 = p2
  b0, c0 = divmod(p1*p2, prod)
  print("      l = i*(p+", p2, ")+", b0, sep="")
  print("      prime", c0, "[l::p] = [0]*(1+(n-l)//p)", sep="")

  for j in range(i+1, i+lng):
    p2 = mod[j]
    a, a0 = p2-a0, p2
    b, c = divmod(p1*p2, prod)
    b, b0 = b-b0, b
    #print("      #l = i*(p+", p2, ")+", b, sep="")
    print("      l += i*", a, "+", b, sep="")
    print("      prime", c, "[l::p] = [0]*(1+(n-l)//p)", sep="")
    #print("      assert (", prod, "*l+", c, ")%p==0, (l, p)", sep="")
print("    i+=1")

print("  res =", prem[:q])
print("  for i in range(n):")
for i in range(lng):
  print("    if prime", mod[i], "[i]: res.append(", prod, "*i+", mod[i], ")", sep="")
for i in range(lng):
  print("  if prime", mod[i], "[n] and (", prod, "*n+", mod[i], ")<=lim: res.append(", prod, "*n+", mod[i], ")", sep="")
print("  return res")
print()
print("erat", prod, "()", sep="")

En pratique, je vais utiliser la roue à 30 dents (erat30), ça fait moins de 200 lignes, et c'est le crible le plus rapide à ma connaissance en Python.

Je vous le laisse en partage, have fun.

--
J'ai écrit le crible segmenté dont j'ai parlé, il est redoutable, mais je ne le publie pas encore.
J'attends de le peaufiner encore.

Francky__laH

20 janvier 2012 à 16:21:35

Il y a un petit bug dans la version précédente, voici une version corrigée et améliorée de 20%.
Je donne la liste des premiers inférieurs à 10⁸ en 1.6s, et ceux inférieurs à 10⁹ en 16s (et avec environ 2,3Go de RAM utilisée).

Version boostée nettement avec l'utilisation de bytearray des tableaux typés (un byte pour True ou False).

J'utilise aussi un très vilain try: while 1: [myJob] except: pass mais qui booste le chrono.
Dès qu'on obtient un premier supérieur à sqrt(lim), on crée un tableau de taille négative, ce qui provoque une erreur et on s'échappe ainsi du while 1: à la vitesse de l'éclair.

def primes(lim):
  """Renvoie la liste des nombres premiers
     inférieurs ou égaux à lim."""
  # cas particuliers
  if lim<7:
    if lim<2: return []
    if lim<3: return [2]
    if lim<5: return [2, 3]
    return [2, 3, 5]
  #  Crible
  n = (lim-1)//30
  m = n+1
  BA = bytearray
  prime1 = BA([1])*m
  prime7 = BA([1])*m
  prime11 = BA([1])*m
  prime13 = BA([1])*m
  prime17 = BA([1])*m
  prime19 = BA([1])*m
  prime23 = BA([1])*m
  prime29 = BA([1])*m
  prime1[0] = 0
  i = 0
  try:
    while 1:
      if prime1[i]:
        p = 30*i+1
        l = i*(p+1)
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*6
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*4
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*2
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*4
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*2
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*4
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*6
        prime29[l::p] = BA(1+(n-l)//p)
      if prime7[i]:
        p = 30*i+7
        l = i*(p+7)+1
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*4+1
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*4
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*4+1
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*6+1
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime7[l::p] = BA(1+(n-l)//p)
      if prime11[i]:
        p = 30*i+11
        l = i*(p+11)+4
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*2
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*2
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*6+2
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*6+2
        prime17[l::p] = BA(1+(n-l)//p)
      if prime13[i]:
        p = 30*i+13
        l = i*(p+13)+5
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*4+1
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*6+3
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*6+3
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*4+1
        prime23[l::p] = BA(1+(n-l)//p)
      if prime17[i]:
        p = 30*i+17
        l = i*(p+17)+9
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*4+3
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*6+3
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*6+3
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*4+3
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime11[l::p] = BA(1+(n-l)//p)
      if prime19[i]:
        p = 30*i+19
        l = i*(p+19)+12
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*6+4
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*6+4
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*2+2
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime23[l::p] = BA(1+(n-l)//p)
      if prime23[i]:
        p = 30*i+23
        l = i*(p+23)+17
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*6+5
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*6+5
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*4+3
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*4+4
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime17[l::p] = BA(1+(n-l)//p)
      if prime29[i]:
        p = 30*i+29
        l = i*(p+29)+28
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*6+6
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*4+4
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*2+2
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*4+4
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*2+2
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*4+4
        prime7[l::p] = BA(1+(n-l)//p)
      i+=1
  except:
    pass
  #  Génération
  RES = [2, 3, 5]
  A = RES.append
  ti=0
  try:
    for i in range(n):
      if prime1[i]: A(ti+1)
      if prime7[i]: A(ti+7)
      if prime11[i]: A(ti+11)
      if prime13[i]: A(ti+13)
      if prime17[i]: A(ti+17)
      if prime19[i]: A(ti+19)
      if prime23[i]: A(ti+23)
      if prime29[i]: A(ti+29)
      ti+=30
  except:
    pass
  if prime1[n] and (30*n+1)<=lim: A(30*n+1)
  if prime7[n] and (30*n+7)<=lim: A(30*n+7)
  if prime11[n] and (30*n+11)<=lim: A(30*n+11)
  if prime13[n] and (30*n+13)<=lim: A(30*n+13)
  if prime17[n] and (30*n+17)<=lim: A(30*n+17)
  if prime19[n] and (30*n+19)<=lim: A(30*n+19)
  if prime23[n] and (30*n+23)<=lim: A(30*n+23)
  if prime29[n] and (30*n+29)<=lim: A(30*n+29)
  return RES



from time import time
primes(10**7)
for e in range(7, 10):
  st=time()
  primes(10**e)
  print("10^", e, " --> done in ", time()-st, "s", sep="")

Une version en forme de générateur est plus pythonique, la voici pour avoir la somme, la consommation mémoire devient bien plus faible (300Mo pour lim=10⁹):

def primes(lim):
  """Renvoie la liste des nombres premiers
     inférieurs ou égaux à lim."""
  # cas particuliers
  if lim<7:
    if lim<2: return
    yield 2
    if lim<3: return
    yield 3
    if lim<5: return
    yield 5
    return
  # cas général
  yield 2
  yield 3
  yield 5
  #  Crible
  n = (lim-1)//30
  m = n+1
  BA = bytearray
  prime1 = BA([1])*m
  prime7 = BA([1])*m
  prime11 = BA([1])*m
  prime13 = BA([1])*m
  prime17 = BA([1])*m
  prime19 = BA([1])*m
  prime23 = BA([1])*m
  prime29 = BA([1])*m
  prime1[0] = 0
  i = 0
  try:
    while 1:
      if prime1[i]:
        p = 30*i+1
        l = i*(p+1)
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*6
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*4
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*2
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*4
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*2
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*4
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*6
        prime29[l::p] = BA(1+(n-l)//p)
      if prime7[i]:
        p = 30*i+7
        l = i*(p+7)+1
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*4+1
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*4
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*4+1
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*6+1
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime7[l::p] = BA(1+(n-l)//p)
      if prime11[i]:
        p = 30*i+11
        l = i*(p+11)+4
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*2
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*2
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*6+2
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*6+2
        prime17[l::p] = BA(1+(n-l)//p)
      if prime13[i]:
        p = 30*i+13
        l = i*(p+13)+5
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*4+1
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*6+3
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*6+3
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*4+1
        prime23[l::p] = BA(1+(n-l)//p)
      if prime17[i]:
        p = 30*i+17
        l = i*(p+17)+9
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*4+3
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*6+3
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*6+3
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*4+3
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime11[l::p] = BA(1+(n-l)//p)
      if prime19[i]:
        p = 30*i+19
        l = i*(p+19)+12
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*6+4
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*6+4
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*2+2
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*4+2
        prime23[l::p] = BA(1+(n-l)//p)
      if prime23[i]:
        p = 30*i+23
        l = i*(p+23)+17
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*6+5
        prime7[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*6+5
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*4+3
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*4+4
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime17[l::p] = BA(1+(n-l)//p)
      if prime29[i]:
        p = 30*i+29
        l = i*(p+29)+28
        prime1[l::p] = BA(1+(n-l)//p)
        l += i*2+1
        prime29[l::p] = BA(1+(n-l)//p)
        l += i*6+6
        prime23[l::p] = BA(1+(n-l)//p)
        l += i*4+4
        prime19[l::p] = BA(1+(n-l)//p)
        l += i*2+2
        prime17[l::p] = BA(1+(n-l)//p)
        l += i*4+4
        prime13[l::p] = BA(1+(n-l)//p)
        l += i*2+2
        prime11[l::p] = BA(1+(n-l)//p)
        l += i*4+4
        prime7[l::p] = BA(1+(n-l)//p)
      i+=1
  except:
    pass
  #  Génération
  ti=0
  try:
    for i in range(n):
      if prime1[i]: yield ti+1
      if prime7[i]: yield ti+7
      if prime11[i]: yield ti+11
      if prime13[i]: yield ti+13
      if prime17[i]: yield ti+17
      if prime19[i]: yield ti+19
      if prime23[i]: yield ti+23
      if prime29[i]: yield ti+29
      ti+=30
  except:
    pass
  if prime1[n] and (30*n+1)<=lim: yield 30*n+1
  if prime7[n] and (30*n+7)<=lim: yield 30*n+7
  if prime11[n] and (30*n+11)<=lim: yield 30*n+11
  if prime13[n] and (30*n+13)<=lim: yield 30*n+13
  if prime17[n] and (30*n+17)<=lim: yield 30*n+17
  if prime19[n] and (30*n+19)<=lim: yield 30*n+19
  if prime23[n] and (30*n+23)<=lim: yield 30*n+23
  if prime29[n] and (30*n+29)<=lim: yield 30*n+29




from time import time
sum(primes(10**7))
for e in range(7, 10):
  st=time()
  S = sum(primes(10**e))
  print("10^", e, " --> ", S, " done in ", time()-st, "s", sep="")

====
Je possède une bonne version du célèbre et cryptique crible d'Atkin, mais une bonne version (comme celle ci) du crible d'Ératosthène est très nettement supérieure.

Remarque : la première (et très longue) partie du code (le crible) est en fait la partie rapide. La génération des nombres, ensuite, à partir des booléens est longue.
Le code tronqué peut servir de test de primalité très efficace dans la situation suivante par exemple : on vous donne 100 millions de nombre inférieurs à 10⁹, lesquels sont premiers ?
C'est bien mieux que de faire 100 millions de tests de primalité même genre Miller-Rabin amélioré.

===

Il se pourrait qu'un jour je mette tout mon travail, repris depuis le début, progressif, pour en faire un tuto du style crible efficace en Python.
Si il y a des remarques , des conseils, des demandes, ... ---> envoyer moi un MP, merci.

Lord Casque Noir

20 janvier 2012 à 23:09:24

Hop, je prends ce code :

def erata(lim):
    fin = (lim-1)//2
    prime = [True]*(fin+1)
    prime[0] = False
    i = 1
    p = 3
    l = 4
    while l<=fin:
        if prime[i]:
            prime[l::p]=[False]*(1+(fin-l)//p)
        i += 1
        p += 2
        l += 4*i
    
    yield 2
    for i in xrange(fin+1):
        if prime[i]: 
            yield 2*i+1

Cette chose (prime[l::p]=[False]*(1+(fin-l)//p)), bien que rusée, n'est pas de l' "algorithmique" mais un remplacement d'une boucle python par une boucle en C.

Tant qu'à faire, autant utiliser numpy (on va prendre pour excuse d'utiliser 1 octet par nombre au lieu de 8) :

import numpy as np

def eratb(lim):
    fin = (lim-1)//2
    prime = np.ones( fin+1, dtype='b1' )
    prime[0] = False
    i = 1
    p = 3
    l = 4
    while l<=fin:
        if prime[i]:
            prime[l::p] = False
        i += 1
        p += 2
        l += 4*i
    
    return np.concatenate(( [2], prime.nonzero()[0]*2+1 ))

Pour eratb( 10**8 )
Version de base : 8.78 s
Version numpy : 2.65 s

Un crible différent, qui ne nécessite qu'une quantité de mémoire proportionnelle à la quantité de nombre premiers trouvés (avec un heap). Bien sûr, il est excessivement lent (ça ne se prête pas terriblement au python, faudrait tester en C pour voir) :

def eratq(lim):
    yield 2
    q = [(9,3)]
    for i in xrange( 3,lim,2 ):
        prime = True
        while i == q[0][0]:
            multiple,base = q[0]
            heapreplace( q, (multiple+base*2,base) )
            prime = False
            
        if prime:
            yield i
            heappush( q, (i*i, i ))

Sinon, pour éviter la partie finale du traitement, tu peux simplement stocker le nombre lui-même dans le tableau primes (mettre 0 si il n'est pas premier) et utiliser filter() pour faire le tri...

Francky__laH

21 janvier 2012 à 10:23:05

Citation : Lord Casque Noir

Cette chose (prime[l::p]=[False]*(1+(fin-l)//p)), bien que rusée, n'est pas de l' "algorithmique" mais un remplacement d'une boucle python par une boucle en C.

Tant qu'à faire, autant utiliser numpy (on va prendre pour excuse d'utiliser 1 octet par nombre au lieu de 8) :

1) On va dire que pour certains défis (SPOJ, FranceIOI, Prologin, ... je m'y débrouille bien

) on n'a pas droit d'avoir numpy. Donc, je le mets hors jeu dans mon cadre. Mais dans l'absolu, numpy est mieux. Justement numpy permet d'éviter ce rusé hacky (prime[l::p]=[False]*(1+(fin-l)//p)) qui n'est que de la compréhension de liste. Je concède volontiers que le calcul de la longueur est plus que pénible, mais on peut pas faire mieux avec Python nature.

2) L'utilisation de bytearray (voir le dernier code) permet en python nature d'utiliser un octet et non 8 par booléen stocké.

3) Pour de vraies performances de fou, c'est le crible segmenté en C, c'est plus le même jeu.

Merci bien pour les variantes sympas.

Francky__laH

1 avril 2012 à 16:53:00

Encore une version qui déchire tout, et qui reste lisible, et moins cryptique que la roue à 30 dents qui doit faire peur.
Le résultat est tout à fait comparable en vitesse.

J'ai repris juste une roue à deux dents, et pour le filtrage, je viens d'opter pour compress d'itertools qui est mieux que filter qui demande une fonction.

Le résultat est fulgurant. 29ms pour la somme des premiers inférieurs à 2 millions. Ma roue à 30 dents utilise 4 fois moins de mémoire et reste un peu plus rapide, mais cette fois le code est digeste.

from itertools import compress
def primes(lim):
  "generate odd primes"
  assert lim>2
  BA = bytearray
  n = (lim-1)//2
  prime = BA([1])*(n+1)
  prime[0] = 0 # 2*0+1 = 1 n'est pas premier
  # Crible
  for i in range((int(lim**0.5)+1)//2):
    if prime[i]:
      p = 2*i+1
      prime[p+i::p]=BA( (n-i)//p )
  return compress(range(1,lim+1,2),prime)

for i in range(3, 30):  # tests
  print(i, list(primes(i)))

st = time() ; print(2+sum(primes(2*10**6)), time()-st)

Je cherche encore a rendre l'avant dernière ligne 100% fonctionnelle en style.
Du genre avec un islice combiné avec un takewhile. (le fonctionnel, c'est cool)
Mais là n'est pas le bottleneck, m'enfin c'est celle là la plus rusée/moche/à-tunner-donc.

===

Pour la roue à 6 dents, j'ai étudié la possibilité de faire un merge avec le module heapq, mais le résultat est plus long.
Mais si l'ordre ne compte pas, alors un compress sur chaque liste et on peut obtenir la somme encore plus rapidement. On va dire que je tiens à garder l'ordre, donc la solution merge, n'est pas pour le moment satisfaisante.

===

Enjoy.

cerium50

1 avril 2012 à 17:15:09

Citation

je viens d'opter pour compress d'itertools qui est mieux que filter qui demande une fonction.

Et par rapport à une simple compréhension de liste/générateur? (si ça t'inquiète ça reste fonctionnel, cf Haskell) Voire yield-er chaque nombre premier dès que t'en trouves un?

Sinon ton code n'est pas vraiment ce que j'appelle fonctionnel (la boucle for)... Ceci dit vouloir faire du fonctionnel pour avoir des bonnes perfs, en Python tout du moins c'est une mauvaise idée: tu payes l'appel de la fonction à chaque fois comme il n'y a pas d'optimisation tail rec.

Francky__laH

1 avril 2012 à 17:47:33

Ben oui, il est pas fonctionnel, mais là j'essaye d'en mettre.

Le compress est plus élégant et bien plus rapide.

return compress(range(lim+1),prime)
#return (p for p in range(2, lim+1) if prime[p]) # plus lent !!!

==

Pour le for qui reste, je sais pas encore si je pourrais le virer.
Déjà, je m'attaque à
prime[i2::p]=BA(1+ (n-i2)//p )
qui est un for déguisé.
C'est très rapide, mais je cherche voir si on peut pas faire mieux.
À la fin, on verra si je peux virer le dernier for.
Mais là, déjà les perfs sont terribles, même si on reste encore assez loin du C.

Merci pour les remarques.

1 2 >>