Pour la définition de ce qu'est un nombre premier, voir ici: http://fr.wikipedia.org/wiki/Nombre_premier
Les fonctions présentées ici renverront la liste des nombres premiers inférieurs à un nombre donné. Par exemple, pour ceux inférieurs à 100:
premiers(100) ⇒ [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]
Combien y a-t-il de nombres premiers? Une infinité. Voilà quelques exemples:
nombres de nombres premiers inférieurs ou égaux à 10**n: 10**1 4 10**2 25 10**3 168 10**4 1229 10**5 9592 10**6 78498 10**7 664579 10**8 5761455 10**9 50847534 10**10 455052511 <= ok 10**11 4118054813 10**12 37607912018 10**13 346065536839 10**14 3204941750802 10**15 29844570422669 10**16 279238341033925 10**17 2623557157654233 10**18 24739954287740860 10**19 234057667276344607 10**20 2220819602560918840
Avec le code le plus rapide parmi ceux de cette page, on peut trouver les 455052511 nombres premiers inférieurs ou égaux à 10 milliards en une heure environ.
Au delà, il faudra une grosse mémoire et beaucoup de patience…
Le principe du calcul est indiqué ici: http://fr.wikipedia.org/wiki/Crible_d'%C3%89ratosth%C3%A8ne
Commençons par un code simple et très proche du principe de base:
def eratosthene(n): """retourne la tableau des nombres premiers <= n (crible d'eratosthene)""" n += 1 tableau = [0,0] + [i for i in xrange(2, n)] for i in xrange(2, n): if tableau[i] != 0: # c'est un nombre 1er: on garde, mais on neutralise ses multiples for j in xrange(i*2, n, i): tableau[j] = 0 return [p for p in tableau if p!=0]
Le code se lit bien:
C'est déjà un code rapide, mais on peut l'améliorer:
Voilà un tel code amélioré:
def eratosthene(n): """retourne la liste des nombres premiers <= n (crible d'eratosthene)""" if n<2: return [] n += 1 tableau = [False,False] + [True]*(n-2) tableau[2::2] = [False]*((n-2)//2 + n%2) # sup. des nb pairs premiers = [2] # initialisation de la tableau des nb 1ers (2 est 1er) racine = int(n**0.5) racine = racine + [1,0][racine%2] # pour que racine soit impair for i in xrange(3, racine+1, 2): if tableau[i]: # on enregistre le nouveau nb 1er premiers.append(i) # on élimine i et ses multiples tableau[i::i] = [False]*((n-i)//i + int((n-i)%i>0)) for i in xrange(racine, n, 2): if tableau[i]: # on enregistre le nouveau nb 1er (pas de multiples à supprimer) premiers.append(i) return premiers
Cette version est de plus très rapide. Par exemple, on trouve les 664579 nombres premiers inférieurs à 10 millions en une seule seconde!!! Et les 5761455 nombres premiers inférieurs à 100 millions en une dizaine de secondes.
Mais au delà, malheureusement, cette méthode consomme beaucoup de mémoire et génère une exception “MemoryError”.
Pour contourner ce problème de mémoire, commençons par créer un itérateur pour éviter d'avoir à contenir en mémoire toute la liste des nombres premiers:
def eratosthene_iter(n): """Itérateur retourne tous les nb premiers <= n (crible d'Eratosthene)""" if n<2: pass # il n'y a aucun nb 1er en dessous de 2! else: n += 1 # pour avoir les nb 1ers <=n et pas seulement <n tableau = [False, False] + [True]*(n-2) tableau[2::2] = [False]*((n-2)//2 + n%2) # supr. des nb pairs yield 2 # 2 est un nombre premier racine = int(n**0.5) racine = racine + [1,0][racine%2] # pour que racine soit impair i, fin, pas = 3, racine+1, 2 while i<fin: # on ne traite que les nb impairs if tableau[i]: yield i # i est un nombre premier # on élimine i et ses multiples tableau[i::i] = [False]*((n-i)//i + int((n-i)%i>0)) i += pas i, fin, pas = racine, n, 2 while i<fin: # on ne traite que les nb impairs if tableau[i]: yield i # i est un nombre premier i += pas
Comme c'est un itérateur, le plus naturel est de l'utiliser dans une boucle for:
for p in eratosthene_iter(100): print p
Mais on peut, à la rigueur, l'utiliser comme suit:
gen = eratosthene_iter(100) while True: try: print gen.next() except StopIteration: break
Comme c'est un itérateur, on peut toujours retrouver la liste complète des nombres renvoyés:
premiers = [p for p in eratosthene_iter(1000)]
On peut même créer une fonction lambda:
premiers = lambda n: [p for p in eratosthene_iter(n)]
Ce code est très rapide, mais le tableau consomme encore beaucoup de mémoire.
Alors, on peut utiliser le module “bitarray' qui ne prend qu'un seul bit pour stocker chaque valeur booléenne.
On trouve ce module ici: http://pypi.python.org/pypi/bitarray. Sous Windows, il faut l'installer avec setuptools.
Voilà le nouveau et dernier code, toujours sous forme d'itérateur:
def eratosthene_ba_iter(n): """Itérateur retourne tous les nb premiers <= n (crible d'Eratosthene) on utilise ici le module 'bitarray' pour stocker les booléens """ if n<2: pass # il n'y a aucun nb 1er en dessous de 2! else: n += 1 # pour avoir les nb 1ers <=n et pas seulement <n tableau = bitarray(n) tableau.setall(True) tableau[0], tableau[1] = False, False # 1 n'est pas un nb 1er yield 2 # 2 est un nombre premier tableau[2::2] = False # on élimine tous les nombres pairs racine = int(n**0.5) racine = racine + [1,0][racine%2] # pour que racine de n soit impair i, fin, pas = 3, racine+1, 2 while i<fin: # on ne traite que les nb impairs if tableau[i]: yield i # on a trouvé un nouveau nb 1er: on le renvoie tableau[i::i] = False # on élimine i et ses multiples i += pas i, fin, pas = racine, n, 2 while i<fin: # on ne traite que les nb impairs if tableau[i]: yield i # on a trouvé un nouveau nb 1er: on le renvoie i += pas
Comme c'est un itérateur, on peut toujours retrouver la liste complète des nombres renvoyés:
premiers = [p for p in eratosthene_ba_iter(1000)]
On peut même créer une fonction lambda:
premiers = lambda n: [p for p in eratosthene_ba_iter(n)]
Cette dernière fonction est vraiment rapide et permet d'obtenir, par exemple, les 455052511 nombres premiers inférieurs à 10 milliard en une heure environ.
Il faut bien se rendre compte de ce que représentent une telle quantité de nombres: un fichier de 5,03 Go (5 403 267 048 octets)!
Le principe du calcul est assez simple. Pour le nombre n donné, on teste chacun des nombres de 2 à n pour savoir s'il est premier. Pour chaque test, le nombre est premier si on ne lui a pas trouvé de diviseur. Quand le nombre est premier, on le conserve dans une liste. A la fin, on renvoie cette liste.
Pour gagner du temps, on utilise plusieurs astuces:
- A part “2” qui est premier, on ne teste pas les nombres pairs qui ne sont forcément pas premiers.
- Pour chaque nombre k à tester, On arrête les divisions à racine de k parce que si on n'a pas encore trouvé de diviseur, on n'en trouvera plus après, et le nombre est donc premier.
- Comme on doit conserver tous les nombres premiers trouvés, on n'utilise qu'eux comme diviseur dans les tests de division!
Pour vérifier, vous pouvez voir ici pas mal de nombres premiers déjà calculés: http://nombrespremiersliste.free.fr/
Voilà le code proposé:
def premiers(n, p=[2,3,5]): """Retourne la liste des nombres premiers <= n (méthode=division)""" k = p[-1]+2 if n < k: return [x for x in p if x<=n] while k <= n: i = 0 while i < len(p): if p[i]*p[i] > k: p.append(k) break if (k % p[i]) == 0: break i += 1 k += 2 return p # Exemple d'utilisation print premiers(100) # affiche [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]
Performance: Cette fonction est très rapide. Par exemple, toujours avec l'accélérateur psyco, la fonction trouve les 5.761.455 nombres premiers inférieurs à 100 millions en une centaine de secondes. A noter que l'affichage de cette très longue liste sera plus longue que son calcul…
Avec cette fonction, on peut avoir n=1 milliard, et on trouve alors les 50.847.534 nombres premiers en à peu près 40 minutes. Mais il ne faut pas rêver: le temps de calcul grimpe très très vite, et ce n'est pas avec ça que vous allez casser le cryptage RSA: le dernier cryptage factorisé, le RSA-640, a demandé plus de 5 mois de calcul avec 30 CPU … Et on en est maintenant au RSA-2048 qui correspond à un nombre d'environ 600 chiffres!
Cas particulier: si n=1, la fonction renvoie une liste vide, parce que “1” n'est pas un nombre premier.
Vous pouvez bien entendu ajouter un test pour empêcher que la fonction ne soit appelée avec n<=0.
Vous pouvez tester cette fonction avec la Calculext ici: http://calculext.jpvweb.com, mais soyez raisonnable: au delà de premiers(100000), vous allez dépasser le temps maxi (8s) de calcul autorisé sur le serveur (qui, de plus, n'a pas psyco).