Pi with Machin's formula (Java)

From LiteratePrograms

Jump to: navigation, search
Other implementations: Erlang | Haskell | Java | Lisp | Python

Machin's formula

A simple way to compute the mathematical constant π ≈ 3.14159 with any desired precision is due to John Machin. In 1706, he found the formula

which he used along with the Taylor series expansion of the arc cotangent function,

to calculate 100 decimals by hand. The formula is well suited for computer implementation, both to compute π with little coding effort (adding up the series term by term) and using more advanced strategies (such as binary splitting) for better speed.

Contents


In order to obtain n digits, we will use fixed-point arithmetic to compute π × 10n as a Java BigDecimal.

Implementation

<<Pi.java>>=
import java.math.BigDecimal;
import java.math.RoundingMode;
public final class Pi {
  private static final BigDecimal TWO = new BigDecimal("2");
  private static final BigDecimal FOUR = new BigDecimal("4");
  private static final BigDecimal FIVE = new BigDecimal("5");
  private static final BigDecimal TWO_THIRTY_NINE = new BigDecimal("239");
  private Pi() {}
  public static BigDecimal pi(int numDigits) {
    int calcDigits = numDigits + 10;
    return FOUR.multiply((FOUR.multiply(arccot(FIVE, calcDigits)))
      .subtract(arccot(TWO_THIRTY_NINE, calcDigits)))
      .setScale(numDigits, RoundingMode.DOWN);
  }
  private static BigDecimal arccot(BigDecimal x, int numDigits) {
    BigDecimal unity = BigDecimal.ONE.setScale(numDigits,
      RoundingMode.DOWN);
    BigDecimal sum = unity.divide(x, RoundingMode.DOWN);
    BigDecimal xpower = new BigDecimal(sum.toString());
    BigDecimal term = null;
    boolean add = false;
    for (BigDecimal n = new BigDecimal("3"); term == null ||
      !term.equals(BigDecimal.ZERO); n = n.add(TWO)) {
      xpower = xpower.divide(x.pow(2), RoundingMode.DOWN);
      term = xpower.divide(n, RoundingMode.DOWN);
      sum = add ? sum.add(term) : sum.subtract(term);
      add = ! add;
    }
    return sum;
  }
}

Applying Machin's formula

The method pi above uses Machin's formula to compute π using the necessary level of precision, which in turn invokes the private arccot method for arccot computation.

To avoid rounding errors in the result, we set the scale of the BigDecimal to 10 more than the requested number of digits, and later round the return value.

High-precision arccot computation

To calculate arccot of an argument value, we start setting the scale of the number 1 to a value which will provide sufficient precision and avoid rounding errors (number of digits requested + 10) and use this as our unity value. This unity value is divided by x to obtain the first term. Note that all BigDecimal division uses the DOWN rounding mode to emulate integer division. We then repeatedly divide by x2 and a counter value that runs over 3, 5, 7, ..., to obtain each next term. The summation is stopped at the first zero term, which in this fixed-point representation corresponds to a real value less than 10-n.

Example usage: 

System.out.println(Pi.pi(100));

Result:

31415926535897932384626433832795028841971693993751058209749445923078164062862089 986280348253421170679


The program can be used to compute tens of thousands of digits in just a few seconds on a modern computer. (More sophisticated techniques are necessary to calculate millions or more digits in reasonable time, although in principle this program will also work.)

Download code
Retrieved from "http://literateprograms.org/Pi_with_Machin%27s_formula_%28Java%29"
Views