Estimador de Monte Carlo de Pi

25

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

keegan
fuente
2
¿están permitidas las funciones trigonométricas? Le sugiero que los prohíba explícitamente.
Level River St
((0..4e9).map{rand**2+rand**2<1}.to_s.sub(/./,"$1.")
John Dvorak
@ JanDvorak ¿Cómo se supone que funciona? ¿No mapte da una serie de truey false?
Martin Ender
@ MartinBüttner Ah, vaya, lo siento. .filter{...}.sizeSin embargo, debería funcionar.
John Dvorak
@ JanDvorak De hecho. Eso es realmente genial :)
Martin Ender

Respuestas:

17

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.
anatolyg
fuente
El comentario debería leer "Calcular x ^ 2% 2 ^ 32"
Cole Johnson
@ColeJohnson No: el número aleatorio está en eax; el mulcomando lo multiplica por sí mismo y pone la parte alta en edx; la parte baja eaxse descarta.
anatolyg
11

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))
falla
fuente
¡Un algoritmo diferente siempre es genial! Además, más eficiente!
anatolyg
7

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
Mate
fuente
6

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)]&
Martin Ender
fuente
99
¿Tienes tres soluciones para la vida, el universo y todo y decides arruinarlo? Herejía.
seequ
6

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í

Optimizador
fuente
2
También probé una respuesta de CJam y logré ingresar 2 bytes debajo de su puntaje. Pero mi código salió tan similar al tuyo, me sentiría sucio publicarlo como una respuesta separada. Todo era igual, excepto la opción variable y estas dos optimizaciones: obtener un número aleatorio de 0 a 1 con en 1dmrlugar de KmrK/, y verificar si la suma de los cuadrados es mayor que 1 con en ilugar de 1>(pensé que esta era particularmente inteligente) .
Runer112
@ Runer112 Gracias. ¡el itruco es realmente genial! Y maldita sea la falta de documentación para1dmr
Optimizer
5

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.

Sp3000
fuente
55
Puede obtener un poco más de precisión sin costo con ...*8000;print a/2e3.
Logic Knight
5

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

marca
fuente
¿Qué (1)hacer?
echristopherson
@echristopherson, lo estás leyendo mal. /no es el símbolo de división, es el personaje producido al escribir SHIFT+Nen un teclado Commodore 64. R/(1)es la forma de acceso directo para RND(1), es decir. "producir un número aleatorio entre 0 y 1 utilizando la semilla RNG actual".
Mark
Oh tienes razon! Buenos viejos personajes gráficos PETSCII.
echristopherson
5

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
randomra
fuente
4

> <> (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
cirpis
fuente
4

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
Steffen
fuente
1
Muy buena idea! Puede guardar dos bytes adicionales con mean(sum(rand(2,4e6).^2)<1)*4.
flawr
4

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

Geobits
fuente
¿Por qué no 9999 iteraciones?
Optimizador
1
@Optimizer Hace que la suma sea más corta. Para 9999 iteraciones, tendría que agregar un número más preciso cada vez, lo que me cuesta dígitos.
Geobits
1
Puede guardar otro byte y mejorar la precisión utilizando "4e5" y "1e-5" para los números.
Vilmantas Baranauskas
@VilmantasBaranauskas Gracias! Siempre me olvido de eso :) Es tentador usar 4e9 y 1e-9 en su lugar, pero eso lleva bastante tiempo ...
Geobits
Protip: cuando juegues al golf, deberías reducir los bytes, no aumentarlos artificialmente
Destructible Lemon
3

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.

Guzman Tierno
fuente
Puede vincular a la respuesta de cfern para dar crédito adecuadamente.
Jonathan Frech
Su respuesta parece ser un fragmento que se anuló E / S . Por favor, arregla o elimina tu publicación.
Jonathan Frech
Lo siento, soy nuevo, no sabía cómo poner el enlace a la solución anterior que ahora parece haber sido eliminado. Con respecto al fragmento: estoy completamente de acuerdo, pero ese fue el código de la solución de JavaScript anterior, que también creo que no era válido por ese motivo. Modifiqué el mío para que sea un programa.
Guzman Tierno
Sí; la respuesta anterior se eliminó porque no era válida: vi su respuesta antes de recomendar la eliminación, de ahí el comentario. +1 por enviar una respuesta válida; bienvenido a PPCG!
Jonathan Frech
2

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.

Peter Taylor
fuente
2

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.

Caballero Lógico
fuente
Puede guardar 3 reemplazando range(9**7)con [0]*9**7o algo así, ya que no lo usa i. Y la lista no es demasiado larga para encontrarse con problemas de memoria.
Sp3000
Gracias. Quería deshacerme de él, range()pero había olvidado por completo ese truco.
Logic Knight
Tengo la sensación de [0]9**7que no es una sintaxis válida.
seequ
Tienes razón. He vuelto a conectar el asterisco perdido (estaba debajo de mi escritorio).
Logic Knight
2

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.

Gato gris
fuente
1

Scala, 87 77 66 bytes

def s=math.pow(math.random,2);Seq.fill(1000)(s+s).count(_<1)/250d
keegan
fuente
Si reemplaza 1000con 8000y 250dcon 2e4ambos, guarde un byte y aumente el número de muestras en un factor de 8.
Dave Swartz
1

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
$ 
Trauma digital
fuente
1

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
seequ
fuente
1

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.
  4. Establezca la escala a 5 puntos decimales.

  5. Lea el número de puntos para muestrear desde la primera línea de entrada.
  6. Ejecutar registro z.
  7. 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.
  8. 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
progrider42
fuente
1

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.

isaacg
fuente
1

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

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);
}
Karl Napf
fuente
por supuesto, probablemente debería publicar la versión sin el espacio en blanco (mantenga la versión actual con el encabezado "sin golf / con espacio en blanco" o algo así)
Destructible Lemon
@DestructibleWatermelon hecho!
Karl Napf
La solución no funciona en GCC sin una nueva línea antes f(). ¿Qué compilador usaste? Ver tio.run/##Pc49C4JAHIDx3U9xGMG9ZdYgwWkgtNbQ1BZ6L/UHO8M07hA/…
eush77
101 bytes
ceilingcat
1

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
Mego
fuente
2
¿Por qué el voto negativo?
Destructible Lemon
1

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
rnso
fuente
1

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

rafa11111
fuente
1
En cada versión puede guardar un byte cambiando 'DO I = 1,1E3' a 'DO I = 1,4E3', cambiando 'A = A + 1' a 'A = A + 1E-3' y cambiando ' PRINT *, A / 250 'a' PRINT *, A '
Penguino
Si, estas seguro! Gracias por la sugerencia!
rafa11111
1

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
Liendre
fuente
1
Puede guardar algunos bytes dejando Mhcalcular la hipotenusa en lugar de hacerlo usted mismo ;-) Además, puede usar xpara tomar la suma de una matriz, en lugar de reducirla mediante la adición:o x@MhMr Mr)<1Ã*4/U
ETHproductions
@ETHproductions Genial, no sabía que puedes usar Mhasí, ¡gracias! Su respuesta aleatoria es casi tan corta como mi respuesta con solo una aleatoria, eso es genial. Recordaré xque tiendo a usar mucho la reducción cuando intento jugar golf, por lo que esto será muy útil.
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

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

Ciaran_McCarthy
fuente
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

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!

Angs
fuente
1

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!

primo
fuente