TP2 -- correction

Eratosthène

Première version :

In [2]:
def Eratosthene(N):  
    # renvoie la liste des premiers entre 2 et N-1
    liste= [x for x in range(2,N)]
    premiers= []
    
    while len(liste) > 0 :
        p= liste[0]
        premiers.append(p)
        liste= [x for x in liste if x%p != 0]
        
    return premiers
In [4]:
print(Eratosthene(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 [5]:
%time L= Eratosthene(20000)

CPU times: user 292 ms, sys: 7.74 ms, total: 300 ms
Wall time: 529 ms

Tentative pour accélérer :

In [6]:
def Eratosthene2(N):
    liste= [x for x in range(2,N)]
    premiers= []
    
    while len(liste) > 0 :
        p= liste[0]
        premiers.append(p)
        x= p
        while x < N:
            if x in liste:
                liste.remove(x)
            x= x+p
        
    return premiers
In [8]:
%time L= Eratosthene2(20000)

CPU times: user 4.2 s, sys: 66.4 ms, total: 4.27 s
Wall time: 6.71 s

En fait c'est bien pire -- l'intention est bonne, mais la construction des listes en Python est trop lente.

Voici la version ultime :

In [9]:
def Eratosthene3(N):
    table= [True for x in range(N)]   # table = [True, True, True, ...]

    p= 2
    
    while p < N:
        # on met "False" pour tous les multiples de p
        i= p + p    # i= 2*p
        while i < N:
            table[i]= False
            i= i+p
        # on cherche le prochain nombre premier
        p= p+1
        while p < N and not(table[p]):
            p= p+1

            
    return [p for p in range(2,N) if table[p] ]
In [10]:
%time L= Eratosthene3(20000)

CPU times: user 14.6 ms, sys: 1.77 ms, total: 16.4 ms
Wall time: 26.6 ms
In [11]:
292/14.0
Out[11]:
20.857142857142858

On a accéléré 20 fois environ.

Bézout

On fait une fonction auxiliaire :

In [61]:
def bezout_aux(x, y, ux, vx, uy, vy):
    if  y == 0:
        return [x, ux, vx]
    else:
        print("On va maintenant calculer le pgcd de", x, "et", y)
        q= x//y
        r= x%y
        print("On a fait la division euclidienne", x, "=", y, "*", q, "+", r)
        print("et le reste s'exprime comme", r, "= a*",ux - q*uy , "+ b*", vx - q*vy )
        return bezout_aux(y, r, uy, vy, ux - q*uy, vx - q*vy)

dans laquelle $x, y, u_x, v_x, u_y, v_y$ sont des entiers. On suppose qu'il existe $a$ et $b$ non donnés tels que

$x = a u_x + b v_x$

$y = a u_y + b v_y$

La fonction bezout_aux(x, y, ux, vx, uy, vy) renvoie $[d, u, v]$ où $d$ est le pgcd de $x$ et $y$, et

$d = au + bv$.

Ensuite on fait:

In [62]:
def bezout(a, b):
    return bezout_aux(a, b, 1, 0, 0, 1)

Notes sur le fonctionnement de bezout_aux.

On note que si $y=0$, alors $pgcd(x,0) = x$ et $x = au_x + bv_x$ donc on renvoie $(x, u_x, v_x)$ dans ce cas.

Sinon on écrit $x= yq + r$ (division euclidienne)

Alors $pgcd(x, y) = pgcd(y, r)$ et $$r = x - yq = au_x + bv_x - q(au_y + bv_y) = a(u_x - q u_y) + b(v_x - q v_y)$$

Donc dans ce cas on renvoie bezout_aux(y, r, uy, vy, ux - quy, vx - qvy)

In [63]:
bezout(24, 15)

On va maintenant calculer le pgcd de 24 et 15
On a fait la division euclidienne 24 = 15 * 1 + 9
et le reste s'exprime comme 9 = a* 1 + b* -1
On va maintenant calculer le pgcd de 15 et 9
On a fait la division euclidienne 15 = 9 * 1 + 6
et le reste s'exprime comme 6 = a* -1 + b* 2
On va maintenant calculer le pgcd de 9 et 6
On a fait la division euclidienne 9 = 6 * 1 + 3
et le reste s'exprime comme 3 = a* 2 + b* -3
On va maintenant calculer le pgcd de 6 et 3
On a fait la division euclidienne 6 = 3 * 2 + 0
et le reste s'exprime comme 0 = a* -5 + b* 8
Out[63]:
[3, 2, -3]

Stirling

Première version :

In [40]:
def stirling(n,k):
    if n==0 or k==0:
        if n==k:
            return 1
        else:
            return 0
    return (n-1)*stirling(n-1,k) + stirling(n-1,k-1)
In [41]:
%time stirling(50, 6)

CPU times: user 8.77 s, sys: 159 ms, total: 8.93 s
Wall time: 12.7 s
Out[41]:
4272912246011833569653034858152559064299307729135569500897280000

Voyons la version avec un dictionnaire. On commence par un dictionnaire vide :

In [53]:
table_stirling= {}

Et on stocke les résultats dedans au fur et à mesure, comme ceci :

In [43]:
def stirling2(n,k):
    # petites valeurs :
    if n==0 or k==0:
        if n==k:
            return 1
        else:
            return 0
    # sinon :
    if (n,k) in table_stirling:
        return table_stirling[(n,k)]
    else:
        valeur= (n-1)*stirling2(n-1,k) + stirling2(n-1,k-1)
        table_stirling[(n,k)]= valeur
        return valeur

Une petite fonction pour afficher joliment le dictionnaire :

In [54]:
def triangle():
    ns= [n for n,k in table_stirling]
    if len(ns) == 0:
        print("<table vide>")
        return
    N= max(ns)
    if N > 50:
        print("on affiche seulement les 50 premières lignes")
        N= 50
    for n in range(N+1):
        if n not in ns:
            continue
        print(n, "|", end=" ")
        ks= [k for a,k in table_stirling if a==n]
        h= max(ks)
        for k in range(1,h+1):
            if k in ks:
                print(table_stirling[(n,k)], end=" ")
            else:
                print("?", end= " ")
        print()

Et maintenant essayons :

In [56]:
triangle()

<table vide>
In [57]:
stirling2(7,4)
Out[57]:
735
In [58]:
%time stirling2(50,6)

CPU times: user 224 µs, sys: 1 µs, total: 225 µs
Wall time: 227 µs
Out[58]:
4272912246011833569653034858152559064299307729135569500897280000