Счастливого Пи Дня всем! Без всякой причины я пытаюсь построить оценку Пи по методу Монте-Карло, которая будет максимально короткой. Можем ли мы построить тот, который может вписаться в твит?

Чтобы уточнить, что я имею в виду, это типичный подход вытягивания случайных точек из единичного квадрата и вычисления отношения, попадающего в единичный круг. Количество образцов может быть жестко закодировано или нет. Если вы жестко закодировали их, вы должны использовать не менее 1000 образцов. Результат может быть возвращен или напечатан как число с плавающей запятой, фиксированная точка или рациональное число.

Никакие функции триггера или константы Пи не должны быть подходом Монте-Карло.

Это код гольф, поэтому выигрывает самое короткое представление (в байтах).

80386 машинный код, 40 38 байт

Hexdump кода:

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

Как получить этот код (с языка ассемблера):

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

Это функция, использующая fastcallсоглашение о вызовах MS (количество итераций передается в регистре ecx). Возвращает результат в stрегистр.

Забавные вещи об этом коде:

• rdrand - всего 3 байта для генерации случайного числа!
• Он использует (без знака) целочисленную арифметику до последнего деления.
• Сравнение квадрата расстояния ( D) с квадратом радиуса ( 2^32) выполняется автоматически - флаг переноса содержит результат.
• Чтобы умножить счет на 4, он считает выборки с шагом 4.

Matlab / Octave, 27 байт

Я знаю, что уже есть ответ Matlab / Octave, но я попробовал свой собственный подход. Я использовал тот факт, что интеграл от 4/(1+x^2)0 до 1 является пи.

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

R, 40 (или 28 или 24 с использованием других методов)

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

Еще один Python, если NumPy разрешен, но очень похож на Matlab / Octave:

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

Mathematica, 42 40 39 байт (или 31/29?)

У меня есть три решения по 42 байта:

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

Все они являются неназванными функциями, которые берут количество выборок nи возвращают рациональное приближение π. Сначала все они генерируют nточки в единичном квадрате в положительном квадранте. Затем они определяют количество тех выборок, которые лежат в единичном круге, а затем делят их на количество выборок и умножают на 4. Разница лишь в том, как они определяют количество образцов в круге единицы:

• Первый использует Countс условием, что Norm[p] < 1.
• Второй вычитает норму каждой точки из, 1а затем округляет. Это превращает числа внутри круга устройства в 1и те, что снаружи 0. После этого я просто подвожу их всех Tr.
• Третий делает то же самое, но вычитает из 1.5, так что я могу использовать Roundвместо Ceiling.

Во время написания этой статьи мне пришло в голову, что действительно есть более короткое решение, если я просто вычту 2и затем использую Floor:

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

или сохранение другого байта с использованием операторов Unicode для настилов или потолков:

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

Обратите внимание, что три решения на основе округления также могут быть записаны с использованием Meanвместо Trи без /#, опять же для тех же байтов.

Если другие подходы, основанные на Монте-Карло, хороши (в частности, тот, который выбрал Питер), я могу сделать 31 байт, оценивая интеграл или 29, используя интеграл , на этот раз в виде числа с плавающей запятой:√(1-x2)1/(1+x2)

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

CJam, 27 23 22 или 20 байтов

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

2 байта сохранены благодаря Runner112, 1 байт сохранен благодаря Sp3000

Он принимает количество итераций из STDIN в качестве входных данных.

Это так же просто, как и получается. Это основные этапы:

• Прочитайте ввод и выполните многократные итерации Монте-Карло
• В каждой итерации получайте сумму квадратов двух случайных чисел с 0 по 1 и смотрите, будет ли она меньше 1
• Получите отношение количества раз, которое мы получили меньше 1 к общему количеству итераций, и умножьте его на 4, чтобы получить PI

Расширение кода :

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

Попробуйте онлайн здесь

Если среднее значение 1/(1+x^2)также считается Монте-Карло, то это можно сделать в 20 байтах:

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

Попробуй здесь

Python 2, 77 75 байт

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

Использует 4000 образцов для сохранения байтов 1e3.

Commodore 64 Basic, 45 байт

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

PETSCII замены: ─= SHIFT+E, /= SHIFT+N, ┌=SHIFT+O

Создает 1000 очков в первом квадранте; для каждого добавляет правдивость «x ^ 2 + y ^ 2 <1» к текущему счету, а затем делит счет на 250, чтобы получить pi. (Наличие знака минус объясняется тем, что на C64 "true" = -1.)

J, 17 байт

Вычисляет среднее значение 40000значений выборки функции 4*sqrt(1-sqr(x))в диапазоне [0,1].

Легко 0 o.xвозвращается sqrt(1-sqr(x)).

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

> <> (Рыба) , 114 байт

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

Теперь> <> не имеет встроенного генератора случайных чисел. Однако у него есть функция, которая отправляет указатель в случайном направлении. Генератор случайных чисел в моем коде:

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

Он в основном генерирует случайные биты, которые составляют двоичное число, а затем преобразует это случайное двоичное число в десятичное.

Все остальное - просто правильные точки в квадрате.

Использование: когда вы запускаете код, вы должны предварительно заполнить стек (-v в интерпретаторе Python), например, количеством выборок.

pi.fish -v 1000

возвращается

3.164

Matlab или Octave 29 байт (благодаря flawr!)

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

(Я не совсем уверен, что <1 в порядке. Я прочитал, что должно быть <= 1. Но насколько велика вероятность нарисовать ровно 1 ...)

Matlab или Octave 31 байт

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

Java, 108 байт

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

Четыре тысячи итераций, добавляя 0,001 каждый раз, когда точка находится внутри единичного круга. Довольно простые вещи.

Примечание. Да, я знаю, что могу сбросить четыре байта, изменив πоднобайтовый символ. Мне нравится это так.

Javascript: 62 байта

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

Я использовал предыдущий (теперь удаленный) ответ javascript и побрил 5 байтов.

GolfScript (34 символа)

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

Онлайн демо

При этом используется фиксированная точка, потому что GS на самом деле не имеет плавающей точки. Это немного злоупотребляет использованием фиксированной точки, поэтому, если вы хотите изменить количество итераций, убедитесь, что оно в два раза больше десяти.

Кредит XNOR для конкретного метода Монте - Карло , используемого.

Python 2, 90 85 81 байт

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

возвращается 3.14120037157к примеру. Количество образцов составляет 4782969 (9 ^ 7). Вы можете получить лучший пи с 9 ^ 9, но вам придется набраться терпения.

Рубин, 39 байт

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

Одним из основных моментов является то, что этот может использовать 8e5нотацию, что делает его расширяемым до ~ 8e9 выборок при одном и том же количестве байтов программы.

Perl 6 , 33 байта

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

Попробуйте онлайн!

Это функция, которая принимает количество выборок в качестве аргумента.

Scala, 87 77 66 байт

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

Чистый Баш, 65 байт

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

Принимает один параметр командной строки, который умножается на 4, чтобы получить количество выборок. Арифметика Bash только целочисленная, поэтому рациональным является вывод. Это может быть передано bc -lдля окончательного разделения:

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

Джо , 20 19 байт

Примечание: этот ответ не является конкурирующим, потому что версия 0.1.2, которая добавила случайность, была выпущена после этого испытания.

Именованная функция F:

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

Безымянная функция:

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

Они оба принимают количество отсчетов в качестве аргумента и возвращают результат. Как они работают?

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

Пример работы:

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

постоянный ток, 59 символов (пробел игнорируется)

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

Я проверил это на Plan 9 и OpenBSD, поэтому я думаю, что он будет работать на Linux (GNU?) dc.

Пояснения к строке:

1. Сохраняет код для [чтения и возведения в квадрат двух чисел; выполнить регистр, iесли 1 больше, чем сумма их квадратов] в регистре u.
2. Сохраняет код в [увеличение регистра xна 1] в регистре i.
3. Сохраняет код [выполнить регистр u, увеличить регистр m, а затем выполнить регистр, zесли регистр mбольше регистра n] в регистре z.

5. Установите масштаб до 5 десятичных знаков.

6. Считайте количество точек для выборки из первой строки ввода.
7. Выполнить регистрацию z.
8. Разделите регистр x(количество попаданий) на регистр n(количество баллов), умножьте результат на 4 и напечатайте.
9. Уволиться.

Однако я обманул

Программе нужен набор случайных чисел от 0 до 1.

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

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

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

Использование:

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

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

Тестовый забег:

% runpi 10000
3.14840

Теперь с меньшим количеством читов (100 байт)

Кто-то указал, что я мог бы включить простой 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

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

Тестовый забег:

$ echo 10000 | dc pigolf.dc
3.13640

Пиф, 19

c*4sm<sm^OQ2 2*QQQQ

Введите желаемое количество итераций в качестве входных данных.

демонстрация

Поскольку у Пита нет функции «Случайное число с плавающей точкой», мне пришлось импровизировать. Программа выбирает два случайных натуральных числа меньше входных, квадратов, сумм и сравнивает их с квадратом ввода. Это выполняется количество раз, равное входному значению, затем результат умножается на 4 и делится на входное значение.

В связанных новостях я скоро добавлю случайную операцию с плавающей запятой в Pyth. Однако эта программа не использует эту функцию.

Если мы интерпретируем «Результат может быть возвращен или напечатан как число с плавающей запятой, фиксированная точка или рациональное число». Обильно, тогда печати числителя и знаменателя полученной дроби должно быть достаточно. В таком случае:

Пиф, 18

*4sm<sm^OQ2 2*QQQQ

Это идентичная программа, с cудаленной операцией деления с плавающей запятой ( ).

Юлия, 37 байт

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

Количество образцов составляет 65536 (= 4 ^ 8).

Немного более длинный вариант: функция с числом выборок sв качестве единственного аргумента:

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

C, 130 байтов

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

На самом деле , 14 байтов (не конкурирующих)

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

Попробуйте онлайн!

Это решение неконкурентное, потому что язык ставит после проблемы. Количество выборок дается как ввод (а не жестко запрограммировано).

Объяснение:

`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

Ракетка 63 байта

Используя метод ответа языка R @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))

Тестирование:

(f 10000)

Выход (дробь):

3 31491308966059784/243801776017028125

В десятичном виде:

(exact->inexact(f 10000))

3.13583200307849

Фортран (GFortran) , 84 83 байта

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

Попробуйте онлайн!

Этот код очень плохо написал. Он потерпит неудачу, если gfortran решит инициализировать переменную Aс другим значением, отличным от 0 (что происходит примерно в 50% сборок ), и, если Aинициализируется как 0, он всегда будет генерировать одну и ту же случайную последовательность для данного начального числа. Затем всегда выводится одно и то же значение для Pi.

Это намного лучшая программа:

Фортран (GFortran) , 100 99 байт

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

Попробуйте онлайн!

(Один байт сохранен в каждой версии; спасибо Penguino).

Japt , 26 или 18 байт

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

Попробуйте онлайн!

Аналогично ответу Оптимизатора , в основном просто пытаюсь выучить Джапта.
Принимает количество итераций для неявного ввода 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.

Если 1/(1+x^2)разрешено (вместо двух отдельных случайностей), то мы можем получить 18 байтов с той же логикой.

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

F #, 149 байт

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

Попробуйте онлайн!

Насколько я могу понять, для выполнения такого рода промежуточного итога в F # создать массив чисел и использовать Seq.sumByметод короче, чем использоватьfor..to..do блок.

Что делает этот код, так это то, что он создает коллекцию чисел с плавающей запятой от 1 до s, выполняет функцию fun x->...для количества элементов в коллекции и суммирует результат. В sколлекции есть элементы, поэтому случайный тест проводится sраз. Фактические числа в коллекции игнорируются (fun x-> но xне используются).

Это также означает, что приложение должно сначала создать и заполнить массив, а затем выполнить итерацию по нему. Так что, вероятно, в два раза медленнее, чемfor..to..do цикл. А при использовании массива создание памяти находится в области O (f ** k)!

Для самого теста вместо использования if then elseоператора он вычисляет расстояние ( q()+q()|>Math.Sqrt) и округляет его Math.Floor. Если расстояние находится внутри круга, оно будет округлено до 0. Если расстояние находится за пределами круга, оно будет округлено до 1.Seq.sumBy метод суммирует эти результаты.

Обратите внимание, что Seq.sumByв итоге мы имеем не точки внутри круга, а точки за его пределами . Так что для результата требуетсяs (наш размер выборки) и вычитает из него сумму.

Также представляется, что выбор размера выборки в качестве параметра короче, чем жесткое кодирование значения. Так что я немного изменяю ...

Хаскелл, 116 114 110 96 байт

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

Потому что иметь дело с import System.Random; r=randoms(mkStdGen 2) заняла бы слишком много драгоценных байтов, я генерирую бесконечный список случайных чисел с линейным конгруэнтным генератором, который, по мнению некоторых, почти криптографически силен: x↦x*9+1 mod 8^9по теореме Халла-Добелла полный период равен 8^9.

gдает, 4если точка случайного числа находится внутри круга для пар случайных чисел в[0..8^9-1] потому что это исключает умножение в используемой формуле.

Использование:

> p 100000
3.14208

Попробуйте онлайн!

Perl 5, 34 байта

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

Количество образцов берется из стандартного ввода. Требуется -p.

Работает потому что:

Попробуйте онлайн!