Finding 3.14 in JavaScript
Tags: javascript, algorithms
Published: 2021-07-18 14:59:57

## Chudnovsky algorithm

To date, one of the fastest and most efficient algorithms for calculating **PI** is considered to be the Chudnovsky algorithm

The principle of this algorithm is the basis of the [2019 record](https://cloud.google.com/blog/products/compute/calculating-31-4-trillion-digits-of-archimedes-constant-on-google-cloud) for the calculation of PI - *31.4 trillion* digits

{% katex %}
 \frac{1}{\pi}=12\sum_{n=0}^{\propto}\frac{(-1)^n(6n)!(A+B_n)}{(n!)^3(3n)!C^{3_n+3/2}}
{% endkatex %}

Skipping all [mathematical transformations](https://bellard.org/pi/pi2700e9/pipcrecord.pdf)

**We get**

{% katex %}
\pi \approx \frac{Q(0,N)}{12T(0,N)+12AQ(0,N)}C^{3/2}
{% endkatex %}

To translate this formula into code, we need to know what Q and T are


**Q and T - Mathematical functions which are expressed as**

{% katex %}
P(n_1, n_2) = P(n_1, m)P(m, n_2)
{% endkatex %}

{% katex %}
Q(n_1, n_2) = Q(n_1, m)P(m, n_2)
{% endkatex %}

{% katex %}
T(n_1, n_2) = T(n_1, m)Q(m, n_2)c^{n_2-m}+P(n_1,m)T(m,n_2)
{% endkatex %}

This looks a little confusing, but let's go over it step by step

**Define the constants**

```js
const A = 13591409
const B = 545140134
const C = 640320
```

**Implementing the algorithm for calculating P, Q and T**

```js
function computePQT(n1, n2) {
  let m = 0
  let PQT = {
    P: 0,
    Q: 0,
    T: 0,
  }

  if (n1 + 1 === n2) {
    PQT.P = n2 * 2 - 1
    PQT.P = PQT.P * (n2 * 6 - 1)
    PQT.P = PQT.P * (n2 * 6 - 5)
    PQT.Q = Math.floor((C * C * C) / 24) * n2 * n2 * n2
    PQT.T = (A + B * n2) * PQT.P
    if (n2 % 2 === 1) {
      PQT.T = -PQT.T
    }
  } else {
    m = Math.floor((n1 + n2) / 2)
    let res1 = computePQT(n1, m)
    let res2 = computePQT(m, n2)
    PQT.P = res1.P * res2.P
    PQT.Q = res1.Q * res2.Q
    PQT.T = res1.T * res2.Q + res1.P * res2.T
  }

  return PQT
}
```

### Find PI

**We need to decide to what decimal place we will count to.** This algorithm at each iteration allows us to find *14.1816474627...* significant digits

*You can try to calculate it yourself*

{% katex %}
\frac{1}{\pi}=\sum_{k=0}^{\propto}c_k
{% endkatex %}

{% katex %}
10^d=lim_{k \mapsto \propto }|c_k/ck+1|
{% endkatex %}

{% katex %}
d=log_{10}151931373056000
{% endkatex %}

After calculating the value, let's **put it in a constant**

```js
const DIGITS_PER_TERM = 14.1816474627
```

**Write a function to calculate PI**

```js
function computePI(digits) {
  if (digits <= 0) {
    return '0'
  }

  const N = Math.floor(digits / DIGITS_PER_TERM) + 1
  const PQT = computePQT(0, N)
  const PI = (PQT.Q / (12 * PQT.T + 12 * A * PQT.Q)) * Math.pow(C, 3 / 2)

  return PI.toFixed(digits)
}
```

Finally, we are ready to **count the decimal places**

```js
const hrstart = process.hrtime()
const PI = computePI(28)
const hrend = process.hrtime(hrstart)

console.log(PI.toString())
console.info(`Execution time (hr): ${hrend[0]}s ${hrend[1] / 1000000}ms`)
```


**Checking the result**

```bash
> node index.js
3.1415926535897935600871733186
Execution time (hr): 0s 0.139102ms
```
### Yeah? Mistake!

We were able to find the number of characters we are interested in, now we can breathe easy and ~~apply the obtained value~~ in practice

But if you look closely, you can find *an error*

**Compare**

    3.1415926535897935600871733186
    3.1415926535897932384626433832

The first value was obtained by us, the second was taken from the [Internet](https://mathshistory.st-andrews.ac.uk/HistTopics/1000_places/)

The divergence starts after 15 characters. **That is how many significant characters the double type has in JavaScript**

### Working on the bugs

To calculate more characters, we need to understand how to work with large numbers in JS

The [BigNumber.js](https://github.com/MikeMcl/bignumber.js/) library for working with large numbers might be suitable for this purpose

But before that we need to simplify the formula a bit by **removing the fractional degree** from it

{% katex %}
\pi=\frac{D\sqrt{E}Q}{AQ+T}
{% endkatex %}

**Rewrite the old constant definitions and add new ones.** At the same time, we **take unnecessary calculations out** of the compute_PQT method

```js
const A = new BigNumber('13591409')
const B = new BigNumber('545140134')
const C = new BigNumber('640320')

const D = new BigNumber('426880')
const E = new BigNumber('10005')

const DIGITS_PER_TERM = new BigNumber('14.1816474627254776555')

const C3_24 = C.multipliedBy(C).multipliedBy(C).dividedToIntegerBy(24)
```

**Rewrite our calculation functions**

```js
function computePI(digits) {
  if (digits <= 0) {
    return '0'
  }

  const DIGITS = new BigNumber(digits)
  const N = DIGITS.dividedToIntegerBy(DIGITS_PER_TERM).plus(1)
  const PREC = DIGITS.multipliedBy(Math.log2(10))

  BigNumber.config({
    DECIMAL_PLACES: Math.ceil(PREC.toNumber()),
    POW_PRECISION: Math.ceil(PREC.toNumber()),
  })

  const PQT = computePQT(new BigNumber(0), N)

  let PI = D.multipliedBy(E.sqrt()).multipliedBy(PQT.Q)
  PI = PI.dividedBy(A.multipliedBy(PQT.Q).plus(PQT.T))

  return PI.toFixed(digits)
}

function computePQT(n1, n2) {
  let m = new BigNumber(0)
  let PQT = {
    P: new BigNumber(0),
    Q: new BigNumber(0),
    T: new BigNumber(0),
  }

  if (n1.plus(1).isEqualTo(n2)) {
    PQT.P = n2.multipliedBy(2).minus(1)
    PQT.P = PQT.P.multipliedBy(n2.multipliedBy(6).minus(1))
    PQT.P = PQT.P.multipliedBy(n2.multipliedBy(6).minus(5))
    PQT.Q = C3_24.multipliedBy(n2).multipliedBy(n2).multipliedBy(n2)
    PQT.T = A.plus(B.multipliedBy(n2)).multipliedBy(PQT.P)
    if (n2.modulo(2).isEqualTo(1)) {
      PQT.T = PQT.T.negated()
    }
  } else {
    m = n1.plus(n2).dividedToIntegerBy(2)

    let res1 = computePQT(n1, m)
    let res2 = computePQT(m, n2)

    PQT.P = res1.P.multipliedBy(res2.P)
    PQT.Q = res1.Q.multipliedBy(res2.Q)
    PQT.T = res1.T.multipliedBy(res2.Q).plus(res1.P.multipliedBy(res2.T))
  }

  return PQT
}
```

**Second try**
```bash
> node index.js
3.1415926535897932384626433833
Execution time (hr): 0s 3.432017ms
```

Note that the running **time of the algorithm is longer**, this is a consequence of storing numbers in strings

**Compare**

    3.1415926535897935600871733186
    3.1415926535897932384626433833
    3.1415926535897932384626433832

### Good!

Only the last digit is different, and that's because we use *toFixed*, which **rounds up the number** when converting it to a string

*Another problem is*

    RangeError: Maximum call stack size exceeded

This error occurs when the node.js runtime has a **call stack overflow**

It can be avoided by giving the runtime the ability to clear the stack

```js
let res1 = await new Promise((resolve) =>
  process.nextTick(async () => resolve(await computePQT(n1, m)))
)
let res2 = await new Promise((resolve) =>
  process.nextTick(async () => resolve(await computePQT(m, n2)))
)
```

**Complete code** can be found on [GitHub](https://github.com/Gerpea/pi)