# Monte Carlo estimator of Pi

Happy Pi Day everyone! For no reason at all, I'm trying to construct a Monte Carlo estimator of Pi that is as short as possible. Can we construct one that can fit in a tweet?

To clarify, what I have in mind is the typical approach of drawing random points from the unit square and calculating the ratio that fall within the unit circle. The number of samples can be hard coded or not. If you hardcode them, you must use at least 1000 samples. The result may be returned or printed as a floating point, fixed point or rational number.

No trig functions or Pi constants, must be a Monte Carlo approach.

This is code golf, so the shortest submission (in bytes) wins.

• are trig functions allowed? I suggest you explicitly ban them. Commented Mar 14, 2015 at 16:08
3.14915

## ><> (Fish), 114 bytes

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

Now, ><> doesn't have a built-in random number generator. It does however have a function that sends the pointer in a random direction. The random number generator in my code:

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

It basically generates random bits that make up a binary number and then converts that random binary number to decimal.

The rest is just the regular points in the square approach.

Usage: when you run the code you must make sure to prepopulate the stack (-v in python interpreter) with the number of samples, for example

pi.fish -v 1000

returns

3.164

Matlab or Octave 29 bytes (thanks to flawr!)

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

(I am no quite sure if <1 is OK. I read it should be <=1. But how big is the probability to draw exactly 1...)

Matlab or Octave 31 bytes

sum(sum(rand(2,4e3).^2)<=1)/1e3
• Very nice idea! You could save two additional bytes with mean(sum(rand(2,4e6).^2)<1)*4. Commented Mar 15, 2015 at 22:20

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

Four thousand iterations, adding 0.001 each time the point is inside the unit circle. Pretty basic stuff.

Note: Yes, I know I can shed four bytes by changing π to a single-byte character. I like it this way.

• why not 9999 iterations ? Commented Mar 14, 2015 at 22:28
• @Optimizer It makes the sum shorter. For 9999 iterations, I'd have to add a more precise number each time instead, which costs me digits. Commented Mar 14, 2015 at 22:34
• You could save another byte and improve precision by using "4e5" and "1e-5" for the numbers. Commented Mar 16, 2015 at 9:12
• @VilmantasBaranauskas Thanks! I always forget about that :) It's tempting to use 4e9 and 1e-9 instead, but that takes quite a while... Commented Mar 16, 2015 at 13:38
• Protip: when golfing, you should actually reduce the bytes, not artificially increase them Commented Oct 28, 2016 at 2:13

## Javascript: 62 Bytes

I used the previous (now deleted) javascript answer and shaved 5 bytes.

• You could link to cfern's answer to properly give credit. Commented Sep 30, 2018 at 23:57
• Your answer appears to be a snippet which is disallowed I/O. Please either fix or delete your post. Commented Oct 1, 2018 at 3:08
• Sorry, I'm new I didn't know how to put the link to the previous solution which now appears have been deleted. Regarding the snippet: I completely agree, but that was the code of the previous javascript solution which I also think it was invalid that reason. I modified mine to be a program. Commented Oct 2, 2018 at 12:35
• Yes; the previous answer was deleted as it was invalid -- I saw your answer before I recommended deletion, thus the comment. +1 for submitting a valid answer; welcome to PPCG! Commented Oct 2, 2018 at 15:54

## GolfScript (34 chars)

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

Online demo

This uses fixed point because GS doesn't really have floating point. It slightly abuses the use of fixed point, so if you want to change the iteration count make sure that it's twice a power of ten.

Credit to xnor for the particular Monte Carlo method employed.

# Python 2, 9085 81 bytes

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

returns 3.14120037157 for example. The sample count is 4782969 (9^7). You can get a better pi with 9^9 but you will have to be patient.

• You can save 3 by replacing range(9**7) with [0]*9**7 or something, since you don't use i. And the list isn't too long to run into memory issues. Commented Mar 15, 2015 at 0:20
• Thanks. I wanted to get rid of range() but I had completely forgotten that trick. Commented Mar 15, 2015 at 6:26
• I have a feeling [0]9**7 isn't valid syntax. Commented Mar 15, 2015 at 13:57
• You are right. I have re-attached the lost asterix (it was under my desk). Commented Mar 15, 2015 at 14:51

# Ruby, 39 bytes

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

One of the highlight is that this one is able to use 8e5 notation, which makes it extendable up to ~8e9 samples in the same program byte count.

# Perl 6, 33 bytes

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

Try it online!

This is a function that takes the number of samples as an argument.

# Scala, 8777 66 bytes

def s=math.pow(math.random,2);Seq.fill(1000)(s+s).count(_<1)/250d
• If you replace 1000 with 8000 and 250d with 2e4 you both save a byte and increase the number of samples by a factor of 8. Commented Mar 15, 2015 at 13:35

# Pure Bash, 65 bytes

for((;i++<$1*4;a+=RANDOM**2+RANDOM**2<32767**2));{ :;} echo$a/$1 Takes a single command-line parameter which is multiplied by 4 to give the number of samples. Bash arithmetic is integer-only, so a rational is output. This may be piped to bc -l for the final division:$ ./montepi.sh 10000
31477/10000
$./montepi.sh 10000|bc -l 3.13410000000000000000$

# Joe, 20 19 bytes

Note: this answer is non-competing, because version 0.1.2, which added randomness, was released after this challenge.

Named function F:

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

These both take the sample count as an argument and return the result. How do they work?

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

Example runs:

% runpi 10000
3.14840

# Now with less cheating (100 bytes)

Someone pointed out that I could include a 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

## 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 # Pyth, 19 c*4sm<sm^OQ2 2*QQQQ Give the desired number of iterations as input. Demonstration Since Pyth doesn't have a "Random floating number" function, I had to improvise. The program choose two random positive integers less than the input, squares, sums, and compared to the input squared. This performed a number of times equal to the input, then the result is multiplied by 4 and divided by the input. In related news, I will be adding a random floating point number operation to Pyth shortly. This program does not use that feature, however. If we interpret "The result may be returned or printed as a floating point, fixed point or rational number." liberally, then printing the numerator and denominator of the resulting fraction should be sufficient. In that case: # Pyth, 18 *4sm<sm^OQ2 2*QQQQ This is an identical program, with the floating point division operation (c) removed. # Julia, 37 bytes 4mean(1-floor(sum(rand(4^8,2).^2,2))) The number of samples is 65536 (=4^8). A sligthly longer variant: a function with number of samples s as the only argument: 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);} 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); } • of course, you should still probably post the version without the whitespace (keep the current version with the header "ungolfed/with whitespace" or something) Commented Oct 28, 2016 at 0:46 • @DestructibleWatermelon done! Commented Oct 28, 2016 at 1:33 • The solution doesn't work in GCC without a newline before f(). What compiler did you use? See tio.run/##Pc49C4JAHIDx3U9xGMG9ZdYgwWkgtNbQ1BZ6L/UHO8M07hA/… Commented Oct 18, 2018 at 4:04 • 101 bytes Commented Mar 21, 2019 at 16:59 ## Racket 63 bytes Using method of R language answer by @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)) Testing: (f 10000) Output (fraction): 3 31491308966059784/243801776017028125 As 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 Try it online! This code is very bad wrote. It will fail if gfortran decide to initialize variable A with other value then 0 (which occurs 50% of the compilations, approximately) and, if A is initialized as 0, it will always generate the same random sequence for the given seed. Then, the same value for Pi is printed always. This is a much better program: # 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 Try it online! (One byte saved in each version; thanks Penguino). • In each version you can save a byte by changing 'DO I=1,1E3' to 'DO I=1,4E3', changing 'A=A+1' to 'A=A+1E-3', and changing 'PRINT*,A/250' to 'PRINT*,A' Commented Apr 2, 2018 at 22:14 • Yes, you are sure! Thank for the suggestion! Commented Apr 2, 2018 at 22:24 # Japt, 26 or 18 bytes o r_+ÂMhMr p +Mr p <1Ã*4/U Try it online! Analogous to Optimizer's answer, mainly just trying to learn Japt. Takes the number of iterations to run for as the implicit input 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. If 1/(1+x^2) is allowed (instead of two separate randoms), then we can achieve 18 bytes with the same logic. o r_Ä/(1+Mr pÃ*4/U • You can save a few bytes by letting Mh calculate the hypotenuse rather than doing it yourself ;-) Also, you can use x to take the sum of an array, rather than reducing by addition: o x@MhMr Mr)<1Ã*4/U Commented Apr 6, 2018 at 14:32 • @ETHproductions Neat, I didn't know you can use Mh like that, thanks! Your two-random answer is nearly as short as my answer with only one random, that's pretty cool. I'll keep x in mind, I tend to use reduction a lot when trying to golf stuff, so this will come in very handy. Commented Apr 6, 2018 at 14:44 # 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 Try it online! As far as I can make out, to do this kind of running total in F# it's shorter to create an array of numbers and use the Seq.sumBy method than use a for..to..do block. What this code does it that it creates a collection of floating-point numbers from 1 to s, performs the function fun x->... for the number of elements in the collection, and sums the result. There are s elements in the collection, so the random test is done s times. The actual numbers in the collection are ignored (fun x->, but x is not used). It also means that the application has to first create and fill the array, and then iterate over it. So it's probably twice as slow as a for..to..do loop. And with the array creation memory usage is in the region of O(f**k)! For the actual test itself, instead of using an if then else statement what it does is it calculates the distance (q()+q()|>Math.Sqrt) and rounds it down with Math.Floor. If the distance is within the circle, it'll be rounded down to 0. If the distance is outside the circle it'll be rounded down to 1. The Seq.sumBy method will then total these results. Note then that what Seq.sumBy has totalled is not the points inside the circle, but the points outside it. So for the result it takes s (our sample size) and subtracts the total from it. It also appears that taking a sample size as a parameter is shorter than hard-coding the value. So I am cheating a little bit... # Haskell, 116114110 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

Because dealing with import System.Random; r=randoms(mkStdGen 2) would take too many precious bytes, I generate an infinite list of random numbers with the linear congruential generator that some say is almost cryptographically strong: x↦x*9+1 mod 8^9, which by the Hull-Dobell Theorem has the full period of 8^9.

g yields 4 if the random number point is inside the circle for pairs of random numbers in [0..8^9-1] because this eliminates a multiplication in the formula used.

Usage:

> p 100000
3.14208

Try it online!

# Perl 5, 34 bytes

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

The number of samples is taken from stdin. Requires -p.

Works because:

Try it online!

# MathGolf, 11 bytes

0♫ô4ƒ²)/+♫/

Try it online!

## Explanation

0            Push 0 (used for the sum)
♫           Push 10000
ô          Start block of length 6
4         Push 4
ƒ²       Push random float in range [0,1) and square
)      Increment by 1
/     Divide
+    Add result to the total sum (last instruction of block)
♫/  Divide by 10000

I'll add an explanation tomorrow. It performs 10000 iterations using the same integral as many other languages.

Assuming number of iterations as input:

# MathGolf, 10 bytes

ô4ƒ²)/+\/(

Try it online!

## Explanation

ô           Start block of length 6
4          Push 4
ƒ²        Push random float in range [0,1) and square
)       Increment by 1
/      Divide
+     Add result to the total sum (last instruction of block)
The first iteration, this adds the result to the implicit input, meaning that if the input is 10000 and the result is 3.1, 10003.1 is pushed on the stack. The end result is 10000 too large.
\    Swap the top two elements on the stack (pushes the implicit input to the top)
/   Divide (Since the sum is wonky, this is around 4.14159...)
(  Decrement by one, resulting in 3.14159...)