// Big_math.java  generate Pi, e to 100 fraction bits (error<1/2^100)
//                and a few other constants and functions
//                sqrt, sin, cos, atan, exp 
import java.math.BigDecimal;

public class Big_math
{
  BigDecimal epsilon; // based on desired precision
  BigDecimal natural_e;
  BigDecimal pi, npi;
  BigDecimal twopi;
  BigDecimal halfpi;
  BigDecimal fourthpi;
  BigDecimal rfact[] = new BigDecimal[50]; // reciprocal factorial about prec/2
  BigDecimal zero = new BigDecimal(0.0);
  BigDecimal one = new BigDecimal(1.0);
  BigDecimal two = new BigDecimal(2.0);
  BigDecimal half = new BigDecimal("0.5");
  int prec  = 100; // digits
                   // precision in bits = about 3.32 precision in digits

  public Big_math() // constructor
  {
    BigDecimal tmp;
    System.out.println("Big_math constructing");
    // "constants" needed by other functions and users
    natural_e = naturalE(prec); /* precision */
    System.out.println("natural_e="+natural_e);
    epsilon = half.pow(prec);
    System.out.println("epsilon="+epsilon);
    // precompute reciprocal factorials for for future use */
    rfact[0] = new BigDecimal(1.0);
    rfact[1] = new BigDecimal(1.0);
    rfact[2] = half;
    for(int i=3; i<50; i++)
    {
      rfact[i] = rfact[i-1].divide(new BigDecimal((double)i),prec,
                                   BigDecimal.ROUND_DOWN);
    }
    System.out.println("rfact[49]="+rfact[49]);
    // precompute pi, pi/4, pi/2, 2pi for future use */
    pi = Npi();
    System.out.println("pi="+pi);
    halfpi = pi.multiply(half);
    System.out.println("halfpi="+halfpi);
    fourthpi = halfpi.multiply(half);
    System.out.println("fourthpi="+fourthpi);
    System.out.println("Big_math constructed");
  }

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

  public 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(half);
      yn = yn.setScale(prec,BigDecimal.ROUND_DOWN);
      if(((yn.subtract(y)).abs()).compareTo(epsilon)<=0) return yn;
      y = yn;
    }
    return yn;
  }

  public 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, prec,BigDecimal.ROUND_DOWN));
    }
    y = y.setScale(prec,BigDecimal.ROUND_DOWN);
    System.out.println("exp_series x="+x);
    System.out.println("           y="+y);
    return y;
  }

  public 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 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, prec,BigDecimal.ROUND_DOWN);
    }
  } // end exp

  public 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


  public BigDecimal cos(BigDecimal x)
  {
      return sin(x.add(halfpi));
  } // end cos

  public 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

  public BigDecimal atan(BigDecimal x)
  {
    // atan(1/x)=Pi/2-atan(x)   atan(-x)=-atan(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, prec,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;
  }

  public BigDecimal atan2(BigDecimal yb, BigDecimal xb) // -pi to pi
  {
    // x = yb/xb  Q1++ sum,  Q2+- pi-sum,  Q3-- -pi+sum, Q4-+ -sum   
    // atan(1/x)=Pi/2-atan(x)   atan(-x)=-atan(x)
    // atan(x)=x - x^3/3 + x^5/5 - x^7/7 + x^9/9 ... -x^27/27
    // x2 = x1/(1+sqrt(1+x1*x1))  2^n*atan(x2) = atan(x1)  repeat
    BigDecimal sum = new BigDecimal(0.0);
    BigDecimal x, x2;
    BigDecimal tenm4 = new BigDecimal(0.0001);
    BigDecimal fct = one;
    if(xb.compareTo(zero)==0)
    {
      if(yb.compareTo(zero)>=0) return halfpi; 
      else return halfpi.negate();
    } 
    x = yb.divide(xb, prec, BigDecimal.ROUND_DOWN);
    x = x.abs();
    if(x.compareTo(one)>0) x = xb.divide(yb, prec, BigDecimal.ROUND_DOWN);
    x = x.abs();
    while(x.compareTo(tenm4)>0)
    {
      x = x.divide(sqrt((x.multiply(x).add(one)).add(one)),
                   prec, BigDecimal.ROUND_DOWN);
      fct = fct.multiply(two);
    }
    x2 = x.multiply(x);
    sum = x2.divide(new BigDecimal(27.0),prec,BigDecimal.ROUND_DOWN);
    for(int i=25; i>4; i=i-4)
    {
      sum = (sum.add(one.divide(new BigDecimal((double)i),
	     prec,BigDecimal.ROUND_DOWN))).multiply(x2);
      sum = (sum.subtract(one.divide(new BigDecimal((double)(i-2)),
	     prec,BigDecimal.ROUND_DOWN))).multiply(x2);
    }
    sum = (sum.add(one)).multiply(x);
    sum = sum.multiply(fct);
    sum = sum.setScale(prec, BigDecimal.ROUND_DOWN);
    if(yb.compareTo(zero)>=0)
    {
      if(xb.compareTo(zero)>=0) return sum; 
      else return pi.subtract(sum);
    }
    else
    {
	if(xb.compareTo(zero)>=0) return sum.negate(); 
	else return (pi.subtract(sum)).negate();
    }
    // never gets here
  } // end atan2
 
  BigDecimal Npi()
  {
    BigDecimal sum = 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 one  = new BigDecimal("1");
    BigDecimal lim  = new BigDecimal("10");

    lim = lim.pow(prec_dig);
    for(int i=1; i<prec_dig; i++)
    {
      fact = fact.multiply(new BigDecimal(i));
      if(fact.compareTo(lim)>0) break;
      sum = sum.add(one.divide(fact, prec_dig, BigDecimal.ROUND_DOWN));
    }
    return sum;
  }
} // end class Big_math