Wszystkiego najlepszego z okazji Dnia Pi! Bez żadnego powodu staram się zbudować estymator Monte Pi Carlo tak krótki, jak to możliwe. Czy możemy zbudować taki, który zmieści się w tweecie?

Aby wyjaśnić, mam na myśli typowe podejście do rysowania losowych punktów z kwadratu jednostkowego i obliczania stosunku, który mieści się w okręgu jednostkowym. Liczba próbek może być zakodowana na stałe lub nie. Jeśli je kodujesz na stałe, musisz użyć co najmniej 1000 próbek. Wynik może zostać zwrócony lub wydrukowany jako zmiennoprzecinkowy, stały punkt lub liczba wymierna.

Żadnych funkcji triggera ani stałych Pi, musi być podejście Monte Carlo.

To jest kod golfowy, więc wygrywa najkrótsze przesłanie (w bajtach).

Kod maszynowy 80386, 40 38 bajtów

Hexdump kodu:

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

Jak zdobyć ten kod (z języka asemblera):

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

Jest to funkcja korzystająca z fastcallkonwencji wywoływania MS (liczba iteracji jest przekazywana do rejestru ecx). Zwraca wynik w strejestrze.

Zabawne rzeczy na temat tego kodu:

• rdrand - zaledwie 3 bajty, aby wygenerować losową liczbę!
• Używa (niepodpisanej) arytmetyki liczb całkowitych do ostatecznego podziału.
• Porównanie kwadratowej odległości ( D) z kwadratowym promieniem ( 2^32) odbywa się automatycznie - flaga przenoszenia zawiera wynik.
• Aby pomnożyć liczbę przez 4, zlicza próbki w krokach co 4.

Matlab / Octave, 27 bajtów

Wiem, że już istnieje odpowiedź Matlab / Octave, ale wypróbowałem własne podejście. Wykorzystałem fakt, że całka 4/(1+x^2)między 0 a 1 to pi.

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

R, 40 (lub 28 lub 24 przy użyciu innych metod)

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

Kolejny Python, jeśli dozwolone jest numpy, ale dość podobny do Matlab / Octave:

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

Mathematica, 42 40 39 bajtów (lub 31/29?)

Mam trzy rozwiązania o wielkości 42 bajtów:

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

Wszystkie są nienazwanymi funkcjami, które pobierają liczbę próbek ni zwracają racjonalne przybliżenie π. Najpierw wszystkie generują npunkty na kwadracie jednostki w kwadrancie dodatnim. Następnie określają liczbę próbek, które znajdują się w okręgu jednostkowym, a następnie dzielą przez liczbę próbek i mnożą przez 4. Jedyna różnica polega na tym, w jaki sposób określają liczbę próbek w okręgu jednostkowym:

• Pierwszy używa Countpod warunkiem, że Norm[p] < 1.
• Drugi odejmuje normę każdego punktu, 1a następnie zaokrągla w górę. To zmienia liczby wewnątrz okręgu jednostki na 1i na zewnątrz na 0. Potem podsumowuję je wszystkie Tr.
• Trzeci robi w zasadzie to samo, ale odejmuje od 1.5, więc mogę użyć Roundzamiast Ceiling.

Aaaaaand podczas pisania tej góry, to dotarło do mnie, że rzeczywiście istnieje krócej rozwiązanie, jeśli po prostu odjąć od 2a następnie użyć Floor:

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

lub zapisywanie innego bajtu za pomocą operatorów podłogowych lub sufitowych Unicode:

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

Należy zauważyć, że trzy rozwiązania oparte zaokrąglenia można również pisane Meanzamiast Tri bez /#, znowu dla tych samych bajtów.

Jeśli inne podejścia oparte na Monte Carlo są w porządku (konkretnie ten, który wybrał Peter), mogę zrobić 31 bajtów poprzez oszacowanie całki lub 29 za pomocą całki , tym razem podanej jako liczba zmiennoprzecinkowa:√(1-x2)1/(1+x2)

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

CJam, 27 23 22 lub 20 bajtów

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

2 bajty zapisane dzięki Runner112, 1 bajt zapisany dzięki Sp3000

Pobiera na wejściu liczbę iteracji ze STDIN.

Jest to tak proste, jak to możliwe. Oto najważniejsze kroki:

• Przeczytaj dane wejściowe i uruchom iteracje Monte Carlo tyle razy
• W każdej iteracji uzyskaj sumę kwadratu dwóch losowych liczb zmiennoprzecinkowych od 0 do 1 i sprawdź, czy jest mniejsza niż 1
• Uzyskaj stosunek tego, ile razy otrzymaliśmy mniej niż 1 przez sumę iteracji i pomnóż go przez 4, aby uzyskać PI

Rozszerzenie kodu :

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

Wypróbuj online tutaj

Jeśli średnia wartość 1/(1+x^2)jest również uważana za Monte Carlo, można to zrobić w 20 bajtach:

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

Wypróbuj tutaj

Python 2, 77 75 bajtów

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

Używa 4000 próbek do zapisania bajtów 1e3.

Commodore 64 Basic, 45 bajtów

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

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

Generuje 1000 punktów w pierwszej ćwiartce; dla każdego dodaje prawdziwość „x ^ 2 + y ^ 2 <1” do bieżącej liczby, a następnie dzieli liczbę przez 250, aby uzyskać pi. (Obecność znaku minus wynika z tego, że na C64 „true” = -1.)

J, 17 bajtów

Oblicza średnią wartość 40000przykładowych wartości funkcji 4*sqrt(1-sqr(x))w zakresie [0,1].

Łatwo 0 o.xwraca sqrt(1-sqr(x)).

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

> <> (Ryba) , 114 bajtów

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

Teraz> <> nie ma wbudowanego generatora liczb losowych. Ma jednak funkcję, która wysyła wskaźnik w losowym kierunku. Generator liczb losowych w moim kodzie:

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

Zasadniczo generuje losowe bity, które tworzą liczbę binarną, a następnie konwertuje tę losową liczbę binarną na dziesiętną.

Reszta to zwykłe punkty w podejściu kwadratowym.

Sposób użycia: po uruchomieniu kodu należy upewnić się, że stos został wstępnie wypełniony (-v w interpreterie w języku Python), na przykład liczbą próbek, na przykład

pi.fish -v 1000

zwraca

3.164

Matlab lub Octave 29 bajtów (dzięki flawr!)

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

(Nie jestem pewien, czy <1 jest w porządku. Przeczytałem, że powinno być <= 1. Ale jak duże jest prawdopodobieństwo narysowania dokładnie 1 ...)

Matlab lub Octave 31 bajtów

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

Java, 108 bajtów

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

Cztery tysiące iteracji, dodając 0,001 za każdym razem, gdy punkt znajduje się wewnątrz okręgu jednostki. Całkiem podstawowe rzeczy.

Uwaga: Tak, wiem, że mogę zrzucić cztery bajty, zmieniając πznak na jednobajtowy. Podoba mi się to w ten sposób.

JavaScript: 62 bajtów

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

Użyłem poprzedniej (teraz usuniętej) odpowiedzi javascript i ogoliłem 5 bajtów.

GolfScript (34 znaki)

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

Demo online

Wykorzystuje stały punkt, ponieważ GS tak naprawdę nie ma zmiennoprzecinkowego. Lekko nadużywa użycia punktu stałego, więc jeśli chcesz zmienić liczbę iteracji, upewnij się, że jest to dwukrotność potęgi dziesięciu.

Kredyt dla XNOR dla danej metody Monte Carlo zatrudniony.

Python 2, 90 85 81 bajtów

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

zwraca 3.14120037157na przykład. Liczba próbek wynosi 4782969 (9 ^ 7). Możesz uzyskać lepsze pi z 9 ^ 9, ale będziesz musiał uzbroić się w cierpliwość.

Rubin, 39 bajtów

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

Jedną z głównych zalet jest to, że można w nim stosować 8e5notację, dzięki czemu można go rozszerzyć do ~ 8e9 próbek w tej samej liczbie bajtów programu.

Perl 6 , 33 bajtów

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

Wypróbuj online!

Jest to funkcja, która przyjmuje liczbę próbek jako argument.

Scala, 87 77 66 bajtów

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

Pure Bash, 65 bajtów

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

Pobiera pojedynczy parametr wiersza polecenia, który jest mnożony przez 4, aby podać liczbę próbek. Arytmetyka Bash jest tylko liczbą całkowitą, więc wynik jest racjonalny. Można to przesłać do bc -lkońcowego podziału:

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

Joe , 20 lat 19 bajtów

Uwaga: ta odpowiedź nie jest konkurencyjna, ponieważ wersja 0.1.2, która dodała losowość, została wydana po tym wyzwaniu.

Nazwana funkcja F:

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

Funkcja bez nazwy:

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

Oba biorą próbkę jako argument i zwracają wynik. Jak oni pracują?

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

Przykładowe przebiegi:

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

dc, 59 znaków (białe znaki są ignorowane)

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

Przetestowałem to na Planie 9 i OpenBSD, więc wyobrażam sobie, że będzie działać na Linuksie (GNU?) dc .

Objaśnienie według wiersza:

1. Przechowuje kod w [odczyt i kwadrat dwóch pływaków; wykonaj rejestr, ijeśli 1 jest większy niż suma ich kwadratów] w rejestrzeu .
2. Przechowuje kod w [inkrementacji rejestru xo 1] w rejestrzei .
3. Przechowuje kod w [wykonaj rejestr u, rejestr przyrostowy m, a następnie uruchom rejestr, zjeśli rejestr mjest większy niż rejestr n] w rejestrze z.

5. Ustaw skalę na 5 miejsc po przecinku.

6. Odczytaj liczbę punktów do próbkowania z pierwszego wiersza wprowadzania.
7. Wykonaj rejestr z .
8. Podziel rejestr x(liczbę trafień) według rejestrun (liczbę punktów), pomnóż wynik przez 4 i wydrukuj.
9. Porzucić.

Oszukałem jednak:

Program potrzebuje zapasów liczb losowych od 0 do 1.

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

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

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

Stosowanie:

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

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

Testowe uruchomienie:

% runpi 10000
3.14840

Teraz z mniejszym oszustwem (100 bajtów)

Ktoś zwrócił uwagę, że mogę załączyć prosty 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

Nie golfił

[
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

Testowe uruchomienie:

$ echo 10000 | dc pigolf.dc
3.13640

Pyth, 19 lat

c*4sm<sm^OQ2 2*QQQQ

Podaj żądaną liczbę iteracji jako dane wejściowe.

Demonstracja

Ponieważ Pyth nie ma funkcji „Losowa liczba zmiennoprzecinkowa”, musiałem improwizować. Program wybiera dwie losowe liczby całkowite dodatnie mniejsze niż wartość wejściowa, kwadraty, sumy i porównywane z kwadratem wejściowym. Wykonano to kilka razy równe wartości wejściowej, a następnie wynik pomnożono przez 4 i podzielono przez wartość wejściową.

W pokrewnych wiadomościach wkrótce dodam do Pytha losową liczbę zmiennoprzecinkową. Ten program nie korzysta jednak z tej funkcji.

Jeśli interpretujemy „Wynik może zostać zwrócony lub wydrukowany jako zmiennoprzecinkowy, stały punkt lub liczba wymierna”. swobodnie, wystarczy wydrukować licznik i mianownik powstałej frakcji. W tym wypadku:

Pyth, 18 lat

*4sm<sm^OQ2 2*QQQQ

Jest to identyczny program, z cusuniętą operacją dzielenia zmiennoprzecinkowego ( ).

Julia, 37 bajtów

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

Liczba próbek wynosi 65536 (= 4 ^ 8).

Odrobinę dłuższy wariant: funkcja z liczbą próbek sjako jedynym argumentem:

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

C, 130 bajtów

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

Nie golfowany:

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

Właściwie 14 bajtów (niekonkurujące)

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

Wypróbuj online!

To rozwiązanie nie konkuruje, ponieważ język jest późniejszy niż wyzwanie. Liczba próbek jest podawana jako dane wejściowe (a nie na stałe).

Wyjaśnienie:

`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

Rakieta 63 bajty

Przy użyciu metody odpowiedzi w języku R autorstwa @Matt:

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

Nie golfowany:

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

Testowanie:

(f 10000)

Wyjście (ułamek):

3 31491308966059784/243801776017028125

W systemie dziesiętnym:

(exact->inexact(f 10000))

3.13583200307849

Fortran (GFortran) , 84 83 bajtów

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

Wypróbuj online!

Ten kod jest bardzo źle napisany. Nie powiedzie się, jeśli gfortran zdecyduje się zainicjować zmienną Ainną wartością niż 0 (co występuje w przybliżeniu w 50% kompilacji) i, jeśliA zostanie zainicjowany jako 0, zawsze wygeneruje tę samą losową sekwencję dla danego ziarna. Następnie zawsze drukowana jest ta sama wartość dla Pi.

To jest znacznie lepszy program:

Fortran (GFortran) , 100 99 bajtów

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

Wypróbuj online!

(Jeden bajt zapisany w każdej wersji; dzięki Penguino).

Japt , 26 lub 18 bajtów

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

Wypróbuj online!

Analogicznie do odpowiedzi Optimizera , głównie po prostu próbując nauczyć się Japt.
Pobiera domyślną liczbę iteracji do uruchomienia 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.

Jeśli 1/(1+x^2)jest dozwolone (zamiast dwóch oddzielnych losowych), możemy osiągnąć 18 bajtów przy użyciu tej samej logiki.

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

F #, 149 bajtów

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

Wypróbuj online!

O ile mi wiadomo, aby wykonać tego rodzaju sumę całkowitą w F #, krótsze jest utworzenie tablicy liczb i użycie Seq.sumBymetody niż użyciefor..to..do bloku.

Co robi ten kod, że tworzy zbiór liczb zmiennoprzecinkowych od 1 do s , wykonuje funkcję fun x->...dla liczby elementów w kolekcji i sumuje wynik. W skolekcji znajdują się elementy, więc losowy test jest wykonywany srazy. Rzeczywiste liczby w kolekcji są ignorowane ( fun x->ale xnie są używane).

Oznacza to również, że aplikacja musi najpierw utworzyć i wypełnić tablicę, a następnie wykonać iterację. Jest więc prawdopodobnie dwa razy wolniejszy niż for..to..dopętla. A przy tworzeniu tablicy wykorzystanie pamięci jest w obszarze O (f ** k)!

Dla samego testu, zamiast używać if then elseinstrukcji, robi to, oblicza odległość ( q()+q()|>Math.Sqrt) i zaokrągla ją w dół Math.Floor. Jeśli odległość mieści się w okręgu, zostanie zaokrąglona w dół do 0. Jeśli odległość znajduje się poza okręgiem, zostanie zaokrąglona w dół do 1. Seq.sumByMetoda następnie sumuje te wyniki.

Zauważ, że Seq.sumBysumą nie są punkty wewnątrz koła, ale punkty poza nim. Więc dla wyniku potrzeba s(nasz rozmiar próbki) i odejmuje od niego sumę.

Wydaje się również, że przyjęcie wielkości próbki jako parametru jest krótsze niż zakodowanie na stałe wartości. Więc trochę oszukuję ...

Haskell, 116 114 110 96 bajtów

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

Ponieważ zajmowanie się import System.Random; r=randoms(mkStdGen 2)zajęłoby zbyt wiele cennych bajtów, generuję nieskończoną listę liczb losowych za pomocą liniowego generatora kongruencjalnego, który według niektórych jest prawie kryptograficznie silny: x↦x*9+1 mod 8^9który według twierdzenia Hulla-Dobella ma pełny okres 8^9.

gdaje, 4jeśli punkt liczb losowych znajduje się wewnątrz koła dla par liczb losowych w[0..8^9-1] ponieważ eliminuje to mnożenie w stosowanej formule.

Stosowanie:

> p 100000
3.14208

Wypróbuj online!

Perl 5, 34 bajtów

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

Liczba próbek jest pobierana ze standardowego wejścia. Wymaga-p .

Działa, ponieważ:

Wypróbuj online!