Estimador de Monte Carlo de Pi

25

Feliz dia Pi todos! Por nenhuma razão, estou tentando construir um estimador de Pi de Monte Carlo o mais curto possível. Podemos construir um que possa caber em um tweet?

Para esclarecer, o que tenho em mente é a abordagem típica de desenhar pontos aleatórios a partir do quadrado unitário e calcular a proporção que se enquadra no círculo unitário. O número de amostras pode ser codificado ou não. Se você codificá-los, use pelo menos 1000 amostras. O resultado pode ser retornado ou impresso como um ponto flutuante, ponto fixo ou número racional.

Nenhuma função trigonométrica ou constante Pi deve ser uma abordagem de Monte Carlo.

Isso é código de golfe, então a submissão mais curta (em bytes) vence.

keegan
fonte
2
são permitidas funções trigonométricas? Eu sugiro que você os bane explicitamente.
Level River St
((0..4e9).map{rand**2+rand**2<1}.to_s.sub(/./,"$1.")
John Dvorak
@JanDvorak Como isso deve funcionar? O não maplhe dá uma matriz de truee false?
Martin Ender
@ MartinBüttner Ah, opa, desculpe. .filter{...}.sizedeve funcionar, no entanto.
John Dvorak
@JanDvorak Indeed. Isso é realmente puro :)
Martin Ender

Respostas:

17

Código da máquina 80386, 40 38 bytes

Hexdump do 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

Como obter este código (da linguagem assembly):

    // 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 é uma função que usa a fastcallconvenção de chamada da Microsoft (o número de iterações é passado no registro ecx). Retorna o resultado no stregistro.

Coisas divertidas sobre este código:

  • rdrand - apenas 3 bytes para gerar um número aleatório!
  • Ele usa aritmética de número inteiro (não assinado) até a divisão final.
  • A comparação da distância ao quadrado ( D) com o raio ao quadrado ( 2^32) é realizada automaticamente - a bandeira de transporte contém o resultado.
  • Para multiplicar a contagem por 4, conta as amostras nas etapas de 4.
anatolyg
fonte
O comentário deve ser "Calcular x ^ 2% 2 ^ 32"
Cole Johnson
@ColeJohnson Não - o número aleatório está dentro eax; o mulcomando multiplica por si mesmo e coloca a parte alta edx; a parte baixa eaxé descartada.
Anatolyg 22/03
11

Matlab / Octave, 27 bytes

Eu sei que já existe uma resposta do Matlab / Octave, mas tentei minha própria abordagem. Eu usei o fato de que a integral 4/(1+x^2)entre 0 e 1 é pi.

mean(4./(1+rand(1,1e5).^2))
flawr
fonte
Um algoritmo diferente é sempre ótimo! Além disso, mais eficiente!
anatolyg
7

R, 40 (ou 28 ou 24 usando outros métodos)

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

Outro Python, se numpy for permitido, mas bem parecido com o Matlab / Octave:

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

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

Eu tenho três soluções, todas em 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 são funções sem nome que recebem o número de amostras e nretornam uma aproximação racional π. Primeiro, todos eles geram npontos no quadrado da unidade no quadrante positivo. Em seguida, eles determinam o número daquelas amostras que estão dentro do círculo unitário e depois dividem pelo número de amostras e multiplicam por 4. A única diferença está em como eles determinam o número de lamelas dentro do círculo unitário:

  • O primeiro usa Countcom a condição de que Norm[p] < 1.
  • O segundo subtrai a norma de cada ponto 1e depois arredonda para cima. Isso transforma os números dentro do círculo da unidade 1e os de fora para 0. Depois, apenas resumo todos eles Tr.
  • O terceiro faz essencialmente o mesmo, mas subtrai o de 1.5, para que eu possa usar em Roundvez de Ceiling.

Aaaaaand, enquanto escrevia isso, ocorreu-me que há realmente uma solução mais curta, se eu apenas subtrair 2e depois usar Floor:

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

ou salvar outro byte usando os operadores de piso ou teto Unicode:

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

Observe que as três soluções baseadas em arredondamento também podem ser escritas com, em Meanvez de Tre sem o /#, novamente para os mesmos bytes.


Se outras abordagens baseadas em Monte Carlo forem boas (especificamente a escolhida por Peter), eu posso fazer 31 bytes estimando a integral de ou 29 usando a integral de , desta vez dado como um número de ponto flutuante:√(1-x2)1/(1+x2)

4Mean@Sqrt[1-1~RandomReal~#^2]&
Mean[4/(1+1~RandomReal~#^2)]&
Martin Ender
fonte
9
Você tem três soluções para a vida, o universo e tudo e decide arruiná-lo? Heresia.
seequ
6

CJam, 27 23 22 ou 20 bytes

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

2 bytes economizados graças ao Runner112, 1 byte economizado graças ao Sp3000

Leva a contagem de iterações de STDIN como uma entrada.

Isso é o mais direto possível. Estes são os principais passos envolvidos:

  • Leia a entrada e execute as iterações de Monte Carlo que muitas vezes
  • Em cada iteração, obtenha a soma do quadrado de dois flutuadores aleatórios de 0 a 1 e veja se é menor que 1
  • Obtenha a proporção de quantas vezes obtivemos menos de 1 pelo total de iterações e multiplique por 4 para obter PI

Expansão do 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";

Experimente online aqui


Se o valor médio de 1/(1+x^2)também for considerado como Monte Carlo, isso poderá ser feito em 20 bytes:

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

Experimente aqui

Optimizer
fonte
2
Eu tentei uma resposta CJam também e consegui entrar em 2 bytes abaixo da sua pontuação. Mas meu código saiu tão parecido com o seu que eu me sentiria sujo ao publicá-lo como uma resposta separada. Tudo era o mesmo, exceto para a escolha da variável e essas duas otimizações: obter um número aleatório de 0 a 1 com em 1dmrvez de KmrK/e verificar se a soma dos quadrados é maior que 1 com em ivez de 1>(eu pensei que este era particularmente inteligente) .
Runer112
@ Runer112 Obrigado. o itruque é realmente legal! E maldita falta de documentação 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

Usa 4000 amostras para salvar bytes 1e3.

Sp3000
fonte
5
Você pode obter um pouco mais de precisão sem nenhum custo ...*8000;print a/2e3.
Logic Knight
5

Commodore 64 Básico, 45 bytes

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

Substituições PETSCII: = SHIFT+E, /= SHIFT+N, =SHIFT+O

Gera 1000 pontos no primeiro quadrante; para cada um, adiciona a veracidade de "x ^ 2 + y ^ 2 <1" a uma contagem contínua e depois divide a contagem por 250 para obter pi. (A presença de um sinal de menos é porque no C64, "true" = -1.)

Marca
fonte
O que (1)faz?
echristopherson
@echristopherson, você está lendo errado. /não é o símbolo de divisão, é o personagem produzido digitando SHIFT+Nem um teclado Commodore 64. R/(1)é o formulário de atalho para RND(1), ie. "produz um número aleatório entre 0 e 1 usando a semente RNG atual".
Mark
Oh, você está certo! Bons e velhos personagens gráficos PETSCII.
echristopherson
5

J, 17 bytes

Calcula o valor médio dos 40000valores de amostra da função 4*sqrt(1-sqr(x))no intervalo[0,1] .

Com folga 0 o.xretornos sqrt(1-sqr(x)).

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

> <> (Peixe) , 114 bytes

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

Agora,> <> não possui um gerador de números aleatórios embutido. No entanto, possui uma função que envia o ponteiro em uma direção aleatória. O gerador de números aleatórios no meu código:

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

Basicamente, gera bits aleatórios que compõem um número binário e, em seguida, converte esse número binário aleatório em decimal.

O resto são apenas os pontos regulares na abordagem quadrada.

Uso: ao executar o código, você deve preencher previamente a pilha (-v no interpretador python) com o número de amostras, por exemplo

pi.fish -v 1000

retorna

3.164
cirpis
fonte
4

Matlab ou Octave 29 bytes (graças a flawr!)

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

(Não tenho certeza se <1 está OK. Li que deveria ser <= 1. Mas qual é a probabilidade de desenhar exatamente 1 ...)

Matlab ou Octave 31 bytes

sum(sum(rand(2,4e3).^2)<=1)/1e3
Steffen
fonte
1
Idéia muito boa! Você pode salvar dois bytes adicionais com 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 π;}

Quatro mil iterações, adicionando 0,001 cada vez que o ponto está dentro do círculo unitário. Coisas bem básicas.

Nota: Sim, sei que posso eliminar quatro bytes mudando πpara um caractere de byte único. Eu gosto assim.

Geobits
fonte
por que não 9999 iterações?
Optimizer
1
@Optimizer Reduz a soma. Para 9999 iterações, eu teria que adicionar um número mais preciso de cada vez, o que me custa dígitos.
Geobits 14/03
1
Você pode salvar outro byte e melhorar a precisão usando "4e5" e "1e-5" para os números.
Vilmantas Baranauskas 16/03/2015
@VilmantasBaranauskas Thanks! Eu sempre esquecer que :) É tentador usar 4E9 e 1e-9 em vez, mas isso leva muito tempo ...
Geobits
Protip: quando golfe, você deve realmente reduzir os bytes, não artificialmente aumentá-los
Destrutível Lemon
3

Javascript: 62 bytes

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

Usei a resposta javascript anterior (agora excluída) e depilei 5 bytes.

Guzman Tierno
fonte
Você pode vincular a resposta da cfern para fornecer o crédito corretamente.
Jonathan Frech
Sua resposta parece ser um trecho de E / S não permitido . Corrija ou exclua sua postagem.
Jonathan Frech
Desculpe, sou novo, não sabia como colocar o link para a solução anterior que agora aparece foi excluída. Em relação ao snippet: concordo plenamente, mas esse foi o código da solução javascript anterior, que também acho inválida por esse motivo. Eu modifiquei o meu para ser um programa.
Guzman Tierno #
Sim; a resposta anterior foi excluída por ser inválida. Vi sua resposta antes de recomendar a exclusão, portanto, o comentário. +1 para enviar uma resposta válida; bem-vindo ao PPCG!
Jonathan Frech
2

GolfScript (34 caracteres)

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

Demonstração online

Isso usa ponto fixo porque o GS realmente não tem ponto flutuante. Abusa levemente o uso do ponto fixo; portanto, se você quiser alterar a contagem de iterações, verifique se é duas vezes a potência de dez.

Crédito para xnor pelo método particular de Monte Carlo empregado.

Peter Taylor
fonte
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

retorna 3.14120037157por exemplo. A contagem da amostra é 4782969 (9 ^ 7). Você pode obter um pi melhor com 9 ^ 9, mas terá que ser paciente.

Cavaleiro Lógico
fonte
Você pode salvar 3 substituindo range(9**7)por [0]*9**7ou algo assim, já que não o usa i. E a lista não é muito longa para encontrar problemas de memória.
Sp3000 15/03/2015
Obrigado. Eu queria me livrar, range()mas havia esquecido completamente esse truque.
Logic Knight
Sinto que a [0]9**7sintaxe não é válida.
seequ
Você está certo. Anexei novamente o asterisco perdido (estava embaixo da minha mesa).
Logic Knight
2

Ruby, 39 bytes

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

Um dos destaques é que este é capaz de usar 8e5notação, o que torna extensível até ~ 8e9 amostras na mesma contagem de bytes de programa.

GreyCat
fonte
2

Perl 6 , 33 bytes

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

Experimente online!

Essa é uma função que recebe o número de amostras como argumento.

Sean
fonte
1

Scala, 87 77 66 bytes

def s=math.pow(math.random,2);Seq.fill(1000)(s+s).count(_<1)/250d
keegan
fonte
Se você substituir 1000com 8000e 250dcom 2e4ambos salvar um byte e aumentar o número de amostras por um fator de 8.
Dave Swartz
1

Pure Bash, 65 bytes

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

Toma um único parâmetro de linha de comando multiplicado por 4 para fornecer o número de amostras. A aritmética do Bash é apenas um número inteiro; portanto, é racional um resultado. Isso pode ser canalizado bc -lpara a divisão final:

$ ./montepi.sh 10000
31477/10000
$ ./montepi.sh 10000|bc -l
3.13410000000000000000
$ 
Trauma Digital
fonte
1

Joe , 20 19 bytes

Nota: esta resposta não é competitiva, porque a versão 0.1.2, que adicionou aleatoriedade, foi lançada após esse desafio.

Função nomeada F:

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

Função sem nome:

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

Ambos pegam a contagem da amostra como argumento e retornam o resultado. Como eles funcionam?

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

Exemplo é executado:

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

dc, 59 caracteres (o espaço em branco é ignorado)

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

Eu testei isso no Plan 9 e no OpenBSD, então imagino que funcione no Linux (GNU?) dc .

Explicação por linha:

  1. Armazena o código para [ler e quadrado dois carros alegóricos; executar registro ise 1 for maior que a soma dos quadrados] no registrou .
  2. Armazena o código para [incrementar o registro xem 1] no registroi .
  3. Armazena o código para [executar o registro u, incrementar o registro me depois executar o registro zse o registro mfor maior que o registro n] no registro z.
  4. Defina a escala para 5 pontos decimais.

  5. Leia o número de pontos a serem amostrados na primeira linha de entrada.
  6. Executar registro z .
  7. Divida o registro x(o número de ocorrências) pelo registron (o número de pontos), multiplique o resultado por 4 e imprima.
  8. Sair.

No entanto, eu trapacei:

O programa precisa de um fornecimento de carros alegóricos aleatórios entre 0 e 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

Execução de teste:

% runpi 10000
3.14840

Agora com menos trapaça (100 bytes)

Alguém apontou que eu poderia incluir uma simples impressão.
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

Ungolfed

[
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

Execução de teste:

$ echo 10000 | dc pigolf.dc
3.13640
progrider42
fonte
1

Pyth, 19

c*4sm<sm^OQ2 2*QQQQ

Forneça o número desejado de iterações como entrada.

Demonstração

Como Pyth não tem uma função "Número flutuante aleatório", tive que improvisar. O programa escolhe dois números inteiros positivos aleatórios menores que a entrada, quadrados, somas e comparados com a entrada ao quadrado. Isso executou um número de vezes igual à entrada e, em seguida, o resultado é multiplicado por 4 e dividido pela entrada.

Em notícias relacionadas, estarei adicionando uma operação aleatória de número de ponto flutuante ao Pyth em breve. Este programa não usa esse recurso, no entanto.


Se interpretarmos "O resultado pode ser retornado ou impresso como um ponto flutuante, ponto fixo ou número racional". liberalmente, a impressão do numerador e denominador da fração resultante deve ser suficiente. Nesse caso:

Pyth, 18

*4sm<sm^OQ2 2*QQQQ

Este é um programa idêntico, com a operação de divisão de ponto flutuante ( c) removida.

isaacg
fonte
1

Julia, 37 bytes

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

O número de amostras é 65536 (= 4 ^ 8).

Uma variante ligeiramente mais longa: uma função com o número de amostras scomo o único argumento:

s->4mean(1-floor(sum(rand(s,2).^2,2)))
pawel.boczarski
fonte
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);}

Ungolfed:

#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
fonte
é claro, você provavelmente ainda deve publicar a versão sem o espaço em branco (mantenha a versão atual com o cabeçalho "ungolfed / with white space" ou algo assim) #
Destructible Lemon
@DestructibleWatermelon done!
Karl Napf
A solução não funciona no GCC sem uma nova linha antes f(). Qual compilador você usou? Veja tio.run/##Pc49C4JAHIDx3U9xGMG9ZdYgwWkgtNbQ1BZ6L/UHO8M07hA/…
eush77
101 bytes
ceilingcat 21/03
1

Na verdade , 14 bytes (não concorrentes)

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

Experimente online!

Esta solução não é competitiva porque o idioma pós-desafio. O número de amostras é dado como entrada (em vez de codificado).

Explicação:

`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
fonte
2
Por que o voto negativo?
Destructible Lemon
1

Raquete 63 bytes

Usando o método da linguagem R, responda por @Matt:

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

Ungolfed:

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

Teste:

(f 10000)

Saída (fração):

3 31491308966059784/243801776017028125

Como decimal:

(exact->inexact(f 10000))

3.13583200307849
rnso
fonte
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

Experimente online!

Este código é muito ruim escreveu. Ele falhará se o gfortran decidir inicializar a variável Acom outro valor, então 0 (que ocorre aproximadamente 50% das compilações) e, seA for inicializado como 0, sempre gerará a mesma sequência aleatória para a semente especificada. Em seguida, o mesmo valor para Pi é impresso sempre.

Este é um programa muito melhor:

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

Experimente online!

(Um byte salvo em cada versão; obrigado Penguino).

rafa11111
fonte
1
Em cada versão, você pode salvar um byte alterando 'DO I = 1,1E3' para 'DO I = 1,4E3', alterando 'A = A + 1' para 'A = A + 1E-3' e alterando ' PRINT *, A / 250 'para' PRINT *, A '
Penguino 02/02
Sim, você tem certeza! Obrigado pela sugestão!
Rafa11111
1

Japt , 26 ou 18 bytes

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

Experimente online!

Análogo à resposta do Optimizer , principalmente apenas tentando aprender japonês.
Leva o número de iterações a serem executadas 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.

Se 1/(1+x^2)for permitido (em vez de dois randoms separados), podemos obter 18 bytes com a mesma lógica.

o r_Ä/(1+Mr pÃ*4/U
Nit
fonte
1
Você pode salvar alguns bytes, deixando Mhcalcular a hipotenusa em vez de fazer você mesmo ;-) Além disso, você pode usar xa soma de uma matriz, em vez de reduzir por adição:o x@MhMr Mr)<1Ã*4/U
ETHproductions
@ETHproductions Legal, eu não sabia que você pode usar Mhassim, obrigado! Sua resposta aleatória é quase tão curta quanto a minha resposta com apenas uma aleatória, isso é bem legal. Eu mantereix em mente que costumo usar muito a redução ao tentar jogar coisas de golfe, então isso será muito ú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

Experimente online!

Tanto quanto eu posso entender, para fazer esse tipo de execução total em F #, é mais curto criar uma matriz de números e usar o método Seq.sumBy método do que usar um for..to..dobloco.

O que esse código faz é que ele cria uma coleção de números de ponto flutuante de 1 a s, executa a função fun x->...para o número de elementos na coleção e soma o resultado. tems elementos na coleção, portanto o teste aleatório é realizado várias svezes. Os números reais na coleção são ignorados ( fun x->, masx não são usados).

Isso também significa que o aplicativo deve primeiro criar e preencher a matriz e, em seguida, iterar sobre ela. Portanto, é provavelmente duas vezes mais lento que umfor..to..do loop. E com a criação da matriz, o uso da memória fica na região de O (f ** k)!

Para o teste em si, em vez de usar uma if then elseinstrução, ele calcula a distância (q()+q()|>Math.Sqrt ) e arredondá-la para baixo Math.Floor. Se a distância estiver dentro do círculo, será arredondado para 0. Se a distância estiver fora do círculo, será arredondado para 1. O Seq.sumBymétodo totalizará esses resultados.

Note então que Seq.sumBy totalizou não são os pontos dentro do círculo, mas os pontos fora dele. Portanto, para o resultado obtido s(nosso tamanho de amostra) e subtrai o total dele.

Parece também que tomar um tamanho de amostra como parâmetro é mais curto do que codificar o valor. Então, eu estou trapaceando um pouco ...

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

Como lidar com isso import System.Random; r=randoms(mkStdGen 2)exigiria muitos bytes preciosos, eu gero uma lista infinita de números aleatórios com o gerador congruencial linear que alguns dizem ser quase criptograficamente forte:x↦x*9+1 mod 8^9 qual, pelo Teorema de Hull-Dobell, tem o período completo 8^9.

g rendimentos 4 se o ponto do número aleatório estiver dentro do círculo para pares de números aleatórios [0..8^9-1]porque isso elimina uma multiplicação na fórmula usada.

Uso:

> p 100000
3.14208

Experimente online!

Angs
fonte
1

Perl 5, 34 bytes

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

O número de amostras é retirado do stdin. Requer-p .

Funciona porque:

Experimente online!

primo
fonte