Définition

Commençons par rappeler la définition d’un nombre premier :

Un nombre premier est un entier naturel qui admet exactement deux diviseurs distincts entiers et positifs : 1 et lui-même.

Cf. Nombre_premier

Ainsi :
– 0 n’est pas premier car divisible par tout nombre non nul
– 1 n’est pas premier car il n’a qu’un seul diviseur (lui-même)
– Au delà de 2 (qui est premier), un nombre premier est forçément impaire
– Autre curiosité, tous les nombres premiers (hormis 2 et 5) se terminent tous par 1, 3, 7 ou 9

Approche naïve

La première approche « naïve » pour déterminer la liste des nombres premiers <= n, est alors :
– une fois le nombre 2 identifié comme premier,
– de parcourir tout les nombres impaires,
– et de tester leurs primalité.

def prime_numbers(num, use_numba=False): """Using a naïve method of determining prime numbers.""" res = [ 2 ] # init list with the prime number 2 for n in range(3, num, 2): # test only odd numbers if is_prime(n): res.append(n) return res

Test de primalité

L’approche « naïve » est alors de tester la primalité par divisions successives :

Pour un nombre num donné :
– si c’est 0 ou 1 ce n’est pas un nombre premier
– si c’est 2 c’est un nombre premier
– s’il est paire ce n’est pas un nombre premier
– sinon tester tous les diviseurs de 2 à num-1 :
– si l’on trouve un nombre qui divise num sans reste alors il n’est pas premier
– si on a parcouru toute la plage sans trouver un diviseur sans reste alors il est premier

def is_prime(num): """Determining whether an input {num} is prime. Return the prime number itself on None. Assertion: {num} is a natural number (positive integers).""" if (num % 2 == 0): return None for n in range(3, num, 2): if (num % n == 0): return None else: return num

On peut optimiser un peu le parcours en s’arrêtant à n/2, en effet si on n’a pas trouvé de diviseur avant n/2 on en trouvera pas au-delà. Mais la complexité restera en o(n) !

L’autre optimisation est de ne tester que les nombres de 2 à sqrt(n), en effet :

si num = pq alors soit p <= sqrt(n) soit q <= sqrt(n)

On passe ici à une complexité de o(sqrt(n)).

P.-S. on peut encore améliorer le parcours :
– On observe que tout nombre premier supérieur à 3 peut être écrit sous la forme 6k +/- 1 (avec k > 0).
– En effet tout entier peut être écrit sous la forme 6k + i (avec i = -1, 0, 1, 2, 3 ou 4).
– Comme 2 divise (6k + 0), (6k + 2) et (6k + 4) et 3 divise (6k + 3), il est plus efficient de tester les divisions par 2 et 3 puis tous les nombres de la forme 6k +/- 1. C’est 3x plus rapide !

À noter qu’en Python3 range est un itérateur [1].

Parallélisme

L’application de la méthode précédente permet d’obtenir les 5 761 455 premiers nombres premiers inférieurs à 100 000 000 en 21 minutes sur mon PC i7-8750H CPU @ 2.20GHz
Mais cet algo n’exploite qu’1 seule unité de calcul, mon PC ayant 12 cœurs l’idée est ici de distribuer les différents tests de primalité sur plusieurs cœurs.

La première approche est d’utiliser le multithreading. Pour rappel le multithreading consiste à exécuter plusieurs threads simultanément au sein d’un même processus.

La création de threads est un processus remarquablement léger, rapide et a une faible empreinte mémoire (ils partagent l’espace mémoire avec leur processus parent).

Malheureusement, en raison du GIL (ou Global Interpreter Lock) en CPython (l’implémentation la plus répendue), l’exécution du bytecode Python dans plusieurs threads est exécuté séquentiellement et non parallèlement.

Le GIL est une exclusion mutuelle (mutex) qui rend les choses plus saines au niveau des tâches permettant d’intégrer facilement des librairies externes qui ne sont pas à fil sécurisé, et exécute le code non-parallèle plus rapidement.

Il y a quelques exceptions :
– Les autres implémentations de Python tels que Jython ou IronPython ne sont pas concernés par le GIL
– Tout ce qui se passe en dehors du GIL (l’appel à des lib C via Cython par exemple)
– Les tâches qui font des I/O (ouvrir/écrire des fichiers, des requêtes http,…)
– Ou l’usage de la librairies comme NumPy.

Pour les autres cas, si on veut exécuter des tâches en parallèle en Python, il faudra se tourner vers le multitraitement.

Multiprocessing

Nous allons ici paralléliser la boucle for de test de primalité.

Paralléliser la boucle signifie répartir tous les processus en parallèle en utilisant plusieurs cœurs, chaque calcul n’attendant pas la fin du précédent.

Pour rappel un processus est une abstraction simple du système d’exploitation, c’est un ensenble d’instructions qui s’exécute indépendamment.
Ainsi, par rapport aux threads (appelé aussi processus léger), le lancement d’un nouveau processus demande l’allocations de ressources (I/O, ports, espace d’adressage) ce qui peut représenter un certain coût et notamment un temps au démarrage.

En Python le package multiprocessing prend en charge la création d’un processus enfant à la demande d’un autre processus en cours.
Assez similaire au module threading, il offre une interface quasi identique à celui-ci, mais privilégie les processus au lieu des threads.

from multiprocessing import Pool def prime_numbers_multi(num, chunksize=1000, use_numba=False): """Using the multiprocessing module to parallelize tasks.""" with Pool(PROCESSES) as p: res = p.imap(is_prime, range(3, num, 2), chunksize=chunksize) return [ 2 ] + [ n for n in res if n ]

On tombe à 6’ pour déterminer les 5 761 455 premiers nombres premiers.

Joblib

Un autre module, joblib utilise le multitraitement pour exécuter sur de multiples cœurs en parallèle un traitement dans une boucle for.
Il fournit un pipeline léger, simple et direct.

Il convient de mentionner que joblib n’est pas magique, le backend ’threading’ souffre du même goulot d’étranglement de GIL et le backend ’multiprocessing’ génère une surcharge importante en raison de la sérialisation de tous les paramètres et valeurs de retour.

import joblib def prime_numbers_joblib(num, jobs=2, use_numba=False): """Using the joblib module to parallelize tasks.""" res = joblib.Parallel(n_jobs=jobs, prefer='processes')(joblib.delayed(is_prime)(n) for n in range(3, num, 2)) return [ 2 ] + [ n for n in res if n ]

Avec la joblib, les 5 761 455 premiers nombres premiers sont déterminé en 6’ sur 10 cœurs.

Numba

Numba utilise un compilation JIT (Just-In-Time) à l’exécution : lors de la première exécution d’une fonction, Numba va générer un code optimisé pour la machine et c’est ce code qui va être utilisé lors des appels successifs de cette fonction.

Particulièrement efficace sur les fonctions, boucles et tableaux NumPy, Numba promet une exécution x100 avec seulement quelques lignes de code !

Appliqué à notre premier algo « naïf » :

import math from numba import njit, prange @njit def is_prime_jit(num): if (num % 2 == 0): return None for n in prange(3, int(math.sqrt(num))+1, 2): if (num % n == 0): return None else: return num

Il réduit le temps de calcul à 1’36" au lieu des 21’ !!!

Crible d’Eratosthène

Le crible d’Ératosthène procède par élimination : il s’agit de supprimer d’une table des entiers de 2 à N tous les multiples d’un entier (autres que lui-même).

À la fin il ne restera que les entiers qui ne sont multiples d’aucun entier à part 1 et eux-mêmes, et qui sont donc les nombres premiers.

On commence par rayer les multiples de 2, puis les multiples de 3 restants, puis les multiples de 5 restants, et ainsi de suite en rayant à chaque fois tous les multiples du plus petit entier restant.

À la fin du processus, tous les entiers qui n’ont pas été rayés sont les nombres premiers inférieurs à N.

def eratosthene(num): """Using the sieve of Eratosthene.""" numbers = [ False ] * (num+1) # => [0..num] numbers[0] = False # out of scope numbers[1] = False for n in range(2, num+1): if numbers[n]: m = 2 while True: try: numbers[n*m] = False # remove all multiples m += 1 except IndexError: break return [ n for n,f in enumerate(numbers) if f ]

Très efficace, on trouve les 5 761 455 premiers nombres premiers en seulement 41".
Cependant, au-delà cet algo requiert beaucoup de mémoire (ici 2 Gio) !

P.-S. quelques optimisations peuvent encore être faites :
– Une fois le nombre 2 identifié comme premier, on ne testera plus que les nombres impaires.
– Dans la boucle de test, une fois qu’on a atteint un nombre premier supérieur à racine de n, il n’y a plus de multiple à supprimer.
– Pour contourner le problème de mémoire, on peut créer un itérateur pour éviter de contenir en mémoire toute la liste des nombres premiers.
– On peut même utiliser un bitarray qui ne prend qu’un seul bit pour stocker chaque valeur booléenne.

Perfs

Un résumé ci-dessous des temps de calcul sur un Dell G5 5587 :
– i7-8750H CPU @ 2.20GHz x 12 ; 16Gio Mem
– Ubuntu 22.04 LTS / Python 3.10 (CPython)

méthode \ limite	1 000 000	jit	10 000 000	jit	100 000 000	jit
simple	1.67"		49"	5"	21’	1’36"
multi(1000)	0.31"				6’	2’
multi(10000)	0.28"				6’	44"
multi(100000)	0.14"		11"	2"	6’	19"
multi(1000000)	0.13"				6’	14"
joblib(10)	3"		20"	18"	6’	4’
Ératosthène	0.33"		4"		41"	5"