Structures répétitives

a) 3 2 1
b) 0 3 6 9
c) (rien : range vide, pas d'erreur)
d) 10
e) 1 0 | 2 0 | 2 1 |
f) 3

Le cas (e) merite attention : pour i = 0, range(0) est vide, donc aucun affichage ; le corps interne est execute \(0+1+2 = 3\) fois au total.

def somme_diviseurs(n):
    """Version naive, en Theta(n)."""
    s = 0
    for d in range(1, n):
        if n % d == 0:
            s += d
    return s

for n in range(2, 10000):
    if somme_diviseurs(n) == n:
        print(n, end=" ")        # 6 28 496 8128

3. Les diviseurs vont par paires \((d, n/d)\) et l'un des deux est toujours \(\le \sqrt n\). Il suffit donc de parcourir \(d\) jusqu'a \(\lfloor\sqrt n\rfloor\) en ajoutant les deux membres de chaque paire, en prenant garde a ne pas compter deux fois \(\sqrt n\) quand \(n\) est un carre parfait, ni \(n\) lui-meme.

def somme_diviseurs(n):
    """Somme des diviseurs stricts de n >= 1, en Theta(sqrt(n))."""
    if n == 1:
        return 0
    s = 1                        # 1 divise n, et n/1 = n est exclu
    d = 2
    while d * d <= n:
        if n % d == 0:
            s += d
            q = n // d
            if q != d:           # eviter le double comptage si n = d^2
                s += q
        d += 1
    return s

for n in range(1, 10000):
    assert somme_diviseurs(n) == sum(d for d in range(1, n) if n % d == 0)

Le gain est spectaculaire : la recherche jusqu'a \(10^6\) passe de plusieurs minutes a moins d'une seconde.

def temps_de_vol(n):
    """Nombre d'etapes de la suite de Syracuse pour passer de n a 1."""
    t = 0
    while n != 1:
        n = n // 2 if n % 2 == 0 else 3*n + 1
        t += 1
    return t

def altitude_max(n):
    m = n
    while n != 1:
        n = n // 2 if n % 2 == 0 else 3*n + 1
        if n > m:
            m = n
    return m

meilleur, tmax = 1, 0
for n in range(1, 10001):
    t = temps_de_vol(n)
    if t > tmax:
        meilleur, tmax = n, t
print(meilleur, tmax)          # 6171 261
print(altitude_max(27))        # 9232, alors que 27 a un temps de vol de 111

4. Non. C'est la conjecture de Collatz : on ignore encore aujourd'hui si la suite atteint 1 pour tout \(n\), et aucun variant n'est connu. La verification a ete menee par ordinateur jusqu'a environ \(2^{68}\) sans contre-exemple, mais aucune preuve n'existe. Cet exemple montre qu'une boucle while parfaitement anodine peut cacher un probleme ouvert : la terminaison n'est jamais gratuite.

def sablier(n):
    """Sablier d'etoiles de hauteur 2n-1."""
    for i in range(n):                 # partie haute, i = 0..n-1
        print(" " * i + "*" * (2*(n-i) - 1))
    for i in range(n - 2, -1, -1):     # partie basse, symetrique
        print(" " * i + "*" * (2*(n-i) - 1))

sablier(4)

Version sans multiplication de chaines :

def sablier_boucles(n):
    for i in range(2*n - 1):
        d = i if i < n else 2*n - 2 - i     # distance au sommet du sablier
        for _ in range(d):
            print(" ", end="")
        for _ in range(2*(n - d) - 1):
            print("*", end="")
        print()

3. La ligne de rang \(d\) (avec \(d = 0,\dots,n-1\)) contient \(2(n-d)-1\) etoiles et apparait deux fois, sauf la ligne \(d = n-1\) (une seule etoile) qui n'apparait qu'une fois. Le total vaut donc \[ 2\sum_{d=0}^{n-1}\bigl(2(n-d)-1\bigr) - 1 = 2\sum_{j=1}^{n}(2j-1) - 1 = 2n^2 - 1, \] en utilisant \(\sum_{j=1}^{n}(2j-1) = n^2\). Pour \(n=4\) : \(2\cdot16-1 = 31\) etoiles.

def crible(N):
    """Liste des nombres premiers <= N (crible d'Eratosthene)."""
    if N < 2:
        return []
    est_premier = [True] * (N + 1)
    est_premier[0] = est_premier[1] = False
    p = 2
    while p * p <= N:
        if est_premier[p]:
            for m in range(p*p, N + 1, p):     # on barre p^2, p^2+p, ...
                est_premier[m] = False
        p += 1
    return [k for k in range(N + 1) if est_premier[k]]

assert crible(30) == [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
assert len(crible(10**6)) == 78498

2. Soit \(m = kp\) un multiple de \(p\) avec \(k < p\). Alors \(m\) admet le facteur \(k < p\), donc \(m\) a deja ete barre lors du traitement d'un nombre premier inferieur a \(p\) (a savoir un facteur premier de \(k\)). Il est donc inutile de recommencer : le premier multiple non encore barre est \(p \cdot p\).

3. Si un entier \(m \le N\) est compose, il admet un facteur premier \(p \le \sqrt m \le \sqrt N\). Il aura donc ete barre par un \(p \le \sqrt N\). Passe cette valeur, tous les composes sont deja elimines.

4. Le crible est en \(O(N \log\log N)\) operations, contre \(O(N\sqrt N)\) pour le test naif applique a chaque entier. Sur \(N = 10^6\), le crible prend moins d'une seconde, le test naif plusieurs dizaines de secondes.

import time
t0 = time.perf_counter(); crible(10**6); print(time.perf_counter() - t0)

1. \((1071,462) \to (462,147) \to (147,21) \to (21,0)\). Le pgcd vaut 21.

2. Terminaison. Variant : \(v = b\). C'est un entier naturel (un reste est toujours \(\ge 0\)). A chaque tour, \(b_{\text{nouveau}} = a \bmod b < b\) car \(b \neq 0\) : le variant decroit strictement. Donc la boucle termine.

3. Correction. Invariant : \(\operatorname{pgcd}(a,b) = \operatorname{pgcd}(a_0, b_0)\), ou \((a_0,b_0)\) sont les valeurs initiales. Il est vrai a l'initialisation. Il est conserve car, pour \(b \neq 0\), en ecrivant \(a = bq + r\), tout diviseur commun de \(a\) et \(b\) divise \(r = a - bq\), et reciproquement tout diviseur commun de \(b\) et \(r\) divise \(a = bq+r\) : les deux couples ont donc les memes diviseurs communs, donc le meme pgcd. A la sortie, \(b = 0\) et l'invariant donne \(\operatorname{pgcd}(a_0,b_0) = \operatorname{pgcd}(a,0) = a\), qui est bien la valeur renvoyee.

4. Soient \((a,b)\) puis \((b, r)\) puis \((r, r')\) trois etats consecutifs, avec \(r = a \bmod b\) et \(r' = b \bmod r\). Deux cas :

• si \(r \le a/2\), c'est termine ;
• si \(r > a/2\), alors comme \(r < b \le a\), la division \(a = bq + r\) impose \(q = 1\) et \(r = a - b\); or \(r > a/2\) donne \(b < a/2\), et \(r' < r\) avec \(r' = b \bmod r < b < a/2\).

Dans tous les cas, le premier terme est divise au moins par 2 toutes les deux iterations. Le nombre de tours est donc majore par \(2\log_2 a = O(\log a)\) : l'algorithme d'Euclide est logarithmique, ce qui explique son efficacite meme sur de tres grands entiers.

def leibniz(n):
    """Somme partielle d'ordre n de la serie de Leibniz, multipliee par 4."""
    s = 0.0
    for k in range(n + 1):
        s += (-1)**k / (2*k + 1)
    return 4 * s

print(leibniz(1000))      # 3.1425916543395442

2. La serie est alternee a termes decroissants en valeur absolue vers 0 : le critere des series alternees donne \(\bigl|\frac\pi4 - S_n\bigr| \le \frac{1}{2n+3}\), donc pour \(\pi\) l'erreur est majoree par \(\frac{4}{2n+3}\). Pour obtenir 3 decimales, il faut \(\frac{4}{2n+3} < 5\cdot10^{-4}\), soit \(n > 3999\). C'est une convergence tres lente : chaque decimale supplementaire coute un facteur 10 sur le nombre de termes.

def leibniz_seuil(eps):
    """Somme jusqu'a ce que le terme soit < eps. Boucle non bornee : while."""
    s, k, terme = 0.0, 0, 1.0
    while terme >= eps:
        s += (-1)**k * terme
        k += 1
        terme = 1 / (2*k + 1)
    return 4*s, k

print(leibniz_seuil(1e-6))    # (3.1415946..., 500000)

On utilise un while car le nombre de termes n'est pas connu a l'avance : c'est precisement le critere de choix entre les deux boucles.

from math import atan, pi
def machin(n):
    """arctan(x) = somme (-1)^k x^(2k+1)/(2k+1)."""
    def arctan(x, n):
        s = 0.0
        for k in range(n + 1):
            s += (-1)**k * x**(2*k+1) / (2*k + 1)
        return s
    return 4 * (4*arctan(1/5, n) - arctan(1/239, n))

print(machin(10) - pi)     # de l'ordre de 1e-16 avec seulement 11 termes !

La formule de Machin converge geometriquement (raison \(1/25\)) : 11 termes suffisent la ou Leibniz en demande des millions. Le choix de la formule compte plus que le choix du langage.

1. Version \(\Theta(n^2)\) (la version \(\Theta(n^3)\) recalcule la somme interne a chaque fois) :

def max_segment_n2(t):
    n = len(t)
    M = t[0]
    for i in range(n):
        s = 0
        for j in range(i, n):
            s += t[j]              # somme du segment t[i..j], calculee de proche en proche
            if s > M:
                M = s
    return M

2. Un segment finissant en \(j\) est soit le singleton \(\{j\}\), soit un segment finissant en \(j-1\) prolonge par \(t_j\). D'ou \[ m_j = \max\bigl(t_j,\ m_{j-1} + t_j\bigr) = t_j + \max(0, m_{j-1}), \qquad m_0 = t_0. \] Et \(M = \max_{0 \le j < n} m_j\).

def kadane(t):
    """Somme maximale d'un segment non vide de t. Theta(n), Theta(1) memoire."""
    m = M = t[0]
    for j in range(1, len(t)):
        m = t[j] + max(0, m)       # meilleur segment finissant en j
        if m > M:
            M = m
    return M

assert kadane([-2, 1, -3, 4, -1, 2, 1, -5, 4]) == 6      # segment [4,-1,2,1]
assert kadane([-5, -2, -9]) == -2                        # que des negatifs
import random
for _ in range(300):
    t = [random.randint(-9, 9) for _ in range(random.randint(1, 12))]
    assert kadane(t) == max_segment_n2(t)

4. Invariant. Avant chaque tour d'indice \(j\), on a :

\[ \texttt{m} = m_{j-1} = \max_{i \le j-1}\sum_{k=i}^{j-1} t_k \quad\text{et}\quad \texttt{M} = \max_{0 \le l \le j-1} m_l = \max_{i \le l \le j-1}\sum_{k=i}^{l} t_k. \]

• Initialisation : pour \(j = 1\), m = M = t[0], ce qui est bien \(m_0\) et le maximum sur les segments inclus dans \([0,0]\).
• Conservation : le tour calcule \(m \leftarrow t_j + \max(0, m)\), ce qui vaut \(m_j\) d'apres la question 2, puis \(M \leftarrow \max(M, m)\), ce qui etend le maximum aux segments finissant en \(j\). L'invariant est retabli pour \(j+1\).
• Conclusion : a la sortie, \(j = n\) et l'invariant donne \(\texttt{M} = \max_{0 \le i \le l < n} \sum_{k=i}^{l} t_k = M\). L'algorithme est correct.

Complexite. Un seul parcours, un nombre constant d'operations par tour : \(\Theta(n)\) en temps et \(\Theta(1)\) en espace. C'est optimal, puisqu'il faut au minimum lire toutes les cases.