// Big_math.java  generate Pi, e to 100 fraction bits (error<1/2^100)
//              and a few other sample computations and functions

import java.math.BigDecimal;

public class Big_math
{
  BigDecimal epsilon; // based on desired precision
  BigDecimal natural_e;
    BigDecimal pi, npi;
  BigDecimal twopi;
  BigDecimal halfpi;
  int prec  = 100; // digits
  int nbits = 332; // precision in bits = about 3.32 precision in digits

  public Big_math() // constructor
  {
    // System.out.println("Big_math");

    // "constants" needed by other functions and users
    natural_e = naturalE(nbits); /* precision */
    epsilon = new BigDecimal("1");
    for(int i=0; i<nbits; i++)
    {
      epsilon = epsilon.multiply(new BigDecimal("0.5"));
    }
    epsilon = epsilon.setScale(nbits,BigDecimal.ROUND_DOWN);
    System.out.println("epsilon="+epsilon);

    System.out.println("test exp");
    BigDecimal e1 = exp(new BigDecimal("1"));
    BigDecimal e2 = exp(new BigDecimal("2"));
    BigDecimal diff = e1.subtract(sqrt(e2));
    System.out.println("     e="+natural_e);
    System.out.println("exp(1)="+e1);
    System.out.println("     err e-exp(1)="+natural_e.subtract(e1));
    System.out.println("exp(2)="+e2);
    System.out.println("exp(1)-sqrt(exp(2))="+diff);
    System.out.println("exp(.001)="+exp(new BigDecimal("0.001")));
    System.out.println("exp(100)="+exp(new BigDecimal("100")));
    System.out.println();

    // test atan and compute Pi
    System.out.println("test atan, compute Pi");
    BigDecimal a1 = atan(new BigDecimal("1"));
    BigDecimal ep1 = epsilon.add(new BigDecimal("1"));
    BigDecimal b1  = (new BigDecimal("2")).subtract((sqrt(new BigDecimal("3"))));
    b1 = b1.setScale(prec,BigDecimal.ROUND_DOWN);
    BigDecimal a2 = atan(ep1.divide(sqrt(new BigDecimal("3")),BigDecimal.ROUND_DOWN));
    BigDecimal a3 = atan(b1);
    BigDecimal b2 = ((new BigDecimal("2")).multiply(sqrt(b1))).subtract(new BigDecimal("1"));
    b2 = b2.setScale(prec,BigDecimal.ROUND_DOWN);
    b2 = b2.divide(b1, BigDecimal.ROUND_DOWN);
    BigDecimal a4 = atan(b2);
    System.out.println("atan(1)="+a1);
    System.out.println("atan(1/sqrt(3))="+a2);
    System.out.println("atan(2-sqrt(3))="+a3);
    System.out.println("atan(b2)="+a4);
    System.out.println("Pi  4*a1="+a1.multiply(new BigDecimal("4")));
    System.out.println("Pi  6*a2="+a2.multiply(new BigDecimal("6")));
    System.out.println("Pi 12*a3="+a3.multiply(new BigDecimal("12")));
    BigDecimal pi4 = a4.multiply(new BigDecimal("24"));
    System.out.println("Pi 24*a4="+pi4);
    npi = Npi();
    System.out.println("Pi Npi  ="+npi);
    System.out.println("     err npi-Pi24="+npi.subtract(pi4));
    System.out.println();

    // test sin
    System.out.println("test sin");
    pi = pi4; // best estimate
    BigDecimal pid2 = pi.divide(new BigDecimal("2"));
    BigDecimal sin50 = sin((pi.multiply(new BigDecimal("50")).add(pid2)));
    System.out.println("sin(50pi+pi/2)= 1 ="+sin50);
    BigDecimal sin50n = sin((pi.multiply(new BigDecimal("-50")).subtract(pid2)));
    System.out.println("sin(-50pi-pi/2)= -1 ="+sin50n);
    System.out.println(" ");
  
    System.out.println("fraction of Pi as bits");
    pi4 = pi4.subtract(new BigDecimal("3")); // fraction
    for(int j=0; j<20; j++)
    {
      for(int i=0; i<50; i++)
      {
        pi4 = pi4.multiply(new BigDecimal("2"));
        if(pi4.compareTo(new BigDecimal("1"))>=0)
        {
          System.out.print("1");
          pi4 = pi4.subtract(new BigDecimal("1"));
        }
        else
        {
          System.out.print("0");
        }
      }
      System.out.println();
    }
    System.out.println();

    System.out.println("fraction of e as bits");
    e1 = natural_e.subtract(new BigDecimal("2")); // make fraction
    for(int j=0; j<20; j++)
    {
      for(int i=0; i<50; i++)
      {
        e1 = e1.multiply(new BigDecimal("2"));
        if(e1.compareTo(new BigDecimal("1"))>=0)
        {
          System.out.print("1");
          e1 = e1.subtract(new BigDecimal("1"));
        }
        else
        {
          System.out.print("0");
        }
      }
      System.out.println();
    }
    System.out.println();

    // exp(Pi*sqrt(163)) integer ?
    BigDecimal psqrt = npi.multiply(sqrt(new BigDecimal("163")));
    psqrt = psqrt.setScale(prec,BigDecimal.ROUND_DOWN);
    BigDecimal iq = exp(psqrt);
    iq = iq.setScale(prec,BigDecimal.ROUND_DOWN);
    System.out.println("e^Pi*sqrt(163)  ="+iq);
    System.out.println();

    System.out.println("end Big_math");
  }

  BigDecimal factorial(BigDecimal n)
  {
    if(n.compareTo(new BigDecimal("1"))<=0) return new BigDecimal("1");
    return n.multiply(factorial(n.subtract(new BigDecimal("1"))));
  }

  BigDecimal sqrt(BigDecimal x)
  {
    BigDecimal y = new BigDecimal("1");
    BigDecimal yn = y;
    BigDecimal xs = x.add(epsilon);
    for(int i=0; i<25; i++)
    {
      yn = (y.add(xs.divide(y,BigDecimal.ROUND_DOWN))).multiply(new BigDecimal("0.5"));
      yn = yn.setScale(nbits,BigDecimal.ROUND_DOWN);
      if(((yn.subtract(y)).abs()).compareTo(epsilon)<=0) return yn;
      y = yn;
    }
    return yn;
  }

  BigDecimal exp_series(BigDecimal x) // abs(x)<=0.5, prec digits 
  {                                   // prec digits returned 
    BigDecimal fact = new BigDecimal("1"); // factorial
    BigDecimal xp   = new BigDecimal("1"); // power of x
    BigDecimal y    = new BigDecimal("1"); // sum of series on x
    int n;

    n = (2*prec)/3;
    for(int i=1; i<n; i++)
    {
      fact = fact.multiply(new BigDecimal(i));
      fact = fact.setScale(prec,BigDecimal.ROUND_DOWN);
      xp   = xp.multiply(x);
      xp = xp.setScale(prec,BigDecimal.ROUND_DOWN);
      y = y.add(xp.divide(fact, BigDecimal.ROUND_DOWN));
    }
    y = y.setScale(prec,BigDecimal.ROUND_DOWN);
    System.out.println("exp_series x="+x);
    System.out.println("           y="+y);
    return y;
  }

  BigDecimal exp(BigDecimal x)
  {
    BigDecimal y = new BigDecimal("1.0");  // sum of series on xc
    BigDecimal xc;                         // x - j
    BigDecimal ep = natural_e;             // e^j
    BigDecimal j  = new BigDecimal("1");
    BigDecimal one = new BigDecimal("1.0");
    BigDecimal half = new BigDecimal("0.5");
    BigDecimal xp;                         // positive, then invert
    
    if(x.abs().compareTo(half) <0) return exp_series(x);
    if(x.compareTo(new BigDecimal("0"))>0) // positive
    {
      while(x.compareTo(j.add(one))>0)
      {
        ep = ep.multiply(natural_e);
        j = j.add(one);
      }
      xc = x.subtract(j);
      y = ep.multiply(exp_series(xc));
      y = y.setScale(prec,BigDecimal.ROUND_DOWN);
      return y;
    }
    else // negative
    {
      xp = x.negate();
      while(xp.compareTo(j.add(one))>0)
      {
        ep = ep.multiply(natural_e);
        j = j.add(one);
      }
      xc = xp.subtract(j);
      y = ep.multiply(exp_series(xc));
      y = y.setScale(prec,BigDecimal.ROUND_DOWN);
      return (one.add(epsilon)).divide(y, BigDecimal.ROUND_DOWN);
    }
  } // end exp

  BigDecimal sin(BigDecimal x)
  {
    BigDecimal y;
    BigDecimal xc;
    BigDecimal tpi = pi.multiply(new BigDecimal("2"));
    
    if(x.abs().compareTo(tpi) <0) return sin_series(x);
    xc = x;
    if(xc.compareTo(new BigDecimal("0"))>0) // positive
    {
      while(xc.compareTo(tpi)>0)
      {
        xc = xc.subtract(tpi);
      }
      y = sin_series(xc);
      y = y.setScale(prec,BigDecimal.ROUND_DOWN);
      return y;
    }
    else // negative
    {
      while(xc.compareTo(tpi)<0)
      {
	xc = xc.add(tpi);
      }
      y = sin_series(xc);
      return y.setScale(prec,BigDecimal.ROUND_DOWN);
    }
  } // end sin

  BigDecimal sin_series(BigDecimal x) // abs(x)<=0.5, prec digits 
  {                                   // prec digits returned 
    BigDecimal fact = new BigDecimal("1"); // factorial
    BigDecimal xp   = x; // power of x
    BigDecimal y    = x; // sum of series on x
    int n;

    n = (2*prec)/3;
    for(int i=3; i<n; i=i+2)
    {
      fact = factorial(new BigDecimal(i));
      fact = fact.setScale(prec,BigDecimal.ROUND_DOWN);
      xp   = ((xp.multiply(x)).multiply(x)).negate();
      xp = xp.setScale(prec,BigDecimal.ROUND_DOWN);
      y = y.add(xp.divide(fact, BigDecimal.ROUND_DOWN));
    }
    y = y.setScale(prec,BigDecimal.ROUND_DOWN);
    return y;
  } // end sin_series

  BigDecimal atan(BigDecimal x)
  {
    // atan(x)=x - x^3/3 + x^5/5 - x^7/7 + x^9/9 ...
    BigDecimal f = new BigDecimal("1");
    BigDecimal y = x;
    y = y.setScale(prec,BigDecimal.ROUND_DOWN);
    BigDecimal xp = y; // x to power
    xp = xp.setScale(prec,BigDecimal.ROUND_DOWN);
    
    BigDecimal xs = xp.add(epsilon); // extended precision for divide
    xs = xs.setScale(prec,BigDecimal.ROUND_DOWN);
    BigDecimal t = epsilon;
    t = t.setScale(prec,BigDecimal.ROUND_DOWN);

    System.out.println("atan("+x+")");
    for(int i=3; i<prec; i=i+2)
    {
      xp = ((xp.multiply(x)).multiply(x)).negate();
      xp = xp.setScale(prec,BigDecimal.ROUND_DOWN);
      xs = xp.add(epsilon);
      xs = xs.setScale(prec,BigDecimal.ROUND_DOWN);
      f = new BigDecimal(i);
      t = xs.divide(f, BigDecimal.ROUND_DOWN);
      t = t.setScale(prec,BigDecimal.ROUND_DOWN);
      y = y.add(t);
      y = y.setScale(prec,BigDecimal.ROUND_DOWN);
      // System.out.println("i="+i+", t="+t);
      // System.out.println("i="+i+", y="+y);
      if((t.abs()).compareTo(epsilon)<=0) return y;
    }
    return y;
  }

  BigDecimal Npi()
  {
    BigDecimal sum = new BigDecimal("1");
    BigDecimal one = new BigDecimal("1");
    BigDecimal two = new BigDecimal("2");
    BigDecimal three = new BigDecimal("3");
    BigDecimal threen = new BigDecimal("3");
    BigDecimal denom = new BigDecimal("3");
    BigDecimal frac;
    BigDecimal twoRootThree;

    twoRootThree = two.multiply(sqrt(threen));
    one = one.setScale(prec,BigDecimal.ROUND_DOWN); // for divide

    for(int i=1; i<prec; i++) // prec digits, subtract, add  each iteration
    {
      frac = denom.multiply(threen);
      frac = frac.setScale(prec,BigDecimal.ROUND_DOWN);
      frac = one.divide(frac, BigDecimal.ROUND_DOWN);
      // System.out.println("1/frac="+frac);
      sum  = sum.subtract(frac);
      threen = threen.multiply(three);
      denom = denom.add(two);
      frac = denom.multiply(threen);
      frac = frac.setScale(prec,BigDecimal.ROUND_DOWN);
      frac = one.divide(frac, BigDecimal.ROUND_DOWN);
      sum = sum.add(frac);
      sum = sum.setScale(prec,BigDecimal.ROUND_DOWN);
      threen = threen.multiply(three);
      denom = denom.add(two);
    }
    sum = twoRootThree.multiply(sum);
    sum = sum.setScale(prec,BigDecimal.ROUND_DOWN);
    return sum;
  }

  BigDecimal naturalE(int prec_dig)
  {
    BigDecimal sum  = new BigDecimal("1");
    BigDecimal fact = new BigDecimal("1");
    BigDecimal del  = new BigDecimal("1");
    BigDecimal one  = new BigDecimal("1");
    BigDecimal ten  = new BigDecimal("10");
    int prec_bits = (prec_dig*332)/100;

    one = one.setScale(prec_dig,BigDecimal.ROUND_DOWN);
    for(int i=0; i<prec_dig; i++) del = del.multiply(ten);
    for(int i=1; i<prec_bits; i++)
    {
      fact = fact.multiply(new BigDecimal(i));
      fact = fact.setScale(prec_dig,BigDecimal.ROUND_DOWN);
      sum = sum.add(one.divide(fact, BigDecimal.ROUND_DOWN));
      if(del.compareTo(fact) <0) break;
    }
    return sum.setScale(prec_dig,BigDecimal.ROUND_DOWN);
  }

  public static void main (String[] args)
  {
    new Big_math();
  }

} // end Big_math