Mailing-List: contact commits-help@commons.apache.org; run by ezmlm
Precedence: bulk
Reply-To: dev@commons.apache.org
Content-Type: text/plain; charset="utf-8"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Subject: svn commit: r995986 -
 /commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java
Date: Fri, 10 Sep 2010 21:23:26 -0000
To: commits@commons.apache.org
From: luc@apache.org
Message-Id: <20100910212326.DB41A23889ED@eris.apache.org>

Author: luc
Date: Fri Sep 10 21:23:25 2010
New Revision: 995986

URL: http://svn.apache.org/viewvc?rev=995986&view=rev
Log:
Use Jonathan and Peter Borwein quartic formula to compute PI,
it is MUCH faster than the previous one especially for large
numbers of digits and allows quicker loading of the class.
It was tested to compute about 10000 decimal digits, just for the fun

Modified:
    commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java

Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java?rev=995986&r1=995985&r2=995986&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java Fri Sep 10 21:23:25 2010
@@ -592,7 +592,7 @@ public class DfpField implements Field<D
         }
     }
 
-    /** Compute &pi; by atan(1/&radic;(3)) = &pi;/6.
+    /** Compute &pi; using Jonathan and Peter Borwein quartic formula.
      * @param one constant with value 1 at desired precision
      * @param two constant with value 2 at desired precision
      * @param three constant with value 3 at desired precision
@@ -600,30 +600,38 @@ public class DfpField implements Field<D
      */
     private static Dfp computePi(final Dfp one, final Dfp two, final Dfp three) {
 
-        Dfp x = three;
-        x = x.sqrt();
-        x = one.divide(x);
-
-        Dfp denom = one;
-
-        Dfp py = new Dfp(x);
-        Dfp y  = new Dfp(x);
-
-        for (int i = 1; i < 10000; i++) {
-            x = x.divide(three);
-            denom = denom.add(two);
-            if ((i&1) != 0) {
-                y = y.subtract(x.divide(denom));
-            } else {
-                y = y.add(x.divide(denom));
-            }
-            if (y.equals(py)) {
+        Dfp sqrt2   = two.sqrt();
+        Dfp yk      = sqrt2.subtract(one);
+        Dfp four    = two.add(two);
+        Dfp two2kp3 = two;
+        Dfp ak      = two.multiply(three.subtract(two.multiply(sqrt2)));
+
+        // The formula converges quartically. This means the number of correct
+        // digits is multiplied by 4 at each iteration! Five iterations are
+        // sufficient for about 160 digits, eight iterations give about
+        // 10000 digits (this has been checked) and 20 iterations more than
+        // 160 billions of digits (this has NOT been checked).
+        // So the limit here is considered sufficient for most purposes ...
+        for (int i = 1; i < 20; i++) {
+            final Dfp ykM1 = yk;
+
+            final Dfp y2         = yk.multiply(yk);
+            final Dfp oneMinusY4 = one.subtract(y2.multiply(y2));
+            final Dfp s          = oneMinusY4.sqrt().sqrt();
+            yk = one.subtract(s).divide(one.add(s));
+
+            two2kp3 = two2kp3.multiply(four);
+
+            final Dfp p = one.add(yk);
+            final Dfp p2 = p.multiply(p);
+            ak = ak.multiply(p2.multiply(p2)).subtract(two2kp3.multiply(yk).multiply(one.add(yk).add(yk.multiply(yk))));
+
+            if (yk.equals(ykM1)) {
                 break;
             }
-            py = new Dfp(y);
         }
 
-        return y.multiply(new Dfp(one.getField(), 6));
+        return one.divide(ak);
 
     }