Bonne journée Pi tout le monde! Pour aucune raison, j'essaie de construire un estimateur de Monte Carlo de Pi aussi court que possible. Pouvons-nous en construire un qui peut tenir dans un tweet?

Pour clarifier, ce que j'ai à l'esprit est l'approche typique consistant à dessiner des points aléatoires à partir du carré unitaire et à calculer le rapport qui se situe dans le cercle unitaire. Le nombre d'échantillons peut être codé en dur ou non. Si vous les codez en dur, vous devez utiliser au moins 1000 échantillons. Le résultat peut être renvoyé ou imprimé sous forme de virgule flottante, de virgule fixe ou de nombre rationnel.

Aucune fonction trig ou constantes Pi, doit être une approche Monte Carlo.

Il s'agit du code golf, donc la soumission la plus courte (en octets) l'emporte.

80386 code machine, 40 38 octets

Hexdump du code:

60 33 db 51 0f c7 f0 f7 e0 52 0f c7 f0 f7 e0 58
03 d0 72 03 83 c3 04 e2 eb 53 db 04 24 58 db 04
24 58 de f9 61 c3

Comment obtenir ce code (à partir du langage d'assemblage):

    // ecx = n (number of iterations)
    pushad;
    xor ebx, ebx; // counter
    push ecx; // save n for later
myloop:
    rdrand eax; // make a random number x (range 0...2^32)
    mul eax; // calculate x^2 / 2^32
    push edx;
    rdrand eax; // make another random number y
    mul eax; // calculate y^2 / 2^32
    pop eax;
    add edx, eax; // calculate D = x^2+y^2 / 2^32 (range 0...2^33)
    jc skip; // skip the following if outside the circle
    add ebx, 4; // accumulate the result multiplied by 4
skip:
    loop myloop;
    push ebx; // convert the result
    fild dword ptr [esp]; // to floating-point
    pop eax;
    fild dword ptr [esp]; // convert n to floating-point
    pop eax;
    fdivp st(1), st; // divide

    popad;
    ret;

Il s'agit d'une fonction utilisant la fastcallconvention d'appel MS (le nombre d'itérations est passé dans le registre ecx). Il renvoie le résultat dans le stregistre.

Choses amusantes sur ce code:

• rdrand - seulement 3 octets pour générer un nombre aléatoire!
• Il utilise l'arithmétique entière (non signée) jusqu'à la division finale.
• La comparaison de la distance au carré ( D) avec le rayon au carré (2^32 ) est effectuée automatiquement - le drapeau de report contient le résultat.
• Pour multiplier le nombre par 4, il compte les échantillons par étapes de 4.

Matlab / Octave, 27 octets

Je sais qu'il existe déjà une réponse Matlab / Octave, mais j'ai essayé ma propre approche. J'ai utilisé le fait que l'intégrale 4/(1+x^2)entre 0 et 1 est pi.

mean(4./(1+rand(1,1e5).^2))

R, 40 (ou 28 ou 24 en utilisant d'autres méthodes)

mean(4*replicate(1e5,sum(runif(2)^2)<1))

mean(4*sqrt(1-runif(1e7)^2))

mean(4/(1+runif(1e7)^2))

Python 2, 56

Un autre Python, si numpy est autorisé, mais assez similaire à Matlab / Octave:

import numpy;sum(sum(numpy.random.rand(2,8e5)**2)<1)/2e5

Mathematica, 42 40 39 octets (ou 31/29?)

J'ai trois solutions à 42 octets:

4Count[1~RandomReal~{#,2},p_/;Norm@p<1]/#&
4Tr@Ceiling[1-Norm/@1~RandomReal~{#,2}]/#&
4Tr@Round[1.5-Norm/@1~RandomReal~{#,2}]/#&

Ce sont toutes des fonctions sans nom qui prennent le nombre d'échantillons net renvoient une approximation rationnelle π. D'abord, ils génèrent tous des npoints dans le carré unitaire du quadrant positif. Ensuite, ils déterminent le nombre d'échantillons qui se trouvent dans le cercle unitaire, puis ils divisent par le nombre d'échantillons et multiplient par 4. La seule différence réside dans la façon dont ils déterminent le nombre d'échantillons à l'intérieur du cercle unitaire:

• Le premier utilise Countà la condition que Norm[p] < 1.
• Le second soustrait la norme de chaque point de 1puis arrondit. Cela transforme les nombres à l'intérieur du cercle d'unité en 1et ceux à l'extérieur vers 0. Après je les résume tous Tr.
• Le troisième fait essentiellement la même chose, mais soustrait le de 1.5, donc je peux utiliser à la Roundplace de Ceiling.

Aaaaaand en écrivant ceci, il m'est venu à l' esprit qu'il existe en effet une solution plus courte, si je soustrais 2et utilise ensuite Floor:

4Tr@Floor[2-Norm/@1~RandomReal~{#,2}]/#&

ou enregistrer un autre octet en utilisant les opérateurs de plancher ou de plafond Unicode:

4Tr@⌊2-Norm/@1~RandomReal~{#,2}⌋/#&
4Tr@⌈1-Norm/@1~RandomReal~{#,2}⌉/#&

Notez que les trois solutions basées sur l'arrondi peuvent également être écrites avec Meanau lieu de Tret sans /#, à nouveau pour les mêmes octets.

Si d'autres approches basées sur Monte Carlo sont correctes (en particulier, celle que Peter a choisie), je peux faire 31 octets en estimant l'intégrale de ou 29 en utilisant l'intégrale de , cette fois donnée comme un nombre à virgule flottante:√(1-x2)1/(1+x2)

4Mean@Sqrt[1-1~RandomReal~#^2]&
Mean[4/(1+1~RandomReal~#^2)]&

CJam, 27 23 22 ou 20 octets

4rd__{{1dmr}2*mhi-}*//

2 octets enregistrés grâce à Runner112, 1 octet enregistré grâce à Sp3000

Il prend le nombre d'itérations de STDIN en entrée.

C'est aussi simple que possible. Ce sont les étapes principales impliquées:

• Lisez l'entrée et exécutez les itérations Monte Carlo plusieurs fois
• À chaque itération, obtenez la somme du carré de deux flottants aléatoires de 0 à 1 et voyez s'il est inférieur à 1
• Obtenez le ratio du nombre de fois où nous avons obtenu moins de 1 par total d'itérations et multipliez-le par 4 pour obtenir PI

Expansion du code :

4rd                     "Put 4 on stack, read input and convert it to a double";
   __{            }*    "Take two copies, one of them determines the iteration"
                        "count for this code block";
      {1dmr}2*          "Generate 2 random doubles from 0 to 1 and put them on stack";
              mh        "Take hypot (sqrt(x^2 + y^2)) where x & y are the above two numbers";
                i       "Convert the hypot to 0 if its less than 1, 1 otherwise";
                 -      "Subtract it from the total sum of input (the first copy of input)";
                    //  "This is essentially taking the ratio of iterations where hypot";
                        "is less than 1 by total iterations and then multiplying by 4";

Essayez-le en ligne ici

Si la valeur moyenne de 1/(1+x^2)est également considérée comme Monte Carlo, cela peut être fait en 20 octets:

Urd:K{4Xdmr_*)/+}*K/

Essayez-le ici

Python 2, 77 75 octets

from random import*;r=random;a=0;exec"a+=r()**2+r()**2<1;"*4000;print a/1e3

Utilise 4000 échantillons pour enregistrer des octets 1e3.

Commodore 64 Basic, 45 octets

1F┌I=1TO1E3:C=C-(R/(1)↑2+R/(1)↑2<1):N─:?C/250

Substitutions PETSCII: ─= SHIFT+E, /= SHIFT+N, ┌=SHIFT+O

Génère 1000 points dans le premier quadrant; pour chacun, ajoute la justesse de "x ^ 2 + y ^ 2 <1" à un décompte en cours, puis divise le décompte par 250 pour obtenir pi. (La présence d'un signe moins est due au fait que sur le C64, "true" = -1.)

J, 17 octets

Calcule la valeur moyenne des valeurs d' 40000échantillon de la fonction 4*sqrt(1-sqr(x))dans la plage [0,1].

Maniablement 0 o.xretours sqrt(1-sqr(x)).

   1e4%~+/0 o.?4e4$0
3.14915

> <> (Poisson) , 114 octets

:00[2>d1[   01v
1-:?!vr:@>x|  >r
c]~$~< |+!/$2*^.3
 .41~/?:-1r
|]:*!r$:*+! \
r+)*: *:*8 8/v?:-1
;n*4, $-{:~ /\r10.

Maintenant,> <> n'a pas de générateur de nombres aléatoires intégré. Il a cependant une fonction qui envoie le pointeur dans une direction aléatoire. Le générateur de nombres aléatoires dans mon code:

______d1[   01v
1-:?!vr:@>x|  >r
_]~$~< |+!/$2*^__
 __________
___________ _
_____ ____ _______
_____ ____~ ______

Il génère essentiellement des bits aléatoires qui composent un nombre binaire, puis convertit ce nombre binaire aléatoire en décimal.

Le reste n'est que les points réguliers de l'approche carrée.

Utilisation: lorsque vous exécutez le code, vous devez vous assurer de préremplir la pile (-v dans l'interpréteur python) avec le nombre d'échantillons, par exemple

pi.fish -v 1000

résultats

3.164

Matlab ou Octave 29 octets (merci à flawr!)

mean(sum(rand(2,4e6).^2)<1)*4

(Je ne sais pas trop si <1 est OK. J'ai lu que cela devrait être <= 1. Mais quelle est la probabilité de dessiner exactement 1 ...)

Matlab ou Octave 31 octets

sum(sum(rand(2,4e3).^2)<=1)/1e3

Java, 108 octets

double π(){double π=0,x,i=0;for(;i++<4e5;)π+=(x=Math.random())*x+(x=Math.random())*x<1?1e-5:0;return π;}

Quatre mille itérations, ajoutant 0,001 à chaque fois que le point se trouve à l'intérieur du cercle unitaire. Des trucs assez basiques.

Remarque: Oui, je sais que je peux perdre quatre octets en passant πà un caractère à un octet. J'aime ça comme ça.

Javascript: 62 octets

for(r=Math.random,t=c=8e4;t--;c-=r()**2+r()**2|0);alert(c/2e4)

J'ai utilisé la réponse javascript précédente (maintenant supprimée) et rasé 5 octets.

GolfScript (34 caractères)

0{^3?^rand.*^.*+/+}2000:^*`1/('.'@

Démo en ligne

Cela utilise un point fixe car GS n'a pas vraiment de virgule flottante. Il abuse légèrement de l'utilisation du point fixe, donc si vous voulez changer le nombre d'itérations, assurez-vous que c'est deux fois une puissance de dix.

Crédit xnor pour la méthode particulière Monte Carlo employé.

Python 2, 90 85 81 octets

from random import*;r=random;print sum(4.for i in[0]*9**7if r()**2+r()**2<1)/9**7

renvoie 3.14120037157par exemple. Le nombre d'échantillons est 4782969 (9 ^ 7). Vous pouvez obtenir un meilleur pi avec 9 ^ 9 mais vous devrez être patient.

Rubis, 39 octets

p (1..8e5).count{rand**2+rand**2<1}/2e5

L'un des points forts est que celui-ci est capable d'utiliser la 8e5notation, ce qui le rend extensible jusqu'à ~ 8e9 échantillons dans le même nombre d'octets de programme.

Perl 6 , 33 octets

{4*sum((1>rand²+rand²)xx$_)/$_}

Essayez-le en ligne!

Il s'agit d'une fonction qui prend le nombre d'échantillons comme argument.

Scala, 87 77 66 octets

def s=math.pow(math.random,2);Seq.fill(1000)(s+s).count(_<1)/250d

Pure Bash, 65 octets

for((;i++<$1*4;a+=RANDOM**2+RANDOM**2<32767**2));{ :;}
echo $a/$1

Prend un seul paramètre de ligne de commande qui est multiplié par 4 pour donner le nombre d'échantillons. L'arithmétique de Bash est uniquement un entier, donc un rationnel est produit. Cela peut être canalisé bc -lpour la division finale:

$ ./montepi.sh 10000
31477/10000
$ ./montepi.sh 10000|bc -l
3.13410000000000000000
$ 

Joe , 20 19 octets

Remarque: cette réponse n'est pas concurrente, car la version 0.1.2, qui ajoutait de l'aléatoire, a été publiée après ce défi.

Fonction nommée F:

:%$,(4*/+1>/+*,?2~;

Fonction sans nom:

%$,(4*/+1>/+*,?2~;)

Ces deux éléments prennent le nombre d'échantillons comme argument et renvoient le résultat. Comment travaillent-ils?

%$,(4*/+1>/+*,?2~;)
   (4*/+1>/+*,?2~;) defines a chain, where functions are called right-to-left
               2~;  appends 2 to the argument, giving [x, 2]
              ?     create a table of random values from 0 to 1 with that shape
            *,      take square of every value
          /+         sum rows, giving a list of (x**2+y**2) values
        1>           check if a value is less than 1, per atom
      /+             sum the results
    4*               multiply by four
%$,                  divide the result by the original parameter

L'exemple s'exécute:

   :%$,(4*/+1>/+*,?2~;
   F400000
3.14154
   F400000
3.14302

dc, 59 caractères (les espaces sont ignorés)

[? 2^ ? 2^ + 1>i]su
[lx 1+ sx]si
[lu x lm 1+ d sm ln>z]sz

5k
?sn
lzx
lx ln / 4* p
q

J'ai testé cela sur Plan 9 et OpenBSD, donc j'imagine que cela fonctionnera sur Linux (GNU?) dc.

Explication par ligne:

1. Stocke le code pour [lire et mettre en carré deux flottants; exécuter le registre isi 1 est supérieur à la somme de leurs carrés] dans le registre u.
2. Stocke le code dans [incrémenter le registre xde 1] dans le registre i.
3. Stocke le code dans [exécuter le registre u, incrémenter le registre m, puis exécuter le registre zsi le registre mest supérieur au registre n] dans le registre z.

5. Réglez l'échelle à 5 décimales.

6. Lisez le nombre de points à échantillonner à partir de la première ligne d'entrée.
7. Exécuter le registre z.
8. Diviser le registre x(le nombre de visites) par registren (le nombre de points), multipliez le résultat par 4 et imprimez.
9. Quitter.

Cependant, j'ai triché:

Le programme a besoin d'une réserve de flotteurs aléatoires entre 0 et 1.

/* frand.c */
#include <u.h>
#include <libc.h>

void
main(void)
{
    srand(time(0));

    for(;;)
        print("%f\n", frand());
}

Usage:

#!/bin/rc
# runpi <number of samples>

{ echo $1; frand } | dc pi.dc

Essai:

% runpi 10000
3.14840

Maintenant avec moins de triche (100 octets)

Quelqu'un a souligné que je pouvais inclure un simple prng.
http://en.wikipedia.org/wiki/RANDU

[lrx2^lrx2^+1>i]su[lx1+sx]si[luxlm1+dsmln>z]sz[0kls65539*2 31^%dsslkk2 31^/]sr?sn5dksk1sslzxlxlm/4*p

Non golfé

[
Registers:
u - routine : execute i if sum of squares less than 1
i - routine : increment register x
z - routine : iterator - execute u while n > m++
r - routine : RANDU PRNG
m - variable: number of samples
x - variable: number of samples inside circle
s - variable: seed for r
k - variable: scale for division
n - variable: number of iterations (user input)
]c
[lrx 2^ lrx 2^ + 1>i]su
[lx 1+ sx]si
[lu x lm 1+ d sm ln>z]sz
[0k ls 65539 * 2 31^ % d ss lkk 2 31 ^ /]sr
? sn
5dksk
1 ss
lzx
lx lm / 4*
p

Essai:

$ echo 10000 | dc pigolf.dc
3.13640

Pyth, 19

c*4sm<sm^OQ2 2*QQQQ

Donnez le nombre d'itérations souhaité en entrée.

Manifestation

Comme Pyth n'a pas de fonction "Nombre flottant aléatoire", j'ai dû improviser. Le programme choisit deux entiers positifs aléatoires inférieurs à l'entrée, les carrés, les sommes et comparés à l'entrée au carré. Cela a été effectué un nombre de fois égal à l'entrée, puis le résultat est multiplié par 4 et divisé par l'entrée.

Dans les nouvelles connexes, j'ajouterai prochainement une opération de nombre à virgule flottante aléatoire à Pyth. Ce programme n'utilise cependant pas cette fonctionnalité.

Si nous interprétons "Le résultat peut être renvoyé ou imprimé sous forme de virgule flottante, de virgule fixe ou de nombre rationnel". libéralement, alors l'impression du numérateur et du dénominateur de la fraction résultante devrait être suffisante. Dans ce cas:

Pyth, 18

*4sm<sm^OQ2 2*QQQQ

Il s'agit d'un programme identique, avec l'opération de division en virgule flottante ( c) supprimée.

Julia, 37 octets

4mean(1-floor(sum(rand(4^8,2).^2,2)))

Le nombre d'échantillons est de 65536 (= 4 ^ 8).

Une variante légèrement plus longue: une fonction avec le nombre d'échantillons scomme seul argument:

s->4mean(1-floor(sum(rand(s,2).^2,2)))

C, 130 octets

#include<stdlib.h>f(){double x,y,c=0;for(int i=0;i<8e6;++i)x=rand(),y=rand(),c+=x*x+y*y<1.0*RAND_MAX*RAND_MAX;printf("%f",c/2e6);}

Non golfé:

#include <stdlib.h>
f(){
 double x,y,c=0;
 for(int i=0; i<8e6; ++i) x=rand(), y=rand(), c+=x*x+y*y<1.0*RAND_MAX*RAND_MAX;
 printf("%f",c/2e6);
}

En fait , 14 octets (non concurrents)

`G²G²+1>`nkæ4*

Essayez-le en ligne!

Cette solution n'est pas concurrente car la langue est postérieure au challenge. Le nombre d'échantillons est donné en entrée (plutôt que codé en dur).

Explication:

`G²G²+1>`nkæ4*
`G²G²+1>`n      do the following N times:
 G²G²+            rand()**2 + rand()**2
      1>          is 1 greater?
          kæ    mean of results
            4*  multiply by 4

Raquette 63 octets

Utilisation de la méthode de réponse en langage R par @Matt:

(/(for/sum((i n))(define a(/(random 11)10))(/ 4(+ 1(* a a))))n)

Non golfé:

(define(f n)
   (/
    (for/sum ((i n))
      (define a (/(random 11)10))
      (/ 4(+ 1(* a a))))
    n))

Essai:

(f 10000)

Sortie (fraction):

3 31491308966059784/243801776017028125

Comme décimal:

(exact->inexact(f 10000))

3.13583200307849

Fortran (GFortran) , 84 83 octets

CALL SRAND(0)
DO I=1,4E3
X=RAND()
Y=RAND()
IF(X*X+Y*Y<1)A=A+1E-3
ENDDO
PRINT*,A
END

Essayez-le en ligne!

This code is very bad wrote. It will fail if gfortran decide to initialize variable A with other value then 0 (which occurs 50% of the compilations, approximately) and, if A is initialized as 0, it will always generate the same random sequence for the given seed. Then, the same value for Pi is printed always.

This is a much better program:

Fortran (GFortran), 100 99 bytes

A=0
DO I=1,4E3
CALL RANDOM_NUMBER(X)
CALL RANDOM_NUMBER(Y)
IF(X*X+Y*Y<1)A=A+1E-3
ENDDO
PRINT*,A
END

Try it online!

(One byte saved in each version; thanks Penguino).

Japt, 26 or 18 bytes

o r_+ÂMhMr p +Mr p <1Ã*4/U

Try it online!

Analogous to Optimizer's answer, mainly just trying to learn Japt.
Takes the number of iterations to run for as the implicit input U.

o                           Take the input and turn it into a range [0, U),
                            essentially a cheap way to get a large array.
  r_                        Reduce it with the default initial value of 0.
    +Â                      On each iteration, add one if
      MhMr p +Mr p          the hypotenuse of a random [0,1)x[0,1) right triangle
                   <1       is smaller than one.
                     Ã*4/U  Multiply the whole result by four and divide by input.

If 1/(1+x^2) is allowed (instead of two separate randoms), then we can achieve 18 bytes with the same logic.

o r_Ä/(1+Mr pÃ*4/U

F#, 149 bytes

open System;
let r=new Random()
let q()=
 let b=r.NextDouble()
 b*b
let m(s:float)=(s-Seq.sumBy(fun x->q()+q()|>Math.Sqrt|>Math.Floor)[1.0..s])*4.0/s

Try it online!

As far as I can make out, to do this kind of running total in F# it's shorter to create an array of numbers and use the Seq.sumBy method than use a for..to..do block.

What this code does it that it creates a collection of floating-point numbers from 1 to s, performs the function fun x->... for the number of elements in the collection, and sums the result. There are s elements in the collection, so the random test is done s times. The actual numbers in the collection are ignored (fun x->, but x is not used).

It also means that the application has to first create and fill the array, and then iterate over it. So it's probably twice as slow as a for..to..do loop. And with the array creation memory usage is in the region of O(f**k)!

For the actual test itself, instead of using an if then else statement what it does is it calculates the distance (q()+q()|>Math.Sqrt) and rounds it down with Math.Floor. If the distance is within the circle, it'll be rounded down to 0. If the distance is outside the circle it'll be rounded down to 1. The Seq.sumBy method will then total these results.

Note then that what Seq.sumBy has totalled is not the points inside the circle, but the points outside it. So for the result it takes s (our sample size) and subtracts the total from it.

It also appears that taking a sample size as a parameter is shorter than hard-coding the value. So I am cheating a little bit...

Haskell, 116 114 110 96 bytes

d=8^9
g[a,b]=sum[4|a*a+b*b<d*d]
p n=(sum.take(floor n)$g<$>iterate((\x->mod(9*x+1)d)<$>)[0,6])/n

Because dealing with import System.Random; r=randoms(mkStdGen 2) would take too many precious bytes, I generate an infinite list of random numbers with the linear congruential generator that some say is almost cryptographically strong: x↦x*9+1 mod 8^9, which by the Hull-Dobell Theorem has the full period of 8^9.

g yields 4 if the random number point is inside the circle for pairs of random numbers in [0..8^9-1] because this eliminates a multiplication in the formula used.

Usage:

> p 100000
3.14208

Try it online!

Perl 5, 34 bytes

$_=$a/map$a+=4/(1+(rand)**2),1..$_

The number of samples is taken from stdin. Requires -p.

Works because:

Try it online!