¡Feliz Día de Pi a todos! Sin ninguna razón, estoy tratando de construir un estimador de Pi de Monte Carlo que sea lo más corto posible. ¿Podemos construir uno que pueda caber en un tweet?

Para aclarar, lo que tengo en mente es el enfoque típico de dibujar puntos aleatorios del cuadrado de la unidad y calcular la relación que cae dentro del círculo de la unidad. El número de muestras puede estar codificado o no. Si los codifica, debe usar al menos 1000 muestras. El resultado puede ser devuelto o impreso como punto flotante, punto fijo o número racional.

Sin funciones trigonométricas o constantes Pi, debe ser un enfoque de Monte Carlo.

Este es el código de golf, por lo que gana el envío más corto (en bytes).

80386 código máquina, 40 38 bytes

Hexdump del código:

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

Cómo obtener este código (del lenguaje ensamblador):

    // 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;

Esta es una función que utiliza la fastcallconvención de llamada de MS (el número de iteraciones se pasa en el registro ecx). Devuelve el resultado en el stregistro.

Cosas divertidas sobre este código:

• rdrand - ¡solo 3 bytes para generar un número aleatorio!
• Utiliza aritmética de enteros (sin signo) hasta la división final.
• La comparación de la distancia al cuadrado ( D) con el radio al cuadrado ( 2^32) se realiza automáticamente: la bandera de acarreo contiene el resultado.
• Para multiplicar el recuento por 4, cuenta las muestras en pasos de 4.

Matlab / Octave, 27 bytes

Sé que ya hay una respuesta Matlab / Octave, pero probé mi propio enfoque. Usé el hecho de que la integral de 4/(1+x^2)entre 0 y 1 es pi.

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

R, 40 (o 28 o 24 usando otros métodos)

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

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

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

Pitón 2, 56

Otro Python, si está permitido numpy, pero bastante similar a Matlab / Octave:

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

Mathematica, 42 40 39 bytes (¿o 31/29?)

Tengo tres soluciones, todas a 42 bytes:

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

Todas son funciones sin nombre que toman el número de muestras nyd devuelven una aproximación racional π. Primero todos generan npuntos en el cuadrado de la unidad en el cuadrante positivo. Luego determinan el número de esas muestras que se encuentran dentro del círculo unitario, y luego se dividen por el número de muestras y se multiplican por 4. La única diferencia está en cómo determinan el número de muestras dentro del círculo unitario:

• El primero usa Countcon la condición de que Norm[p] < 1.
• El segundo resta la norma de cada punto 1y luego lo redondea. Esto convierte a los números dentro del círculo unidad a 1y a los que están fuera 0. Después solo los sumo a todos Tr.
• El tercero hace esencialmente lo mismo, pero resta el de 1.5, por lo que puedo usarlo en Roundlugar de Ceiling.

Aaaaa, y mientras escribía esto, se me ocurrió que existe una solución más corta, si solo resta 2y luego uso Floor:

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

o guardar otro byte utilizando los operadores de piso o techo Unicode:

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

Tenga en cuenta que las tres soluciones basadas en redondeo también se pueden escribir con, en Meanlugar de Try sin /#, nuevamente, para los mismos bytes.

Si otros enfoques basados ​​en Monte Carlo están bien (específicamente, el que Peter ha elegido), puedo hacer 31 bytes estimando la integral de o 29 usando la integral de , esta vez dada como un número de coma flotante:√(1-x2)1/(1+x2)

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

CJam, 27 23 22 o 20 bytes

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

2 bytes guardados gracias a Runner112, 1 byte guardado gracias a Sp3000

Toma el recuento de iteraciones de STDIN como entrada.

Esto es tan sencillo como parece. Estos son los principales pasos involucrados:

• Lea la entrada y ejecute las iteraciones de Monte Carlo que muchas veces
• En cada iteración, obtenga la suma del cuadrado de dos flotantes aleatorios de 0 a 1 y vea si es menor que 1
• Obtenga la proporción de cuántas veces obtuvimos menos de 1 por iteraciones totales y multiplíquelo por 4 para obtener PI

Expansión de código :

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";

Pruébalo en línea aquí

Si el valor promedio de 1/(1+x^2)también se considera como Monte Carlo, entonces esto se puede hacer en 20 bytes:

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

Pruébalo aquí

Python 2, 77 75 bytes

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

Utiliza 4000 muestras para guardar bytes con 1e3.

Commodore 64 Basic, 45 bytes

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

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

Genera 1000 puntos en el primer cuadrante; para cada uno, agrega la veracidad de "x ^ 2 + y ^ 2 <1" a una cuenta corriente, luego divide la cuenta por 250 para obtener pi. (La presencia de un signo menos se debe a que en el C64, "verdadero" = -1.)

J, 17 bytes

Calcula el valor medio de los valores de 40000muestra de la función 4*sqrt(1-sqr(x))en el rango [0,1].

Prácticamente 0 o.xregresa sqrt(1-sqr(x)).

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

> <> (Pescado) , 114 bytes

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

Ahora,> <> no tiene un generador de números aleatorios incorporado. Sin embargo, tiene una función que envía el puntero en una dirección aleatoria. El generador de números aleatorios en mi código:

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

Básicamente genera bits aleatorios que forman un número binario y luego convierte ese número binario aleatorio a decimal.

El resto son solo los puntos regulares en el enfoque cuadrado.

Uso: cuando ejecuta el código, debe asegurarse de rellenar previamente la pila (-v en el intérprete de Python) con el número de muestras, por ejemplo

pi.fish -v 1000

devoluciones

3.164

Matlab u Octave 29 bytes (¡gracias a flawr!)

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

(No estoy muy seguro si <1 está bien. Leí que debería ser <= 1. Pero, ¿qué tan grande es la probabilidad de sacar exactamente 1 ...?)

Matlab u Octave 31 bytes

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

Java, 108 bytes

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

Cuatro mil iteraciones, agregando 0.001 cada vez que el punto está dentro del círculo unitario. Cosas bastante básicas.

Nota: Sí, sé que puedo eliminar cuatro bytes cambiando πa un carácter de un solo byte. Me gusta así.

Javascript: 62 bytes

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

Utilicé la respuesta de JavaScript anterior (ahora eliminada) y afeité 5 bytes.

GolfScript (34 caracteres)

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

Demostración en línea

Esto usa punto fijo porque GS realmente no tiene punto flotante. Abusa un poco del uso del punto fijo, por lo que si desea cambiar el recuento de iteraciones asegúrese de que sea el doble de una potencia de diez.

Crédito a xnor por el método particular de Monte Carlo empleado.

Python 2, 90 85 81 bytes

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

vuelve 3.14120037157por ejemplo. El recuento de la muestra es 4782969 (9 ^ 7). Puede obtener una mejor pi con 9 ^ 9 pero tendrá que ser paciente.

Ruby, 39 bytes

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

Uno de los aspectos más destacados es que este puede usar 8e5notación, lo que lo hace extensible hasta ~ 8e9 muestras en el mismo conteo de bytes del programa.

Perl 6 , 33 bytes

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

Pruébalo en línea!

Esta es una función que toma el número de muestras como argumento.

Scala, 87 77 66 bytes

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

Pure Bash, 65 bytes

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

Toma un único parámetro de línea de comandos que se multiplica por 4 para dar el número de muestras. La aritmética de Bash es solo de enteros, por lo que se genera un racional. Esto se puede canalizar bc -lpara la división final:

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

Joe , 20 19 bytes

Nota: esta respuesta no es competitiva, porque la versión 0.1.2, que agregó aleatoriedad, se lanzó después de este desafío.

Función denominada F:

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

Función sin nombre:

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

Ambos toman el recuento de muestra como argumento y devuelven el resultado. ¿Cómo trabajan?

%$,(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

Ejecuciones de ejemplo:

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

dc, 59 caracteres (se ignora el espacio en blanco)

[? 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

Probé esto en Plan 9 y OpenBSD, así que imagino que funcionará en Linux (GNU?) dc.

Explicación por línea:

1. Almacena el código para [leer y cuadrar dos carrozas; ejecutar el registro isi 1 es mayor que la suma de sus cuadrados] en el registro u.
2. Almacena el código para [incrementar el registro xen 1] en el registroi .
3. Almacena código para [ejecutar registro u, incrementar registro my luego ejecutar registro zsi el registro mes mayor que el registro n] en el registro z.

5. Establezca la escala a 5 puntos decimales.

6. Lea el número de puntos para muestrear desde la primera línea de entrada.
7. Ejecutar registro z.
8. Divida el registro x(el número de visitas) por el registro n(el número de puntos), multiplique el resultado por 4 e imprima.
9. Dejar.

Sin embargo, hice trampa:

El programa necesita un suministro de flotadores aleatorios entre 0 y 1.

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

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

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

Uso:

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

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

Prueba de funcionamiento:

% runpi 10000
3.14840

Ahora con menos trampa (100 bytes)

Alguien señaló que podría incluir 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

Sin 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

Prueba de funcionamiento:

$ echo 10000 | dc pigolf.dc
3.13640

Pyth, 19

c*4sm<sm^OQ2 2*QQQQ

Dé el número deseado de iteraciones como entrada.

Demostración

Como Pyth no tiene una función de "número flotante aleatorio", tuve que improvisar. El programa elige dos enteros positivos aleatorios menores que la entrada, cuadrados, sumas y en comparación con la entrada al cuadrado. Esto se realizó varias veces igual a la entrada, luego el resultado se multiplica por 4 y se divide por la entrada.

En noticias relacionadas, agregaré una operación de números aleatorios de coma flotante a Pyth en breve. Sin embargo, este programa no usa esa función.

Si interpretamos "El resultado puede ser devuelto o impreso como punto flotante, punto fijo o número racional". generosamente, luego imprimir el numerador y el denominador de la fracción resultante debería ser suficiente. En ese caso:

Pyth, 18

*4sm<sm^OQ2 2*QQQQ

Este es un programa idéntico, con la operación de división de coma flotante ( c) eliminada.

Julia, 37 bytes

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

El número de muestras es 65536 (= 4 ^ 8).

Una variante ligeramente más larga: una función con número de muestras scomo único argumento:

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

C, 130 bytes

#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);}

Sin 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 realidad , 14 bytes (no competitivos)

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

Pruébalo en línea!

Esta solución no es competitiva porque el lenguaje es posterior al desafío. La cantidad de muestras se proporciona como entrada (en lugar de codificada).

Explicación:

`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

Raqueta 63 bytes

Usando el método de respuesta en lenguaje R por @Matt:

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

Sin golf:

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

Pruebas:

(f 10000)

Salida (fracción):

3 31491308966059784/243801776017028125

Como decimal:

(exact->inexact(f 10000))

3.13583200307849

Fortran (GFortran) , 84 83 bytes

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

Pruébalo en línea!

Este código está muy mal escrito. Fallará si gfortran decide inicializar la variable Acon otro valor que no sea 0 (que ocurre aproximadamente en el 50% de las compilaciones) y, si Ase inicializa como 0, siempre generará la misma secuencia aleatoria para la semilla dada. Entonces, siempre se imprime el mismo valor para Pi.

Este es un programa mucho mejor:

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

Pruébalo en línea!

(Un byte guardado en cada versión; gracias Penguino).

Japt , 26 o 18 bytes

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

Pruébalo en línea!

Análogo a la respuesta de Optimizer , principalmente tratando de aprender Japt.
Toma el número de iteraciones para ejecutar como entrada implícita 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.

Si 1/(1+x^2)está permitido (en lugar de dos randoms separados), entonces podemos lograr 18 bytes con la misma lógica.

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

Pruébalo en línea!

Por lo que puedo entender, para hacer este tipo de total acumulado en F # es más corto crear una matriz de números y usar el Seq.sumBy método que usar un for..to..dobloque.

Lo que hace este código es que crea una colección de números de punto flotante del 1 al s, realiza la función fun x->...para el número de elementos en la colección y suma el resultado. Hay selementos en la colección, por lo que la prueba aleatoria se realiza sveces. Los números reales en la colección se ignoran ( fun x->, perox no se usan).

También significa que la aplicación primero debe crear y completar la matriz, y luego iterar sobre ella. Por lo tanto, es probablemente el doble de lento que un for..to..dobucle. ¡Y con la creación de la matriz, el uso de la memoria está en la región de O (f ** k)!

Para la prueba en sí misma, en lugar de usar una if then elsedeclaración, lo que hace es calcular la distancia ( q()+q()|>Math.Sqrt) y redondearla hacia abajo conMath.Floor . Si la distancia está dentro del círculo, se redondeará hacia abajo a 0. Si la distancia está fuera del círculo, se redondeará hacia abajo a 1. El Seq.sumBymétodo totalizará estos resultados.

Tenga en cuenta entonces que lo que Seq.sumByha totalizado no son los puntos dentro del círculo, sino los puntos fuera de él. Entonces, para el resultado que se necesitas (nuestro tamaño de muestra) y resta el total de él.

También parece que tomar un tamaño de muestra como parámetro es más corto que codificar el valor. Así que estoy haciendo trampa un poco ...

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

Debido a que tratar import System.Random; r=randoms(mkStdGen 2)tomaría demasiados bytes preciosos, genero una lista infinita de números aleatorios con el generador lineal congruencial que algunos dicen que es casi criptográficamente fuerte: x↦x*9+1 mod 8^9que según el Teorema de Hull-Dobell tiene el período completo 8^9.

gproduce 4si el punto de número aleatorio está dentro del círculo para pares de números aleatorios en[0..8^9-1] porque esto elimina una multiplicación en la fórmula utilizada.

Uso:

> p 100000
3.14208

Pruébalo en línea!

Perl 5, 34 bytes

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

El número de muestras se toma de stdin. Requiere-p .

Funciona porque:

Pruébalo en línea!