Estimateur de Monte Carlo de Pi

25

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.

keegan
la source
2
les fonctions trigonométriques sont-elles autorisées? Je vous suggère de les interdire explicitement.
Level River St
((0..4e9).map{rand**2+rand**2<1}.to_s.sub(/./,"$1.")
John Dvorak
@JanDvorak Comment est-ce censé fonctionner? Le ne mapvous donne- t-il pas un tableau de trueet false?
Martin Ender
@ MartinBüttner Ah, oups, désolé. .filter{...}.sizedevrait fonctionner, cependant.
John Dvorak
@JanDvorak Effectivement. C'est vraiment bien :)
Martin Ender

Réponses:

17

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.
anatolyg
la source
Le commentaire doit se lire "Calculer x ^ 2% 2 ^ 32"
Cole Johnson
@ColeJohnson Non - le nombre aléatoire est dedans eax; la mulcommande la multiplie par elle-même et y met la partie haute edx; la partie basse en eaxest jetée.
anatolyg
11

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))
flawr
la source
Un algorithme différent est toujours génial! Aussi, plus efficace!
anatolyg
7

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
Mat
la source
6

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)]&
Martin Ender
la source
9
Vous avez trois solutions à la vie, à l'univers et à tout et vous décidez de le ruiner? Hérésie.
seequ
6

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

Optimiseur
la source
2
J'ai également essayé une réponse CJam et j'ai réussi à entrer 2 octets sous votre score. Mais mon code est sorti de manière similaire au vôtre, je me sentirais sale de le poster comme une réponse distincte. Tout était le même sauf le choix de la variable et ces deux optimisations: obtenir un nombre aléatoire de 0 à 1 avec 1dmrau lieu de KmrK/, et vérifier si la somme des carrés est supérieure à 1 avec iau lieu de 1>(je pensais que celui-ci était particulièrement intelligent) .
Runer112
@ Runer112 Merci. l' iastuce est vraiment soignée! Et sacrément le manque de documentation pour1dmr
Optimizer
5

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.

Sp3000
la source
5
Vous pouvez obtenir un peu plus de précision sans frais avec ...*8000;print a/2e3.
Logic Knight
5

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.)

marque
la source
Que fait (1)-il?
echristopherson
@echristopherson, vous l'avez mal lu. /n'est pas le symbole de division, c'est le caractère produit en tapant SHIFT+Nsur un clavier Commodore 64. R/(1)est la forme de raccourci pour RND(1), ie. msgstr "produire un nombre aléatoire entre 0 et 1 en utilisant la graine RNG actuelle".
Mark
Oh, tu as raison! Bons vieux personnages graphiques PETSCII.
echristopherson
5

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
randomra
la source
4

> <> (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
cirpis
la source
4

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
Steffen
la source
1
Très belle idée! Vous pouvez enregistrer deux octets supplémentaires avec mean(sum(rand(2,4e6).^2)<1)*4.
flawr
4

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.

Géobits
la source
pourquoi pas 9999 itérations?
Optimizer
1
@Optimizer Il rend la somme plus courte. Pour 9999 itérations, je devrais plutôt ajouter un nombre plus précis à chaque fois, ce qui me coûte des chiffres.
Geobits
1
Vous pouvez enregistrer un autre octet et améliorer la précision en utilisant "4e5" et "1e-5" pour les nombres.
Vilmantas Baranauskas
@VilmantasBaranauskas Merci! J'oublie toujours ça :) C'est tentant d'utiliser à la place 4e9 et 1e-9, mais cela prend un certain temps ...
Geobits
Protip: lorsque vous jouez au golf, vous devez en fait réduire les octets, et non les augmenter artificiellement
Destructible Lemon
3

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.

Guzman Tierno
la source
Vous pouvez créer un lien vers la réponse de cfern pour attribuer correctement le crédit.
Jonathan Frech
Votre réponse semble être un extrait dont les E / S sont interdites . Veuillez corriger ou supprimer votre message.
Jonathan Frech
Désolé, je suis nouveau, je ne savais pas comment mettre le lien vers la solution précédente qui semble maintenant avoir été supprimée. Concernant l'extrait de code: je suis tout à fait d'accord, mais c'était le code de la solution javascript précédente qui, je pense aussi, n'était pas valide pour cette raison. J'ai modifié le mien pour en faire un programme.
Guzman Tierno
Oui; la réponse précédente a été supprimée car elle n'était pas valide - j'ai vu votre réponse avant de recommander la suppression, d'où le commentaire. +1 pour avoir soumis une réponse valide; bienvenue chez PPCG!
Jonathan Frech
2

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é.

Peter Taylor
la source
2

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.

Logic Knight
la source
Vous pouvez enregistrer 3 en remplaçant range(9**7)par [0]*9**7ou quelque chose, puisque vous n'utilisez pas i. Et la liste n'est pas trop longue pour rencontrer des problèmes de mémoire.
Sp3000
Merci. Je voulais m'en débarrasser range()mais j'avais complètement oublié cette astuce.
Logic Knight
J'ai l'impression [0]9**7que la syntaxe n'est pas valide.
seequ
Vous avez raison. J'ai ré-attaché l'astérisque perdu (il était sous mon bureau).
Logic Knight
2

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.

GreyCat
la source
2

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.

Sean
la source
1

Scala, 87 77 66 octets

def s=math.pow(math.random,2);Seq.fill(1000)(s+s).count(_<1)/250d
keegan
la source
Si vous remplacez 1000par 8000et 250davec 2e4vous enregistrez un octet et augmentez le nombre d'échantillons d'un facteur 8.
Dave Swartz
1

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
$ 
Traumatisme numérique
la source
1

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
seequ
la source
1

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.
  4. Réglez l'échelle à 5 décimales.

  5. Lisez le nombre de points à échantillonner à partir de la première ligne d'entrée.
  6. Exécuter le registre z.
  7. Diviser le registre x(le nombre de visites) par registren (le nombre de points), multipliez le résultat par 4 et imprimez.
  8. 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
progrider42
la source
1

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.

isaacg
la source
1

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)))
pawel.boczarski
la source
1

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);
}
Karl Napf
la source
bien sûr, vous devriez toujours poster la version sans les espaces (gardez la version actuelle avec l'en-tête "ungolfed / with white" ou quelque chose comme ça)
Destructible Lemon
@DestructibleWatermelon done!
Karl Napf
La solution ne fonctionne pas dans GCC sans une nouvelle ligne avant f(). Quel compilateur avez-vous utilisé? Voir tio.run/##Pc49C4JAHIDx3U9xGMG9ZdYgwWkgtNbQ1BZ6L/UHO8M07hA/…
eush77
101 octets
plafondcat
1

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
Mego
la source
2
Pourquoi le downvote?
Destructible Lemon
1

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
rnso
la source
1

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).

rafa11111
la source
1
In each version you can save a byte by changing 'DO I=1,1E3' to 'DO I=1,4E3', changing 'A=A+1' to 'A=A+1E-3', and changing 'PRINT*,A/250' to 'PRINT*,A'
Penguino
Yes, you are sure! Thank for the suggestion!
rafa11111
1

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
Nit
la source
1
You can save a few bytes by letting Mh calculate the hypotenuse rather than doing it yourself ;-) Also, you can use x to take the sum of an array, rather than reducing by addition: o x@MhMr Mr)<1Ã*4/U
ETHproductions
@ETHproductions Neat, I didn't know you can use Mh like that, thanks! Your two-random answer is nearly as short as my answer with only one random, that's pretty cool. I'll keep x in mind, I tend to use reduction a lot when trying to golf stuff, so this will come in very handy.
Nit
1

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...

Ciaran_McCarthy
la source
1

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!

Angs
la source
1

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!

primo
la source