[Question]
  [
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.
  
[Answer]
# 80386 machine code, 40 38 bytes
Hexdump of the code:
```
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
```
How to get this code (from assembly language):
```
    // 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;
```
This is a function using MS `fastcall` calling convention (number of iterations is passed in register `ecx`). It returns the result in the `st` register.
Fun things about this code:
* `rdrand` - just 3 bytes to generate a random number!
* It uses (unsigned) integer arithmetic until the final division.
* The comparison of squared distance (`D`) with squared radius (`2^32`) is performed automatically - the carry flag contains the result.
* To multiply the count by 4, it counts the samples in steps of 4.
[Answer]
# Matlab/Octave, 27 bytes
I know there already is a Matlab/Octave answer, but I tried my own approach. I used the fact that the integral of `4/(1+x^2)` between 0 and 1 is pi.
```
mean(4./(1+rand(1,1e5).^2))
```
[Answer]
# CJam, ~~27 23~~ 22 or 20 bytes
```
4rd__{{1dmr}2*mhi-}*//
```
*2 bytes saved thanks to Runner112, 1 byte saved thanks to Sp3000*
It takes the iteration count from STDIN as an input.
This is as straight forward as it gets. These are the major steps involved:
* Read the input and run the Monte Carlo iterations that many times
* In each iteration, get sum of square of two random floats from 0 to 1 and see if it is less than 1
* Get the ratio of how many times we got less than 1 by total iterations and multiply it by 4 to get PI
**Code expansion**:
```
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";
```
[Try it online here](http://cjam.aditsu.net/#code=4rd__%7B%7B1dmr%7D2*mhi-%7D*%2F%2F&input=100000)
---
If average value of `1/(1+x^2)` is also considered as Monte Carlo, then this can be done in 20 bytes:
```
Urd:K{4Xdmr_*)/+}*K/
```
[Try it here](http://cjam.aditsu.net/#code=Urd%3AK%7B4Xdmr_*)%2F%2B%7D*K%2F&input=1000000)
[Answer]
# R, 40 (or 28 or 24 using other methods)
```
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
Another Python one, if numpy is allowed, but pretty similar to Matlab/Octave:
```
import numpy;sum(sum(numpy.random.rand(2,8e5)**2)<1)/2e5
```
[Answer]
# Commodore 64 Basic, 45 bytes
```
1F┌I=1TO1E3:C=C-(R/(1)↑2+R/(1)↑2<1):N─:?C/250
```
PETSCII substitutions: `─` = `SHIFT+E`, `/` = `SHIFT+N`, `┌` = `SHIFT+O`
Generates 1000 points in the first quadrant; for each, adds the truthness of "x^2+y^2<1" to a running count, then divides the count by 250 to get `pi`. (The presence of a minus sign is because on the C64, "true" = -1.)
[Answer]
# Mathematica, ~~42~~ ~~40~~ 39 bytes (or 31/29?)
I've got three solutions all at 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}]/#&
```
They are all unnamed functions that take the number of samples `n` andd return a rational approximating π. First they all generate `n` points in the unit square in the positive quadrant. Then they determine the number of those samples that lie within the unit circle, and then they divide by the number of samples and multiply by `4`. The only difference is in how they determine the number of sampples inside the unit circle:
* The first one uses `Count` with the condition that `Norm[p] < 1`.
* The second one subtracts the norm of each point from `1` and then rounds up. This turns numbers inside the unit circle to `1` and those outside to `0`. Afterwards I just sum them all up with `Tr`.
* The third one does essentially the same, but subtracts the from `1.5`, so I can use `Round` instead of `Ceiling`.
Aaaaaand *while* writing this up, it occurred to me that there is indeed a shorter solution, if I just subtract from `2` and then use `Floor`:
```
4Tr@Floor[2-Norm/@1~RandomReal~{#,2}]/#&
```
or saving another byte by using the Unicode flooring or ceiling operators:
```
4Tr@⌊2-Norm/@1~RandomReal~{#,2}⌋/#&
4Tr@⌈1-Norm/@1~RandomReal~{#,2}⌉/#&
```
Note that the three rounding-based solutions can also be written with `Mean` instead of `Tr` and without the `/#`, again for the same bytes.
---
If other Monte Carlo based approaches are fine (specifically, the one Peter has chosen), I can do 31 bytes by estimating the integral of `√(1-x2)` or 29 using the integral of `1/(1+x2)`, this time given as a floating point number:
```
4Mean@Sqrt[1-1~RandomReal~#^2]&
Mean[4/(1+1~RandomReal~#^2)]&
```
[Answer]
# Python 2, ~~77~~ 75 bytes
```
from random import*;r=random;a=0;exec"a+=r()**2+r()**2<1;"*4000;print a/1e3
```
Uses 4000 samples to save bytes with `1e3`.
[Answer]
# J, 17 bytes
Computes the mean value of `40000` sample values of the function `4*sqrt(1-sqr(x))` in the range `[0,1]`.
Handily `0 o.x` returns `sqrt(1-sqr(x))`.
```
   1e4%~+/0 o.?4e4$0
3.14915
```
[Answer]
## [><> (Fish)](http://esolangs.org/wiki/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
```
[Answer]
**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
```
[Answer]
# 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.
[Answer]
## **Javascript: 62 Bytes**
```
for(r=Math.random,t=c=8e4;t--;c-=r()**2+r()**2|0);alert(c/2e4)
```
I used the previous (now deleted) javascript answer and shaved 5 bytes.
[Answer]
## GolfScript (34 chars)
```
0{^3?^rand.*^.*+/+}2000:^*`1/('.'@
```
[Online demo](http://golfscript.apphb.com/?c=MHteMz9ecmFuZC4qXi4qKy8rfTIwMDA6XipgMS8oJy4nQA%3D%3D)
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](https://codegolf.stackexchange.com/users/20260/xnor) for the particular Monte Carlo method employed.
[Answer]
# 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
```
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.
[Answer]
# 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.
[Answer]
# [Perl 6](http://perl6.org/), 33 bytes
```
{4*sum((1>rand²+rand²)xx$_)/$_}
```
[Try it online!](https://tio.run/##K0gtyjH7n1upoJamYPu/2kSruDRXQ8PQrigxL@XQJm0IpVlRoRKvqa8SX/tfrzixUiEtv0ghTcPQwMBAU6GiQsHQ4D8A "Perl 6 – Try It Online")
This is a function that takes the number of samples as an argument.
[Answer]
# Scala, 87 77 66 bytes
```
def s=math.pow(math.random,2);Seq.fill(1000)(s+s).count(_<1)/250d
```
[Answer]
# 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
$ 
```
[Answer]
# [Joe](https://github.com/JaniM/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:
```
   :%$,(4*/+1>/+*,?2~;
   F400000
3.14154
   F400000
3.14302
```
[Answer]
# dc, 59 characters (whitespace is ignored)
```
[? 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
```
I tested this on Plan 9 and OpenBSD, so I imagine it will work on Linux (GNU?) `dc`.
## Explanation by line:
1. Stores code to [read and square two floats; execute register `i` if 1 is greater than the sum of their squares] in register `u`.
2. Stores code to [increment register `x` by 1] in register `i`.
3. Stores code to [execute register `u`, increment register `m`, and then execute register `z` if register `m` is greater than register `n`] in register `z`.
4. 
5. Set the scale to 5 decimal points.
6. Read the number of points to sample from the first line of input.
7. Execute register `z`.
8. Divide register `x` (the number of hits) by register `n` (the number of points), multiply the result by 4, and print.
9. Quit.
## However, I cheated:
The program needs a supply of random floats between 0 and 1. 
```
/* frand.c */
#include <u.h>
#include <libc.h>
void
main(void)
{
    srand(time(0));
    for(;;)
        print("%f\n", frand());
}
```
## Usage:
```
#!/bin/rc
# runpi <number of samples>
{ echo $1; frand } | dc pi.dc
```
## Test run:
```
% 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
```
## Test run:
```
$ echo 10000 | dc pigolf.dc
3.13640
```
[Answer]
# Pyth, 19
```
c*4sm<sm^OQ2 2*QQQQ
```
Give the desired number of iterations as input.
[Demonstration](https://pyth.herokuapp.com/?code=c*4sm%3Csm%5EOQ2+2*QQQQ&input=10000)
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.
[Answer]
# 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)))
```
[Answer]
# 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);
}
```
[Answer]
## 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
```
[Answer]
# [Fortran (GFortran)](https://gcc.gnu.org/fortran/), ~~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!](https://tio.run/##S8svKilKzNNNT4Mw/v93dvTxUQgOcvRz0TDQ5HLxV/C0NdQxdDXmirAFC2pyRcIYnm4aEVoR2pFakTaGmo62jtqGXK5@Li7@XAFBnn4hWjqO@kamBiCh//8B "Fortran (GFortran) – 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)](https://gcc.gnu.org/fortran/), ~~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!](https://tio.run/##S8svKilKzNNNT4Mw/v93tDXgcvFX8LQ11DF0NeZydvTxUQhy9HPx9433C/V1cg3SiNDEJhqpyeXpphGhFaEdqRVpY6jpaOuobcjl6ufi4s8VEOTpF6Kl46hvZGoAEvr/HwA "Fortran (GFortran) – Try It Online")
(One byte saved in each version; thanks Penguino).
[Answer]
# [Japt](https://github.com/ETHproductions/japt), 26 or 18 bytes
```
o r_+ÂMhMr p +Mr p <1Ã*4/U
```
[Try it online!](https://tio.run/##y0osKPn/P1@hKF77cJNvhm@RQoGCNpi0MTzcrGWiH/r/v6EBEAAA)
Analogous to [Optimizer's answer](https://codegolf.stackexchange.com/a/47767/16484), 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
```
[Answer]
# 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!](https://tio.run/##VU49a8MwEN39K45AQUqxokCnNvJQ2m5JTD2GDHJzag2WZEvnxIb8d8d1s/Tg4PE@7p2J6ZcPOI6@QQfFEAntS1IjQVAOL/Cp3clbxmeqZVwl8ItKFcQOe3rzXVnjJEO5LGePZfHZ1F4TVyymBbYidvZ1YKZz0KfZdOJx2mu21fQjijbQHX7U3gd@WAspRDzy5ZOQqzgeNu@OwpD7ylF2/CvQlQMdvs9w/wUjVVYT7k1egQILazmPkAk0YQoaB4sHs/hnTECONw)
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...
[Answer]
# 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
```
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!](https://tio.run/##FclLDoIwFEDROat4gw7aggqJMZJQNlJr8ppWJNBPSk0YuPeKd3juG7fFrmspRtyffTVJbLQS28fJ6xc51prrwXCjqghe0MPPGRdLX2sICTwj00DGOduE2VL62E@jC4b2fK87ZtjxmGybm2IXXxzOHgTENPsMBCJ07b/yAw "Haskell – Try It Online")
[Answer]
# Perl 5, 34 bytes
```
$_=$a/map$a+=4/(1+(rand)**2),1..$_
```
The number of samples is taken from stdin. Requires `-p`.
Works because:
[![](https://chart.googleapis.com/chart?cht=tx&chf=bg,s,FFFFFF00&chl=%5Cdisplaystyle%5Cint_0%5E1%5C+%7B1%5Cover+1%2Bx%5E2%7D%5C+%5Cmathrm%7Bd%7Dx%3D%7B%5Cpi%5Cover+4%7D)](https://www.wolframalpha.com/input/?i=integral+%7B+1%2F(1%2Bx%5E2)+%7D+for+x+%3D+0..1)
[Try it online!](https://tio.run/##K0gtyjH9/18l3lYlUT83sUAlUdvWRF/DUFujKDEvRVNLy0hTx1BPTyX@/39DAxD4l19QkpmfV/xftwAA "Perl 5 – Try It Online")
[Answer]
# [MathGolf](https://github.com/maxbergmark/mathgolf), 11 bytes
```
0♫ô4ƒ²)/+♫/
```
[Try it online!](https://tio.run/##y00syUjPz0n7/9/g0czVh7eYHJt0aJOmvjaQo///PwA "MathGolf – 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](https://github.com/maxbergmark/mathgolf), 10 bytes
```
ô4ƒ²)/+\/(
```
[Try it online!](https://tio.run/##y00syUjPz0n7///wFpNjkw5t0tTXjtHX@P/f0MDAAAA "MathGolf – 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...)
```
]