Finding 3.14 in JavaScript

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 for the calculation of PI - 31.4 trillion digits

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

Skipping all mathematical transformations

We get

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

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

Q and T - Mathematical functions which are expressed as

P(n_1, n_2) = P(n_1, m)P(m, n_2)

Q(n_1, n_2) = Q(n_1, m)P(m, n_2)

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)

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

Define the constants

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

Implementing the algorithm for calculating P, Q and T

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

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

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

d=log_{10}151931373056000

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

const DIGITS_PER_TERM = 14.1816474627

Write a function to calculate PI

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

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

> 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

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 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

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

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

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

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

> 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

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)))
)