Let me start out by warning that this post does not deal strictly with
Forte, so you purists can go ahead and click the delete button now.
OK, I got to thinking about PI. I just knew there was an algorithm that
could reach a reasonably close value at something less that several
million iterations. So, by searching through some really old notes --
and Google (not a tool that was available when I was going through
college, sigh) -- I found one.
Below is the best algorithm for calculating PI I have yet come across.
It will generate pi to the limits of double floating point in only five
iterations.
Well, almost. Running the app generates the following output:
Generated PI in 5 iterations.
The generated value of pi is 3.141592653589794
The Math.PI value of pi is 3.141592653589793
Actual PI to 53 places is
3.14159265358979323846264338327950288419716939937510582
As you can see, the answer in the last digit is off by one. But this is
to be expected. We are trying to reach PI by generating a sequence of
values, each value closer than the one before. By the end of the fifth
iteration, we have reached the limit of our precision to adjust the
value. That is, the next adjustment would be generated by using an
amount smaller than the smallest amount we can represent with double.
So no matter what precision is available to us, our calculation can
only get us to within the last digit or two. (Using float arrives at an
"exact" value, because it happens to round to the correct digit.)
To illustrate this, I rewrote the procedure using BigDecimal. I did
this because 1) I had never used BigDecimal before and want to at least
play with it and because 2) being "between jobs," I have gobs of time so
why not.
I will attach the BigDecimal source separately, as it is a bit longer (I
had to include a sqrt function for BigDecimal!). But the output is this:
Generated PI in 6 iterations.
Actual PI to 53 places is
3.14159265358979323846264338327950288419716939937510582
The generated value of pi is
3.14159265358979323846264338327950288419716939937461
The Math.PI value of pi is 3.141592653589793
So this tells me I have found a good algorithm. Going from 15 decimal
places to 50 required only one additional iteration to arrive at the
answer within the limits of that precision. Not bad. This algorithm
can't be used to generate millions of places of PI, but it can probably
whip out a few hundred thousand in no time at all (maybe take all of 8
loops).
If anyone has any critiques of my use of BigDecimal, any suggestions for
improvement, please send them to me. I now see there are more uses for
BigDecimal than I had thought originally. For example, in this instance
of calculating PI, I could have calculated it using BigDecimal at, say
16 places. It would still have taken only 5 iterations and could have
returned the value as a double accurate all the way out. So if you are
going to be doing calculations where the 14th and 15th places of your
double value are going to be significant, do the calculations in
BigDecimal and convert the final answer back to double.
Another benefit. In the original algorithm, the loop is exited when C
(capC in the code) reaches zero. When I wrote it that way in the double
version, the loop never exited. The value would reach
1.232595164407831E-32, something very close to zero, and would stay
there. As double is supposed to be able to handle as small a number as
2.0e-149, the problem is obviously in the calculation of C rather than
the representation of C. So I compare the result of each loop to the
previous result and exit when they are the same. This actually causes
the loop to iterate one extra time. Big deal.
This problem does not occur using BigDecimal. The interim calculations
actually increase the scale as needed. The value C ultimately reaches
zero (or at least zero out to whatever scale is chosen) and the loop
exits. As minor as this problem was in this particular instance, it
could be more significant in other applications (such as balancing very
large amounts of money to the penny).
I have been converted.
Tomm
=================== source ======================
/*
* BrentSalamin.java
*
* Created on September 23, 2002, 1:54 PM
*/
/**
*
* @author Tomm Carr
*/
public class BrentSalamin {
private static final double INITIAL_A = 1.0,
INITIAL_B = 1.0 / Math.sqrt( 2.0 ),
INITIAL_R = 1.0,
INITIAL_S = 0.5;
public double getPi() {
double retVal = 0,
lastPi,
smallA = INITIAL_A,
smallB = INITIAL_B,
smallR = INITIAL_R,
smallS = INITIAL_S,
capA,
capB,
capC,
capR,
capS,
capP;
int iterations = 0;
while (true) {
iterations++;
lastPi = retVal;
capA = (smallA + smallB) / 2;
capB = Math.sqrt( smallA * smallB );
capC = Math.pow( capA - smallB, 2.0 );
capR = 2 * smallR;
capS = smallS - capC * capR;
retVal = (2 * capA * capA) / capS;
if (retVal == lastPi)
// Keep looping until we have calculated pi to as much
precision as we can
break;
smallA = capA;
smallB = capB;
smallR = capR;
smallS = capS;
}
while (lastPi != retVal);
// track loop iterations and print for illustration only
System.out.println( "Generated PI in " + iterations + "
iterations." );
return retVal;
}
/**
* @param args the command line arguments
*/
public static void main(String[] args) {
BrentSalamin piCalc = new BrentSalamin();
double pi = piCalc.getPi();
System.out.println( "The value of pi is " + pi );
}
}
/*
* BrentSalamin.java
*
* Created on September 23, 2002, 1:54 PM
*/
import java.math.BigDecimal;
/**
*
* @author Tomm Carr
*/
public class BrentSalaminBD {
private static final int DB_SCALE = 50,
DB_MODE =
BigDecimal.ROUND_HALF_DOWN;
// BDFUDGE is used by the sqrt function. So, if you move these declarations
around,
// make sure BDFUDGE is declared before BDSQRT2.
private static final BigDecimal BDFUDGE // A value that is 'close enough' to
zero
= new BigDecimal(
"0.00000000000000000000000000000000000000000000000001" );
private static final BigDecimal BD0 // Zero to DB_SCALE decimal
places
= new BigDecimal( 0.0 ).setScale( DB_SCALE );
private static final BigDecimal BD1 // One to DB_SCALE decimal
places
= new BigDecimal( 1.0 ).setScale( DB_SCALE );
private static final BigDecimal BD2 // Two to DB_SCALE decimal
places
= new BigDecimal( 2.0 ).setScale( DB_SCALE );
private static final BigDecimal BDSQRT2 // sqrt(2) to DB_SCALE decimal places
= sqrt( BD2 );
private static final BigDecimal INITIAL_A = BD1;
private static final BigDecimal INITIAL_B = BD1.divide( BDSQRT2, DB_MODE );
private static final BigDecimal INITIAL_R = BD1;
private static final BigDecimal INITIAL_S = BD1.divide( BD2, DB_MODE );
public BigDecimal getPi() {
BigDecimal retVal = new BigDecimal(0.0),
smallA = INITIAL_A,
smallB = INITIAL_B,
smallR = INITIAL_R,
smallS = INITIAL_S,
capA,
capB,
capC,
capR,
capS,
capP;
int iterations = 0;
while (true) {
iterations++;
capA = smallA.add(smallB).divide( BD2, DB_MODE );
capB = sqrt( smallA.multiply( smallB ) );
capC = sqr( capA.subtract( smallB ));
capR = smallR.multiply( BD2 ).setScale( DB_SCALE, DB_MODE );
capS = smallS.subtract( capC.multiply( capR ) ).setScale(
DB_SCALE, DB_MODE );
retVal = BD2.multiply(( capA.multiply( capA )).divide( capS,
DB_MODE )).setScale( DB_SCALE, DB_MODE );
if (capC.compareTo( BD0 ) == 0) {
// Keep looping until we have calculated pi to as much
precision as we can.
// In the case of BigDecimals, precision is limited
only by the scale.
// The algorithm will loop until 'C' (capC) goes to
zero (BD0). Goes to
// zero, that is, beyond the scale we are using.
break;
}
smallA = capA;
smallB = capB;
smallR = capR;
smallS = capS;
}
// track loop iterations and print for illustration only
System.out.println( "Generated PI in " + iterations + " iterations." );
return retVal;
}
/**
* @param args the command line arguments
*/
public static void main(String[] args) {
BrentSalaminBD piCalc = new BrentSalaminBD();
BigDecimal pi = piCalc.getPi();
System.out.println( "Actual PI to 53 places
is\t3.14159265358979323846264338327950288419716939937510582" );
System.out.println( "The generated value of pi is\t" + pi );
System.out.println( "The Math.PI value of pi is\t" + Math.PI );
}
// A few necessary BigDecimalMath routines
// Returns the square of the value
private static BigDecimal sqr( BigDecimal a ) {
return a.multiply( a ).setScale( DB_SCALE, DB_MODE );
}
// Returns the absolute value of the value
private static BigDecimal abs( BigDecimal a ) {
BigDecimal ans = a;
if (a.compareTo( BD0 ) < 0) {
ans = a.negate();
}
return ans;
}
// Returns the square root of the value. Uses the old "guess and adjust"
method.
private static BigDecimal sqrt( BigDecimal w ) {
return sqrt2( w, BD1 );
}
private static BigDecimal sqrt2( BigDecimal w, BigDecimal g ) {
BigDecimal newGuess = g.subtract( guess( w, g ).divide( g.multiply( BD2 ),
DB_MODE ));
if (closeEnough( newGuess, g ))
return newGuess;
else
return sqrt2( w, newGuess ); // Recurse with a better guess
}
private static BigDecimal guess( BigDecimal w, BigDecimal g ) {
return sqr( g ).subtract( w ).setScale( DB_SCALE, DB_MODE );
}
private static boolean closeEnough( BigDecimal a, BigDecimal b ) {
return (BDFUDGE.compareTo( abs( a.subtract( b ))) >= 0);
}
}
To change your JDJList options, please visit: http://www.sys-con.com/java/list.cfm
