13.3 Generating Billions of Digits

	Java Number Cruncher: The Java Programmer's Guide to Numerical Computing By Ronald Mak
	Table of Contents

	Chapter 13. Computing p

In 1985, two brothers who are mathematicians, Jonathan Borwein and Peter Borwein, published an algorithm that can generate billions of digits of p assuming the computer has sufficient memory! Program 13 C3 implements this algorithm, which is simple to describe.

Set and . Then iterate

graphics/13equ09.gif

and the a _i will converge quartically toward .

Iteration #1 produces 8 correct decimal digits, iteration #2 produces 41 digits, iteration #3 produces 171 digits, iteration #4 produces 694 digits, and each subsequent iteration increases the number of correct digits by more than a factor of 4. Because it needs to test for convergence, the program will do one more iteration than necessary. But it stops partway through the last iteration when it realizes that's the last one.

This program computes the digits of p in several phases, and each phase can consist of several tasks. Table 13-3 shows the phases and tasks .

Listing 13-3a shows class PiBorweinAlgorithm , which implements the algorithm.

Listing 13-3a Class `PiBorweinAlgorithm.`

 package numbercruncher.program13_3; import java.math.BigDecimal; import numbercruncher.mathutils.BigFunctions; /**  * Implement the Borwein algorithm for pi.  */ class PiBorweinAlgorithm implements PiBorweinConstants {     /** number of pi digits */  private int digits;     /** computation scale */    private int scale;     /** iteration counter */    private int iterations;     /** value of pi */          private BigDecimal      pi;     /** parent object */        private PiBorweinParent parent;     /**      * Constructor.      * @param digits the number of digits to compute      * @param parent the parent object      */     PiBorweinAlgorithm(int digits, PiBorweinParent parent)     {         this.digits = digits;         this.parent = parent;         this.scale  = Math.max(((int) (1.005f*digits)), (digits + 5));     }     /**      * Get the number of iterations.      * @return the number of iterations      */     int getIterations() { return iterations; }     /**      * Get the value of pi as a string.      * @return the string      */     String getPi() { return pi.toString(); }     /**      * Compute the digits of pi.      * Notify the parent of each phase and task.      */     void compute()     {         parent.notifyScale(scale);         // Initialization phase.         parent.notifyPhase(INITIALIZING);         parent.notifyTask("constants");         BigDecimal big1   = BigDecimal.valueOf(1);         BigDecimal big2   = BigDecimal.valueOf(2);         BigDecimal big4   = BigDecimal.valueOf(4);         BigDecimal big6   = BigDecimal.valueOf(6);         BigDecimal power2 = big2;         parent.notifyTask("sqrt2");         BigDecimal sqrt2 = BigFunctions.sqrt(big2, scale);         parent.notifyTask("y");         BigDecimal y = sqrt2.subtract(BigDecimal.valueOf(1));         parent.notifyTask("a");         BigDecimal a =             big6.subtract(                 big4.multiply(sqrt2)                     .setScale(scale, BigDecimal.ROUND_HALF_EVEN));         BigDecimal y2, y4, yRoot4, yNumerator, yDenominator;         BigDecimal aTerm4, aTerm;         // Loop once per iteration.         for (;;) {             parent.notifyPhase(++iterations);             parent.notifyTask("y4");             y4 = y.multiply(y)                     .setScale(scale, BigDecimal.ROUND_HALF_EVEN);             y4 = y4.multiply(y4)                     .setScale(scale, BigDecimal.ROUND_HALF_EVEN);             parent.notifyTask("yRoot4");             yRoot4 = BigFunctions.sqrt(big1.subtract(y4), scale);             yRoot4 = BigFunctions.sqrt(yRoot4, scale);             parent.notifyTask("y");             yNumerator   = big1.subtract(yRoot4);             yDenominator = big1.add(yRoot4);             y = yNumerator.divide(yDenominator, scale,                                   BigDecimal.ROUND_HALF_EVEN);             if (y.signum() == 0) break;             parent.notifyTask("aTerm");             aTerm4 = big1.add(y);             aTerm4 = aTerm4.multiply(aTerm4)                         .setScale(scale, BigDecimal.ROUND_HALF_EVEN);             aTerm4 = aTerm4.multiply(aTerm4)                         .setScale(scale, BigDecimal.ROUND_HALF_EVEN);             a = a.multiply(aTerm4)                         .setScale(scale, BigDecimal.ROUND_HALF_EVEN);             parent.notifyTask("power2");             power2 = power2.multiply(big4)                         .setScale(scale, BigDecimal.ROUND_HALF_EVEN);             parent.notifyTask("y2");             y2 = y.multiply(y)                         .setScale(scale, BigDecimal.ROUND_HALF_EVEN);             parent.notifyTask("a");             aTerm = big1.add(y).add(y2);             aTerm = power2.multiply(y).multiply(aTerm)                         .setScale(scale, BigDecimal.ROUND_HALF_EVEN);             a = a.subtract(aTerm)                         .setScale(scale, BigDecimal.ROUND_HALF_EVEN);         }         // Inversion phase.         parent.notifyPhase(INVERTING);         pi = big1.divide(a, digits, BigDecimal.ROUND_DOWN);         parent.notifyPhase(DONE);     } }

Table 13-3. Phases and tasks of the Borwein algorithm for p .

Phase	Task	Computation
Initialization	constants	`BigDecimal` 1, 2, 4, 6, and 5 x 10 ^-( ^scale ⁺¹⁾
	sqrt2
	y
	a
Iteration i	y ⁴
	yRoot4
	y
	aTerm	a _i _- ₁ (1 + y _i ) ⁴
	power2	2 ² ⁱ ⁺ ¹
	y2
	a	a _i = a _i _- ₁ (1 + y _i ) ⁴ - 2 ² ⁱ ⁺ ¹ y _i (1 + y _i + )
Inverting		1/ a _i

Method compute() executes the algorithm. It notifies the parent object of the current phase and task. The parent object must implement interface PiBorweinParent , shown in Listing 13-3b.

To ensure the accuracy of the last few computed digits of p , method compute() does its computations with a scale that is 0.5% greater than the desired number of digits (at least five digits more).

Listing 13-3b Interface `PiBorweinParent` .

 package numbercruncher.program13_3; /**  * Interface for the parent of a PiBorwein object.  */ interface PiBorweinParent {     /**      * Scale notification.      * @param scale the scale being used      */     void notifyScale(int scale);     /**      * Phase notification.      * @param phase the current phase      */     void notifyPhase(int phase);     /**      * Task notification.      * @param task the current computation task      */     void notifyTask(String task); }

Phase notification is done with a set of constants in interface PiBorweinConstants , which is shown in Listing 13-3c.

Listing 13-3c Constants for the Borwein algorithm in interface `PiBorweinConstants` .

 package numbercruncher.program13_3; /**  * Constants for the Borwein pi algorithm.  */ interface PiBorweinConstants {     static final int INITIALIZING = -1;     static final int INVERTING    = -2;     static final int DONE         = -3;     static final int STOPPED      = -4;     static final String STOPPED_EXCEPTION = "STOPPED"; }

Finally, Listing 13-3d shows Program 13 C3 and its output. The noninteractive version of the program prints timestamped lines, each of which contains the current phase and task. As in Program 13 C2, each timestamp includes the elapsed time (hh:mm:ss) of the previous phase.

Listing 13-3d Program 13 C3 and its output.

 package numbercruncher.program13_3; import java.math.BigDecimal; import numbercruncher.mathutils.BigFunctions; import numbercruncher.piutils.PiFormula; /**  * PROGRAM 13-3: The Borwein Pi Algorithm  *  * Compute digits of pi by the Borwein algorithm.  */ public class PiBorwein extends PiFormula                        implements PiBorweinParent, PiBorweinConstants {     /** number of digits to compute */  private int digits;     /** the Borwein algorithm */     private PiBorweinAlgorithm algorithm;     /**      * Compute the digits of pi using the Borwein algorithm.      * @param digits the number of digits of pi to compute      */     private void compute(int digits)     {         this.digits = digits;         startTime = System.currentTimeMillis();         System.out.println(timestamp(startTime) + " START TIME\n");         algorithm = new PiBorweinAlgorithm(digits, this);         algorithm.compute();     }     /**      * Main.      * @param args the array of program arguments      */     public static void main(String args[])     {         PiBorwein pi = new PiBorwein();         try {             int digits = Integer.parseInt(args[0]);             pi.compute(digits);         }         catch(Exception ex) {             System.out.println("ERROR: " + ex.getMessage());         }     }     //-----------------------------------//     // Implementation of PiBorweinParent //     //-----------------------------------//     /**      * Scale notification.      * @param scale the scale being used      */     public void notifyScale(int scale)     {         System.out.println("digits = " + digits);         System.out.println("scale  = " + scale);     }     /**      * Phase notification.      * @param phase the current phase      */     public void notifyPhase(int phase)     {         switch (phase) {             case INITIALIZING: {                 System.out.print("\n" + timestamp(startTime) +                                  " Initialization:");                 break;             }             case INVERTING: {                 System.out.println("\n" + timestamp(markTime) +                                    " Inverting");                 break;             }             case DONE: {                 System.out.println(timestamp(markTime) +                                    " pi computed");                 String totalTime = timestamp(startTime);                 printPi(algorithm.getPi());                 System.out.println("\n" + algorithm.getIterations() +                                    " iterations");                 System.out.println(totalTime + " TOTAL COMPUTE TIME");                 break;             }             default: {                 System.out.print("\n" + timestamp(markTime) +                                  " Iteration " + phase + ":");                 break;             }         }         markTime = System.currentTimeMillis();     }     /**      * Task notification.      * @param task the current computation task      */     public void notifyTask(String task)     {         System.out.print(" " + task);     } }

Output:

 01:29:30.720 (00:00:00) START TIME digits = 700 scale  = 705 01:29:30.720 (00:00:00) Initialization: constants sqrt2 y a 01:29:31.320 (00:00:01) Iteration 1: y4 yRoot4 y aTerm power2 y2 a 01:29:32.370 (00:00:01) Iteration 2: y4 yRoot4 y aTerm power2 y2 a 01:29:33.580 (00:00:01) Iteration 3: y4 yRoot4 y aTerm power2 y2 a 01:29:35.170 (00:00:02) Iteration 4: y4 yRoot4 y aTerm power2 y2 a 01:29:36.870 (00:00:02) Iteration 5: y4 yRoot4 y aTerm power2 y2 a 01:29:38.520 (00:00:02) Iteration 6: y4 yRoot4 y 01:29:39.230 (00:00:01) Inverting 01:29:39.620 (00:00:00) pi computed pi = 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510        58209 74944 59230 78164 06286 20899 86280 34825 34211 70679        82148 08651 32823 06647 09384 46095 50582 23172 53594 08128        48111 74502 84102 70193 85211 05559 64462 29489 54930 38196        44288 10975 66593 34461 28475 64823 37867 83165 27120 19091        45648 56692 34603 48610 45432 66482 13393 60726 02491 41273        72458 70066 06315 58817 48815 20920 96282 92540 91715 36436        78925 90360 01133 05305 48820 46652 13841 46951 94151 16094        33057 27036 57595 91953 09218 61173 81932 61179 31051 18548        07446 23799 62749 56735 18857 52724 89122 79381 83011 94912        98336 73362 44065 66430 86021 39494 63952 24737 19070 21798        60943 70277 05392 17176 29317 67523 84674 81846 76694 05132        00056 81271 45263 56082 77857 71342 75778 96091 73637 17872        14684 40901 22495 34301 46549 58537 10507 92279 68925 89235 6 iterations 01:29:39.620 (00:00:09) TOTAL COMPUTE TIME

Screen 13-3 shows the interactive ^[5] version of Program 13 C3 after it has computed 100,200 digits of p . This applet was run in Microsoft's Internet Explorer version 6.0 (on a 266 MHz Pentium II), which evidently contains a very fast Java virtual machine.

^[5] Each interactive program in this book has two versions, an applet and a standalone application. You can download all the Java source code. See the downloading instructions in the preface of this book.

Screen 13-3. The interactive version of Program 13 C3 after it computed 100,200 digits of p .

graphics/13scr03.jpg

The first 100,200 digits of p were computed in 1961 on an IBM 7090 computer using the arctangent formula

It required 8 hours and 43 minutes.

Top

Listing 13-3a Class PiBorweinAlgorithm.