1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import numpy as np
 
class Polynome_lagrange:
    ... 
    def __matmul__(self, other):
        # méthode permettant de réaliser le produit scalaire entre les polynômes self et other
        # (y0, ..., yn)⋅(z0, ..., zn)  = y0⋅z0 + ... + yn⋅zn
 
        # si other est un objet de la classe Polynome_lagrange
        if isinstance(other, Polynome_lagrange):
 
            # si les 2 vecteurs self et other sont définis dans la même base (l0, l1, ..., ln)
            if self.x==other.x:
 
                # calcul du produit scalaire des vecteurs self.y et other.y
                ps = self.y @ other.y
 
                # renvoie la valeur du produit scalaire
                return ps
 
            else: # sinon
 
                print("Les deux vecteurs ne sont pas définis dans la même base !")
 
        else: # sinon
 
            print("Le produit scalaire doit porter sur deux objets de la classe Polynome_lagrange !")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# création du 1er objet de la classe Polynome_lagrange
p1 = Polynome_lagrange(x=[0, 1, 2], y=[1, 2, 1])
 
# création du 2e objet de la classe Polynome_lagrange
p2 = Polynome_lagrange(x=[0, 1, 2], y=[-1, 1, -1])
 
# affiche les expressions des polynômes p1 et p2
print("p1 =", p1)
print("p2 =", p2)
 
print()
 
# affiche l'appartenance de p1 et p2 à leur espace vectoriel : R2[X], Base = [L0, L1, L2], x = [0, 1, 2]
print("p1 ∈ " + p1.espace() + ", Base = " + p1.base())
print("p2 ∈ " + p2.espace() + ", Base = " + p2.base())
 
print()
 
# calcul du produit scalaire de p1 et p2
ps = (p1 @ p2)
 
# affiche la valeur du produit scalaire de p1 et p2
print("p1⋅p2 = {0}".format(ps))

1
2
3
4
5
6
7
8
9
10
11
12
13
import numpy as np
 
class Polynome_lagrange:
    ... 
    def norme(self):
        # méthode permettant d'évaluer la norme du vecteur self défini dans la base (l0, l1, ..., ln)
        # norme((y0, ..., yn)) = (y0*y0 + ... + yn*yn)^0.5
 
        # calcul de la norme du vecteur self.y
        n = np.linalg.norm(self.y)
 
        # renvoie la valeur de la norme du vecteur
        return n

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# création de l'objet de la classe Polynome_lagrange
p = Polynome_lagrange(x=[0, 1, 2], y=[1, 1, 1])
 
# affiche l'expression du polynôme p
print("p =", p)
 
print()
 
# affiche l'appartenance de p à leur espace vectoriel : R2[X], Base = [L0, L1, L2], x = [0, 1, 2]
print("p ∈ " + p.espace() + ", Base = " + p.base())
 
print()
 
# calcul la norme de p
n = p.norme()
 
# affiche la valeur de la norme de p
print("norme(p) = {0}".format(n))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import numpy as np
 
class Polynome_lagrange:
    ... 
    def coef_corr(self, other):
        # méthode permettant d'évaluer le coefficient de corrélation entre self et other
 
        # si other est un objet de la classe Polynome_lagrange
        if isinstance(other, Polynome_lagrange):
 
            # si les 2 vecteurs self et other sont définis dans la même base (l0, l1, ..., ln)
            if self.x==other.x:
 
                # calcul du coefficient de corrélation entre self.y et other.y
                r = round(np.corrcoef(self.y, other.y)[0, 1], 12)
 
                # renvoi du coefficient de corrélation
                return r
 
            else: # sinon
 
                print("Les deux vecteurs ne sont pas définis dans la même base !")
 
        else: # sinon
 
            print("Le coefficient de corrélation doit porter sur deux objets de la classe Polynome_lagrange !")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# création du 1er objet de la classe Polynome_lagrange
p1 = Polynome_lagrange(x=[0, 1, 2], y=[1, 2, 5])
 
# création du 2e objet de la classe Polynome_lagrange
p2 = Polynome_lagrange(x=[0, 1, 2], y=[1, 3, 9])
 
# affiche les expressions des polynômes p1 et p2
print("p1 =", p1)
print("p2 =", p2)
 
print()
 
# affiche l'appartenance de p1 et p2 à leur espace vectoriel : R2[X], Base = (L0, L1, L2), x = [0, 1, 2]
print("p1 ∈ " + p1.espace() + ", Base = " + p1.base())
print("p2 ∈ " + p2.espace() + ", Base = " + p2.base())
 
print()
 
# calcul du coefficient de corrélation entre p1 et p2
r = p1.coef_corr(p2)
 
# affiche la valeur du coefficient de corrélation
print("coef_corr(p1,⋅p2) = {0}".format(r))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
import numpy as np
 
class Polynome_lagrange():
 
    def __init__(self, x, y):
        # méthode constructeur de la classe
 
        # si les deux listes ont la  même taille
        if len(x)==len(y):
 
            # on définit la liste de valeurs en x permettant d'obtenir la famille de polynômes (l0, l1, ..., ln) représentant une base dans ℝn[X]
            self.x = x
 
            # on définit le tableau de valeurs en y représentant les coordonnées du vecteur (tableau numpy)
            self.y = np.array(y)
 
        else: # sinon
 
            print("Les deux listes de valeurs x et y doivent avoir la même dimension !")
 
 
    def __str__(self):
        # permet d'afficher le polynôme sous la forme y0*L0(X) + ... + yn*Ln(X)
 
        Lx = ''
        # parcours des y
        for j in range(len(self.y)):
            Lx +=  str(self.y[j]) + "∙L" + str(j) + "(X)" + " + "
 
        return Lx[:-3]    
 
    def espace(self):
        # retourne l'espace vectoriel du polynôme self : R2[X]
 
        n = len(self.x)-1
 
        # retourne l'espace vectoriel du polynôme d'interpolation de Lagrange
        return "ℝ{0}[X]".format(n)
 
 
    def base(self):
        # retourne la base de polynômes de self : (L0, L1, L2), x = [0, 1, 2]
 
        n = len(self.x)-1
        base = "" 
 
        # parcours des indices des polynômes de lagrange Li
        for i in range(n+1):
            base += "L" + str(i) + ", "
 
        # création de la famille : (L0, L1, L2)
        base = "(" + base[:-2] + ")"
 
        # retourne la base du polynôme d'interpolation de Lagrange
        return "{0}, x = {1}".format(base, self.x)
 
 
 
    def __add__(self, other):
        # méthode permettant de redéfinir l'opérateur « + » pour 2 polynômes d'interpolation de Lagrange ou 2 vecteurs :
        # (y0, ..., yn) + (z0, ..., zn) = (y0 + z0, ..., yn + zn)
 
        # si other est un objet de la classe Polynome_lagrange
        if isinstance(other, Polynome_lagrange):
 
            # si les 2 vecteurs self et other sont définis dans la même base (l0, l1, ..., ln)
            if self.x==other.x:
 
                # addition des 2 vecteurs self.y et other.y (2 tableaux numpy)
                y = self.y + other.y
 
                # renvoie le polynômes résultat de l'addition des 2 polynômes self et other
                return Polynome_lagrange(self.x, y)
 
            else: # sinon
 
                print("Les deux vecteurs ne sont pas définis dans la même base !")
 
        else: # sinon
 
            print("Le produit scalaire doit porter sur deux objets de la classe Polynome_lagrange !")
 
 
    def __mul__(self, s):
        # méthode permettant de multiplier other par un scalaire s.
        # s*(y0, ..., yn) = (s*y0, ..., s*yn)
 
        # multiplication du vecteur self.y par le scalaire s
        y = s*self.y
 
        # renvoie le polynôme d'interpolation de Lagrange résultat de la multiplication de self par le scalaire s
        return Polynome_lagrange(self.x, y) 
 
 
    def __matmul__(self, other):
        # méthode permettant de réaliser le produit scalaire entre les polynômes self et other
        # (y0, ..., yn)⋅(z0, ..., zn)  = y0⋅z0 + ... + yn⋅zn
 
        # si other est un objet de la classe Polynome_lagrange
        if isinstance(other, Polynome_lagrange):
 
            # si les 2 vecteurs self et other sont définis dans la même base (l0, l1, ..., ln)
            if self.x==other.x:
 
                # calcul du produit scalaire des vecteurs self.y et other.y
                ps = self.y @ other.y
 
                # renvoie la valeur du produit scalaire
                return ps
 
            else: # sinon
 
                print("Les deux vecteurs ne sont pas définis dans la même base !")
 
        else: # sinon
 
            print("Le produit scalaire doit porter sur deux objets de la classe Polynome_lagrange !")
 
 
    def norme(self):
        # méthode permettant d'évaluer la norme du vecteur self défini dans la base (l0, l1, ..., ln)
        # norme((y0, ..., yn)) = (y0*y0 + ... + yn*yn)^0.5
 
        # calcul de la norme du vecteur self.y
        n = np.linalg.norm(self.y)
 
        # renvoie la valeur de la norme du vecteur
        return n
 
 
    def coef_corr(self, other):
        # méthode permettant d'évaluer le coefficient de corrélation entre self et other
 
        # si other est un objet de la classe Polynome_lagrange
        if isinstance(other, Polynome_lagrange):
 
            # si les 2 vecteurs self et other sont définis dans la même base (l0, l1, ..., ln)
            if self.x==other.x:
 
                # calcul du coefficient de corrélation entre self.y et other.y
                r = round(np.corrcoef(self.y, other.y)[0, 1], 12)
 
                # renvoi du coefficient de corrélation
                return r
 
            else: # sinon
 
                print("Les deux vecteurs ne sont pas définis dans la même base !")
 
        else: # sinon
 
            print("Le coefficient de corrélation doit porter sur deux objets de la classe Polynome_lagrange !")
 
 
    def __eq__(self, other):
        # méthode permettant de redéfinir l'opérateur « == » pour 2 polynômes d'interpolation de Lagrange
 
        # renvoie True si les bases et les listes de coordonnées des vecteurs sont identiques
        return (self.x==other.x) and np.array_equal(self.y,other.y)
 
 
    def __call__(self, X):
        # permet d'évaluer le polynôme d'interpolation de Lagrange en fonction de X et des valeurs de x et y
 
        # initialisation de la variable Lx
        Lx = 0
 
        # parcours des y
        for j in range(len(self.y)):
            lj = 1
            # parcours des x
            for i in range(len(self.x)):
                if i!=j:
                    lj *= (X - self.x[i])/(self.x[j] - self.x[i])
            # ajour de la valeur de lj*y[j] à Lx
            Lx +=self.y[j]*lj
 
        # renvoie la valeur de Lx
        return Lx
 
 
# fonction lambda permettant de changer le type d'appel pour évaluer le coefficient de corrélation
# r = p1.coef_corr(p2) -> r = coef_corr(p1, p2)
coef_corr = lambda p1,p2 : p1.coef_corr(p2)
 
print("IV-A. Évaluation du produit scalaire de deux polynômes\n")
 
# création du 1er objet de la classe Polynome_lagrange
p1 = Polynome_lagrange(x=[0, 1, 2], y=[1, 2, 1])
 
# création du 2e objet de la classe Polynome_lagrange
p2 = Polynome_lagrange(x=[0, 1, 2], y=[-1, 1, -1])
 
# affiche les expressions des polynômes p1 et p2
print("p1 =", p1)
print("p2 =", p2)
 
print()
 
# affiche l'appartenance de p1 et p2 à leur espace vectoriel : R2[X], Base = (L0, L1, L2), x = [0, 1, 2]
print("p1 ∈ " + p1.espace() + ", Base = " + p1.base())
print("p2 ∈ " + p2.espace() + ", Base = " + p2.base())
 
print()
 
# calcul du produit scalaire de p1 et p2
ps = (p1 @ p2)
 
# affiche la valeur du produit scalaire de p1 et p2
print("p1⋅p2 = {0}".format(ps))
 
print();print()
 
 
print("IV-B. Évaluation de la norme d'un polynôme\n")
 
# création de l'objet de la classe Polynome_lagrange
p = Polynome_lagrange(x=[0, 1, 2], y=[1, 1, 1])
 
# affiche l'expression du polynôme p
print("p =", p)
 
print()
 
# affiche l'appartenance de p à leur espace vectoriel : R2[X], Base = (L0, L1, L2), x = [0, 1, 2]
print("p ∈ " + p.espace() + ", Base = " + p.base())
 
print()
 
# calcul la norme de p
n = p.norme()
 
# affiche la valeur de la norme de p
print("norme(p) = {0}".format(n))
 
print();print()
 
 
print("IV-C. Évaluation du coefficient de corrélation entre deux vecteurs\n")
 
# création du 1er objet de la classe Polynome_lagrange
p1 = Polynome_lagrange(x=[0, 1, 2], y=[1, 2, 5])
 
# création du 2e objet de la classe Polynome_lagrange
p2 = Polynome_lagrange(x=[0, 1, 2], y=[1, 3, 9])
 
# affiche les expressions des polynômes p1 et p2
print("p1 =", p1)
print("p2 =", p2)
 
print()
 
# affiche l'appartenance de p1 et p2 à leur espace vectoriel : R2[X], Base = (L0, L1, L2), x = [0, 1, 2]
print("p1 ∈ " + p1.espace() + ", Base = " + p1.base())
print("p2 ∈ " + p2.espace() + ", Base = " + p2.base())
 
print()
 
# calcul du coefficient de corrélation entre p1 et p2
#r = p1.coef_corr(p2)
r = coef_corr(p1, p2)
 
# affiche la valeur du coefficient de corrélation
print("coef_corr(p1,⋅p2) = {0}".format(r))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
def partitionnement_k_means(points, k, k_means=None):
    # fonction de partitionnement en k-moyennes
    # sources : https://en.wikipedia.org/wiki/K-means_clustering#Standard_algorithm_(naive_k-means)
 
    # arguments :
    # points : liste des points ;
    # k : nombre de groupes dans le résultat du découpage ;
    # k_means : k moyennes initiales, argument optionnel.
 
    if not k_means: # si l'argument k_means est ommis
        # tirage au hasard des k points représentant les moyennes initiales
        k_means = random.sample(points, k)
 
    # convergence à False au début
    convergence = False
 
    # tant qu'il n'y a pas convergence (que les k moyennes ont changé)
    while not convergence:
 
        # initialisation de la liste des clusters
        clusters = [[] for _ in range(k)] 
 
        # parcours des points
        for xj in points:
 
            # génération de la liste des distances du point xj aux moyennes des clusters
            distances2_k_means = [distance2(xj, mi) for mi in k_means]
 
            # indice du cluster dont la moyenne est la plus proche de xj
            indice_cluster = np.argmin(distances2_k_means)
 
            # ajout de xj au cluster repéré par indice_cluster
            clusters[indice_cluster].append(xj)
 
 
        # la liste des moyennes devient la liste des moyennes précédentes
        prev_k_means = k_means
 
        # génération de la liste des nouvelles moyennes
        k_means = [centroide(ci) for ci in clusters]
 
        convergence = (k_means==prev_k_means)
 
    return clusters

1
2
3
>>> import numpy as np
>>> np.argmin([3, 2, 4, 2, 5])
1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
def partitionnement_k_means(points, k, k_means=None):
    # fonction de partitionnement en k-moyennes
    # sources : https://en.wikipedia.org/wiki/K-means_clustering#Standard_algorithm_(naive_k-means)
 
    # arguments :
    # points : liste des points ;
    # k : nombre de groupes dans le résultat du découpage ;
    # k_means : k moyennes initiales, argument optionnel.
 
    if not k_means: # si l'argument k_means est ommis
        # tirage au hasard des k points représentant les moyennes initiales
        k_means = random.sample(points, k)
 
    # convergence à False au début
    convergence = False 
 
    no_iteration = 1
 
    # tant qu'il n'y a pas convergence (que les k moyennes ont changé)
    while not convergence:
 
        print(f"Itération n°{no_iteration}\n")
        print("k-means =", k_means)
 
        affectations = ""
 
        # initialisation de la liste des clusters
        clusters = [[] for _ in range(k)]
 
        # parcours des points
        for xj in points:
 
            # génération de la liste des distances du point xj aux moyennes des clusters
            distances2_k_means = [distance2(xj, mi) for mi in k_means]
 
            # indice du cluster dont la moyenne est la plus proche de xj
            indice_cluster = np.argmin(distances2_k_means)
 
            # ajout de xj au cluster repéré par indice_cluster
            clusters[indice_cluster].append(xj)
 
            affectations += str(xj) + " → clusters[{0}]; ".format(indice_cluster)
 
        print("affectations :", affectations[0:-2])        
        print("clusters =", clusters)
        print("somme_carrés_distances =", somme_distances2_k_means(clusters)); print()
 
        # la liste des moyennes devient la liste précédente
        prev_k_means = k_means
 
        # génération de la liste des nouvelles moyennes
        k_means = [centroide(ci) for ci in clusters]
 
        convergence = (k_means==prev_k_means)
 
        no_iteration+=1
 
    return clusters

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
# liste des points 
points = [(2,10),(2,5),(8,4),(5,8),(7,5),(6,4),(1,2),(4,9)]
 
print("Liste des points :")
print(points)
 
print()
 
# nombre de groupes
k = 3
 
print("k =", k)
print()
 
# moyennes initiales
k_means = [(8,4),(5,8),(7,5)]
 
print("k moyennes initiales :")
print(k_means); print() 
 
# partitionnement des n points en k groupes : algorithme des k-means
resultat_k_means = partitionnement_k_means(points, k, k_means)
 
print()
 
print("Résultat :")
print(resultat_k_means)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
# liste des points 
points = [(2,10),(5,0),(8,4),(5,8),(7,5),(6,4),(1,2),(4,9)]
 
print("Liste des points :")
print(points)
 
print()
 
# nombre de groupes
k = 3
 
print("k =", k)
print()
 
# variable représentant la somme des distances des points à la moyenne de leur cluster
somme_distances2_mini = 10e10
 
# nombre de tirages aléatoires des k moyennes initiales : nombre d'exécutions de la fonction de partitionnement
nb_tirages = 5
 
print("nombre de tirages =", nb_tirages)
print()
 
# initialisation du générateur de nombres aléatoires
random.seed()
 
# on réitère le partitionnement pour minimiser la somme des distances des points à la moyenne de leur cluster
for no_tirage in range(1, nb_tirages+1):
 
    print("\n---- Tirage n°{0} ----\n".format(no_tirage))
 
    # partitionnement des n points en k groupes : algorithme des k-means
    clusters = partitionnement_k_means(points, k)
 
    # calcul de la somme des carrés des distances de chacun des points à la moyenne de leur cluster
    somme_distances2 = somme_distances2_k_means(clusters)
 
    # si la somme des distances obtenue est inférieure à somme_distances2_mini
    if somme_distances2<somme_distances2_mini:
        # mémorisation des clusters et de la somme minimum des carrés des distances
        resultat_k_means, somme_distances2_mini = clusters, somme_distances2
        print("<<<-- Résultat retenu -->>>\n")
 
print()
 
print("Résultat :")
print(resultat_k_means)
 
print()
 
print("Somme des carrés des distances :")
print(somme_distances2_mini)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
import random
import numpy as np
 
def centroide(points):
    # renvoie le centroïde des points
 
    x = [p[0] for p in points]
    y = [p[1] for p in points]
 
    return (sum(x)/len(x), sum(y)/len(y))
 
 
def distance2(x1, x2):
    # calcul du carré de la distance entre x1 et x2
 
    return pow((x1[0] - x2[0]),2) + pow((x1[1] - x2[1]),2)
 
 
def somme_distances2_k_means(clusters):
    # somme des carrés des distances des points à la moyenne de leur cluster
    # fonction à minimiser
 
    # initialisation de la somme des distances
    somme_distances2 = 0
 
    # parcours des clusters
    for ci in clusters:
        # calcul du centroïde de ci
        mi = centroide(ci)
 
        # ajout de la somme des carrés des distances de chaque point à la moyenne de son cluster 
        somme_distances2 +=  sum([distance2(xj, mi) for xj in ci])
 
    # retourne la somme des carrés des distances        
    return somme_distances2
 
 
def partitionnement_k_means(points, k, k_means=None):
    # fonction de partitionnement en k-moyennes
    # sources : https://en.wikipedia.org/wiki/K-means_clustering#Standard_algorithm_(naive_k-means)
 
    # arguments :
    # points : liste des points ;
    # k : nombre de groupes dans le résultat du découpage ;
    # k_means : k moyennes initiales, argument optionnel.
 
    if not k_means: # si l'argument k_means est ommis
        # tirage au hasard des k points représentant les moyennes initiales
        k_means = random.sample(points, k)
 
    # convergence à False au début
    convergence = False
 
    no_iteration = 1
 
    # tant qu'il n'y a pas convergence (que les k moyennes ont changé)
    while not convergence:
 
        print(f"Itération n°{no_iteration}\n")
        print("k-means =", k_means)
 
        affectations = ""
 
        # initialisation de la liste des clusters
        clusters = [[] for _ in range(k)]
 
        # parcours des points
        for xj in points:
 
            # génération de la liste des distances du point xj aux moyennes des clusters
            distances2_k_means = [distance2(xj, mi) for mi in k_means]
 
            # indice du cluster dont la moyenne est la plus proche de xj
            indice_cluster = np.argmin(distances2_k_means)
 
            # ajout de xj au cluster repéré par indice_cluster
            clusters[indice_cluster].append(xj)
 
            affectations += str(xj) + " → clusters[{0}]; ".format(indice_cluster)
 
        print("affectations :", affectations[0:-2])        
        print("clusters =", clusters)
        print("somme_carrés_distances =", somme_distances2_k_means(clusters)); print()
 
        # la liste des moyennes devient la liste précédente
        prev_k_means = k_means
 
        # génération de la liste des nouvelles moyennes
        k_means = [centroide(ci) for ci in clusters]
 
        convergence = (k_means==prev_k_means)
 
        no_iteration+=1
 
    return clusters 
 
print("--- Algorithme des k-moyennes ---\n")
 
print("I. Partitionnement en k-moyennes : choix de k moyennes initiales :\n")
 
# liste des points 
points = [(2,10),(2,5),(8,4),(5,8),(7,5),(6,4),(1,2),(4,9)]
 
print("Liste des points :")
print(points)
 
print()
 
# nombre de groupes
k = 3
 
print("k =", k)
print()
 
# moyennes initiales
k_means = [(8,4),(5,8),(7,5)]
 
print("k moyennes initiales :")
print(k_means); print() 
 
# partitionnement des n points en k groupes : algorithme des k-means
resultat_k_means = partitionnement_k_means(points, k, k_means)
 
print()
 
print("Résultat :")
print(resultat_k_means)
 
print();print()
 
wait = input("Appuyer su la touche Entrée pour continuer..")
 
print();print()
 
print("II. Partitionnement en k-moyennes : tirage au hasard des k moyennes initiales :\n")
 
# liste des points 
points = [(2,10),(5,0),(8,4),(5,8),(7,5),(6,4),(1,2),(4,9)]
 
print("Liste des points :")
print(points)
 
print()
 
# nombre de groupes
k = 3
 
print("k =", k)
print()
 
# variable représentant la somme des distances des points à la moyenne de leur cluster
somme_distances2_mini = 10e10
 
# nombre de tirages aléatoires des k moyennes initiales : nombre d'exécutions de la fonction de partitionnement
nb_tirages = 5
 
print("nombre de tirages =", nb_tirages)
print()
 
# initialisation du générateur de nombres aléatoires
random.seed()
 
# on réitère le partitionnement pour minimiser la somme des distances des points à la moyenne de leur cluster
for no_tirage in range(1, nb_tirages+1):
 
    print("\n---- Tirage n°{0} ----\n".format(no_tirage))
 
    # partitionnement des n points en k groupes : algorithme des k-means
    clusters = partitionnement_k_means(points, k)
 
    # calcul de la somme des carrés des distances de chacun des points à la moyenne de leur cluster
    somme_distances2 = somme_distances2_k_means(clusters)
 
    # si la somme des distances obtenue est inférieure à somme_distances2_mini
    if somme_distances2<somme_distances2_mini:
        # mémorisation des clusters et de la somme minimum des carrés des distances
        resultat_k_means, somme_distances2_mini = clusters, somme_distances2
        print("<<<-- Résultat retenu -->>>\n")
 
print()
 
print("Résultat :")
print(resultat_k_means)
 
print()
 
print("Somme des carrés des distances :")
print(somme_distances2_mini)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def proportion_premiers(x1, xn, k):
    # x1 : borne inférieure de l'intervalle des nombres entiers
    # xn : borne supérieure de l'intervalle des entiers
    # k : nombre de valeurs tirées au hasard dans l'intervalle [x1, xn] 
 
    # séquence de nombres entiers compris dans l'intervalle [x1, xn]
    nombres_entiers = range(x1, xn + 1)
 
    # initialisation du compteur de nombres premiers
    qte_premiers = 0
 
    # tirage au hasard de k entiers parmi la séquence de n entiers
    for j in range(k):
 
        # choix d'un nombre entier au hasard parmi la liste de nombres entiers
        N = random.choice(nombres_entiers)
 
        # si N est premier
        if test_primalite(N):
            # incrémentation du compteur de nombres premiers
            qte_premiers+=1
 
    # proportion de nombre premiers choisis au hasard dans l'intervalle [x1, xn]
    f = qte_premiers/k
 
    # renvoi de la proportion de nombres premiers compris dans l'intervalle
    return f
 
 
def quantite_premiers(x1, xn, k):
    # x1 : borne inférieure de l'intervalle des nombres entiers
    # xn : borne supérieure de l'intervalle des entiers
    # k : nombre de valeurs tirées au hasard dans l'intervalle [x1, xn]
 
    return proportion_premiers(x1, xn, k)*(xn-x1+1)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
# nombre d'entiers par intervalle
n = 4000000
 
# séquence des x avec un pas de longueur n
x_values = range(2000000, 14000001, n)
 
# parcours des x avec un pas de longueur n
for x in x_values:
 
    # proportion de nombres premiers dans cet intervalle suivant la fonction Li(x)
    p = (Li(x+n-1)-Li(x-1))/n
 
    print("Intervalle des valeurs : " + str([x, x+n-1]))
    print("Proportion estimée avec Li(x) : " + str(round(p,7)))
    print("Qté estimée avec Li(x) : " + str(round(p*n)))
 
    print()
 
    # initialisation de la liste des résultats
    results = []
 
    # parcours des exposants de 10 : 10^3, .., 10^5
    for i in range(3, 6):
        k = 10**i
 
        # évaluation de la proportion de nombres premiers compris dans [x, x+n-1] : méthode de Monte-Carlo
        f = round(proportion_premiers(x, x+n-1, k),7)
 
        # bornes inférieure et supérieure de l'intervalle de confiance
        borne_inf, borne_sup = round(f - 1/math.sqrt(k),7), round(f + 1/math.sqrt(k),7)
 
	# ligne de résultats : [k, f, [borne_inf, borne_sup], n*f, [n*borne_inf, n*borne_sup]]
        result = [k, f, str([borne_inf, borne_sup]), round(n*f), str([round(n*borne_inf), round(n*borne_sup)])]
 
        # ajout du résultat à la liste
        results.append(result)
 
    # création du tableau numpy à 2 dimensions contennant la liste de résultats
    tb_valeurs = np.array(results)
 
    # création du dataframe
    df = pd.DataFrame(tb_valeurs, columns = ['  k', '  Proportion observée avec MC', '  Intervalle de confiance', '  Qté estimée avec MC', '  Intervalle de confiance'])
 
    # affiche le dataframe
    print(df)
 
    print("------------------------------------------------------------------------------------------------------------------")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
mport matplotlib.pyplot as plt
 
...
 
# valeur du pas dans la séquence de nombres entiers
n = 200000
 
# séquence des x avec un pas de longueur n
x_values = range(n, 40*n+1, n)
 
# liste des quantités de nombres premiers estimées à l'aide de la fonction Li(x) 
Li_values =[Li(x) for x in x_values]
 
# nombre de tirages par intervalle 
k = 1000
 
# initialisation de la liste des qtés de nombres premiers : méthode de Monte-Carlo
y_values = []
 
# initialisation de la quantité de nombres premiers
qte_premiers = 0
 
# borne inférieure de l'intervalle [x1, x]
x1=1 
 
# parcours des x avec un pas de longueur n
for x in x_values:
 
    # évaluation de la quantité de nombres premiers inférieurs ou égaux à x  : méthode de Monte-Carlo
    qte_premiers += quantite_premiers(x1, x, k)
 
    # ajout de la quantité à la liste y_values
    y_values.append(qte_premiers)
 
    x1 = x+1 # borne inférieure prochain intervalle
 
# lissage des points
# y_values = savgol_filter(y_values, len(y_values)//2+1, 3)
 
plt.plot(x_values, Li_values, label = "Li(x)")
plt.plot(x_values, y_values, label = "Quantités évaluées avec MC")
 
plt.ylim(0)
 
plt.legend()
 
# titre du graphique 
plt.title('Qtés de nombres premiers inférieurs ou égaux à x')
 
# affichage du graphique 
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
from interface_graphique import Interface_graphique
 
...
 
# valeur du pas dans la séquence de nombres entiers
n = 200000
 
# nombre de tirages par intervalle 
k = 1000
 
# création de l'objet Interface_graphique avec les paramètres : valeur_mini, valeur_maxi, valeur_pas, nb_tirages, Li, quantite_premiers
interface_graphique = Interface_graphique(n, 40*n, n, k, Li, quantite_premiers)
 
# centre la fenêtre
interface_graphique.eval('tk::PlaceWindow . center')
 
interface_graphique.mainloop()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
import math
import random
import numpy as np
import scipy.integrate as integrate
import pandas as pd
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
from interface_graphique import Interface_graphique
 
 
def Li(x):
    s = integrate.quad(lambda t: 1/math.log(t), 2, x)
    return s[0]
 
def proportion_premiers(x1, xn, k):
    # x1 : borne inférieure de l'intervalle des nombres entiers
    # xn : borne supérieure de l'intervalle des entiers
    # k : nombre de valeurs tirées au hasard dans l'intervalle [x1, xn] 
 
    # séquence de nombres entiers compris dans l'intervalle [x1, xn]
    nombres_entiers = range(x1, xn + 1)
 
    # initialisation du compteur de nombres premiers
    qte_premiers = 0
 
    # tirage au hasard de k entiers parmi la séquence de n entiers
    for j in range(k):
 
        # choix d'un nombre entier au hasard parmi la liste de nombres entiers
        N = random.choice(nombres_entiers)
 
        # si N est premier
        if test_primalite(N):
            # incrémentation du compteur de nombres premiers
            qte_premiers+=1
 
    # proportion de nombre premiers choisis au hasard dans l'intervalle [x1, xn]
    f = qte_premiers/k
 
    # renvoi de la proportion de nombres premiers compris dans l'intervalle
    return f
 
 
def quantite_premiers(x1, xn, k):
    # x1 : borne inférieure de l'intervalle des nombres entiers
    # xn : borne supérieure de l'intervalle des entiers
    # k : nombre de valeurs tirées au hasard dans l'intervalle [x1, xn]
 
    return proportion_premiers(x1, xn, k)*(xn-x1+1)
 
 
def test_primalite(N):
    # fonction permettant de tester si N est premier en utilisant la méthode naïve
    # N : nombre entier à tester
 
    if N==1: return False # 1 n'est pas un nombre premier
 
    # détermination du dernier entier m à tester comme diviseur de N : m est égal à la partie entière de √𝑁
    m = int(math.sqrt(N))
 
    # parcours des entiers : 2 -> m
    for i in range(2, m+1):
        # Si i divise N
        if N % i == 0:
            # N possède un diviseur donc n'est pas premier
            # renvoie le tuple False => est_premier
            return False
 
    # N n'a pas de diviseur autre que 1 et lui-même, il est donc premier.
    # renvoie True => est_premier
    return True
 
 
print("I. Méthode de Monte-Carlo : évaluation des proportions et des qtés de nombres premiers compris dans des intervalles successifs ..\n")
 
pd.set_option('display.max_columns', 500)
pd.set_option('display.width', 1000)
 
# nombre d'entiers par intervalle
n = 4000000
 
# séquence des x avec un pas de longueur n
x_values = range(2000000, 14000001, n)
 
# parcours des x avec un pas de longueur n
for x in x_values:
 
    # proportion de nombres premiers dans cet intervalle suivant la fonction Li(x)
    p = (Li(x+n-1)-Li(x-1))/n
 
    print("Intervalle des valeurs : " + str([x, x+n-1]))
    print("Proportion estimée avec Li(x) : " + str(round(p,7)))
    print("Qté estimée avec Li(x) : " + str(round(p*n)))
 
    print()
 
    # initialisation de la liste des résultats
    results = []
 
    # parcours des exposants de 10 : 10^3, .., 10^5
    for i in range(3, 6):
        k = 10**i
 
        # évaluation de la proportion de nombres premiers compris dans [x, x+n-1] : méthode de Monte-Carlo
        f = round(proportion_premiers(x, x+n-1, k),7)
 
        # bornes inférieure et supérieure de l'intervalle de confiance
        borne_inf, borne_sup = round(f - 1/math.sqrt(k),7), round(f + 1/math.sqrt(k),7)
 
	# ligne de résultats : [k, f, [borne_inf, borne_sup], n*f, [n*borne_inf, n*borne_sup]]
        result = [k, f, str([borne_inf, borne_sup]), round(n*f), str([round(n*borne_inf), round(n*borne_sup)])]
 
        # ajout du résultat à la liste
        results.append(result)
 
    # création du tableau numpy à 2 dimensions contennant la liste de résultats
    tb_valeurs = np.array(results)
 
    # création du dataframe
    df = pd.DataFrame(tb_valeurs, columns = ['  k', '  Proportion observée avec MC', '  Intervalle de confiance', '  Qté estimée avec MC', '  Intervalle de confiance'])
 
    # affiche le dataframe
    print(df)
 
    print("------------------------------------------------------------------------------------------------------------------")
 
print()
 
print("II. Courbes représentatives des qtés de nombres premiers inférieurs où égaux à x : valeurs de Li(x) / qtés estimées avec MC\n")
 
print("II-A. Génération du graphique avec la librairie matplotlib ..\n")
 
# valeur du pas dans la séquence de nombres entiers
n = 200000
 
# séquence des x avec un pas de longueur n
x_values = range(n, 40*n+1, n)
 
# liste des quantités de nombres premiers estimées à l'aide de la fonction Li(x) 
Li_values =[Li(x) for x in x_values]
 
# nombre de tirages par intervalle 
k = 1000
 
# initialisation de la liste des qtés de nombres premiers : méthode de Monte-Carlo
y_values = []
 
# initialisation de la quantité de nombres premiers
qte_premiers = 0
 
# borne inférieure de l'intervalle [x1, x]
x1=1 
 
# parcours des x avec un pas de longueur n
for x in x_values:
 
    # évaluation de la quantité de nombres premiers inférieurs ou égaux à x  : méthode de Monte-Carlo
    qte_premiers += quantite_premiers(x1, x, k)
 
    # ajout de la quantité à la liste y_values
    y_values.append(qte_premiers)
 
    x1 = x+1 # borne inférieure prochain intervalle
 
# lissage des points
# y_values = savgol_filter(y_values, len(y_values)//2+1, 3)
 
plt.plot(x_values, Li_values, label = "Li(x)")
plt.plot(x_values, y_values, label = "Quantités évaluées avec MC")
 
plt.ylim(0)
 
plt.legend()
 
# titre du graphique 
plt.title('Qtés de nombres premiers inférieurs ou égaux à x')
 
# affichage du graphique 
plt.show()
 
 
print("II-B. Paramétrage du graphique sur une interface Tkinter ..") 
 
# création de l'objet Interface_graphique avec les paramètres : valeur_mini, valeur_maxi, valeur_pas, nb_tirages, Li, quantite_premiers
interface_graphique = Interface_graphique(n, 40*n, n, k, Li, quantite_premiers)
 
# centre la fenêtre
interface_graphique.eval('tk::PlaceWindow . center')
 
interface_graphique.mainloop()

1
2
3
4
5
6
7
8
9
10
Fonction testPrimaliteFermat(N)
    choisir un entier positif a < N au hasard   

    Si aN-1 ≡ 1 mod N alors
        renvoyer oui (N est probablement premier)
    Sinon
        renvoyer non (N est composé)
    Fin Si 

Fin Fonction

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Fonction testPrimaliteFermatItere(N)
     
    pour i de 1 jusqu'à k

	choisir un entier positif ai < N au hasard

        Si aiN-1 ≢ 1 mod N alors
            renvoyer non (N est composé : sortie de la fonction)
        Fin Si 

    Fin Pour

    renvoyer oui (N est probablement premier : sortie de la fonction)

Fin fonction

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def test_primalite(N):
    # fonction permettant de tester si N est premier en utilisant la méthode naïve
    # N : nombre entier à tester
 
    # détermination du dernier entier m à tester comme diviseur de N : m est égal à la partie entière de √𝑁
    m = int(math.sqrt(N))
 
    # parcours des entiers : 2 -> m
    for i in range(2, m+1):
        # Si i divise N
        if N % i == 0:
            # N possède un diviseur donc n'est pas premier
            # renvoie le tuple : (False, i-1) => (est_premier, nbre_iterations)
            return (False, i-1)
 
    # N n'a pas de diviseur autre que 1 et lui-même, il est donc premier.
    # renvoie le tuple : (True, m-1) => (est_premier, nbre_iterations)
    return (True, m-1)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
print("Test de primalité classique :\n")
 
# valeur de N
N = 99999980021
 
print("N = " + str(N))
 
print()
 
# Teste si le nombre entier N est premier par la méthode naïve.
est_premier, nbre_iterations = test_primalite(N)
 
# Affiche le résultat
if est_premier:
    print("Résultat du test : {0} est premier.".format(N))
else:
    print("Résultat du test : {0} est composé.".format(N))
 
print("------------------------------------------")
print("Nombre de division(s) : " + str(nbre_iterations))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def test_primalite_fermat(N):
    # fonction permettant de tester si N est premier avec un entier positif a < N choisi au hasard
    # N : nombre entier à tester
 
    # séquence de nombres entiers compris entre 2 et N-1 : 2 -> N-1
    nombres_entiers = range(2, N)
 
    # choix au hasard de l'entier positif a < N
    a = random.choice(nombres_entiers)
 
    # Si a^(N-1) mod N = 1 : a^(N-1) est congru à 1 modulo N
    if pow(a, N-1, N) == 1:
        # N est probablement premier
        return True
 
    else: # sinon : N est composé
        return False

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
def test_primalite_fermat_iter(N, k, n=None):
    # fonction permettant de tester si N est premier en itérant au plus k fois le test de Fermat
    # N : nombre entier à tester
    # k : nombre maxi. d'itérations ou d'entiers positifs a1, ... ak < N choisis au hasard
    # n : valeur maxi de la séquence de nombre entiers inférieurs à N : 2 -> n, si n est omis ou si n ≥ N, alors n = N-1
 
    # si l'argument n est omis ou si n ≥ N
    if (not n) or (n >= N) : n = N-1
 
    # séquence de nombres entiers compris entre 2 et n : 2 -> n, avec n < N
    nombres_entiers = range(2, n+1)
 
    # on itère k fois le test de Fermat : 0 -> k-1
    for i in range(k):
 
        # choix de l'entier positif ai < N au hasard
        ai = random.choice(nombres_entiers)
 
        # Si ai^(N-1) mod N ≠ 1 : ai^(N-1) n'est pas congru à 1 modulo N
        #if not test_primalite_fermat(N):
        if pow(ai, N-1, N) != 1:
            # N n'est pas premier :
            # renvoie le tuple : (False, i+1) => (est_premier, nbre_iterations)
            return (False, i+1)
 
    # N est probablement premier
    # renvoie le tuple : (True, k) => (est_premier, nbre_iterations)
    return (True, k)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
print("Test de primalité de Fermat :\n")
 
# valeur de N
N = 99999980021
 
print("N = " + str(N))
 
# paramètre de fiabilité k : nombres d'entiers positifs a1, ..., ak < N choisis au hasard 
k = 20
 
print("N = " + str(N)) 
print("k = " + str(k))
 
print()
 
# test de primalité de Fermat pour un nombre entier N et au plus k itérations
est_premier, nbre_iterations = test_primalite_fermat_iter(N, k, n=1000)
 
# Affiche le résultat
if est_premier: # si N est premier
    print("Résultat du test : {0} est probablement premier.".format(N))
else: # sinon
    print("Résultat du test : {0} est composé.".format(N))
 
print("------------------------------------------------------------------------------------------")
print("Nombre de tirage(s) : " + str(nbre_iterations))
print("Nombre d'exponentiation(s) modulaire(s) : " + str(nbre_iterations))
if est_premier:
    # risque d'erreur maxi si ce n'est pas un nombre de Carmichael : t = 1/(2^k)
    t = pow(2,-k)
    print("Risque d'erreur si N n'est pas un nombre de Carmichael inférieur à : " + str(t))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
print("Mesure du temps d'exécution des deux fonctions de test sur une série de 20 000 nombres entiers ..\n")
 
# initialisation des compteurs de tests positifs
cpt1 = 0; cpt2 = 0
 
print("Version n°1 - Test de primalité classique :\n")
 
# heure Unix de début (exprimée en secondes)
start = time.time()
 
# parcours des 20 0000 nombres entiers : 99999980000 -> 99999999999
for N in range(99999980000, 100000000000):
 
    # test classique pour N
    est_premier, nbre_iterations = test_primalite(N)
 
    if est_premier: # si le test est positif
        cpt1+=1 # incrémentation du compteur
 
# heure Unix de fin (exprimée en secondes)
end = time.time()
 
# Affiche la durée d'exécution.
print(f"Durée d'exécution : {end - start} sec.\n")
 
 
print("Version n°2 - Test de primalité de Fermat\n")
 
# heure Unix de début (exprimée en secondes)
start = time.time()
 
# parcours des 20 0000 nombres entiers : 99999980000 -> 99999999999
for N in range(99999980000, 100000000000):
 
    # test de primalité de Fermat pour N et k
    est_premier, nbre_iterations = test_primalite_fermat_iter(N, k=10, n=1000)
 
    if est_premier: # si le test est positif
        cpt2+=1 # incrémentation du compteur
 
# heure Unix de fin (exprimée en secondes)
end = time.time()
 
# Affiche la durée d'exécution.
print(f"Durée d'exécution : {end - start} sec.")
print("Taux d'erreur : " + str((cpt2-cpt1)/20000))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
import math
import random
import time
 
def test_primalite(N):
    # fonction permettant de tester si N est premier en utilisant la méthode naïve
    # N : nombre entier à tester
 
    # détermination du dernier entier m à tester comme diviseur de N : m est égal à la partie entière de √𝑁
    m = int(math.sqrt(N))
 
    # parcours des entiers : 2 -> m
    for i in range(2, m+1):
        # Si i divise N
        if N % i == 0:
            # N possède un diviseur donc n'est pas premier
            # renvoie le tuple : (False, i-1) => (est_premier, nbre_iterations)
            return (False, i-1)
 
    # N n'a pas de diviseur autre que 1 et lui-même, il est donc premier.
    # renvoie le tuple : (True, m-1) => (est_premier, nbre_iterations)
    return (True, m-1)
 
 
def test_primalite_fermat(N):
    # fonction permettant de tester si N est premier avec un entier positif a < N choisi au hasard
    # N : nombre entier à tester
 
    # séquence de nombres entiers compris entre 2 et N-1 : 2 -> N-1
    nombres_entiers = range(2, N)
 
    # choix de l'entier positif a < N au hasard
    a = random.choice(nombres_entiers)
 
    # Si a^(N-1) mod N = 1 : a^(N-1) est congru à 1 modulo N
    if pow(a, N-1, N) == 1:
        # N est probablement premier
        return True
 
    else: # sinon : N est composé
        return False
 
 
def test_primalite_fermat_iter(N, k, n=None):
    # fonction permettant de tester si N est premier en itérant au plus k fois le test de Fermat
    # N : nombre entier à tester
    # k : nombre maxi. d'itérations ou d'entiers positifs a1, ... ak < N choisis au hasard
    # n : valeur maxi de la séquence de nombre entiers inférieurs à N : 2 -> n, si n est omis ou si n ≥ N, alors n = N-1
 
    # si l'argument n est omis ou si n ≥ N
    if (not n) or (n >= N) : n = N-1
 
    # séquence de nombres entiers compris entre 2 et n : 2 -> n, avec n < N
    nombres_entiers = range(2, n+1)
 
    # on itère k fois le test de Fermat : 0 -> k-1
    for i in range(k):
 
        # choix de l'entier positif ai < N au hasard
        ai = random.choice(nombres_entiers)
 
        # Si ai^(N-1) mod N ≠ 1 : ai^(N-1) n'est pas congru à 1 modulo N
        #if not test_primalite_fermat(N):
        if pow(ai, N-1, N) != 1:
            # N n'est pas premier :
            # renvoie le tuple : (False, i+1) => (est_premier, nbre_iterations)
            return (False, i+1)
 
    # N est probablement premier
    # renvoie le tuple : (True, k) => (est_premier, nbre_iterations)
    return (True, k) 
 
 
print("I. Test de primalité..\n")
 
print("I-A. Version n°1 - Méthode naïve :\n")
 
# valeur de N
N = 99999980021 # 99999980051, 99999980053, 99999980057, 99999980089, 99999980099, 99999980123, 99999980147, 99999980159, 99999980167
 
print("N = " + str(N))
 
print()
 
# Teste si le nombre entier N est premier par la méthode naïve.
est_premier, nbre_iterations = test_primalite(N)
 
# Affiche le résultat
if est_premier:
    print("Résultat du test : {0} est premier.".format(N))
else:
    print("Résultat du test : {0} est composé.".format(N))
 
print("--------------------------------------------")
print("Nombre de division(s) : " + str(nbre_iterations))
 
print();print()
 
print("I-B. Version n°2 - Test probabiliste de Fermat :\n")
 
# paramètre de fiabilité k : nombres d'entiers positifs a1, ..., ak < N choisis au hasard
k = 20
 
print("N = " + str(N)) 
print("k = " + str(k))
 
print()
 
# test de primalité de Fermat pour un nombre entier N et au plus k itérations
est_premier, nbre_iterations = test_primalite_fermat_iter(N, k, n=1000)
 
# Affiche le résultat
if est_premier: # si N est premier
    print("Résultat du test : {0} est probablement premier.".format(N))
else: # sinon
    print("Résultat du test : {0} est composé.".format(N))
 
print("------------------------------------------------------------------------------------------")
print("Nombre de tirage(s) : " + str(nbre_iterations))
print("Nombre d'exponentiation(s) modulaire(s) : " + str(nbre_iterations))
if est_premier:
    # risque d'erreur maxi si ce n'est pas un nombre de Carmichael : t = 1/(2^k)
    t = pow(2,-k)
    print("Risque d'erreur si N n'est pas un nombre de Carmichael inférieur à : " + str(t))
 
print(); print()
 
print("II. Mesure du temps d'exécution des deux fonctions de test sur une série de 20 000 nombres entiers ..\n")
 
# initialisation des compteurs de tests positifs
cpt1 = 0; cpt2 = 0
 
print("Version n°1 - Test de primalité classique :\n")
 
# heure Unix de début (exprimée en secondes)
start = time.time()
 
# parcours des 20 0000 nombres entiers : 99999980000 -> 99999999999
for N in range(99999980000, 100000000000):
 
    # test classique pour N
    est_premier, nbre_iterations = test_primalite(N)
 
    if est_premier: # si le test est positif
        cpt1+=1 # incrémentation du compteur
 
# heure Unix de fin (exprimée en secondes)
end = time.time()
 
# Affiche la durée d'exécution.
print(f"Durée d'exécution : {end - start} sec.\n")
 
 
print("Version n°2 - Test de primalité de Fermat :\n")
 
# heure Unix de début (exprimée en secondes)
start = time.time()
 
# parcours des 20 0000 nombres entiers : 99999980000 -> 99999999999
for N in range(99999980000, 100000000000):
 
    # test de primalité de Fermat pour N et k
    est_premier, nbre_iterations = test_primalite_fermat_iter(N, k=10, n=1000)
 
    if est_premier: # si le test est positif
        cpt2+=1 # incrémentation du compteur
 
# heure Unix de fin (exprimée en secondes)
end = time.time()
 
# Affiche la durée d'exécution.
print(f"Durée d'exécution : {end - start} sec.")
 
print("Taux d'erreur : " + str((cpt2-cpt1)/20000))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Algo arrangements(elements, k, arr=()):
    Si k=0 alors
        afficher(arr)
    Sinon
        Pour i de 0 jusqu'à longueur(elements)-1
            # nouvelle liste d'éléments sans elements[i] (indice débutant à 0 : comme en Python)
	    new_elements = elements[0:i] + elements[i+1:]

            # nouveau tuple avec elements[i] en plus
	    new_arr = arr + elements[i]

            # appel récursif pour k-1
            arrangements(new_elements, k-1, new_arr) 
	Fin pour
    Fin si   
Fin Algo

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def generer_arrangements(elements, k, arr=()):
    # fonction permettant de générer la liste des arrangements de k éléments pris dans un ensemble de n éléments
    # generer_arrangements(['a','b','c'], 2) → ('a','b'), ('a','c'), ('b','a'), ('b','c'), ('c','a'), ('c','b')
 
    # si k=0
    if k==0:
        # On renvoie une liste contenant l'arrangement arr. 
        return [arr]
 
    else: # sinon
        # initialisation de la liste des arrangements
        arrangements = []
 
        # parcours des éléments et de leur indice
        for i,e in enumerate(elements):
            # appel récursif pour k-1 : generer_arrangements(elements[0:i] + elements[i+1:], k-1, arr + (e,))
            # arguments : elements[0:i] + elements[i+1:] → elements, k-1 → k, arr + (e,) → arr 
            arrangements += generer_arrangements(elements[0:i] + elements[i+1:], k-1, arr + (e,))
 
        # Renvoie la liste des arrangements.
        return arrangements

1
2
3
4
5
6
7
8
9
10
11
# valeur de k
k = 3
 
# création de la liste d'éléments
elements = ['a','b','c','d']
 
# Génère la liste des arrangements de k éléments pris dans un ensemble de n éléments.
arrangements = generer_arrangements(elements, k)
 
# Affiche la liste des arrangements.
print(arrangements)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def generateur_arrangements(elements, k, arr=()):
    # fonction permettant de générer tous les arrangements de k éléments pris dans un ensemble de n éléments
    # generateur_arrangements(['a','b','c'], 2) → ('a','b'), ('a','c'), ('b','a'), ('b','c'), ('c','a'), ('c','b')
 
    # si k=0 
    if k==0:
        # On génère l'arrangement arr. 
        yield arr
 
    else: # sinon
        # parcours des éléments et de leur indice
        for i,e in enumerate(elements):
            # appel récursif pour k-1 : generer_arrangements(elements[0:i] + elements[i+1:], k-1, arr + (e,)) : yield from ..
            # arguments : elements[0:i] + elements[i+1:] → elements, k-1 → k, arr + (e,) → arr 
            yield from generateur_arrangements(elements[0:i] + elements[i+1:], k-1, arr + (e,))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# valeur de k
k = 3
 
# création de la liste d'éléments
elements = ['a','b','c','d']
 
# Crée le générateur permettant de parcourir la liste des arrangements de k éléments pris dans un ensemble de n éléments.
gen_arrangements = generateur_arrangements(elements, k)
 
print("Liste des arrangements :\n")
 
# parcours des arrangements un à un
for arrangement in gen_arrangements:
    # Affiche l'arrangement.
    print(arrangement)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# valeur de k
k = 3
 
print("k = " + str(k))
 
# création de la liste d'éléments : numéros au départ
elements = [1, 2, 3, 4]
 
print("E = " + str(elements).replace("[","{").replace("]","}"))
 
print()
 
# Crée le générateur permettant de parcourir la liste des arrangements de k éléments pris dans un ensemble de n éléments
gen_arrangements = generateur_arrangements(elements, k)
 
print("Liste des tiercés dans l'ordre possibles :\n")
 
# parcours les arrangements un à un : tiercés possibles
for arrangement in gen_arrangements:
    # Affiche l'arrangement : tiercé possible
    print(arrangement)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
def generer_arrangements(elements, k, arr=()):
    # fonction permettant de générer la liste des arrangements de k éléments pris dans un ensemble de n éléments
    # generer_arrangements(['a','b','c'], 2) → ('a','b'), ('a','c'), ('b','a'), ('b','c'), ('c','a'), ('c','b')
 
    # si k=0
    if k==0:
        # On renvoie une liste contenant l'arrangement arr. 
        return [arr]
 
    else: # sinon
        # initialisation de la liste des arrangements
        arrangements = []
 
        # parcours des éléments et de leur indice
        for i,e in enumerate(elements):
            # appel récursif pour k-1 : generer_arrangements(elements[0:i] + elements[i+1:], k-1, arr + (e,))
            # arguments : elements[0:i] + elements[i+1:] → elements, k-1 → k, arr + (e,) → arr 
            arrangements += generer_arrangements(elements[0:i] + elements[i+1:], k-1, arr + (e,))
 
        # Renvoie la liste des arrangements.
        return arrangements
 
 
def generateur_arrangements(elements, k, arr=()):
    # fonction permettant de générer tous les arrangements de k éléments pris dans un ensemble de n éléments
    # generateur_arrangements(['a','b','c'], 2) → ('a','b'), ('a','c'), ('b','a'), ('b','c'), ('c','a'), ('c','b')
 
    # si k=0 
    if k==0:
        # On génère l'arrangement arr.
        yield arr
 
    else: # sinon
        # parcours des éléments et de leur indice
        for i,e in enumerate(elements):
            # appel récursif pour k-1 : generer_arrangements(elements[0:i] + elements[i+1:], k-1, arr + (e,)) : yield from ..
            # arguments : elements[0:i] + elements[i+1:] → elements, k-1 → k, arr + (e,) → arr  
            yield from generateur_arrangements(elements[0:i] + elements[i+1:], k-1, arr + (e,))    
 
 
print("Génération des arrangements de k éléments pris dans un ensemble de n éléments..\n")
 
# valeur de k
k = 3
 
print("k = " + str(k))
 
# création de la liste d'éléments
elements = ['a','b','c','d']
 
print("E = " + str(elements).replace("[","{").replace("]","}"))
 
print()
 
print("I. Version n°1 : fonction récursive classique\n")
 
# Génère la liste des arrangements de k éléments pris dans un ensemble de n éléments.
arrangements = generer_arrangements(elements, k)
 
print("Liste des arrangements :\n")
 
# Affiche la liste des arrangements
print(arrangements)
 
print();print()
 
print("II. Version n°2 : fonction génératrice avec yield et yield from\n")
 
# Crée le générateur permettant de parcourir la liste des arrangements de k éléments pris dans un ensemble de n éléments.
gen_arrangements = generateur_arrangements(elements, k)
 
print("Liste des arrangements :\n")
 
# parcours des arrangements un à un
for arrangement in gen_arrangements:
    # Affiche l'arrangement.
    print(arrangement)
 
print();print()
 
print("III. Application : tiercés dans l'ordre\n")
 
# valeur de k
k = 3
 
print("k = " + str(k))
 
# création de la liste d'éléments : numéros au départ
elements = [1, 2, 3, 4]
 
print("E = " + str(elements).replace("[","{").replace("]","}"))
 
print()
 
# Crée le générateur permettant de parcourir la liste des arrangements de k éléments pris dans un ensemble de n éléments
gen_arrangements = generateur_arrangements(elements, k)
 
print("Liste des tiercés dans l'ordre possibles :\n")
 
# parcours les arrangements un à un : tiercés possibles
for arrangement in gen_arrangements:
    # Affiche l'arrangement : tiercé possible
    print(arrangement)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
def generer_combinaisons(elements, k):
    # fonction permettant de générer la liste des combinaisons avec répétition de k éléments pris dans un ensemble de n éléments
    # generer_combinaisons(['a','b','c'], 2) → ('a','a'), ('a','b'), ('a','c'), ('b','b'), ('b','c'), ('c','c')
 
    # nombre d'éléments de l'ensemble
    n = len(elements)
 
    # initialisation de la liste d'indices correspondant à une k-combinaison : [0, 0, .., 0] → ('a', 'a', .., 'a')
    indices = [0] * k
 
    # On crée la 1re combinaison à partir de la liste d'indices. ex. : [0, 0, 0] → ('a', 'a', 'a')
    k_combinaison = tuple(elements[i] for i in indices)
 
    # initialisation de la liste avec la 1re combinaison
    k_combinaisons = [k_combinaison]
 
    # tant que vrai
    while True:
 
        # parcours des indices de la liste
        for i in reversed(range(k)):
            # Si la valeur de l'indice est différente de n-1.
            if indices[i] != n - 1:
                break # On sort de la boucle : indice trouvé
        else: # sinon
            break # On sort de la boucle.
 
        # copie (k-i) fois la valeur (indices[i] + 1) à partir de l'indice i dans la liste 
        indices[i:] = [indices[i] + 1] * (k - i)
 
        # On génère la combinaison correspondant à la liste d'indices.
        k_combinaison = tuple(elements[i] for i in indices)
 
        # ajout de la combinaison à la liste
        k_combinaisons.append(k_combinaison)
 
    # renvoi la liste des combinaisons
    return k_combinaisons

1
2
3
4
5
6
7
8
9
10
11
# valeur de k
k = 4
 
# création de la liste d'éléments
elements = ['a','b','c']
 
# Génère la liste des combinaisons avec répétition de k éléments pris dans un ensemble de n éléments.
combinaisons = generer_combinaisons(elements, k)
 
# Affiche la liste des combinaisons.
print(combinaisons)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
def generateur_combinaisons(elements, k):
    # fonction permettant de générer la liste des combinaisons avec répétition de k éléments pris dans un ensemble de n éléments
    # generateur_combinaisons(['a','b','c'], 2) → ('a','a'), ('a','b'), ('a','c'), ('b','b'), ('b','c'), ('c','c')
 
    # nombre d'éléments de l'ensemble
    n = len(elements)
 
    # initialisation de la liste d'indices correspondant à une k-combinaison : [0, 0, .., 0] → ('a', 'a', .., 'a')
    indices = [0] * k
 
    # On génère la 1re combinaison. ex. : [0, 0, 0] → ('a', 'a', 'a')
    yield tuple(elements[i] for i in indices)
 
    while True:
 
        # parcours des indices de la liste indices
        for i in reversed(range(k)):
            # Si la valeur de l'indice est différente de n-1
            if indices[i] != n - 1:
                break # On sort de la boucle : indice trouvé
        else:
            return
 
        # copie (k-i) fois la valeur (indices[i] + 1) à partir de l'indice i dans la liste
        indices[i:] = [indices[i] + 1] * (k - i)
 
        # On génère la combinaison correspondant aux indices.
        yield tuple(elements[i] for i in indices)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# valeur de k
k = 4
 
# création de la liste d'éléments
elements = ['a','b','c']
 
# Crée le générateur permettant de parcourir la liste des combinaisons avec répétition de k éléments pris dans un ensemble de n éléments.
gen_combinaisons = generateur_combinaisons(elements, k)
 
print("Liste des combinaisons avec répétitions :\n")
 
# Affiche les combinaisons une à une.
for combinaison in gen_combinaisons:
    print(combinaison)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# valeur de k
k = 4
 
# création de la liste d'éléments
elements = ['X', 'Y', 'Z']
 
# Crée le générateur permettant de parcourir la liste des combinaisons avec répétition de k éléments pris dans un ensemble de n éléments.
gen_combinaisons = generateur_combinaisons(elements, k)
 
print("Liste des combinaisons et leur triplet (x, y, z) tel que x + y + z = 4 :\n")
 
# Affiche les triplets (x, y, z) tels que x + y + z = 4.
for combinaison in gen_combinaisons:
    # détermination du nombre de 'X', de 'Y' et de 'Z' pour chaque combinaison
    x, y, z = combinaison.count("X"), combinaison.count("Y"), combinaison.count("Z")
 
    # Affiche la combinaison et son triplet (x, y, z).
    print(str(combinaison) + " → " + str((x, y, z)))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def generer_combinaisons(elements, k):
    # fonction permettant de générer la liste des combinaisons avec répétition de k éléments pris dans un ensemble de n éléments
    # generer_combinaisons(['a','b','c'], 2) → ('a','a'), ('a','b'), ('a','c'), ('b','b'), ('b','c'), ('c','c')
 
    # nombre d'éléments de l'ensemble
    n = len(elements)
 
    # initialisation de la liste d'indices correspondant à une k-combinaison : [0, 0, .., 0] → ('a', 'a', .., 'a')
    indices = [0] * k
 
    # On crée la 1re combinaison à partir de la liste d'indices. ex. : [0, 0, 0] → ('a', 'a', 'a')
    k_combinaison = tuple(elements[i] for i in indices)
 
    # initialisation de la liste avec la 1re combinaison
    k_combinaisons = [k_combinaison]
 
    # tant que vrai
    while True:
 
        # parcours des indices de la liste
        for i in reversed(range(k)):
            # Si la valeur de l'indice est différente de n-1.
            if indices[i] != n - 1:
                break # On sort de la boucle : indice trouvé
        else: # sinon
            break # On sort de la boucle.
 
        # copie (k-i) fois la valeur (indices[i] + 1) à partir de l'indice i dans la liste 
        indices[i:] = [indices[i] + 1] * (k - i)
 
        # On génère la combinaison correspondant à la liste d'indices.
        k_combinaison = tuple(elements[i] for i in indices)
 
        # ajout de la combinaison à la liste
        k_combinaisons.append(k_combinaison)
 
    # renvoi la liste des combinaisons
    return k_combinaisons
 
 
def generateur_combinaisons(elements, k):
    # fonction permettant de générer la liste des combinaisons avec répétition de k éléments pris dans un ensemble de n éléments
    # generateur_combinaisons(['a','b','c'], 2) → ('a','a'), ('a','b'), ('a','c'), ('b','b'), ('b','c'), ('c','c')
 
    # nombre d'éléments de l'ensemble
    n = len(elements)
 
    # initialisation de la liste d'indices correspondant à une k-combinaison : [0, 0, .., 0] → ('a', 'a', .., 'a')
    indices = [0] * k
 
    # On génère la 1re combinaison. ex. : [0, 0, 0] → ('a', 'a', 'a')
    yield tuple(elements[i] for i in indices)
 
    while True:
 
        # parcours des indices de la liste indices
        for i in reversed(range(k)):
            # Si la valeur de l'indice est différente de n-1
            if indices[i] != n - 1:
                break # On sort de la boucle : indice trouvé
        else:
            return
 
        # copie (k-i) fois la valeur (indices[i] + 1) à partir de l'indice i dans la liste
        indices[i:] = [indices[i] + 1] * (k - i)
 
        # On génère la combinaison correspondant aux indices.
        yield tuple(elements[i] for i in indices)
 
 
print("Génération des combinaisons avec répétition de k éléments pris dans un ensemble de n éléments..\n")
 
# valeur de k
k = 4
 
print("k = " + str(k))
 
# création de la liste d'éléments
elements = ['a','b','c']
 
print("E = " + str(elements).replace("[","{").replace("]","}"))
 
print()
 
print("I. Version n°1 : méthode classique\n")
 
# Génère la liste des combinaisons avec répétition de k éléments pris dans un ensemble de n éléments.
combinaisons = generer_combinaisons(elements, k)
 
print("Liste des combinaisons avec répétition :\n")
 
# Affiche la liste des combinaisons
print(combinaisons)
 
print();print()
 
print("II. Version n°2 : fonction génératrice avec yield\n")
 
# Crée le générateur permettant de parcourir la liste des combinaisons avec répétition de k éléments pris dans un ensemble de n éléments.
gen_combinaisons = generateur_combinaisons(elements, k)
 
print("Liste des combinaisons avec répétition :\n")
 
# Affiche les combinaisons une à une.
for combinaison in gen_combinaisons:
    print(combinaison)
 
print();print()
 
print("III. Application : génération des triplets (x, y, z) tels que x + y + z = 4\n")
 
# valeur de k
k = 4
 
print("k = " + str(k))
 
# création de la liste d'éléments
elements = ['X', 'Y', 'Z']
 
print("E = " + str(elements).replace("[","{").replace("]","}"))
 
print()
 
# Crée le générateur permettant de parcourir la liste des combinaisons avec répétition de k éléments pris dans un ensemble de n éléments
gen_combinaisons = generateur_combinaisons(elements, k)
 
print("Liste des combinaisons et leur triplet (x, y, z) tel que x + y + z = 4 :\n")
 
# Affiche les triplets (x, y, z) tels que x + y + z = 4
for combinaison in gen_combinaisons:
    # détermination du nombre de 'X', de 'Y' et de 'Z' pour chaque combinaison
    x, y, z = combinaison.count("X"), combinaison.count("Y"), combinaison.count("Z")
 
    # Affiche la combinaison et son triplet (x, y, z).
    print(str(combinaison) + " → " + str((x, y, z)))