archive-gh-me/mathjs
1
0
Fork 0
mirror of https://github.com/josdejong/mathjs.git synced 2026-01-25 15:07:57 +00:00
Code Issues Projects Releases Wiki Activity

A faster algorithm for BigNumber pi

Browse Source
This commit is contained in:
jos 2014-07-21 21:32:18 +02:00
parent 8c224d3e71
commit 2a0ccce90c
4 changed files with 156 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
lib/constants.js
View File

@ -22,39 +22,44 @@ module.exports = function (math, config) {
}
/**
* Calculate BigNumber pi
* arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + x^9/9 - ...
* = x - x^2*x^1/3 + x^2*x^3/5 - x^2*x^5/7 + x^2*x^7/9 - ...
* @param {BigNumber} x
* @returns {BigNumber} arc tangent of x
*/
function arctan(x) {
var y = x;
var yPrev = NaN;
var x2 = x.times(x);
var num = x;
var sign = -1;
for (var k = 3; !y.equals(yPrev); k += 2) {
num = num.times(x2);
yPrev = y;
y = (sign > 0) ? y.plus(num.div(k)) : y.minus(num.div(k));
sign = -sign;
}
return y;
}
/**
* Calculate BigNumber pi.
*
* Uses Machin's formula: pi / 4 = 4 * arctan(1 / 5) - arctan(1 / 239)
* http://milan.milanovic.org/math/english/pi/machin.html
* @returns {BigNumber} Returns pi
*/
function bigPi() {
// the Bailey-Borwein-Plouffe formula
// http://stackoverflow.com/questions/4484489/using-basic-arithmetics-for-calculating-pi-with-arbitary-precision
var p16 = new BigNumber(1);
var k8 = new BigNumber(0);
var pi = new BigNumber(0);
// we calculate pi with a few decimal places extra to prevent round off issues
var Big = BigNumber.constructor({precision: BigNumber.config().precision + 4});
var pi4th = new Big(4).times(arctan(new Big(1).div(5)))
.minus(arctan(new Big(1).div(239)));
var one = new BigNumber(1);
var two = new BigNumber(2);
var four = new BigNumber(4);
for(var k = new BigNumber(0); k.lte(config.precision); k = k.plus(1)) {
// pi += 1/p16 * (4/(8*k + 1) - 2/(8*k + 4) - 1/(8*k + 5) - 1/(8*k+6));
// p16 *= 16;
//
// a little simplified (faster):
// pi += p16 * (4/(8*k + 1) - 2/(8*k + 4) - 1/(8*k + 5) - 1/(8*k+6));
// p16 /= 16;
var f = four.div(k8.plus(1))
.minus(two.div(k8.plus(4)))
.minus(one.div(k8.plus(5)))
.minus(one.div(k8.plus(6)));
pi = pi.plus(p16.times(f));
p16 = p16.div(16);
k8 = k8.plus(8);
}
return pi;
// the final pi has the requested number of decimals
return new BigNumber(4).times(pi4th);
}
/**

8
test/contstants.test.js
View File

@ -7,7 +7,7 @@ describe('constants', function() {
describe('number', function () {
it('should have pi', function() {
approx.equal(math.pi, 3.14159265358979);
approx.equal(math.pi, 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664);
approx.equal(math.sin(math.pi / 2), 1);
approx.equal(math.PI, math.pi);
});
@ -60,7 +60,7 @@ describe('constants', function() {
var bigmath = math({number: 'bignumber', precision: 64});
it('should have bignumber pi', function() {
assert.equal(bigmath.pi.toString(), '3.141592653589793238462643383279502884197169399375105820974944591');
assert.equal(bigmath.pi.toString(), '3.141592653589793238462643383279502884197169399375105820974944592');
});
it('should have bignumber tau', function() {
@ -91,8 +91,8 @@ describe('constants', function() {
assert.equal(bigmath.LOG10E.toString(), '0.4342944819032518276511289189166050822943970058036665661144537832');
});
it('should have bignumber PI', function() {
assert.equal(bigmath.PI.toString(), '3.141592653589793238462643383279502884197169399375105820974944591');
it('should have bignumber PI (upper case)', function() {
assert.equal(bigmath.PI.toString(), '3.141592653589793238462643383279502884197169399375105820974944592');
});
it('should have bignumber SQRT1_2', function() {

54
test/pi_bailey-borwein-plouffe.html Normal file
View File

@ -0,0 +1,54 @@
<html>
<head>
<title>Pi Bailey-Borwein-Plouffe</title>
<script src="../node_modules/decimal.js/decimal.js"></script>
</head>
<body>
<script>
Decimal.config({precision: 100});
function pi() {
// the Bailey-Borwein-Plouffe formula
// http://stackoverflow.com/questions/4484489/using-basic-arithmetics-for-calculating-pi-with-arbitary-precision
var zero = new Decimal(0);
var one = new Decimal(1);
var two = new Decimal(2);
var four = new Decimal(4);
var p16 = one;
var pi = zero;
var precision = Decimal.config().precision;
var k8 = zero;
for (var k = zero; k.lte(precision); k = k.plus(one)) {
// pi += 1/p16 * (4/(8*k + 1) - 2/(8*k + 4) - 1/(8*k + 5) - 1/(8*k+6));
// p16 *= 16;
//
// a little simpler:
// pi += p16 * (4/(8*k + 1) - 2/(8*k + 4) - 1/(8*k + 5) - 1/(8*k+6));
// p16 /= 16;
var f = four.div(k8.plus(1))
.minus(two.div(k8.plus(4)))
.minus(one.div(k8.plus(5)))
.minus(one.div(k8.plus(6)));
pi = pi.plus(p16.times(f));
p16 = p16.div(16);
k8 = k8.plus(8);
}
return pi;
}
console.time('estimation');
var calculatedPi = pi();
console.timeEnd('estimation');
document.write('<code>Real: 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664...</code><br>');
document.write('<code>Est:&nbsp;&nbsp;' + calculatedPi.toString() + '</code>');
</script>
</body>
</html>

64
test/pi_machin.html Normal file
View File

@ -0,0 +1,64 @@
<html>
<head>
<title>Pi Machin</title>
<script src="../node_modules/decimal.js/decimal.js"></script>
</head>
<body>
<script language="javascript">
Decimal.config({precision: 100});
// arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + x^9/9 - ...
// = x - x^2*x^1/3 + x^2*x^3/5 - x^2*x^5/7 + x^2*x^7/9 - ...
function arctan(x) {
var y = x;
var yPrev = NaN;
var x2 = x.times(x);
var num = x;
var sign = -1;
for (var k = 3; !y.equals(yPrev); k += 2) {
num = num.times(x2);
yPrev = y;
y = (sign > 0) ? y.plus(num.div(k)) : y.minus(num.div(k));
sign = -sign;
}
return y;
}
function pi() {
// Machin: Pi / 4 = 4 * arctan(1 / 5) - arctan(1/239)
// http://milan.milanovic.org/math/english/pi/machin.html
// we calculate pi with one decimal place extra to prevent round off issues
var DecimalPlus = Decimal.constructor({precision: Decimal.config().precision + 1});
var pi4th = new DecimalPlus(4).times(arctan(new DecimalPlus(1).div(5)))
.minus(arctan(new DecimalPlus(1).div(239)));
// the final pi has the requested number of decimals
return new Decimal(4).times(pi4th);
}
console.time('calculation');
var calculatedPi = pi();
console.timeEnd('calculation');
var mathematicaPi = '3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664...';
document.write('<code>Calculated:&nbsp;&nbsp;' + calculatedPi + '</code><br>');
document.write('<code>Mathematica:&nbsp;' + mathematicaPi + '</code>');
/*
var mathematicaPi = '3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664...';
document.write('<code>Calculated:&nbsp;&nbsp;' + pi + '</code><br>');
document.write('<code>Mathematica:&nbsp;' + mathematicaPi + '</code>');
*/
</script>
</body>
</html>