From 2a0ccce90c26cb20c61bc99e49cc3faf7b26c2f2 Mon Sep 17 00:00:00 2001
From: jos <wjosdejong@gmail.com>
Date: Mon, 21 Jul 2014 21:32:18 +0200
Subject: [PATCH] A faster algorithm for BigNumber pi

---
 lib/constants.js                    | 63 +++++++++++++++-------------
 test/contstants.test.js             |  8 ++--
 test/pi_bailey-borwein-plouffe.html | 54 ++++++++++++++++++++++++
 test/pi_machin.html                 | 64 +++++++++++++++++++++++++++++
 4 files changed, 156 insertions(+), 33 deletions(-)
 create mode 100644 test/pi_bailey-borwein-plouffe.html
 create mode 100644 test/pi_machin.html

diff --git a/lib/constants.js b/lib/constants.js
index fcb6e57e1..834e2b4e2 100644
--- a/lib/constants.js
+++ b/lib/constants.js
@@ -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);
   }
 
   /**
diff --git a/test/contstants.test.js b/test/contstants.test.js
index f7557a746..14ea1035d 100644
--- a/test/contstants.test.js
+++ b/test/contstants.test.js
@@ -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() {
diff --git a/test/pi_bailey-borwein-plouffe.html b/test/pi_bailey-borwein-plouffe.html
new file mode 100644
index 000000000..dff6f7b95
--- /dev/null
+++ b/test/pi_bailey-borwein-plouffe.html
@@ -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>
\ No newline at end of file
diff --git a/test/pi_machin.html b/test/pi_machin.html
new file mode 100644
index 000000000..d28dde4ff
--- /dev/null
+++ b/test/pi_machin.html
@@ -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>
\ No newline at end of file