L3 Informatique MPI4 – TD 8 : 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]:
def primes(n):
    if n<2: return []
    pp = [2]; mm = [m for m in range(3,n) if m%2]
    while mm:
        p = mm[0]
        pp.append(p)
        mm = [m for m in mm[1:] if m%p]
    return pp
In [2]:
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 [3]:
# Miller-Rabin

from random import randrange

def test_base(a,n):
    m = n-1; k = 0
    while not m%2:
        k+=1; m/=2
    b = pow(a,m,n)
    if b==1: return True
    for i in range(k):
        x = b
        b = pow(b,2,n)
        if b==1:
            if x != n-1: return False
            else: return True
    return False

def miller_rabin(n, tests=4):
    if n in [2,3]: return True
    if not n%2: return False
    for i in range(tests):
        a = randrange(2,n-1)
        if not test_base(a,n): return False
    return True

pp = primes(10000)
mr = [p for p in range(2,10000) if miller_rabin(p)]
In [4]:
[q for q in mr if q not in pp]
Out[4]:
[]
In [5]:
miller_rabin(1729, tests=4)
Out[5]:
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≤p2\), 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 [6]:
def trial_division(n):
    if n == 1: return 1
    if n%2 == 0: return 2
    if n%3 == 0: return 3
    if n%5 == 0: return 5
    if n<10: return n
    p = 3
    while p*p <= n:
        if n%p == 0: return p
        p +=2
    return n
    
In [8]:
def F(n):
    return 2**(2**n)+1

print F(6)
print trial_division(F(6))
18446744073709551617
274177

In [9]:
def trial_division(n, bound=None):

    if n == 1: return 1
    for p in [2, 3, 5]:
        if n%p == 0: return p
    if bound == None: bound = n
    dif = [6, 4, 2, 4, 2, 4, 6, 2]
    m = 7; i = 1
    while m <= bound and m*m <= n:
        if n%m == 0:
            return m
        m += dif[i%8]
        i += 1
    return n

def factor(n):
    if n in [-1, 0, 1]: return []
    if n < 0: n = -n
    F = []
    while n != 1:
        p = trial_division(n)
        e = 1
        n /= p
        while n%p == 0:
            e += 1; n /= p
        F.append((p,e))
    F.sort()
    return F
In [10]:
factor(F(5))
Out[10]:
[(641, 1), (6700417, 1)]
In [11]:
factor(F(6))
Out[11]:
[(274177, 1), (67280421310721L, 1)]

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.

In [12]:
def lucas(n, tests=4):
    bases = [randrange(2,n-1) for i in range(tests)]
    for a in bases:
        if pow(a,n-1,n) != 1: return False # n n'est pas permier 
    Q = [x[0] for x in factor(n-1)]
    for a in bases:
        if all([(pow(a,(n-1)/q,n) != 1) for q in Q]):
            return True # n est premier
    return 'FAIL'       # On ne peut pas conclure, il faut d'autres essais
In [13]:
n = 2**107-1
lucas(n)
Out[13]:
True
In [14]:
lucas(F(4))
Out[14]:
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 \(< 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 [15]:
def lehmer(p, verbose=False):
    assert miller_rabin(p)
    M = pow(2,p)-1
    s = 4
    for i in range(p-2):
        u,s = s, (s*s-2)%M
    if s%M == 0:
        if verbose: print 'M(%d) = (%d) est premier' % (p,M)
        return True
    else:
        if verbose: print 'M(%d) = (%d) est composé' % (p,M)
        return False
In [16]:
n=2**107-1;n
Out[16]:
162259276829213363391578010288127L
In [17]:
lucas(n,tests=8)
Out[17]:
True
In [18]:
lehmer(107)
Out[18]:
True
In [19]:
lehmer(107,verbose=True)
M(107) = (162259276829213363391578010288127) est premier

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

In [ ]: