// WallisPi.java
// Andrew Davison, ad@coe.psu.ac.th, April 2023

/* Calculate the Wallis formula for pi:

    PI     2     2     4     4     6     6     8     8     10     10 
   ---- = --- * --- * --- * --- * --- * --- * --- * --- * ---- * ---- * ....
    2      1     3     3     5     5     7     7     9      9     11 

   https://en.wikipedia.org/wiki/Wallis_product

   Calculate it for several different term sizes, 
   using both doubles and BigDecimals. The error between the results
   and Math.PI are reported, along with the execution time. 
   Also the error when using 355/133 is included for comparison.

   In this case, BigDecimal adds no extra accuracy, but is noticably slower.
   
   Usage:
     > compile WallisPi.java
     > run WallisPi
*/

import java.math.*;
import java.util.function.*;


public class WallisPi
{
  private static BigDecimal two = BigDecimal.valueOf(2);


  public static void main(String[] args)
  {
    System.out.println("Math.PI: " + Math.PI);
    double piApprox = 355.0/113.0;
    double errPercent = Math.abs((piApprox - Math.PI)*100/Math.PI);
    System.out.printf("355/133: %.10f; err:%.10f%%\n", piApprox, errPercent);

    int[] ts = { 50, 100, 200, 2000, 20000, 200000 };
    System.out.println("\nwallisPI()");
    for (int numTerms : ts)
      runAndTime(numTerms, WallisPi::wallisPI);

    System.out.println("\nwallisPI_BD()");
    for (int numTerms : ts)
      runAndTime(numTerms, WallisPi::wallisPI_BD);

  } /* main */


  private static void runAndTime(int numTerms, IntFunction<Double> fn)
  {
    long startTime = System.currentTimeMillis();
    double piVal = fn.apply(numTerms);
    double errPercent = Math.abs((piVal - Math.PI)*100.0/Math.PI);

    System.out.printf("Terms %6d) %.8f (%d ms); err: %.5f%%\n",
                 numTerms, piVal, (System.currentTimeMillis() - startTime), 
                 errPercent);
  }  // end of runAndTime()


  private static double wallisPI(int numTerms)
  {
	if (numTerms <= 0)
      return 0;

	double product = 1;
	double numer = 0;
	double denom;
	for (int i=0; i < numTerms; i++){
	  denom = numer+1;
	  if (i%2 == 0)
        numer += 2;
	  product *= numer/denom;
	}
	return product*2;
  }  // end of wallisPI()


  private static double wallisPI_BD(int numTerms)
  // a BigDecimal version of wallisPI()
  {
	if (numTerms <= 0)
      return 0;

	BigDecimal numer = BigDecimal.ZERO;
	BigDecimal denom = BigDecimal.ONE;

    BigDecimal totNumer = BigDecimal.ONE;
    BigDecimal totDenom = BigDecimal.ONE;

	for (int i=0; i < numTerms; i++){
	  denom = numer.add(BigDecimal.ONE);
	  if (i%2 == 0)
        numer = numer.add(two);
      totNumer = totNumer.multiply(numer);
      totDenom = totDenom.multiply(denom);
	}
	return totNumer.multiply(two).
              divide(totDenom, 20, RoundingMode.HALF_UP).doubleValue();
  }  // end of wallisPI_BD()

}  // end of WallisPi class