4. Convergences
La méthode des fractions
La documentation officielle de Python est ici particulièrement éclairante puisque c’est justement une approximation de \(\pi\) qui a été choisie pour illustrer la fonction limit.denominator() du module fractions :
# Réf. https://docs.python.org/3/library/fractions.html#module-fractions
from fractions import Fraction
Fraction('3.1415926535897932').limit_denominator(1_000)
Fraction(355, 113)
Nous retrouvons Milü 密率, la fameuse approximation de Zu Chongzhi.
Il n’est pas difficile de construire, autour de cette fonction, un programme qui aura pour objet de répertorier, pour une plage d’entiers donnée, les fractions donnant les approximations de \(\pi\) les plus intéressantes.
L’idée est d’incrémenter une variable qui va représenter la limite du dénominateur, de façon à obtenir différentes fractions proches de \(\dfrac{\pi}{4}\) pour un rang déterminé d’entiers.
Car c’est bien \(\dfrac{\pi}{4}\) que l’on va chercher à représenter sous forme de fraction, et non pas \(\pi\) lui-même. Ainsi, le numérateur nous donnera directement le nombre de points placés à l’intérieur du cercle, tandis que le dénominateur indiquera le nombre total de points.
On devrait donc rencontrer parmi ces fractions remarquables 355 / 452 et non pas 355 / 113 puisque c’est \(4\times\dfrac{355}{452}\) qui nous donnera l’approximation de \(\pi\). Quant au pas d’incrément (step), il doit pouvoir être modifié de façon à filtrer les résultats.
from fractions import Fraction
from decimal import Decimal as D, getcontext
from tabulate import tabulate # présentation sous forme de tableau
# Cf. https://pypi.org/project/tabulate/
from pi_string import pi # module local (voir plus haut, partie II)
precision = 30 # la précision peut être augmentée si nécessaire
pi_ref = pi(precision) # π sous forme de string avec 30 décimales
getcontext().prec += 4 # + de précision pour les calculs
pi_quart = (D(pi_ref)/4).quantize(D(pi_ref)) # π/4 avec 30 décimales
den_lim = 1
results = []
def get_digits(pi_mc):
"""Renvoie le nombre de décimales exactes pour une valeur approchée de π"""
k = 0
while pi_mc[:k] in pi_ref :
k += 1
return max(k - 3, 0)
# Programme principal
for i in range(1, 1_000_000, 10_000):
"""
Calcule pour une plage d’entiers donnée les fractions qui donnent les
valeurs les plus proches de π.
"""
fpi = Fraction(pi_quart).limit_denominator(i)
num, den = fpi.numerator, fpi.denominator
if den > den_lim:
den_lim = den
pi_mc = str((4 * D(num) / D(den)).quantize(D(pi_ref)))
results.append((den, num, pi_mc, get_digits(pi_mc)))
# tri, formatage et affichage des résultats
results.sort(key= lambda x: x[3]) # tri par nombre de décimales
headers=['Total', 'Cible', 'Valeur approchée de π', 'Nb déc.']
print(tabulate(results,
headers,
disable_numparse=True,
colalign=(4*'right'.split())))
Total Cible Valeur approchée de π Nb déc.
------- ------- -------------------------------- ---------
452 355 3.141592920353982300884955752212 6
19655 15437 3.141592470109386924446705672856 6
29599 23247 3.141592621372343660258792526774 7
33215 26087 3.141592653921421044708715941592 9
198838 156167 3.141592653315764592281153501845 9
232053 182254 3.141592653402455473534063338978 9
265268 208341 3.141592653467436705520454785349 9
298483 234428 3.141592653517955796477521332873 10
331698 260515 3.141592653558357300918305205337 10
364913 286602 3.141592653591403978482542414219 10
Le calcul est très rapide et permet d’obtenir quasi-instantanément des fractions approchant \(\pi\) avec une grande précision pour des valeurs élevées de n. Si l’on calcule les fractions possibles jusqu’à un nombre de points de 1 milliard avec un pas d’incrément de 100 millions (ce qui filtre beaucoup de résultats), on obtient :
Total Cible Valeur approchée de π Nb déc.
--------- --------- -------------------------------- ---------
27235615 21390802 3.141592653589794098646202775299 14
183749173 144316263 3.141592653589793299369026275835 16
340262731 267241724 3.141592653589793235392564929480 17
La valeur donnant le résultat le plus proche de \(\pi\) obtenue pour un nombre total de points inférieur à un milliard (109) est celle qui correspond à 267 241 724 points placés dans la cible sur un total de 340 262 731. La précision atteinte est alors de 17 décimales !
Le plafond des 20 décimales exactes est atteint (et dépassé) pour n = 26 805 949 036 :
Total Cible Valeur approchée de π Nb déc.
----------- ----------- -------------------------------- ---------
26805949036 21053343141 3.141592653589793238462381742774 21
Grâce à la puissance conjuguée des modules fractions et decimal, la précision théorique que l’on peut atteindre est quasiment sans limite.
Pour aller plus loin, il suffit d’augmenter la valeur de la variable precision que nous avons ici arbitrairement fixée à 30.
Ainsi, en choisissant une précision \(> 50 \) et en faisant une recherche sur la plage range(int(1e25), int(1e25),int(1e24)), on trouve facilement que la première fraction, dans l’ordre des n croissants, à donner une approximation de \(\pi\) avec 50 décimales exactes, est obtenue à partir des valeurs suivantes :
Total : 33 297 080 577 553 090 318 600 632 Cible : 26 151 465 932 107 044 561 886 949
Il est facile de le vérifier :
from decimal import Decimal as D, getcontext
getcontext().prec = 60
4 * D(26_151_465_932_107_044_561_886_949) / D(33_297_080_577_553_090_318_600_632)
Decimal('3.14159265358979323846264338327950288419716939937510414085630')
Cela nous donne bien une valeur de \(\pi\) exacte jusqu’à la cinquantième décimale :
from pi_string import pi
print(pi(50))
3.14159265358979323846264338327950288419716939937510
On notera que la fonction range() n’accepte que les entiers, ce qui explique que les nombres utilisant la notation scientifique (1e25 pour 1025, considérés comme des flottants doivent être passés comme arguments sous la forme int(1e25). Il est aussi possible (quoique fastidieux) de noter 1 suivi de 25 zéros.
Nous nous sommes ici appuyés sur la méthode limit_denominator() de la classe Fraction du module fractions qui permet de trouver des approximations rationnelles de nombres flottants.
Une autre approche possible est d’utiliser la puissance de NumPy pour chercher par itération les tirages pour lesquels l’espérance de la loi binomiale \(\mathit{B}\left(n, \dfrac{\pi}{4}\right)\) est proche d’un entier.
La méthode des entiers
Intuitivement, on comprend que puisque la loi binomiale \(\mathit{B}\left(n, \dfrac{\pi}{4}\right)\) a pour espérance \(n\times\dfrac{\pi}{4}\), si cette espérance est proche d’un entier, alors cet entier, une fois multiplié par 4 et divisé par n sera proche de \(\pi\).
Prenons l’exemple de l’approximation de Zu Chongzhi citée plus haut, 355 / 113 qui nous donne 6 décimales exactes. C’est la fraction 355 / 452 que nous chercherons puisque, une fois multipliée par 4, elle donnera le résultat attendu.
L’espérance de la loi binomiale pour cette valeur de n est \(n\times\dfrac{\pi}{4}\), soit :
from math import pi
452 * pi / 4
354.9999698556466
Comme on le voit, ce chiffre est très proche de 355. C’est la raison pour laquelle la fraction 355 / 113 est une aussi bonne approximation de \(\pi\).
L’idée sera donc d’incrémenter une variable, de calculer le produit de cette variable par \(\dfrac{\pi}{4}\) (probabilité de la loi binomiale), et de filtrer les résultats pour ne garder que les valeurs les plus proches d’un entier.
Le programme ci-dessous est fondé sur ce principe. Il permet de calculer les \(\pi_{mc}\) pour des tranches d’entiers données. Comme nous nous intéressons surtout au nombre de décimales exactes, une procédure de comptage a été ajoutée pour n’afficher que les résultats qui présentent un nombre minimal de décimales exactes.
from tabulate import tabulate
import numpy as np
from decimal import Decimal as D, getcontext
from pi_string import pi # module local (voir plus haut)
###############################################################################
# Paramètres modifiables :
alpha = 1 # borne inférieure de la tranche d’entiers
omega = 1_000_000 # borne supérieure
dec_min = 10 # nombre minimum de décimales exactes recherchées
###############################################################################
pi_ref = pi(30) # valeur de π de référence (30 décimales)
getcontext().prec += 4 # augmentation de la précision pour calculs
prob = np.longdouble(pi_ref) / 4 # probabilité de la loi binomiale (π/4)
n_prox = 100 # nombre d’éléments analysés dans une tranche donnée
tabsize = int(1e6) # taille des tableaux
results = []
def get_digits(pi_mc):
"""Renvoie le nombre de décimales exactes pour une valeur approchée de π"""
k = 0
while pi_mc[:k] in pi_ref :
k += 1
return max(k - 3, 0)
def checking_arrays(n_start, n_end):
"""
Récupère les *n_prox* (= 100 par défaut) éléments d’un tableau d’entiers
dont le produit par *prob* (= π/4) est le plus proche (par excès ou
par défaut) d’un entier, calcule le produit de cet entier par π/4 et
vérifie le nombre de décimales exactes de π exactes que compte le résultat,
enfin l’ajoute à la liste *result* si le nombre de décimales exactes est
supérieur ou égal à *n_dec*.
"""
tab = np.arange(n_start, n_end)
esp = tab * prob
s_int = np.rint(esp)
s_id = np.argpartition(abs(s_int - esp), n_prox)
s_tab = tab[s_id[:n_prox]]
s_ins = s_int[s_id[:n_prox]].astype(np.int64)
for ins, nx in sorted(zip(s_ins, s_tab)):
pi_mc = str((4 * ins / D(int(nx))).quantize(D(pi_ref)))
n_dec = get_digits(pi_mc)
if n_dec >= dec_min:
results.append((nx, ins, pi_mc, n_dec))
# programme principal
floor, remain = divmod(omega - alpha, tabsize)
loops = [tabsize for i in range(floor)]+([remain] if remain else [])
for i, loop in enumerate(loops):
start = alpha + i*tabsize
end = start + loop
checking_arrays(start, end)
# tri, formatage et affichage des résultats
if results:
results.sort(key= lambda x: x[3]) # tri par nombre de décimales
headers=['Total', 'Cible', 'Valeur approchée de π', 'Nb déc.']
print(tabulate(results,
headers,
disable_numparse=True,
colalign=(4*'right'.split())))
else:
print(f"Aucun résultat de {dec_min} décimales pour n \u2208 "
f"[{alpha};{omega}].")
Total Cible Valeur approchée de π Nb déc. ------- ------- -------------------------------- --------- 298483 234428 3.141592653517955796477521332873 10 331698 260515 3.141592653558357300918305205337 10 364913 286602 3.141592653591403978482542414219 10 596966 468856 3.141592653517955796477521332873 10 630181 494943 3.141592653539221271349025121354 10 663396 521030 3.141592653558357300918305205337 10 696611 547117 3.141592653575668486429298417624 10 729826 573204 3.141592653591403978482542414219 10 862234 677197 3.141592653502413497959950547067 10 895449 703284 3.141592653517955796477521332873 10 928664 729371 3.141592653532386309795577302447 10 961879 755458 3.141592653545820212313607012940 10 995094 781545 3.141592653558357300918305205337 10
Bien que ce programme soit un peu plus lent que le précédent, il donne bien plus de résultats. Le module fractions réduit automatiquement les fractions, alors qu’ici elles sont présentes sous une forme non simplifiée pour autant qu’elles soient dans la tranche analysée.
Ainsi, en lançant une recherche sur la tranche [1; 10 000] avec un minimum demandé de 6 décimales, on obtient plus d’une vingtaine de résultats, là où le programme fondé sur le module fractions ne renvoie que l’approximation de Zu Chongzhi comme résultat à 6 décimales.
Ces deux programmes complémentaires permettent de trouver des couples d’entiers n et nb_cible tels que \(4\times\dfrac{nb\_cible}{n}\) donne un résultat proche de \(\pi\) pour la précision demandée.
Les couples ainsi obtenus sont dans un rapport idéal, donnant le maximum de précision au calcul de \(\pi\) par la méthode de Monte-Carlo.
⁂
Probabilités
Intuitivement, il semble évident que la probabilité de « tomber » sur l’une de ces paires convoitées est plus grande si le nombre de points est petit.
Pour en faire l’expérience, lançons au hasard le générateur de NumPy pour voir combien de boucles il lui faut pour arriver, en tirant au hasard 452 points, à en placer exactement 355 à l’intérieur du cercle, donnant ainsi le ratio de Zu Chongzhi :
import numpy as np
rng = np.random.default_rng()
count = 1
def get_in_circle(n):
tab = rng.random((n, 2))
return np.sum(tab[:, 0]**2 + tab[:, 1]**2 < 1)
while get_in_circle(452) != 355:
count += 1
print("Nombre de boucles :", count)
Nombre de boucles : 13
Le résultat est différent à chaque fois qu’on lance le programme et peut varier de 1 à plusieurs dizaines. En outre, il n’est pas reproductible puisque le seed du PRNG n’a pas été initialisé.
Mais on sait que la probabilité d’obtenir k succès (avec \(0\le k\le n\)) pour n épreuves est donnée par la formule :
$$\mathbb{P}(X=k) = \binom{n}{k} p^k(1-p)^{n-k}$$
Sachant que p = \(\dfrac{\pi}{4}\) que n = 452, on peut calculer P(X = 355) :
from scipy import stats, pi
B = stats.binom(452, pi/4)
B.pmf(355)
0.045665083662869
À chaque tirage (ou boucle), on a donc environ 4,6 chances sur 100 d’obtenir le résultat espéré : 355 points à l’intérieur du cercle.
Essayons d’en savoir plus…
- Quelle est la probabilité d’obtenir 355 avant d’avoir effectué une dizaine de boucles ?
- Combien faudra-t-il de boucles pour avoir 9 chances sur 10 de parvenir au résultat ?
Pour répondre à ces deux questions, il faut faire intervenir une seconde loi binomiale, dont le premier paramètre, nb (le nombre de boucles), est inconnu et dont la probabilité de succès est cette fois égale à 0,046 (B.pmf(355)).
- Obtenir 355 avant d’avoir effectué nb boucles signifie qu’il faut au moins une fois obtenir 355 avant d’avoir effectué nb tirages.
Si l’on appelle Y la variable qui compte les succès, on a donc :
\(P(Y > 0) = 1 - P(Y = 0)\). Et \(P(Y = 0) = (1 - p)^{nb}\) puisque aucun succès signifie que 355 ne sort pas nb fois de suite.
Cela revient à calculer \(P(Y > 0) = 1 - (1 - p)^{nb}\), où nb est le nombre de boucles et \(p \approx 0,046\).
La probabilité d’obtenir 355 avant d’avoir fait 10 boucles est donc de 1 - (1 - 0,046)10, soit environ 0,38.
Il y a 38 % de chances d’obtenir au moins une fois 355 avant d’avoir parcouru 10 boucles. - Ici, nb est l’inconnue, et l’on connaît \(P(Y > 0) = 0,9\). L’équation à résoudre est donc :
\(1 - (1-p)^{nb} = 0,9 \Leftrightarrow nb = \dfrac{log(1 - 0,9)}{log(1 - p)}\) où \(p \approx 0,046\).
On trouve nb = 50 (en prenant la forme non arrondie de p), puisque nb doit être un entier et qu’il faut ici arrondir le résultat à l’entier supérieur.
Si l’on fait tourner un programme répétant des séquences aléatoires de 452 points, il y aura donc 9 chances sur 10 de placer au moins une fois exactement 355 points à l’intérieur du cercle avant d’avoir lancé 50 simulations.
Généralisons ces résultats et écrivons une fonction qui, pour un nombre total de points donnés (n) et un nombre de points placés dans la cible (nb_cible) nous donnera, si l’on indique combien de fois la séquence est répétée (nb), la probabilité (chance) de succès, sachant qu’on entend par « succès » arriver au moins une fois à placer exactement ce nombre de points dans la cible.
À l’inverse, pour une probabilité donnée, la fonction devra renvoyer le nombre de boucles à parcourir pour avoir cette chance.
from math import pi, log, ceil
from scipy.stats import binom
def prob(n, nb_cible, nb=None, chance=None):
B = binom(n, pi/4)
p = B.pmf(nb_cible)
if nb and chance is None:
return f'chance = {1 - (1 - p)**nb:.6f}'
if chance and nb is None:
return f'nb = {ceil(log(1 - chance) / log(1 - p))}'
Pour savoir combien le programme doit parcourir de boucles pour avoir 95 % de chances d’arriver à placer 355 points dans la cible sur un total de 452, il suffit de demander :
prob(452, 355, chance=95/100)
'nb = 65'
En parcourant 65 boucles — et donc en procédant à 65 tirages — on a donc une quasi-certitude (95 %) de parvenir au résultat.
Autre exemple : le programme checking_arrays, abordé plus haut, a permis de montrer que pour 298 483 points créés au total, si 234 428 sont placés dans le cercle, alors l’approximation de \(\pi\) qui en résulte atteindra une précision de 10 décimales. C’est d’ailleurs le plus petit nombre de points qui permet d’arriver à ce résultat.
On veut savoir combien de fois il faudra répéter l’expérience (= tracer 298 483 points aléatoires) pour avoir 80 % de chances d’arriver au moins une fois à en placer exactement 234 428 dans la cible :
prob(298_483, 234_428, chance=80/100)
'nb = 905'
Il faudra donc parcourir au moins 905 boucles pour arriver à une probabilité de 0,8.
Et si l’on veut connaître les chances de parvenir au but en parcourant seulement 100 boucles :
prob(298_483, 234_428, nb=100)
'chance = 0.163076'
Le calcul des probabilités permet aussi d’estimer le temps de calcul nécessaire pour parvenir à un succès puisqu’il est directement proportionnel au nombre de boucles parcourues.
Mais on ne peut pas se contenter d’obtenir par hasard le résultat espéré. Encore faut-il que l’expérience soit reproductible.
⁂
Chercher la graine (seed )
Le programme recherchant le seed conduisant à un succès peut s’écrire :
import numpy as np
def search_seed():
rng = np.random.default_rng(seed) # le PRNG est réinitialisé à chaque appel
tab = rng.random((n, 2))
return np.sum(tab[:, 0]**2 + tab[:, 1]**2 < 1)
seed = 0 # le seed est initialisé à 0
n = 298_483 # le nombre total de points tracés
# le nombre de points qui doivent être à l’intérieur du cercle
nb_cible = 234_428
while search_seed() != nb_cible:
seed += 1
print('Seed trouvé :', seed)
Seed trouvé : 167
Il a suffi ici d’une poignée de secondes pour parvenir au résultat. On sait désormais qu’il suffira d’initialiser le générateur aléatoire avec la valeur 167 comme seed pour être certain qu’en générant 298 483 points pseudo-aléatoires, 234 428 d’entre eux seront dans le cercle. Cela nous donnera 10 décimales de \(\pi\) exactes :
4 * 234_428 / 298_483
3.141592653517956
Le résultat est bien sûr reproductible :
import numpy as np
rng = np.random.default_rng(167)
tab = rng.random((298_483, 2))
print(np.sum(tab[:, 0]**2 + tab[:, 1]**2 < 1))
234428
Et voilà ! Le générateur de nombres aléatoires n’en est plus tout à fait un puisqu’en le nourrissant avec des graines soigneusement choisies, on peut l’amener au résultat attendu. Toute la difficulté est bien sûr de trouver les bonnes graines.
Plus le nombre de points augmente, plus la probabilité de placer exactement le nombre voulu de points dans la cible diminue, et donc plus le nombre d’itérations nécessaires pour trouver le seed qui mène à ce résultat augmente. Dans le même temps, chaque itération prend plus de temps puisque NumPy doit traiter des tableaux contenant plus d’éléments.
À titre d’illustration, voyons quelles sont les chances d’arriver à trouver les seeds qui correspondent à ces trois résultats remarquables :
Total Cible Valeur approchée de π Nb déc.
------------ ----------- -------------------------------- -------
156513558 122925461 3.141592653589793160283277184204 15
340262731 267241724 3.141592653589793235392564929481 17
26805949036 21053343141 3.141592653589793238462381742775 21
Ces trois résultats ont été obtenus avec la fonction check_arrays() que l’on a vue plus haut. Chacun correspond à un record : le premier est le n le plus bas permettant d’obtenir une précision de 15 décimales, le second est la plus grande précision (17 décimales) qu’il est possible d’atteindre pour n inférieur à 1 milliard. Enfin, le troisième est la première combinaison total-cible qui nous donne une précision atteignant (et même dépassant) 20 décimales exactes.
Il serait intéressant de trouver les seeds initiant les séquences aléatoires qui permettent d’arriver à ces résultats.
Mais est-ce bien réaliste en terme de temps de calcul ?
La fonction prob() que l’on a vue plus haut peut apporter des réponses :
# Variables
p15 = (156_513_558, 122_925_461)
p17 = (340_262_731, 267_241_724)
p21 = (26_805_949_036, 21_053_343_141)
sample = p15, p17, p21
# Quelle est la probabilité d’avoir un succès avant d’avoir
# parcouru 20_000 boucles ?
for params in ((n, nb_cible, 20_000) for n, nb_cible in sample):
print(prob(*params))
chance = 0.788500 chance = 0.651326 chance = 0.111915
On voit donc que pour les deux premiers, il y a de bonnes chances (79 % et 65 %) de parvenir à placer le nombre voulu de points dans la cible avant d’avoir parcouru 20 000 boucles, alors pour le troisième, cette probabilité est beaucoup plus faible (11 %).
# Combien de boucles faut-il pour avoir 2 chances sur 3
# de parvenir au résultat ?
for params in ((n, nb_cible, None, 2/3) for n, nb_cible in sample):
print(prob(*params))
nb = 14144 nb = 20855 nb = 185127
En fonction de ces résultats, on peut décider ou non de se lancer dans l’aventure. Le temps de calcul peut être estimé par exemple en mesurant le temps que met le programme pour parcourir une boucle, et en multipliant par le nombre de boucles qu’il faut pour parvenir à une probabilité raisonnable de succès.
Mais il reste toujours une part d’incertitude. Rien n’empêche a priori le seed convoité de se présenter dès les premières itérations.
Optimisation de l’algorithme
Avant de se lancer dans de longues heures de calcul, il est préférable d’optimiser le programme de recherche.
- Répartir la charge
Plutôt que de tenter de traiter un tableau de 156 513 558 \(\times\) 2 valeurs aléatoires pour modéliser le placement de 156 513 558 points, on pourra parcourir 156 sous-tableaux de 1 million \(\times\) 2 éléments et un dernier de 513 558 \(\times\) 2 éléments. - Traitement asynchrone
On gagnera du temps en exécutant de manière asynchrone l’itération sur les seeds. Les sous-classesProcessPoolExecutoretThreadPoolExecutordu module concurrent.futures permettent, au choix, d’utiliser des threads ou des processus séparés, ces derniers étant répartis sur plusieurs processeurs.
Si l’on veut connaître le nombre de processeurs dont on dispose, on peut utiliser la commande :
import os
os.cpu_count()
4
Mais cette information n’est pas indispensable puisque, par défaut, ProcessPoolExecutor utilisera le maximum de processeurs à sa disposition (avec une limite de 61 sous Windows, d’après la documentation).
L’utilisation du multithreading ou du multiprocessing (l’intérêt d’utiliser l’un plutôt que l’autre peut varier selon la configuration matérielle), permet de raccourcir sensiblement le temps de traitement sur un grand nombre d’itérations. Dans le cas présent, j’ai constaté un gain de temps de l’ordre de 100 % par rapport à une solution non parallélisée, mais cela peut dépendre de la configuration.
Voici le programme de recherche de seed, avec ces deux optimisations que sont la répartition par tableaux et la parallélisation par multiprocessing.
Nous lançons ici une recherche pour trouver le seed associé à p15 = (156 513 558, 122 925 461) qui doit nous donner 15 décimales de \(\pi\) exactes.
import numpy as np
import os
import concurrent.futures as cf
n = 156_513_558 # nombre de points total
objectif = 122_925_461 # nombre idéal de pts à l’intérieur du cercle
inter = range(0, 20_000) # intervalle des valeurs de seeds à tester
# répartition en tableaux
tabsize = int(1e6)
floor, remain = divmod(n, tabsize)
loops = [tabsize for _ in range(floor)]+([remain] if remain else [])
def search(sd):
"""
Calcule le nombre de points placés à l’intérieur du cercle
le PRNG étant initialisé avec un seed donné en argument.
Renvoie cette valeur.
"""
rng = np.random.default_rng(seed=sd)
inside = 0
for ntab in loops:
tab = rng.random((ntab, 2))
inside += np.sum(tab[:, 0]**2 + tab[:, 1]**2 < 1)
return inside
# Traitement principal. Répartition du calcul sur plusieurs processeurs
with cf.ProcessPoolExecutor() as executor: # ou ThreadPoolExecutor
for sd, ins in zip(inter, executor.map(search, inter)):
if ins == objectif:
msg = f'Seed trouvé : {sd}'
break
else:
msg = f'Aucun résultat dans {inter}.'
# Affichage du résultat
print(msg)
# Ou enregistrement dans un ficher (+ extinction ordinateur)
# with open('seed_found.txt', 'a') as f:
# f.write(msg+'\n')
#
# os.system('shutdown')
Seed trouvé : 17823
C’est le résultat qu’on finit par obtenir en laissant tourner ce programme quelques heures (on peut le laisser tourner la nuit, lui demander d’écrire le résultat dans un ficher puis d’éteindre l’ordinateur).
Il est facile de vérifier que ce résultat est correct :
import numpy as np
rng = np.random.default_rng(17_823)
tab = rng.random((156_513_558, 2))
print(np.sum(tab[:, 0]**2 + tab[:, 1]**2 < 1))
122925461
C’est bien le nombre idéal de points que l’on voulait avoir à l’intérieur du cercle, celui qui nous donne une précision de 15 décimales :
4 * 122_925_461 / 156_513_558
3.141592653589793
De la même manière, avec un peu de patience, on finit par trouver pour un objectif de 267 241 724 points placés dans la cible sur un total de 340 262 731, un seed qui nous mène exactement à ce résultat :
Seed trouvé : 16584
from decimal import Decimal
print(4 * Decimal(267_241_724) / Decimal(340_262_731))
3.141592653589793235392564929
On atteint cette fois une précision de 17 décimales.
En revanche, pour n = 26 805 949 036, avec l’objectif de 21 053 343 141 points dans la cible (précision de 20 décimales), la recherche du seed demande un peu plus de temps. Mais avec une bonne puissance de calcul, on devrait y parvenir sans trop de difficulté.
Fort de ces résultats, et grâce à l’animation créée dans la partie II, il est possible d’observer les mouvements de la courbe aux moments où l’approximation atteint des sommets de précision.
