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
Skipping all mathematical transformations
We get
To translate this formula into code, we need to know what Q and T are
Q and T - Mathematical functions which are expressed as
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
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
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)))
)
Complete code can be found on GitHub
Top comments (0)