On a retrouvé π au milieu du hasard : calcul d'une probabilité

Bonjour,

Je cherche à calculer numériquement la probabilité que deux nombres naturels tirés au hasard soient premiers entre eux.
Je me suis basé sur cet article : https://zestedesavoir.com/articles/3452/on-a-retrouve-p-au-beau-milieu-du-hasard/

Ne connaissant pas le Rust, j’ai essayé d’écrire un premier jet de code en PHP. Le PHP est un langage basé sur le C, de haut-niveau, idéal donc pour implémenter une première fois l’algorithme sans se compliquer la tâche avec toutes les subtilités langages du C. On perd en performance, mais mon objectif est d’abord d’écrire l’algorithme correctement avec un langage plus souple.

Idéalement, mon programme doit fournir en sortie la probabilité à calculer : $\ P = \dfrac{6}{\pi^2} \approx 0,607927...$.
Le problème, c’est qu’il renvoie $P = 0,6177$ en réalisant 10 000 fois l’expérience.

Auriez-vous une idée d’où pourrait venir le problème ?
Code ci-dessous.

<?php

/**
  *     Ce programme cherche à calculer numériquement
  *     la probabilité que deux nombres entiers positifs
  *     tirés au hasard soient premiers entre eux.
  *
  *     On s'attend à une valeur proche de 6/π² = 0,607927......
  */




// On décide de tirer les nombres entiers au hasard entre 0 et $_MAX inclus
$_MAX = 100000;


// On liste tous les nombres premiers depuis 2 jusqu'à $_MAX
$premiers = array();

for ($i = 2; $i <= $_MAX; ++$i)
{
    $j = 2;
    $s = sqrt($i);
    
    while ($j < $s && $i%$j != 0)
        ++$j;
    
    if ($j > $s)
        $premiers[] = $i;
}

// Affiche la liste de tous les nombres premiers entre 2 et $_MAX inclus
// var_dump($premiers);


/* --------------------------------------------------------- */



$N      = 0;            // Nombre de succès de l'expérience
$essais = 10000;        // Nombre de fois que l'expérience est réalisée

for ($i = 0; $i < $essais; ++$i)
{
    $a      =   mt_rand(0, $_MAX);
    $b      =   mt_rand(0, $_MAX);
    $min    =   min($a, $b);
    
    $in = true;
    
    for ($j = 0; $premiers[$j] <= $min; ++$j)
    {
        $c1 = $a % ($premiers[$j]);
        $c2 = $b % ($premiers[$j]);
        
        if ($c1 + $c2 == 0) // Les nombres ne sont pas premiers entre eux
        {
            $in = false;
            break;
        }
    }
    
    if ($in)
        $N++;
}

$p = $N/$essais;        // Probabilité expérimentale
echo 'P = '.$p;

Effectivement, j’ai un peu zappé la partie de l’article qui parlait des maths et j’ai zappé ce paragraphe :

La valeur de la probabilité est juste, certes, c’est bien 6/π26/\pi^{26/π2. Mais le fait est que cette probabilité vaut 6/π26/\pi}26/π2 uniquement si on tire un nombre dans l’ensemble des entiers positifs non nuls (cf. la formulation sous forme de limite plus haut, si vous connaissez cette notion).

Problème : la convergence est extrêmement lente. Il faut énormément augmenter la taille de l’échantillon pour avoir une précision de 6/π26/\pi^26/π2 qui augmente un peu. Et quand ce n’est pas le cas, on peut toujours faire des centaine de milliards de tirages : la probabilité ne permet pas d’atteindre π, mais une valeur qui s’en approche.

Ça reste sympa à essayer

Si tu aimes ce genre de trucs, tu as le Project Euler qui pose des problèmes plus ou moins ardus qu’on peut résoudre numériquement. À part les premiers problèmes, ça demande un peu de réflexion plutôt qu’une implémentation bourrin de l’énoncé.

Dans la même idée mais en français, tu as le Défi Turing.