DEV Community

German Cyganov
German Cyganov

Posted on • Originally published at Medium

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

1π=12n=0(1)n(6n)!(A+Bn)(n!)3(3n)!C3n+3/2 \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

πQ(0,N)12T(0,N)+12AQ(0,N)C3/2 \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(n1,n2)=P(n1,m)P(m,n2) P(n_1, n_2) = P(n_1, m)P(m, n_2)
Q(n1,n2)=Q(n1,m)P(m,n2) Q(n_1, n_2) = Q(n_1, m)P(m, n_2)
T(n1,n2)=T(n1,m)Q(m,n2)cn2m+P(n1,m)T(m,n2) 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
Enter fullscreen mode Exit fullscreen mode

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
}
Enter fullscreen mode Exit fullscreen mode

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

1π=k=0ck \frac{1}{\pi}=\sum_{k=0}^{\propto}c_k
10d=limkck/ck+1 10^d=lim_{k \mapsto \propto }|c_k/ck+1|
d=log10151931373056000 d=log_{10}151931373056000

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

const DIGITS_PER_TERM = 14.1816474627
Enter fullscreen mode Exit fullscreen mode

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)
}
Enter fullscreen mode Exit fullscreen mode

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`)
Enter fullscreen mode Exit fullscreen mode

Checking the result

> node index.js
3.1415926535897935600871733186
Execution time (hr): 0s 0.139102ms
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode

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

π=DEQAQ+T \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)
Enter fullscreen mode Exit fullscreen mode

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
}
Enter fullscreen mode Exit fullscreen mode

Second try

> node index.js
3.1415926535897932384626433833
Execution time (hr): 0s 3.432017ms
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode

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)))
)
Enter fullscreen mode Exit fullscreen mode

Complete code can be found on GitHub

Discussion (0)