On a retrouvé π au beau milieu du hasard !

Alors bien sûr on ne va pas tirer des nombres au hasard jusqu’à tomber par (énorme) chance sur π ! On va être un peu plus intelligent que ça… en tirant deux nombres. Attendez, vous allez comprendre…

Une histoire de nombres premiers entre eux

Tout nombre entier peut s’écrire comme le produit d’autres nombres : jusque-là, rien de bien compliqué. Si on prend quelques nombres et qu’on les écrit de la sorte, on obtient quelque chose comme ça :

$\begin{array}{rcl}12 &=& 3 \times 4\\11 &=& 11 \ \ \text{(il n'y a rien de mieux)}\\100 &=& 4 \times 5 \times 5\\40 &=& 2 \times 2 \times 2 \times 5\end{array}$

Si on regarde les nombres deux par deux, on aura dans certains cas des facteurs en commun (par exemple, 100 et 40 partagent 2 et 5). Dans d’autres cas, il n’y en a pas (11 et 12 ne partagent rien) : on dit alors que ces nombres sont premiers entre eux.

Quand on tire deux nombres au hasard, il y a donc deux possibilités : ils peuvent être premiers entre eux, ou ne pas l’être.

Par exemple, j’ai tiré avec deux dés à cent faces les nombres 77 et 87 : 87 s’écrit 3 × 29 et 77 = 7 × 11 — ils sont premiers entre eux, ne partageant aucun facteur commun. J’ai relancé mes dés pour trouver 2 et 14 : 2 est premier et 14 s’écrit 2 × 7 — ils ne sont donc, eux, pas premiers entre eux.

Expérimentons : tirons plein de nombres !

Mais so what, me diriez-vous, où trouve-t-on π là-dedans ? Dans les probabilités ! Car en tirant deux nombres au hasard, on n’a pas autant de chance d’avoir deux nombres premiers entre eux que de nombres qui ne le sont pas. Lancez plein de dés, vous verrez : vous obtiendrez plus souvent des nombres premiers entre eux que le contraire. Cette probabilité se calcule, et devinez quoi ? Elle est liée à π !

En effet (et je le justifierai plus loin),

Comment vérifier ça ? Un moyen simple est de… tirer plein de nombres. En réalisant l’expérience tout plein de fois et en comptant quel pourcentage de couples premiers entre eux on a obtenu, on approche la probabilité expérimentalement .

Oh, on pourrait tirer des milliers de couples de dés et vérifier que les nombres obtenus sont premiers entre eux… mais ce serait un peu longuet.¹ Comme j’ai un peu la flemme, j’ai décidé d’écrire un petit programme qui fait tout le boulot pour moi, ce qui me permet de faire des millions (ou milliards si l’envie m’en prend) de tirages sans effort. Ha !

Commencer doucement…

Comme il faut savoir commencer petit, j’ai démarré avec un millier de tirages. J’ai obtenu 643 couples de nombres premiers entre eux, ce qui fait une probabilité expérimentale, basée sur ce tirage, de :

$P = \frac{643}{1000} = 0,643$

On aurait donc $6 / \pi^2 = 0,643$ ? Si oui, ça donnerait :

$\pi^2 = 6 \times \frac{1}{0,643} \approx 9,331$

et donc :

$\pi \approx \sqrt{6 \times \frac{1}{0,643}} \approx 3,054$

Je vous l’accorde, la précision n’est pas extraordinaire… mais tout de même, pour à peine un millier de paires de nombres tirés au hasard, retrouver π avec une précision de 2 %, c’est pas si mal ! Cela dit, on peut faire mieux.

Comme ce n’est pas moi qui lance les dés, autant se faire plaisir : partons sur un million de tirages. En moins d’une seconde, c’est bon, et on a trouvé 607 780 couples de nombres premiers entre eux. En appliquant le même calcul que plus haut, avec une probabilité expérimentale de $\frac{607\,780}{1\,000\,000} = 0,607\,780$, on trouve :

$\pi \approx 3,141\,972$

Voilà qui n’est pas si mal ! Certes, dès la quatrième décimale on perd la précision (elle devrait être un 5, la vraie valeur commençant par 3,14159) — mais encore une fois, on tire cette valeur du hasard complet, donc une précision de 0,01 % c’est relativement acceptable. Surtout qu’en pratique, prendre 3,14 comme approximation de π, ça suffit à la majorité des usages…

…mais finir gourmand

Dans un brin de folie (et surtout car ça travaille tout seul pendant que je fais autre chose), j’ai demandé dix milliards de tirages (oui, 10 000 000 000 de tirages). Après seulement trois heures et quart de calcul (ahrem)², 6 079 281 509 couples de nombres premiers entre eux avaient été tirés, ce qui nous permet d’estimer que π a pour valeur…

$\pi \approx 3,141\,589\,943$

C’est déjà mieux : on a quatre décimales justes (la vraie valeur étant π = 3,141592…) !

extern crate gcd;
extern crate rand;

use std::env;

use gcd::Gcd;
use rand::Rng;

use std::f64::consts::PI;

const DEFAULT_DRAWS: u64 = 1_000_000;
const MAX_NUMBER: u128 = std::u128::MAX;

#[inline(always)]
fn are_randoms_coprime<R: Rng>(rng: &mut R) -> bool {
    rng.gen_range(1, MAX_NUMBER).gcd(rng.gen_range(1, MAX_NUMBER)) == 1
}

fn main() {
    let args: Vec<String> = env::args().collect();
    let draws: u64 = args.get(1).unwrap_or(&String::new()).parse().unwrap_or(DEFAULT_DRAWS);

    let mut rng = rand::thread_rng();

    let coprimes = (1..draws)
        .filter(|_| are_randoms_coprime(&mut rng))
        .count() as f64;

    let pi = (6_f64 / (coprimes / draws as f64)).sqrt();
    let precision = (pi - PI).abs() / PI * 100_f64;

    println!("{} coprimes in {} draws", coprimes, draws);
    println!("π ≈ {} (precision {} %)", pi, precision);
}

$ cargo run --release --quiet               
607780 coprimes in 1000000 draws
π ≈ 3.141972812647824 (precision 0.01210083864935896 %)

$ cargo run --release --quiet -- 10000000000
6079281509 coprimes in 10000000000 draws
π ≈ 3.14158994300916 (precision 0.0000862804128715215 %)

Bon allez, il est temps de rentrer un peu dans les mathématiques afin de comprendre d’où vient cette probabilité !

Partons de ce qu’on cherche. Avec $a$ et $b$ deux entiers naturels quelconques (ceux qu’on tire au dé), on cherche la probabilité qu’ils soient premiers entre eux, autrement dit que leur plus grand diviseur commun (PGCD) soit égal à 1.

Pour la suite, on va noter $p_1$ cette probabilité (vous comprendrez assez vite pourquoi ce 1 en indice), et on la note ainsi :

$p_1 = P\left(\text{pgcd}(a, b) = 1\right)$

Une première observation que l’on peut faire, c’est qu’on peut avoir des nombres arbitrairement grands ici (et même potentiellement très grands, la seule limite étant l’infini). Le plus grand diviseur commun de deux nombres peut donc être, lui aussi, aussi grand que l’on veut. $\text{pgcd}(a, b)$ peut valoir 1, ou 2, ou 3, ou… ou $k$, ou $k + 1$, ou… jusqu’à l’infini.

Mais alors, il est certain que le PGCD a l’une de ces valeurs. Donc si on somme la probabilité que le PGCD soit 1, et soit 2, et 3, et…, et $k$, et…, jusqu’à l’infini, on doit trouver 1 (la probabilité certaine).

Autrement dit, si on note $p_k$ la probabilité que le PGCD de $a$ et $b$ soit $k$ :

$p_k = P\left(\text{pgcd}(a, b) = k\right)$

…alors on sait que :

$p_1 + p_2 + \cdots + p_k + \cdots = 1$

Ce qu’on peut écrire d’une façon un peu plus compacte :

$\sum\limits_{k = 1}^\infty p_k = 1$

On va donc tenter de déterminer la probabilité que le PGCD de $a$ et $b$ soit $k$, c’est peut-être plus simple…

Alors, dans quels cas est-ce que $\, \text{pgcd}(a, b) = k\,$ ?

Réfléchissons. Pour que le PGCD de $a$ et $b$ soit $k$, il faut respecter deux critères en même temps :

• $k$ doit diviser $a$ et $b$ (c’est mieux… sinon il aura du mal à être le plus grand diviseur) ; et
• il ne doit pas y avoir de diviseur plus grand que $k$.

Pour trouver la probabilité $p_k$, il nous faut donc trouver la probabilité de ces deux événements et les multiplier (car on veut les deux en même temps, et que ces événements sont indépendants les uns des autres).

Probabilité que $k$ divise $a$ et $b$

Quand on y pense,

• 2 divise la moitié des nombres (tous les nombres pairs) ;
• 3 en divise un tiers (3, 6, 9, 12, …) ;
• 4 en divise un quart…

En généralisant, un nombre quelconque $k$ divise un $k$-ième des nombres… la probabilité qu’il en divise un en particulier est donc $1/k$ — par exemple, 2 divise la moitié des nombres, donc en prenant un nombre au hasard, il y a une chance sur deux qu’il soit divisible par deux (autrement dit, pair).

Ici on veut que $k$ divise $a$ et $b$. La probabilité qu’on cherche est donc celle que $k$ divise $a$ multipliée par celle que $k$ divise $b$. Les deux nombres n’ayant pas particulièrement de différence de traitement, cette probabilité est bien sûr la même, et la probabilité que $k$ divise $a$ et $b$ vaut donc $1/k^2$.

Probabilité qu’il n’y ait pas de diviseur plus grand que $k$

On veut qu’il n’y ait pas de diviseur commun plus grand que $k$ : c’est équivalent au fait qu’il n’y ait pas de diviseur plus grand que 1 si on divise par $k$ (autrement dit que le $\text{pgcd}\left(\frac{a}{k}, \frac{b}{k}\right)$ soit 1)¹.

Mais attendez, la probabilité que deux nombres aient un PGCD de 1… on l’a déjà nommée plus haut, c’est $p_1$ !

Revenons donc à ce qu’on cherchait…

Si on rassemble ces deux points, on peut donc affirmer que $p_k$, la probabilité que le PGCD de $a$ et $b$ vaille $k$, vaut le produit de ces deux valeurs trouvées, donc :

$p_k = \frac{1}{k^2} \times p_1 = \frac{p_1}{k^2}$

Or, on avait plus haut une relation avec $p_k$ ! Et si on essayait d’injecter cette nouvelle valeur ? On obtient :

$\sum\limits_{k = 1}^\infty p_k = p_1 + \frac{p_1}{2^2} + \frac{p_1}{3^2}+ \frac{p_1}{4^2} + \cdots = 1$

Or $p_1$ est une constante, on peut donc factoriser pour la sortir de la somme :

$p_1 \times \left(1 + \frac{1}{2^2} + \frac{1}{3^2}+ \frac{1}{4^2} + \cdots\right) = 1$

Et donc en réorganisant :

$p_1 = \frac{1}{1 + \frac{1}{2^2} + \frac{1}{3^2}+ \frac{1}{4^2} + \cdots} = \frac{1}{\sum\limits_{k = 1}^\infty \frac{1}{k^2}}$

On avance ! Il ne nous reste plus qu’à déterminer ce mystérieux terme : la somme des inverses des carrés de tous les nombres entiers positifs non nuls… mais il se trouve que c’est un problème connu, le problème de Bâle, résolu par Euler en 1741. Et il se trouve que justement… on y retrouve… π !

$\sum\limits_{k = 1}^\infty \frac{1}{k^2} = \frac{\pi^2}{6}$

De là, il est trivial de conclure que :

$p_1 = \frac{1}{\pi^2 / 6} = \frac{6}{\pi^2}$

Et voilà.

Alors oui, on peut calculer π avec cette méthode, en tirant plein de nombres au dé… mais est-elle efficace ? En pratique, pas vraiment. Ou plutôt, vraiment pas.

Le problème, c’est que pour que ça fonctionne correctement, il faudrait pouvoir tirer des nombres sur l’ensemble des nombres entiers positifs qui existent, donc jusqu’à l’infini. Et l’infini, c’est grand.

Si on ne tire pas des nombres sur tout ça, la précision s’en ressent beaucoup (vraiment beaucoup). Avec un gros dé on tire entre quoi, un et deux-cent ? C’est ridiculement petit… et on obtient au mieux zéro ou une décimale. Et même avec le logiciel que j’ai écrit pour faire les tirages à ma place, en prenant le plus grand nombre possible facilement, un nombre à trente-huit chiffres tout de même¹, tellement grand qu’il est difficile de se l’imaginer — on est très, très loin de l’infini, et on peine à dépasser quatre ou cinq décimales.

Lorsqu’on calcule π avec une très grande précision, espérant des milliards de chiffres après la virgule, on utilise des méthodes très différentes et bien plus efficaces, donc on pourrait éventuellement parler s’il y a de l’intérêt parmi mes lecteurs. Mais ces méthodes n’ont pas la poésie de tirer π du hasard le plus total…