Bibliothèques mathématiques

a) 4.0            un FLOTTANT, meme pour un carre parfait
b) ValueError: math domain error
c) 4              racine entiere, exacte
d) 1.2246467991473532e-16   et non 0.0 : pi est approche
e) 3.0
f) (-3, -2)       floor arrondit vers -inf, ceil vers +inf
g) False          nan n'est egal a rien, pas meme a lui-meme
h) nan
i) 1              les deux bornes sont incluses
j) None           shuffle modifie en place et ne renvoie rien

Le point (d) est important : \(\sin(\pi)\) ne vaut pas 0 en machine parce que la constante pi n'est qu'une approximation de \(\pi\) sur 53 bits. Le resultat obtenu est de l'ordre de l'erreur sur \(\pi\), soit \(10^{-16}\). Tester sin(pi) == 0 serait faux ; utiliser isclose(sin(pi), 0, abs_tol=1e-12).

import statistics as st

def moyenne(t):
    """Moyenne arithmetique. Precondition : t non vide."""
    return sum(t) / len(t)

def variance_pop(t):
    """Variance de POPULATION (division par n)."""
    m = moyenne(t)
    return sum((x - m)**2 for x in t) / len(t)

def variance_ech(t):
    """Variance d'ECHANTILLON (division par n-1, estimateur sans biais)."""
    m = moyenne(t)
    return sum((x - m)**2 for x in t) / (len(t) - 1)

def ecart_type(t, population=True):
    return (variance_pop(t) if population else variance_ech(t)) ** 0.5

def mediane(t):
    u = sorted(t); n = len(u)
    return u[n//2] if n % 2 else (u[n//2 - 1] + u[n//2]) / 2

def quartiles(t):
    """(Q1, Q2, Q3) par la methode des mediane des demi-listes."""
    u = sorted(t); n = len(u)
    bas = u[:n//2]
    haut = u[(n+1)//2:]            # on exclut la mediane si n est impair
    return mediane(bas), mediane(u), mediane(haut)

def resume(t):
    q1, q2, q3 = quartiles(t)
    return {"n": len(t), "min": min(t), "Q1": q1, "mediane": q2, "Q3": q3,
            "max": max(t), "moyenne": moyenne(t),
            "ecart_type": ecart_type(t), "etendue": max(t) - min(t),
            "interquartile": q3 - q1}

notes = [12, 15, 9, 18, 11, 15, 7, 14]
assert abs(moyenne(notes) - st.mean(notes)) < 1e-12
assert abs(variance_ech(notes) - st.variance(notes)) < 1e-12
assert abs(variance_pop(notes) - st.pvariance(notes)) < 1e-12
assert mediane(notes) == st.median(notes)
print(resume(notes))

{'n': 8, 'min': 7, 'Q1': 10.0, 'mediane': 13.0, 'Q3': 15.0, 'max': 18,
 'moyenne': 12.625, 'ecart_type': 3.339..., 'etendue': 11, 'interquartile': 5.0}

import random
from collections import Counter

random.seed(42)

# 1. uniformite
tirages = [random.randint(1, 6) for _ in range(10000)]
c = Counter(tirages)
for f in range(1, 7):
    print(f, c[f], f"{c[f]/10000:.4f}  (theorie 0.1667)")

# 2. somme de deux des
sommes = Counter(random.randint(1,6) + random.randint(1,6) for _ in range(100000))
print(f"{'s':>3}{'simule':>10}{'theorie':>10}")
for s in range(2, 13):
    theorie = (6 - abs(s - 7)) / 36
    print(f"{s:>3}{sommes[s]/100000:>10.4f}{theorie:>10.4f}")

# 3. au moins un six en quatre lancers
n = 200000
succes = sum(1 for _ in range(n)
             if any(random.randint(1,6) == 6 for _ in range(4)))
print(f"simule {succes/n:.4f}, exact {1 - (5/6)**4:.4f}")   # ~ 0.5177

# 4. paradoxe des anniversaires
def anniversaires(n, essais=20000):
    """Probabilite qu'au moins deux personnes parmi n partagent une date."""
    collisions = 0
    for _ in range(essais):
        dates = [random.randrange(365) for _ in range(n)]
        if len(set(dates)) < n:
            collisions += 1
    return collisions / essais

def exacte(n):
    """1 - 365!/((365-n)! * 365^n)."""
    p = 1.0
    for k in range(n):
        p *= (365 - k) / 365
    return 1 - p

for n in range(20, 26):
    print(f"n = {n} : simule {anniversaires(n):.4f}, exact {exacte(n):.4f}")

n = 20 : simule 0.4102, exact 0.4114
n = 21 : simule 0.4457, exact 0.4437
n = 22 : simule 0.4738, exact 0.4757
n = 23 : simule 0.5086, exact 0.5073   <- premiere fois > 1/2
n = 24 : simule 0.5390, exact 0.5383
n = 25 : simule 0.5670, exact 0.5687

La reponse est \(n = 23\). Le resultat surprend car on confond intuitivement « quelqu'un partage mon anniversaire » (probabilite \(1 - (364/365)^{22} \approx 6\%\)) avec « deux personnes quelconques partagent un anniversaire ». Le nombre de paires est \(\binom{23}{2} = 253\), et non 22 : c'est cette croissance quadratique qui explique le phenomene.

import random
from math import pi, exp, sqrt, gamma

random.seed(1)

def estime_pi(n):
    d = sum(1 for _ in range(n)
            if (lambda x, y: x*x + y*y <= 1)(random.random(), random.random()))
    return 4 * d / n

def integrale_mc(f, a, b, n):
    """Estime l'integrale de f sur [a,b] par la moyenne des valeurs."""
    return (b - a) * sum(f(random.uniform(a, b)) for _ in range(n)) / n

def volume_boule_mc(d, n):
    """Volume de la boule unite en dimension d, par rejet dans le cube [-1,1]^d."""
    dedans = 0
    for _ in range(n):
        s = sum(random.uniform(-1, 1)**2 for _ in range(d))
        if s <= 1:
            dedans += 1
    return (2 ** d) * dedans / n

def volume_exact(d):
    """pi^(d/2) / Gamma(d/2 + 1)."""
    return pi**(d/2) / gamma(d/2 + 1)

print(f"pi ~ {estime_pi(10**6):.5f}")
print(f"integrale ~ {integrale_mc(lambda x: exp(-x*x), 0, 1, 10**6):.6f}"
      f"  (valeur 0.746824)")

for d in (2, 3, 5, 10):
    print(f"d = {d:>2} : simule {volume_boule_mc(d, 200000):.4f}, "
          f"exact {volume_exact(d):.4f}")

d =  2 : simule 3.1424, exact 3.1416
d =  3 : simule 4.1858, exact 4.1888
d =  5 : simule 5.2864, exact 5.2638
d = 10 : simule 2.5395, exact 2.5502

# 4. decroissance de l'erreur
print(f"{'N':>9}{'erreur':>12}{'erreur x sqrt(N)':>20}")
for k in range(2, 8):
    N = 10**k
    err = abs(estime_pi(N) - pi)
    print(f"{N:>9}{err:>12.6f}{err * sqrt(N):>20.4f}")

        N      erreur    erreur x sqrt(N)
      100    0.058407              0.5841
     1000    0.045593              1.4418
    10000    0.010393              1.0393
   100000    0.004433              1.4019
  1000000    0.000546              0.5460
 10000000    0.000390              1.2334

Le produit \(\text{erreur} \times \sqrt N\) reste de l'ordre de l'unite quand \(N\) varie sur cinq ordres de grandeur : cela confirme \(\text{erreur} = \Theta(1/\sqrt N)\), conformement au theoreme central limite (la constante est \(\sigma = \sqrt{\pi(4-\pi)} \approx 1{,}64\)).

from fractions import Fraction
from math import gcd

def harmonique(n):
    """H_n exact."""
    return sum(Fraction(1, k) for k in range(1, n + 1))

for n in range(1, 11):
    h = harmonique(n)
    print(f"H_{n:<2} = {h}  = {float(h):.10f}")

# 2. H_n n'est jamais entier pour n >= 2
for n in range(2, 200):
    assert harmonique(n).denominator != 1
print("H_n n'est jamais entier pour 2 <= n < 200")

# 3.
h20 = harmonique(20)
approche = sum(1/k for k in range(1, 21))
print(h20)                                  # 55835135/15519504
print(float(h20), approche, float(h20) - approche)

2. Preuve (theoreme de Kurschak). Soit \(n \ge 2\) et \(2^m\) la plus grande puissance de 2 inferieure ou egale a \(n\). On montre que \(2^m\) est le seul entier de \([1,n]\) divisible par \(2^m\) : en effet \(2\cdot 2^m = 2^{m+1} > n\) par maximalite. Multiplions \(H_n\) par le ppcm \(L\) de \(1,\dots,n\) : tous les termes \(L/k\) sont pairs sauf \(L/2^m\), qui est impair. La somme \(L\,H_n\) est donc impaire, alors que \(L\) est pair (puisque \(n \ge 2\)). Si \(H_n\) etait entier, \(L H_n\) serait pair : contradiction.

class Frac:
    """Fraction rationnelle reduite."""

    def __init__(self, num, den=1):
        if den == 0:
            raise ZeroDivisionError("denominateur nul")
        if den < 0:
            num, den = -num, -den          # le signe est porte par le numerateur
        d = gcd(abs(num), den)
        self.num, self.den = num // d, den // d

    def __add__(self, autre):
        return Frac(self.num*autre.den + autre.num*self.den, self.den*autre.den)

    def __sub__(self, autre):
        return Frac(self.num*autre.den - autre.num*self.den, self.den*autre.den)

    def __mul__(self, autre):
        return Frac(self.num*autre.num, self.den*autre.den)

    def __truediv__(self, autre):
        if autre.num == 0:
            raise ZeroDivisionError
        return Frac(self.num*autre.den, self.den*autre.num)

    def __eq__(self, autre):
        return self.num == autre.num and self.den == autre.den    # formes reduites

    def __lt__(self, autre):
        return self.num*autre.den < autre.num*self.den            # den > 0

    def __repr__(self):
        return f"{self.num}" if self.den == 1 else f"{self.num}/{self.den}"

    def __float__(self):
        return self.num / self.den

a, b = Frac(1, 3), Frac(1, 6)
assert a + b == Frac(1, 2)
assert a * b == Frac(1, 18)
assert a / b == Frac(2, 1)
assert Frac(2, 4) == Frac(1, 2)
assert Frac(1, 3) < Frac(1, 2)
assert repr(Frac(4, 2)) == "2"
print(sum((Frac(1, k) for k in range(2, 11)), Frac(0)))    # 1627/2520

La reduction systematique par le pgcd dans le constructeur garantit une forme canonique unique, ce qui rend la comparaison == triviale et evite l'explosion de la taille des entiers.

import random
from math import sqrt, pi

random.seed(7)

def marche_1d(n):
    """Position finale apres n pas."""
    return sum(random.choice((-1, 1)) for _ in range(n))

def statistiques_1d(n, essais=20000):
    positions = [marche_1d(n) for _ in range(essais)]
    m = sum(positions) / essais
    m2 = sum(x*x for x in positions) / essais
    return m, m2

for n in (10, 100, 1000):
    m, m2 = statistiques_1d(n, 20000)
    print(f"n = {n:>5} : E[X] ~ {m:>7.3f} (theorie 0), "
          f"E[X^2] ~ {m2:>9.2f} (theorie {n})")

n =    10 : E[X] ~   0.010 (theorie 0), E[X^2] ~      9.94 (theorie 10)
n =   100 : E[X] ~  -0.036 (theorie 0), E[X^2] ~     99.55 (theorie 100)
n =  1000 : E[X] ~   0.163 (theorie 0), E[X^2] ~   1002.68 (theorie 1000)

2. Theorie. Posons \(X_n = \sum_{i=1}^n \varepsilon_i\) avec \(\varepsilon_i\) independantes valant \(\pm1\) avec probabilite \(1/2\). Alors \(\mathbb{E}[\varepsilon_i] = 0\) donc \(\mathbb{E}[X_n] = 0\). Et \[ \mathbb{E}[X_n^2] = \sum_i \mathbb{E}[\varepsilon_i^2] + \sum_{i \ne j} \mathbb{E}[\varepsilon_i]\mathbb{E}[\varepsilon_j] = n + 0 = n, \] en utilisant l'independance. La distance typique a l'origine est donc \(\sqrt n\) : c'est la loi de diffusion, universelle en physique statistique.

def retour_avant(n, essais=20000):
    """Probabilite de repasser par 0 avant n pas."""
    c = 0
    for _ in range(essais):
        x = 0
        for k in range(n):
            x += random.choice((-1, 1))
            if x == 0:
                c += 1
                break
    return c / essais

for n in (10, 100, 1000):
    print(f"n = {n:>5} : retour avant n pas ~ {retour_avant(n):.4f}")

def marche_2d(n):
    x = y = 0
    for _ in range(n):
        dx, dy = random.choice([(1,0), (-1,0), (0,1), (0,-1)])
        x += dx; y += dy
    return sqrt(x*x + y*y)

for n in (100, 400, 1600):
    d = sum(marche_2d(n) for _ in range(5000)) / 5000
    print(f"n = {n:>5} : distance moyenne ~ {d:>7.3f}, sqrt(n) = {sqrt(n):.3f}, "
          f"rapport {d/sqrt(n):.3f}")

n =   100 : distance moyenne ~   8.874, sqrt(n) = 10.000, rapport 0.887
n =   400 : distance moyenne ~  17.712, sqrt(n) = 20.000, rapport 0.886
n =  1600 : distance moyenne ~  35.412, sqrt(n) = 40.000, rapport 0.885

Le rapport est stable autour de \(\sqrt{\pi}/2 \approx 0{,}886\) : c'est la valeur theorique de \(\mathbb{E}[\|X_n\|]/\sqrt n\) pour une marche 2D, la loi limite etant celle de Rayleigh.

from itertools import combinations, permutations
from math import comb, factorial, e, sqrt
import random

# 1. mains de 5 cartes contenant au moins un as
couleurs = "PCKT"
hauteurs = "789XVDRA"
jeu = [h + c for c in couleurs for h in hauteurs]
assert len(jeu) == 32

mains = list(combinations(jeu, 5))
print(len(mains), comb(32, 5))                  # 201376 201376
avec_as = sum(1 for m in mains if any(c[0] == "A" for c in m))
theorie = comb(32, 5) - comb(28, 5)             # total - sans aucun as
print(avec_as, theorie, f"{avec_as/len(mains):.4f}")   # 103096 103096 0.5120

# 2. voyageur de commerce, force brute
random.seed(3)
n = 8
villes = [(random.uniform(0, 100), random.uniform(0, 100)) for _ in range(n)]

def longueur(tour):
    d = 0.0
    for i in range(len(tour)):
        a, b = villes[tour[i]], villes[tour[(i+1) % len(tour)]]
        d += sqrt((a[0]-b[0])**2 + (a[1]-b[1])**2)
    return d

# on fixe la ville 0 : les (n-1)! permutations du reste suffisent
meilleur, dmin = None, float("inf")
for p in permutations(range(1, n)):
    tour = (0,) + p
    d = longueur(tour)
    if d < dmin:
        meilleur, dmin = tour, d
print(f"tour optimal : {meilleur}, longueur {dmin:.2f}")
print(f"{factorial(n-1)} permutations examinees")

# 3. derangements
def derangements(n):
    return sum(1 for p in permutations(range(n))
               if all(p[i] != i for i in range(n)))

print(f"{'n':>3}{'D(n)':>10}{'n!/e':>12}{'ecart':>10}")
for n in range(1, 10):
    d = derangements(n)
    print(f"{n:>3}{d:>10}{factorial(n)/e:>12.3f}{abs(d - factorial(n)/e):>10.3f}")

  n      D(n)        n!/e     ecart
  1         0       0.368     0.368
  2         1       0.736     0.264
  3         2       2.207     0.207
  4         9       8.829     0.171
  5        44      44.146     0.146
  6       265     264.873     0.127
  7      1854    1854.112     0.112
  8     14833   14832.899     0.101
  9    133496  133496.089     0.089

Le nombre de derangements vaut exactement \(D_n = n!\sum_{k=0}^{n}\frac{(-1)^k}{k!}\), qui est l'entier le plus proche de \(n!/e\) (l'ecart tend vers 0). Consequence celebre : la probabilite qu'une permutation aleatoire n'ait aucun point fixe tend vers \(1/e \approx 0{,}368\), independamment de \(n\).

4. Limites de la force brute.

Le voyageur de commerce est NP-difficile : au-dela d'une quinzaine de villes, il faut recourir a la programmation dynamique (\(O(n^2 2^n)\), utilisable jusqu'a \(n \approx 25\)), a la programmation lineaire en nombres entiers, ou a des heuristiques (2-opt, recuit simule) qui donnent de bonnes solutions sans garantie d'optimalite.

import random
from math import gcd

class Congruentiel:
    """Generateur congruentiel lineaire x -> (a x + c) mod m."""

    def __init__(self, a, c, m, graine=1):
        self.a, self.c, self.m, self.x = a, c, m, graine

    def suivant(self):
        self.x = (self.a * self.x + self.c) % self.m
        return self.x

    def uniforme(self):
        return self.suivant() / self.m

randu = Congruentiel(65539, 0, 2**31, 1)
bon   = Congruentiel(1664525, 1013904223, 2**32, 1)   # Numerical Recipes

print([randu.suivant() for _ in range(5)])

2. Theoreme de Hull-Dobell (1962). Le generateur \(x_{n+1} = (a x_n + c) \bmod m\) avec \(c \ne 0\) est de periode maximale \(m\) (pour toute graine) si et seulement si :

1. \(\operatorname{pgcd}(c, m) = 1\) ;
2. \(a - 1\) est divisible par tout facteur premier de \(m\) ;
3. si \(4 \mid m\), alors \(4 \mid (a-1)\).

def periode_maximale(a, c, m):
    """Verifie les trois conditions de Hull-Dobell."""
    if gcd(c, m) != 1:
        return False
    n, p = m, 2
    facteurs = set()
    while p * p <= n:
        while n % p == 0:
            facteurs.add(p); n //= p
        p += 1
    if n > 1:
        facteurs.add(n)
    if any((a - 1) % p != 0 for p in facteurs):
        return False
    return not (m % 4 == 0 and (a - 1) % 4 != 0)

assert periode_maximale(1664525, 1013904223, 2**32)
assert not periode_maximale(65539, 0, 2**31)      # c = 0 : jamais maximal

3. Le defaut de RANDU. On verifie la relation exacte, valable pour tout \(n\) :

\[ x_{n+2} = 6x_{n+1} - 9x_n \pmod{2^{31}}. \]

Preuve : \(a = 65539 = 2^{16}+3\), donc \(a^2 = 2^{32} + 6\cdot2^{16} + 9 \equiv 6\cdot 2^{16} + 9 = 6a - 9 \pmod{2^{31}}\) (car \(2^{32} \equiv 0\) et \(6a = 6\cdot2^{16}+18\)). Ainsi \(x_{n+2} = a^2 x_n \equiv (6a-9)x_n = 6x_{n+1} - 9x_n\).

g = Congruentiel(65539, 0, 2**31, 1)
x = [g.suivant() for _ in range(20)]
for n in range(18):
    assert (x[n+2] - 6*x[n+1] + 9*x[n]) % 2**31 == 0
print("relation x[n+2] = 6x[n+1] - 9x[n] verifiee : RANDU est defectueux")

Consequence geometrique : tous les triplets \((x_n, x_{n+1}, x_{n+2})\) se placent sur seulement 15 plans de l'espace, au lieu d'etre repartis dans le cube. RANDU, diffuse par IBM dans les annees 1960 et utilise pendant deux decennies, a invalide un nombre considerable de simulations scientifiques.

def test_chi2(generateur, n=100000, cases=10):
    """Test du chi2 d'uniformite sur [0,1[. Renvoie la statistique."""
    effectifs = [0] * cases
    for _ in range(n):
        effectifs[min(cases - 1, int(generateur() * cases))] += 1
    attendu = n / cases
    return sum((o - attendu)**2 / attendu for o in effectifs)

random.seed(0)
print(f"random Python : chi2 = {test_chi2(random.random):.2f}")
g1 = Congruentiel(1664525, 1013904223, 2**32, 12345)
print(f"congruentiel   : chi2 = {test_chi2(g1.uniforme):.2f}")
mauvais = Congruentiel(3, 0, 101, 1)
print(f"mauvais (m=101): chi2 = {test_chi2(mauvais.uniforme):.2f}")

random Python : chi2 = 5.82
congruentiel   : chi2 = 8.31
mauvais (m=101): chi2 = 254.13

Sous l'hypothese d'uniformite, la statistique suit approximativement une loi \(\chi^2\) a \(k-1 = 9\) degres de liberte, d'esperance 9 et de quantile 95 % egal a 16,92. Les deux premiers generateurs passent le test ; le troisieme (module 101, donc periode au plus 100) est rejete sans ambiguite.