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

/* Calculate pi using three different Machin-like formulae,
   and compare them with the Nut.PI 100 dp constant.

   Pi Machin formulae:
   https://en.wikipedia.org/wiki/Machin-like_formula

   The Machin formulae require atan() for BigDecimal, 
   which comes from Big-Math.

   Then generate PI using the JML Pi class and load the
   10,000 digit pi stored in pi10000.txt, and count
   their digits. This is to observe if pi is a 
   Normal number (i.e. its digits are distributed
   uniformly).

   The digit spread is also displayed as a vectorgram,
   inspired by John Venn's "The Logic of Chance", 
   1888, 3ed ed.

   https://en.wikipedia.org/wiki/Normal_number
   https://mathworld.wolfram.com/NormalNumber.html


   Usage:
     > compile CalcPi.java
     > run CalcPi
*/


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

import de.tilman_neumann.jml.transcendental.Pi;
import de.tilman_neumann.jml.precision.Precision;
import de.tilman_neumann.jml.precision.Scale;

import ch.obermuhlner.math.big.BigDecimalMath;


public class CalcPi
{

  public static void main(String[] args)
  {
    // built-in PI values
    System.out.println(" Math.PI = " + Math.PI);
    System.out.println("  Nut.PI = " + Nut.PI);   // 100 dp

    MathContext mc = new MathContext(100);

    BigDecimal bd4 = BigDecimal.valueOf(4);

    // Machin 1). pi = 4 * arctan(1)
    BigDecimal piSimple = bd4.multiply(
         BigDecimalMath.atan(BigDecimal.ONE,mc));
    System.out.println("\nSimple pi= " + piSimple);


    BigDecimal bd12 = BigDecimal.valueOf(0.5);
    BigDecimal bd15 = BigDecimal.valueOf(0.2);
    BigDecimal bd18 = BigDecimal.valueOf(1.0/8);

    // Machin 2). pi = 4(arctan(1/2) + arctan(1/5) + arctan(1/8))
    BigDecimal piDase =
       bd4.multiply(
         BigDecimalMath.atan(bd12,mc).
         add(BigDecimalMath.atan(bd15,mc)).
         add(BigDecimalMath.atan(bd18,mc))
                   );
    System.out.println("\n Dase pi = " + piDase);
    // https://en.wikipedia.org/wiki/Zacharias_Dase

    BigDecimal bd16 = BigDecimal.valueOf(16);
    // BigDecimal bd1239 = BigDecimal.valueOf(1.0/239);  // BAD!!
           // Leads to a very poor result, since 1/239
           // is a recurring fraction
    BigDecimal bd1239 = (BigDecimal.ONE).
          divide(BigDecimal.valueOf(239), 100, RoundingMode.HALF_UP);

    // Machin 3). pi = 16*arctan(1/5) - 4*arctan(1/239)
    BigDecimal piENIAC = 
        (bd16.multiply(
            BigDecimalMath.atan(bd15,mc))).
        subtract(
          (bd4.multiply(
            BigDecimalMath.atan(bd1239,mc)))
                );
    System.out.println("\nENIAC pi = " + piENIAC);
     /* The first electronic computation of pi was made in 
        September 1949 on the ENIAC (Electronic Numerical Integrator 
        and Computer) at the Army Ballistic Research Lab in 
        Aberdeen, Maryland. It calculated pi to 2037 decimal places 
        in 70 hours using this formula. */

    // -----------------
    System.out.println("\nDistribution of digits in pi");
    
    BigDecimal piJML= Pi.pi(Scale.valueOf(100));  // 100 dp
    System.out.println("\n  JML pi = " + piJML);
      // https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula

    int[] counts = countDigits(piJML.toString());
    System.out.println("Counts of its digits: " + 
                        Arrays.toString(counts));

    System.out.println("Generating 10,000 digit version...");
    piJML= Pi.pi(Scale.valueOf(10000));   
    counts = countDigits(piJML.toString());
    System.out.println("Counts of its digits:\n  " + 
                        Arrays.toString(counts));

    String pi10000 = Nut.loadNumberString("pi10000.txt");
    System.out.println("\nFirst 30 dp of pi10000: " + 
                         pi10000.substring(0, 32));  // with "3."
    counts = countDigits(pi10000);
    System.out.println("Counts of digits in pi10000:\n  " + 
                                     Arrays.toString(counts));

    // Nut.vectorgram(10, piJML.toString());
    Nut.vectorgram(10, pi10000);   // ten directions
  } /* main */


  private static int[] countDigits(String s)
  {
    int[] counts = new int[10];
    for(int i=0; i < s.length(); i++) {
      char ch = s.charAt(i);
      if (Character.isDigit(ch))
        counts[ch - '0']++;
    }
    return counts;
  }  // end of countDigits()


}  // end of CalcPi class