# On the computation of π

## 1 Asking the math library

My computer tells me that π is *approximatively*

from math import *pi

3.141592653589793

## 2 Buffon's needle

By the method of Buffon's needle we
get the **approximation**

import numpy as np N = 10000 np.random.seed(seed=42) x = np.random.uniform(size=N, low=0, high=1) theta = np.random.uniform(size=N, low=0, high=pi/2) P = sum(x + np.sin(theta) > 1) / N 2/P

3.128911138923655

## 3 Using Monte Carlo experiments

A method that is easier to understand and does not make use of the \(\sin\) function is based on the fact that if \(X \sim U(0,1)\) and \(Y \sim U(0,1)\), then \(P[X^2+Y^2 \leq 1] = \pi/4\) (see "Monte Carlo method" on Wikipedia). The following code uses this approach:

import matplotlib.pyplot as plt np.random.seed(seed=42) x = np.random.uniform(size=N, low=0, high=1) y = np.random.uniform(size=N, low=0, high=1) inside = (x**2 + y**2) <= 1 outside = np.logical_not(inside) fig, ax = plt.subplots(1) ax.scatter(x[inside], y[inside], c='b', alpha=0.2, edgecolor=None) ax.scatter(x[outside], y[outside], c='r', alpha=0.2, edgecolor=None) ax.set_aspect('equal') plt.savefig(matplot_lib_filename) print(matplot_lib_filename)

It is then straightforward to obtain an approximation to π by counting how many times, on average, \(X^2 + Y^2 \le 1\):

4 * np.mean(inside)

3.1556

## 4 Through Wallis' product

John Wallis' product for π states that the infinite product below converges to \(\pi / 2\). Check out the geometrical proof.

prod = 1 for i in range(2, N): prod *= i / (i-1) if i % 2 == 0 else (i-1) / i prod * 2

3.141435562175509

## 5 As the root of a function

As we know, \(sin(\pi) = 0\), thus value for π can be found by approximating the root (one of them, that is) of \(sin(x)\). This can be done using Newton's method, for instace.

from scipy.optimize import fsolve results = fsolve(sin, x0=3) results[0]

3.141592653589793

## Comments