L3 Informatique MPI4 – TD8 : Nombres premiers

Premiers pas.

Écrire une fonction primes(n) retournant la liste des nombres premiers inférieurs à $n$. On pourra utiliser le crible d'Eratosthène.

Calculer primes(10000) et comparer avec ce qu'on obtient en filtrant range(10000) avec 4 tests de Miller-Rabin.

In [1]:
from TD8 import *
print primes(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]

In [2]:
pp = primes(10000)
mr = [p for p in range(2,10000) if miller_rabin(p)]
In [3]:
[q for q in mr if q not in pp]
Out[3]:
[1729, 2465]
In [4]:
miller_rabin(1729, tests=10)
Out[4]:
False

Factorisation

Force brute

Pour factoriser $n$, précalculer une liste des $M$ premiers nombres premiers, tester s'ils divisent $n$, et répéter récursivement sur le quotient. Si $p$ est le dernier de la liste et $n≤p^2$, c'est fini, sinon il faut continuer : tester les facteurs $f_k=p+2k$ jusqu'à ce que $n≤f^2_k$. On pourra écrire une fonction intermédiaire trial_division(n) qui retourne le premier diviseur rencontré. Est-il nécessairement premier ?

Force brute améliorée

On se contente de la liste [2,3,5] et on repart de 7. On a $2\times 3\times 5=30$, et si le facteur $f$ testé n'est pas congru à 1 ou à un nombre premier modulo 30, alors $f\wedge 30\not= 1$ et il ne peut pas être premier. On incrémentera donc $f$ de 4,2,4,2,4,6,2,6 (périodicité 8) de manière à ce que $f\mod 30$ soit successivement 11,13,17,19,23,29,1,7. Comparer les temps d'exécution des deux méthodes.

In [5]:
def F(n):
    return 2**(2**n)+1
F(5)
Out[5]:
4294967297
In [6]:
factor(F(5))
Out[6]:
[(641, 1), (6700417, 1)]
In [7]:
factor(F(6))
Out[7]:
[(274177, 1), (67280421310721L, 1)]
In [8]:
F(7)
Out[8]:
340282366920938463463374607431768211457L

Test de Lucas

C'est la réciproque du théorème de Fermat. On rappelle que ce dernier est un cas particulier du théorème de Lagrange : dans un groupe fini, l'ordre de tout élément est un diviseur de l'ordre du groupe. Ainsi, si $n$ est premier, ${\mathbb Z}/n{\mathbb Z}$ est un corps, et son groupe multiplicatif est d'ordre $n−1$. Mais si $n$ n'est pas premier, l'ordre de ce groupe est strictement inférieur à $n−1$. Donc, si on peut trouver un élément a d'ordre exactement n−1, n sera nécessairement premier.

Pour appliquer cette idée, on factorise $n−1$ pour avoir la liste $Q$ de ses diviseurs premiers. On calcule ensuite pour diverses valeurs de $a$ telles que $a^{n−1}\equiv 1 \mod n$, les $a^{(n−1)/q} \mod n$ pour tous les $q\in Q$. Si pour un certain a aucun n'est égal à 1, on a la preuve que n est premier. Voir ici pour un complément d'information.

In [9]:
def lucas(n, tests=4):
    pass
In [10]:
n = 2**107-1
lucas(n)
Out[10]:
True
In [11]:
def F(n):
    return 2**(2**n)+1
In [12]:
lucas(F(4))
Out[12]:
True

Test de Lucas-Lehmer.

Il permet de savoir si un nombre de Mersenne est premier. Voir pour les détails.

Pour prouver que $M_p=2^p-1$ est premier, on suppose qu'il possède un diviseur premier $q$ avec$q^2\le M_p$. Le groupe $G$ des éléments inversibles de ${\mathbb Z}_q[\sqrt{3}]$ est d'ordre au plus $q^2-1 \le 2^p-2$. Si on peut trouver un élément $a$ d'ordre supérieur à celui de ce groupe, cela prouvera qu'un tel $q$ n'existe pas.

Soient $z=2+\sqrt{3}$ et $\bar z=2-\sqrt{3}$. On a $z\bar z=1$, et la suite $s_n =z^{(2^n)}+\bar z^{(2^n)}$ vérifie $s_0=4$ et $s_n=s_{n-1}^2-2$. Si on suppose que $M_p|s_{p-2}$, on a donc $$ z^{(2^{p-2})}+\bar z^{(2^{p-2})} = kM_p$$ donc $z^{(2^{p-1})}=kM_pz^{(2^{p-2})}-1$, d'où $$z^{(2^{p-1})}\equiv -1\mod q$$ et $$\ z^{(2^{p})}\equiv 1\mod q$$ Ainsi, $z$ est d'ordre $2^p$ dans $G$. Comme l'ordre de $G$ est $\le 2^p-2$, c'est impossible.

Donc, si $M_p|s_{p-2}$, alors $M_p$ est premier.

Le plus grand nombre premier connu est en général un nombre de Mersenne.

In [13]:
    
def lehmer(p, verbose=False):
    pass
In [14]:
n=2**107-1;n
Out[14]:
162259276829213363391578010288127L
In [15]:
lucas(n,tests=8)
Out[15]:
True
In [16]:
lehmer(107)
Out[16]:
True
In [17]:
lehmer(107,verbose=True)
M(107) = (162259276829213363391578010288127) est premier

Out[17]:
True
In [18]:
pp = primes(10000)
In [19]:
mm = [p for p in pp if lehmer(p)]
In [20]:
print mm
[3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941]

In [21]:
len(pp)
Out[21]:
1229
In [21]: