San Vu Ngoc - IRMAR

Initialisations

In [1]:
%matplotlib inline
# inline can be replaced by notebook to get interactive plots
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as img
import matplotlib.cm as cm
import matplotlib.colors as colors
In [2]:
plt.rcParams['figure.figsize'] = (20.0, 12.0) # set figures display bigger

Algorithme d'Euclide pour le pgcd

On fait ici le cas $a\geq 0$ et $b\geq 0$. À adapter pour le cas général !

L'algorithme est récursif. Il suffit de programmer deux étapes, qui correspondent à l'initialisation et l'hérédité dans le principe de récurrence.

  1. l'initialisation: En fait, en programmation récursive on appelle plutôt ça la terminaison: c'est ce qui permet d'assurer que le programme ne va pas tourner indéfiniment: "Si b=0, alors le pgcd(a,b) est a"

  2. l'hérédité: Pour calculer pgcd(a,b), avec $a\geq b$, on se ramène à des entiers "plus petits" grâce à la formule pgcd(a,b) = pgcd(b,r) où r est le reste de la division euclidienne de a par b.

    Si jamais $a<b$ ce n'est pas grave car la méthode ci-dessus va appeler la fonction pgcd(b,a) ! (et ensuite, le problème ne se posera plus car $r<b$).

Important: Pour utiliser une méthode récursive, il faut démontrer que la notion de "plus petits" utilisée dans l'hérédité finit bien par aboutir à l'étape de terminaison en un nombre fini d'itérations...
Ici on sait que $r<b\leq a$. Donc $\min(b,r)\leq \min(a,b) - 1$. Donc au pire au bout de $N$ étapes avec $N=\min(a,b)$, on sait que le "reste" $r$ sera nul.

In [3]:
def pgcd (a,b):
    if b == 0:
        return a  
    else:
        return pgcd(b, a % b)
In [4]:
pgcd(7,123)
Out[4]:
1
In [5]:
pgcd(600,124)
Out[5]:
4
In [6]:
pgcd(12,0), pgcd(0,0), pgcd(0,12)
Out[6]:
(12, 0, 12)
In [7]:
pgcd(1234,2020)
Out[7]:
2

La version 'originale' d'Euclide avec seulement des soustractions

In [8]:
def pgcdEuclide (a,b):
    if b == 0:
        return a  
    elif a < b:
        return pgcdEuclide (b,a)
    else:
        return pgcdEuclide(b, a - b)
In [9]:
pgcdEuclide(600,124)
Out[9]:
4
In [10]:
pgcdEuclide(1234,2020), pgcdEuclide(0,12), pgcdEuclide(0,0)
Out[10]:
(2, 12, 0)

Version non récursive

Les algorithmes récursifs sont très pratiques lorsqu'on peut raisonner par récurrence. Mais, en python, ils ne sont pas optimisés, il demandent davantage de mémoire (mémoire de pile ou stack) et ne sont donc pas adaptés pour un grand nombre d'itérations (en gros supérieur à 1000).

(Pour le pgcd ce n'est pas très grave car l'algorithme d'Euclide est très efficace: peu d'itérations même pour des entiers très grands, cf. ci-dessous. Il faudrait des nombres avec plus de 200 chiffres pour que ça plante...)

In [11]:
pgcd(9896667081063042353430215182058733503146666634384944373401765774819916952466570536757090685909051250971877555701522243102171464782938632257729627142425602168458229597715225909882920542638682441199086679, 16013143712502213357018698163653060025603471927161902106364041769343651616369837753524801567948803360284511232110892962053473043906000050639883022618085781173913128777782344520942246774419401664791597285)
Out[11]:
1
In [12]:
pgcd(61164766314391710035884829815943265224568052927769577329622759945237346639032672167677108820397521093126336764093707189513015791230614183821533954756601790054548991800671186110593262317807192235925106064, 98966670810630423534302151820587335031466666343849443734017657748199169524665705367570906859090512509718775557015222431021714647829386322577296271424256021684582295977152259098829205426386824411990866793)

---------------------------------------------------------------------------
RecursionError                            Traceback (most recent call last)
<ipython-input-12-020b46408f83> in <module>()
----> 1 pgcd(61164766314391710035884829815943265224568052927769577329622759945237346639032672167677108820397521093126336764093707189513015791230614183821533954756601790054548991800671186110593262317807192235925106064, 98966670810630423534302151820587335031466666343849443734017657748199169524665705367570906859090512509718775557015222431021714647829386322577296271424256021684582295977152259098829205426386824411990866793)

<ipython-input-3-0703eb7d39b0> in pgcd(a, b)
      3         return a
      4     else:
----> 5         return pgcd(b, a % b)

... last 1 frames repeated, from the frame below ...

<ipython-input-3-0703eb7d39b0> in pgcd(a, b)
      3         return a
      4     else:
----> 5         return pgcd(b, a % b)

RecursionError: maximum recursion depth exceeded in comparison
In [13]:
def pgcdI (a,b):  # version "itérative"
    while b != 0:
        (a, b) = (b, a % b)
    return a
In [14]:
pgcdI(1234,2020), pgcdI(12,18), pgcdI(7,123), pgcdI(0,12)
Out[14]:
(2, 6, 1, 12)
In [15]:
pgcdI(61164766314391710035884829815943265224568052927769577329622759945237346639032672167677108820397521093126336764093707189513015791230614183821533954756601790054548991800671186110593262317807192235925106064, 98966670810630423534302151820587335031466666343849443734017657748199169524665705367570906859090512509718775557015222431021714647829386322577296271424256021684582295977152259098829205426386824411990866793)
Out[15]:
1

Images

1. pgcd

Rappel: on a tracé la semaine dernière la relation "i divise j". Maintenant on regarde le "graphe" de "pgcd(i,j)".

In [16]:
def GCD (n):
    T = np.zeros((n,n))
    for i in range(0, n):
        for j in range(0, n):
            T[i,j] = pgcdI(j,i)
    return T
In [17]:
plt.imshow(GCD(100), interpolation='none', cmap=cm.Spectral)
plt.colorbar()
#img.imsave('pgcd-log1000.png', GCD(1000), cmap=cm.Spectral)
Out[17]:
<matplotlib.colorbar.Colorbar at 0x7fc69c31d978>

2. Reste de la division de i par j

In [18]:
n=100
R = np.zeros((n,n))
for i in range(1, n):
        for j in range(0, n):
            R[i,j] = j % i + 1  # on ajoute 1 pour pouvoir faire 
                                # une échelle de couleurs logarithmique
In [19]:
plt.imshow(R, interpolation='none', cmap=cm.Spectral, norm=colors.LogNorm(1,n))
plt.colorbar()
Out[19]:
<matplotlib.colorbar.Colorbar at 0x7fc69c21c048>
In [20]:
Top = R[1:20,::]
plt.imshow(Top, interpolation='nearest', cmap=cm.Spectral, norm=colors.LogNorm(1,n))
Out[20]:
<matplotlib.image.AxesImage at 0x7fc69c1a1278>

3. Nombre d'itérations de l'algorithme d'Euclide

On compte le nombre de divisions euclidiennes

In [21]:
def pgcdIter (a,b):
    r = a % b
    if r == 0:
        return 1
    else:
        return 1 + pgcdIter(b,r)
    
def pgcdIterI (a,b):  # version "itérative"
    n = 0
    while b != 0:
        (a, b) = (b, a % b)
        n = n + 1
    return n
In [22]:
pgcdIter(100,1), pgcdIterI(100,1)
Out[22]:
(1, 1)
In [23]:
pgcdIter(1234,2020), pgcdIterI(1234,2020)
Out[23]:
(9, 9)
In [24]:
def ITER (n):
    T = np.zeros((n,n))
    for i in range(1, n-1):
        for j in range(1, n-1):
            T[i,j] = pgcdIter(j,i)
    return T
In [26]:
%matplotlib notebook
plt.rcParams['figure.figsize'] = (10.0, 6.0) # set figures display bigger
plt.imshow(ITER(100), interpolation='none', cmap=cm.gist_earth)
plt.colorbar()
Out[26]:
<matplotlib.colorbar.Colorbar at 0x7fc697ff37b8>
In [27]:
ar=ITER(1000)
%matplotlib inline
plt.rcParams['figure.figsize'] = (20.0, 12.0) # set figures display bigger
plt.imshow(ar, cmap=cm.gist_earth)
plt.colorbar()
#img.imsave('image1000.png', ar, cmap=cm.gist_earth)
Out[27]:
<matplotlib.colorbar.Colorbar at 0x7fc69c1569e8>
In [28]:
plt.imshow(ITER(1000), cmap=cm.Spectral)
plt.colorbar()
#img.imsave('image1000-2.png', ar, cmap=cm.Spectral)
Out[28]:
<matplotlib.colorbar.Colorbar at 0x7fc69508ada0>
In [29]:
Top = ar[1:200,::]
plt.imshow(Top, interpolation='nearest', cmap=cm.gist_earth)
Out[29]:
<matplotlib.image.AxesImage at 0x7fc694ff2e48>
In [30]:
pgcdIter(546,337)
Out[30]:
12
In [31]:
pgcdIter(377,610)
Out[31]:
14
In [32]:
pgcdIter(55,89)
Out[32]:
10

On peut regarder l'image avec le gimp pour voir les pixels... cf (377,610)

Suite de Fibonacci

digression ?

https://youtu.be/DxmFbdp7v9Q?t=474

In [33]:
def fibPaire (n):
    if n == 0:
        return (0,1)
    else:
        a,b = fibPaire (n-1)
        return (b,a+b)
    
def fib (n):
    a,b = fibPaire (n)
    return a
In [34]:
[fib(n) for n in range(0,17)]
Out[34]:
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987]

La complexité la "pire" est obtenue par des paires de nombres de Fibonacci consécutifs. Leur pgcd est toujours 1 (ils sont premiers entre eux). À démontrer par récurrence...

In [35]:
pgcdIter(9227465, 14930352)
Out[35]:
35

Pour tester des grands nombres, il faut une version non récursive de l'algorithme.

In [36]:
def fibI (n):
    a,b = 0,1
    while n > 0:
        (a,b) = (b,a+b)
        n = n-1
    return a
In [37]:
[fibI(n) for n in range(0,15)]
Out[37]:
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377]

On constate que les points où le nombre d'itération est le plus grand correspondent à des couples (i,j) où les i,j sont des éléments consécutifs dans la suite de Fibonacci !

In [38]:
fibI(974), fibI(974+1)
Out[38]:
(160131437125022133570186981636530600256034719271619021063640417693436516163698377535248015679488033602845112321108929620534730439060000506398830226180857811739131287777823445209422467744194016647915972857,
 259098107935652557104489133457117935287501385615468464797658075441635685688364082902818922538578546112563887878124152051556445086889386828976126497605113833423713583754975704308251673170580841059906839650)
In [39]:
float(fibI(974))
Out[39]:
1.6013143712502212e+203
In [40]:
pgcdI(fibI(974), fibI(974+1))
Out[40]:
1
In [41]:
pgcdIterI(fibI(974), fibI(974+1))
Out[41]:
974
In [42]:
fibI(972), fibI(972+1)
Out[42]:
(61164766314391710035884829815943265224568052927769577329622759945237346639032672167677108820397521093126336764093707189513015791230614183821533954756601790054548991800671186110593262317807192235925106064,
 98966670810630423534302151820587335031466666343849443734017657748199169524665705367570906859090512509718775557015222431021714647829386322577296271424256021684582295977152259098829205426386824411990866793)
In [43]:
pgcdIterI(fibI(972), fibI(972+1))
Out[43]:
972

Mauvais (i,j)

On trace les (i,j) qui sont les $(F_k,F_{k+1})$ et leurs multiples. La couleur indique la complexité de l'algorithme d'Euclide pour pgcd(i,j). On commence à $k=k0$.

In [44]:
def BAD (n,k0):
    T = np.zeros((n,n))
    F1, F2 = fibPaire(k0)   # on part les k0-ème et 
    m = 1                   # (k0+1)-ème nombres de Fibonacci.    
    while F2 < n:
        i, j = F1, F2
        while j < n:
            T[i,j] = m
            T[j,i] = m
            i, j = i+F1, j+F2
        m = m+1
        F1, F2 = F2, F1+F2   # on passe au couple de Fibonacci suivant
    return T
In [45]:
ar=BAD(100,4)
plt.imshow(ar, interpolation='none', cmap=cm.gist_earth)
plt.colorbar()
Out[45]:
<matplotlib.colorbar.Colorbar at 0x7fc694f19a90>

Le nombre d'or

Les pentes des droites de l'image ci-dessus sont données par le nombre d'or !

$\phi = \frac{1+\sqrt{5}}2$