User:Hakerh400/How to calculate n-th digit of pi

In this article we explain the simplest and most intuitive (but not the optimal or fastest) way to calculate $n$-th digit of $\pi$ after decimal point for any given non-negative integer $n$. The reason that the author decided to write this article is because the author wants to show that we can compute anything computable as easy as we can describe it. In other words, if we can define what $\pi$ is, we can write a simple program to compute that number $\pi$ (assuming $\pi$ is computable). People usually think that computing digits of $\pi$ with arbitrary precision would require some advanced math formulas and proofs. However, we show in this article that no math formulas are needed, besides the definition of $\pi$ itself. All algorithms are explained in pseudocode.

Definition of $\pi$

Number $\pi$ can be defined in a lot of different ways. In this article we use the definition that $\pi$ is the number $4$ times larger than the number obtained by dividing the area of circle with radius $r$ by the area of square with side length $2r$. Basically:

$\displaystyle \pi =4{\frac {A_{\text{O}}}{A_{\square }}}$

This definition is enough to make a program that can compute any digit of $\pi$.

Approach

First, in order to write such a program, we need to formally specify what we want. Given a non-negative integer $n$, we want $n$-th decimal digit of $\pi$ after decimal point ($1$-indexed). It means that for $n=0$ our function must return $1$, for $n=1$ our function must return $4$, etc. So, we want to calculate the result of dividing the area of circle by the area of square, so that the error is smaller than $10^{-(n+1)}$

There are two problems:

1. How to calculate the area of circle without using $\pi$?
2. How to calculate quotient so that the rounding error is guaranteed to be smaller than some constant?

The first problem can be solved by using the definition of circle. Circle is a set of points in plane and every point of the circle has distance from the center smaller than or equal to the radius of the circle. We can determine whether some point belongs to the circle without using square root - we can simply calculate square of the radius and compare it to the sum of the squares of the coordinates. Point with coordinates $(x,y)$ belongs to the circle iff

$\displaystyle x^{2}+y^{2}\leq r^{2}$

We can approximate one quarter of a circle by allocating a matrix $r\times r$ in memory and for each cell of the matrix we calculate whether the point with coordinates $(x,y)$ belongs to circle with radius $r$. So, we allocate the matrix and fill it with booleans such that the cell has value true iff the pixel belongs to the circle:

matrix[x, y] = x * x + y * y <= r * r

Then we calculate the number of cells in the matrix that have value true, so we obtain the approximate area of circle with radius $r$. Then we can multiply that number by $10^{n}$ and then perform integer division by $r^{2}$. The last digit of that number would be the $n$-th digit of $\pi$.

However, the problem is how to determine how large $r$ we need depending on $n$. We want a function that for given $n$ determines $r$ that will be enough to correctly calculate $n$-th digit of $\pi$. We want the rounding error to be smaller than $10^{-(n+1)}$. The error of approximating a quarter of a circle using a matrix is smaller than $N$, because each boundary cell can have error at most $\textstyle {\frac {1}{2}}$ and there are less than $2N$ such cells. Just to be sure, we will double that number and add $1$:

$\displaystyle r=2\cdot 4\cdot 10^{n+1}+1=8\cdot 10^{n+1}+1$

Optimizations

We do not really need to allocate a matrix. We can use two loops and a single counter. Also, we know that if we walk around the circle quarter, if one coordinate increases, the other must decrease, because the sum of their squares is constant. So, we can keep track of the previous value of the coordinate and we can use bisection to determine the new boundary coordinate.

Pseudocode

Here is the algorithm in pseudocode. All variables are integers with unlimited size and all divisions are integer divisions.

function nthDigitOfPi(n):
  k = 10n
  r = 8 * k + 1
  sum = 0
  x = 0

  for y = r - 1 downto 0:
    d = r2 - y2
    x1 = x
    x2 = r

    while x2 - x1 >= 2:
      x3 = (x1 + x2) / 2
      if x3 * x3 < d:
        x1 = x3
      else:
        x2 = x3

    if x1 = x2 or x2 * x2 < d:
      x = x2
    else:
      x = x1

    sum = sum + x

  return (4 * k * sum / r2) mod 10

Implementation

There is also an implementation in JavaScript (that function is $0$-indexed).

Time complexity

The algorithm has exponential time complexity, because the outer loop iterates $r$ times and $r=8\cdot 10^{n+1}+1$