Chapitre 36 — programmer les algorithmes de l'arithmétique : division euclidienne, PGCD et Bezout, congruences, nombres premiers, puis les chiffrements qui en découlent.
Naviguez avec les flèches ← → du clavier, ou via le Sommaire.
//, %, diviseurs, algorithme d'Euclidepow(a, k, n), restes chinoisLes entiers de Python sont exacts et de taille illimitée ($2^{2026}$ a $610$ chiffres) : tant qu'on reste dans les entiers, avec // et %, tout calcul arithmétique est juste.
Pour $a\in\mathbb{Z}$ et $b\in\mathbb{Z}^*$ : a//b est le quotient, a%b le reste, divmod(a,b) le couple des deux. La division a/b renvoie un flottant : on ne l'utilise jamais en arithmétique.
Si $b>0$, a%b $\in\llbracket 0,b-1\rrbracket$, même pour $a<0$ : c'est le reste du cours. Si $b<0$, le reste prend le signe de $b$.
>>> 17 // 5, 17 % 5 (3, 2) >>> -17 // 5, -17 % 5 (-4, 3) >>> 17 // -5, 17 % -5 (-4, -3) >>> divmod(-37, 11) (-4, 7) >>> 2**100 1267650600228229401496703205376 >>> 3**200 % 7 2
$k$ divise $n$ si et seulement si le reste est nul : n % k == 0.
def diviseurs(n) : L = [] for k in range(1, n+1) : if n % k == 0 : L.append(k) return L >>> 2023 % 7 == 0, 2024 % 7 == 0 (True, False) >>> diviseurs(2025) [1, 3, 5, 9, 15, 25, 27, 45, 75, 81, 135, 225, 405, 675, 2025]
Les diviseurs vont par paires $d$ et $\frac{n}{d}$ : tester les $d$ tels que $d^2\leqslant n$ suffit. Pour $n=10^7$, on passe de dix millions de divisions à $3\,162$.
Si $r$ est le reste de $a$ par $b$ :
def pgcd(a, b) : while b != 0 : r = a % b a = b b = r return a >>> pgcd(7119, 4977), pgcd(2235282, 32718) (63, 42)
b décroît strictement dans $\mathbb{N}$ : la boucle s'arrête. $\operatorname{pgcd}(\mathtt{a},\mathtt{b})$ est un invariant de boucle. Moins de $2\log_2 b$ tours : l'algorithme est en $O(\log b)$.
Avec $r_0=a$, $r_1=b$, $(u_0,v_0)=(1,0)$, $(u_1,v_1)=(0,1)$ et $q_{k+1}$ le quotient de $r_{k-1}$ par $r_k$ :
Au dernier reste non nul : $au+bv=\operatorname{pgcd}(a,b)$.
def euclideetendu(a, b) : u0, v0 = 1, 0 # a = a*u0 + b*v0 u1, v1 = 0, 1 # b = a*u1 + b*v1 while b != 0 : q, r = divmod(a, b) a, b = b, r u0, v0, u1, v1 = u1, v1, u0-q*u1, v0-q*v1 return a, u0, v0 >>> euclideetendu(7119, 4977) (63, 7, -10) >>> euclideetendu(2022, 557) (1, 73, -265)
$7119\times 7-4977\times 10=63$ et $2022\times 73-557\times 265=1$, sans aucun essai.
Si $\operatorname{pgcd}(a,n)=1$, Bezout donne $au+nv=1$, donc $au\equiv 1\ [n]$ : l'inverse de $a$ est le reste de $u$ modulo $n$. Sinon, $a$ n'a pas d'inverse. Python l'obtient aussi par pow(a, -1, n).
def inversemod(a, n) : d, u, v = euclideetendu(a, n) if d != 1 : return None # pas d'inverse return u % n >>> inversemod(9, 26), inversemod(7, 26), inversemod(4, 26) (3, 15, None) >>> pow(9, -1, 26) 3
Avec $d=\operatorname{pgcd}(a,n)$. Si $d\nmid b$ : aucune solution ($6x\equiv 5\ [26]$). Sinon, on divise tout par $d$ et on multiplie par l'inverse : $d$ classes modulo $n$. Ainsi $7x\equiv 3\ [26]\iff x\equiv 19$, et $6x\equiv 4\ [26]\iff x\equiv 5$ ou $18$.
Solutions si et seulement si $d=\operatorname{pgcd}(a,b)$ divise $c$. Euclide étendu fournit une solution $(x_0,y_0)$, et toutes s'écrivent :
def diophante(a, b, c) : d, u, v = euclideetendu(a, b) if c % d != 0 : return "vide" x0, y0 = u*(c//d), v*(c//d) return x0, y0, b//d, -a//d >>> diophante(114, 30, 54) (-9, 36, 5, -19) >>> diophante(114, 30, 55) 'vide' >>> diophante(2022, 557, 1) (73, -265, 557, -2022)
$114x+30y=54$ : solutions $(-9+5k,\ 36-19k)$ ; $k=2$ redonne $(1,-2)$. Une recherche par essais dans $\llbracket -10,9\rrbracket^2$ échouerait sur $2022x+557y=1$.
On réduit modulo $n$ à chaque étape : au plus $2\log_2 k$ multiplications, sur des entiers inférieurs à $n^2$. C'est l'algorithme de pow(a, k, n).
def expomod(a, k, n) : r = 1 a = a % n while k > 0 : if k % 2 == 1 : r = (r*a) % n a = (a*a) % n k = k // 2 return r >>> expomod(3, 13, 7), expomod(7, 2026, 10) (3, 9) >>> pow(3, 13, 7), pow(7, 2026, 10) (3, 9) >>> pow(2, 10**18, 10**9 + 7) 719476260
$13=1101_2=8+4+1$, donc $3^{13}=3^1\times 3^4\times 3^8\equiv 3\ [7]$ : quatre tours au lieu de treize multiplications. $2^{10^{18}}$ a $3\times 10^{17}$ chiffres, son reste est pourtant immédiat.
Si $p$ est premier et ne divise pas $a$, alors $a^{p-1}\equiv 1\ [p]$. Donc si $a^{n-1}\not\equiv 1\ [n]$ avec $a$ premier avec $n$, $n$ n'est pas premier : $853661=41\times 20821$ est démasqué sans chercher de diviseur.
def fermat(n, a) : return pow(a, n-1, n) == 1 >>> fermat(131071, 2), fermat(853661, 2) (True, False) >>> [n for n in range(3, 2000, 2) if fermat(n, 2) and not estpremier(n)] [341, 561, 645, 1105, 1387, 1729, 1905] >>> all(fermat(561, a) for a in range(2, 561) if pgcd(a, 561) == 1) True
$341=11\times 31$ passe le test en base $2$ : c'est un pseudo-premier. $561=3\times 11\times 17$ le passe pour toutes les bases premières avec lui (nombre de Carmichaël). Condition nécessaire, jamais une preuve.
Avec $\operatorname{pgcd}(m,n)=1$, le système $x\equiv a\ [m]$, $x\equiv b\ [n]$ a une unique classe de solutions modulo $mn$. On pose $x=a+mt$ : il reste $mt\equiv b-a\ [n]$, d'où $t\equiv (b-a)u\ [n]$ avec $u$ l'inverse de $m$ modulo $n$.
def chinois(a, m, b, n) : u = inversemod(m, n) # m*u = 1 [n] t = (b - a)*u % n return (a + m*t) % (m*n) >>> chinois(1, 5, 5, 7) 26 >>> x = chinois(2, 3, 3, 5) >>> x = chinois(x, 15, 2, 7) >>> x, x % 3, x % 5, x % 7 (23, 2, 3, 2)
Au IIIe siècle, Sun Zi cherche un nombre de « reste $2$ par $3$, $3$ par $5$ et $2$ par $7$ ». En enchaînant deux appels (modules $15$ puis $105$), on trouve $23$.
Tout entier non premier $n\geqslant 2$ admet un diviseur premier $p$ tel que $p^2\leqslant n$ : il suffit de chercher un diviseur jusqu'à $\sqrt{n}$. La fonction renvoie un booléen, pour servir dans un if ou une liste.
def estpremier(n) : if n < 2 : return False d = 2 while d*d <= n : if n % d == 0 : return False d = d + 1 return True >>> [n for n in range(30) if estpremier(n)] [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] >>> estpremier(131071), estpremier(853661) (True, False) >>> estpremier(2**31 - 1) True
Le pire cas est un nombre premier : $\sqrt{n}$ divisions, soit $46\,341$ pour $2^{31}-1$ (instantané), mais des années pour trente chiffres.
Une liste de $N+1$ booléens, tous True au départ. Pour chaque $d$ non barré, on barre ses multiples à partir de $d^2$ : les plus petits l'ont déjà été par un facteur plus petit. Les non barrés sont premiers.
def crible(N) : premier = [True]*(N+1) premier[0] = False premier[1] = False for d in range(2, N+1) : if premier[d] : for m in range(d*d, N+1, d) : premier[m] = False return [n for n in range(N+1) if premier[n]] >>> crible(50) [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47] >>> len(crible(1000)), len(crible(10**6)) (168, 78498)
$168$ premiers sous $1000$, $78\,498$ sous un million ; $\dfrac{N}{\ln N}$ en donne l'ordre de grandeur. Crible jusqu'à $10^6$ : $0{,}04$ s, contre $1{,}7$ s avec estpremier.
On divise $n$ par $d=2,3,4,\dots$ tant que c'est possible, et l'on s'arrête dès que $d^2>n$ : l'entier restant, s'il est différent de $1$, est premier.
def facteurs(n) : d = 2 L = [] while d*d <= n : while n % d == 0 : L.append(d) n = n//d d = d + 1 if n > 1 : L.append(n) return L >>> facteurs(5632) [2, 2, 2, 2, 2, 2, 2, 2, 2, 11] >>> facteurs(9785873259) [3, 3, 19, 229, 269, 929] >>> facteurs(3**36) == [3]*36 True
n = n/d transforme $n$ en flottant (12/2 vaut 6.0), exact seulement jusqu'à $2^{53}$ : au-delà, la boucle ne se termine plus ($3^{36}\simeq 1{,}5\times 10^{17}$). Toujours n = n//d.
Si $n=p_1^{k_1}\times\dots\times p_r^{k_r}$ :
exposants(n) regroupe d'abord les facteurs de facteurs(n) en couples $[p,k]$.
def diviseursfacto(n) : nb, s = 1, 1 for p, k in exposants(n) : nb = nb*(k+1) s = s*(p**(k+1) - 1)//(p - 1) return nb, s >>> exposants(4116) [[2, 2], [3, 1], [7, 3]] >>> diviseursfacto(4116) (24, 11200) >>> len(diviseurs(4116)), sum(diviseurs(4116)) (24, 11200) >>> [n for n in range(2, 10000) if diviseursfacto(n)[1] == 2*n] [6, 28, 496, 8128]
$4116=2^2\times 3\times 7^3$ : $\tau=3\times 2\times 4=24$ et $\sigma=7\times 4\times 400=11\,200$. La dernière ligne liste les nombres parfaits inférieurs à $10^4$.
A a le rang $0$, Z le rang $25$ : rang(c) vaut ord(c) - ord("A") et lettre(k) fait l'inverse avec chr. César : $y\equiv x+k\ [26]$. Affine : $y\equiv ax+b\ [26]$, déchiffré par $x\equiv u(y-b)\ [26]$ où $u$ est l'inverse de $a$.
def cesar(texte, k) : resultat = "" for c in texte : resultat = resultat + lettre((rang(c) + k) % 26) return resultat def affine(texte, a, b) : resultat = "" for c in texte : resultat = resultat + lettre((a*rang(c) + b) % 26) return resultat >>> cesar("ARITHMETIQUE", 3) 'DULWKPHWLTXH' >>> cesar("DULWKPHWLTXH", -3) 'ARITHMETIQUE' >>> affine("CODE", 5, 8) 'SAXC'
$a$ doit être premier avec $26$ : avec $a=2$, $b=1$, B et O donnent toutes deux D. $12$ valeurs de $a$, $26$ de $b$ : $312$ clés, essayées en un instant.
Si $\delta=ad-bc$ est premier avec $26$, d'inverse $\delta'$ : $x_1\equiv \delta'(dy_1-by_2)$ et $x_2\equiv \delta'(-cy_1+ay_2)\ [26]$. C'est l'inverse d'une matrice $2\times 2$, où $\delta'$ remplace $\frac{1}{\delta}$.
def hill(x1, x2, a, b, c, d) : y1 = (a*x1 + b*x2) % 26 y2 = (c*x1 + d*x2) % 26 return y1, y2 def dehill(y1, y2, a, b, c, d) : delta = inversemod(a*d - b*c, 26) x1 = delta*(d*y1 - b*y2) % 26 x2 = delta*(-c*y1 + a*y2) % 26 return x1, x2 >>> inversemod(5*7 - 13*2, 26) 3 >>> hill(2, 14, 5, 13, 2, 7), dehill(13, 2, 5, 13, 2, 7) ((10, 24), (13, 4))
Clé $(5,13,2,7)$, $\delta'=3$ : le bloc CO $=(2,14)$ devient KY $=(10,24)$ et NC se déchiffre en NE. Avec numpy, ne pas utiliser np.linalg.inv (inverse réel).
def clesrsa(p, q, e) : n = p*q phi = (p-1)*(q-1) d = inversemod(e, phi) return (n, e), (n, d) def rsa(M, cle) : n, e = cle return pow(M, e, n) >>> publique, privee = clesrsa(61, 53, 17) >>> publique, privee ((3233, 17), (3233, 2753)) >>> rsa(2026, publique), rsa(352, privee) (352, 2026) >>> all(rsa(rsa(M, publique), privee) == M for M in range(3233)) True
Fermat donne $M^{ed}\equiv M$ modulo $p$ et modulo $q$, Gauss modulo $pq$. Casser la clé demande de factoriser $n$ : immédiat pour $999\,985\,999\,949=999\,983\times 1\,000\,003$, hors de portée pour $617$ chiffres.
I.N.S.E.E. : clé $97-r$, où $r$ est le reste des treize premiers chiffres modulo $97$ (clefinsee(2620219223078) vaut $89$). I.S.B.N. à treize chiffres : valide si
def verifisbn(code) : s = 0 for i in range(13) : if i % 2 == 0 : s = s + int(code[i]) else : s = s + 3*int(code[i]) return s % 10 == 0 >>> verifisbn("9782100545247") True >>> verifisbn("9782100545274") False
Un chiffre faux : toujours. Deux chiffres voisins échangés : sauf si leur écart vaut $\pm 5$. L'ancien code modulo $11$, premier, détectait toutes les transpositions.
a//b, a%b, divmod(a,b)n % k == 0euclideetendu(a, b) : $(d,u,v)$, $au+bv=d$pow(a, -1, n), si $\operatorname{pgcd}(a,n)=1$pow(a, k, n) : $O(\log k)$/ donne un flottant, exact seulement jusqu'à $2^{53}$ : en arithmétique, toujours //.a**k % n calcule d'abord $a^k$ en entier : pour un grand $k$, utiliser pow(a, k, n).-17 % 5 vaut $3$, pas $-2$ : le reste est dans $\llbracket 0,b-1\rrbracket$ dès que $b>0$.inversemod renvoie None.range(1, n+1) pour aller jusqu'à $n$ : la borne de droite est exclue.-17 // 5, -17 % 5 ?(-4, 3) : $-17=5\times(-4)+3$, avec $0\leqslant 3<5$.pow(a, k, n) plutôt que a**k % n ?inversemod(4, 26) ?None : $\operatorname{pgcd}(4,26)=2\neq 1$, donc $4$ n'a pas d'inverse modulo $26$.fermat(n, 2) vaut True : $n$ est-il premier ?