Compute the value of π

Python already knows the value of π to 15 significant digits

import sys
import math

print(f"math.pi = {math.pi}")
print(f"sys.float_info.dig = {sys.float_info.dig}")

math.pi = 3.141592653589793
sys.float_info.dig = 15

The Monte Carlo method of computing π

The area of a circle of radius r is

A = πr²

A unit circle is a circle of radius 1. Therefore the area of a unit circle is

A = πr² = π

So we can compute π by computing the area of a unit circle. Let’s compute the area of the unit circle by the Monte Carlo method. About 78% of the random points will fall within or on the circle.

pi.py

The error is
6.665358979329739 × 10^–5
= 6.665358979329739 × .00001
= .00006665358979329739
So the e-05 means “move the decimal point 5 places to the right”. The e stands for “exponent”.

estimated pi = 3.141526
   actual pi = 3.141592653589793
       error = 6.665358979329739e-05

Things to try

Print the lowercase Greek letter π. The Unicode code chart says that the code number of this character is hexadecimal 03C0. In binary, this number would be
0000001111000000
To make it easier to read, divide it into four nibbles:
0000 0011 1100 0000
For a \u in a quoted string, see Escape Sequences.
```
print(f"estimated \u03C0 = {estimatedPi}")
print(f"   actual \u03C0 = {math.pi}")
```
```
 estimated π = 3.1413028
    actual π = 3.141592653589793
```

I get a different value for π each time I run the script, because the random numbers are different each time.

estimated π = 3.1414048
   actual π = 3.141592653589793
      error = 0.00018785358979300781

estimated π = 3.1415796
   actual π = 3.141592653589793
      error = 1.305358979308835e-05

To get the same random numbers each time,

#immediately after the imports
random.seed(0)

 estimated π = 3.1413028
    actual π = 3.141592653589793
       error = 0.0002898535897930543

 estimated π = 3.1413028
    actual π = 3.141592653589793
       error = 0.0002898535897930543

You would want the same random numbers each time if you were debugging the script, and you wanted it to behave the same way each time you ran it.

How long did it take to loop ten million times on my MacBook Pro with a 2.8 gigahertz Intel Core i7 processor? To time a piece of Python code, you first have to package it as a user-defined function. You can also compute the average by using statistics.mean.

"""
Time how long a user-defined function takes to execute.
"""

import sys
import timeit

def f():      #Define (i.e., create) a function named f containing 2 statements.
    sum = 10 + 20
    print(f"sum = {sum}")

s = timeit.timeit(f, number = 2)             #This s is a float.
print(f"{s} seconds")
print()                                      #Skip a line.

s = timeit.repeat(f, repeat = 3, number = 2) #This s is a list of 3 floats.
print(f"{s} seconds")
print(f"Average is {sum(s) / len(s)} seconds.")
sys.exit(0)

sum = 30
sum = 30
0.10639075400000003 seconds

sum = 30
sum = 30
sum = 30
sum = 30
sum = 30
sum = 30
[0.099797052, 0.10035522399999997, 0.11670459899999996] seconds
Average is 0.1056189583333333 seconds.

How long did it take to loop ten million times on my MacBook Pro with a 2.8 gigahertz Intel Core i7 processor?

"""
pi.py

Compute the value of pi using the Monte Carlo method.
Generate n random points, each lying in the 1 by 1 square whose lower left
corner is at the origin.  Count how many of these points fall in the upper right
quadrant of the unit circle centered at the origin.  The fraction count/n is the
fraction of the square that is occupied by the quadrant.
"""

import sys
import random
import math
import timeit

n = 10_000_000 #how many random points.  Underscores ignored.
count = 0      #how many random points fell within the circle

def f():
    global count                    #Don't need to say global n
    for i in range(n):
        x = random.random()         #a random number in the range 0 <= x < 1
        y = random.random()         #another random number in the same range
        distance = math.hypot(x, y) #distance from (0, 0) to (x, y)
        if distance <= 1:           #if the point (x, y) is in or on the circle,
            count += 1

repeat = 3
number = 2

s = timeit.repeat(f, repeat = repeat, number = number)
print(f"{s} seconds")
print(f"Average is {sum(s) / len(s)} seconds.")
print()

estimatedPi = 4 * count / (n * repeat * number)
print(f"estimated \u03C0 = {estimatedPi}")
print(f"   actual \u03C0 = {math.pi}")
print(f"      error = {abs(estimatedPi - math.pi)}")
sys.exit(0)

[7.88843082, 7.6932537, 7.683352416] seconds
Average is 7.755012312000001 seconds.

estimated π = 3.1418729333333335
   actual π = 3.141592653589793
      error = 0.0002802797435403903

We currently have

    if math.hypot(x, y) <= 1:   #if the point (x, y) is in or on the circle,
        count += 1

which means

    if math.sqrt(x**2 + y**2) <= 1:
        count += 1

We can eliminate the square root by squaring both sides of the inequality.

    if x**2 + y**2 <= 1:
        count += 1

Does this let the script run faster?

[8.397018482999998, 8.219414262, 8.227936939] seconds
Average is 8.281456561333332 seconds.

estimated π = 3.1415746
   actual π = 3.141592653589793
      error = 1.805358979289906e-05

Will it run faster if we change the exponentiation to plain old multiplication?

    if x*x + y*y <= 1:
        count += 1

[6.834141311, 6.711918756999999, 6.686868610000001] seconds
Average is 6.744309559333334 seconds.

estimated π = 3.141641133333333
   actual π = 3.141592653589793
      error = 4.847974354005302e-05

If we increase the number of iterations from 10 million to 100 million, will the error get smaller? The correct digits are in yellow. We got a lot more of them this time.

n = 100_000_000 #how many random points.  Underscores ignored.

[69.296841571, 67.593704348, 67.090559848] seconds
Average is 67.99370192233333 seconds.

estimated π = 3.1415928866666665
   actual π = 3.141592653589793
      error = 2.330768733571631e-07

[For space cadets only.] Compute the area of one “hump” of a sine curve using the Monte Carlo method.