import com.sgi.math.CR; // See http://www.hpl.hp.com/personal/Hans_Boehm/crcalc/CRCalc.html
import java.util.Random;
class TrueSine {

  static CR ln2 = CR.valueOf(2).ln();
  static CR zero = CR.valueOf(0);
  static CR half = CR.valueOf(0.5);
  static CR one = CR.valueOf(1);
  static CR two = CR.valueOf(2);
  
  /**
     If paranoid, then compare two different ways
     of computing highOrderBit.
   */
  static boolean paranoid = false;

  /*
    This is used in a series of comparisons.
    Taking logarithms and converting to integers
    seems to have a minor rounding bug in it.
  */
  static CR[] negativePowersOfTwo;
  static {
    negativePowersOfTwo = new CR[2000];
    negativePowersOfTwo[0] = one;
    for (int i = 1; i < negativePowersOfTwo.length; i++) {
      negativePowersOfTwo[i] = negativePowersOfTwo[i-1].divide(two); 
    }
  }

  /** 
      Return the power of two equal to the high order bit
      of b.  Assumes |b| < 2.
  */
  static CR highOrderBit(CR b) {
    b = b.abs();

    for (int i = 0; i < negativePowersOfTwo.length; i++) {
      CR test = negativePowersOfTwo[i];
      if (b.compareTo(test, -64, -128) >= 0) {
	return test;
      }
    }
    return zero;
  }

  /** 
      Buggy version of highOrderBit
   */
  static CR buggyHighOrderBit(CR b) {
    long d = b.abs().ln().divide(ln2).subtract(half).longValue();
    CR hob = CR.valueOf(d).multiply(ln2).exp(); 
    return hob;
  }
  
  /**
     Takes the difference of a and b divided by the high
     order bit of b, and returns the negative of the
     logarithm base 2 of that quantity.  It is the "bit"
     in which a and b differ, where the actual mantissa
     bits are numbered from 0 (implicit leading 1) to 52.
  */
  static double relativeDiff(CR a, CR b) {
    CR tb = highOrderBit(b);
    if (paranoid) {
      CR tb2 = buggyHighOrderBit(b);
      if (tb.doubleValue() != tb2.doubleValue()) {
	System.out.println("High order bits differ for " + b.doubleValue() +
			   ", good=" + tb.doubleValue() + 
			   ", bad=" + tb2.doubleValue()); 
      }
    }
    return - (a.subtract(b).divide(tb).abs().ln().divide(ln2).doubleValue());
  }

  /**
     For unexceptional values, next(x) and prev(x) are the closest 
     unequal doubles greater and less than d in magnitude.
  */
  static double next(double d) {
    return Double.longBitsToDouble(Double.doubleToLongBits(d) + 1);
  }

  static double prev(double d) {
    return Double.longBitsToDouble(Double.doubleToLongBits(d) - 1);
  }

  /**
     Do a test of the high order bit function.
  */
  static void testHighOrderBit() {
    double a = 1.0;
    double b = next(a);
    double c = prev(a);
    for (int i = 0; i < 64; i++) {
      double hoba = highOrderBit(CR.valueOf(a)).doubleValue();
      double hobb = highOrderBit(CR.valueOf(b)).doubleValue();
      double hobc = highOrderBit(CR.valueOf(c)).doubleValue();

      System.out.println("a = " + a + ", b = " + b + ", c = " + c +
			 ", hoba = " + hoba + ", hobb = " + hobb + ", hobc = " + hobc);

      a = a / 2.0;
      b = b / 2.0;
      c = c / 2.0;
      }
  }

  /**
     Do a test of the high order bit function.
     This test doesn't seem to reveal the problems
     exposed by setting "paranoid" to true.
  */
  static void testBuggyHighOrderBit() {
    double a = 1.0;
    double b = Math.sqrt(2.0);

    for (int i = 0; i < 64; i++) {
      double a_minus = prev(a);
      double a_plus  = next(a);
      double b_minus = prev(b);
      double b_plus  = next(b);

      double hoba = highOrderBit(CR.valueOf(a)).doubleValue();
      double hobb = highOrderBit(CR.valueOf(b)).doubleValue();

      double hobam = highOrderBit(CR.valueOf(a_minus)).doubleValue();
      double hobbm = highOrderBit(CR.valueOf(b_minus)).doubleValue();

      double hobap = highOrderBit(CR.valueOf(a_plus)).doubleValue();
      double hobbp = highOrderBit(CR.valueOf(b_plus)).doubleValue();

      System.out.println(" a = " +    a + ",  am = " + a_minus + ",  ap = " + a_plus);
      System.out.println("ha = " + hoba + ", ham = " +   hobam + ", hap = " + hobap);

      System.out.println(" b = " +    b + ",  bm = " + b_minus + ",  bp = " + b_plus);
      System.out.println("hb = " + hobb + ", hbm = " +   hobbm + ", hbp = " + hobbp);

      System.out.println();

      a = a / 2.0;
      b = b / 2.0;
      }
  }

  /**
     Do a test of the relative difference function
  */
  static void testRelativeDiff() {
    double a = 0.5;
    double b = 1.0;
    for (int i = 1; i < 60; i++) {
      double c = a + b;
      if (c == b) break;
      System.out.println("i = " + i + ", a = " + a + ", rd = " + relativeDiff(CR.valueOf(c), CR.valueOf(b)));
      a = a / 2.0;
    }
  }

  /**
     Test Math.sin() (repeatably).
     If the result differs in by more than 1/2 LSB,
     report it.
  */
  static void testSine() {
    double scale9 = 1024 * 1024 * 1024; // 10 ** 9
    double scale18 = scale9 * scale9;  // 10 ** 18
    double scale36 = scale18 * scale18;  // 10 ** 36
    double scale288 = scale36 * scale36 * scale36 * scale36 * scale36 * scale36 * scale36 * scale36 ;  // 10 ** 288;
    double scale = scale288 * scale9 * 1024.0; // 10 ** 300;
      
    Random r = new Random(0x123456787654321L);
    for (int i = 0; i < 100000; i++) {
      double d = r.nextDouble() * Math.PI * scale;
      double js = Math.sin(d);
      CR     crs = CR.valueOf(d).sin();
      double dcrs = crs.doubleValue();

      double relative_diff = relativeDiff(crs, CR.valueOf(js));

      if (relative_diff < 53.0) {
	System.out.println("I=" + i + "\tX=" + d + "\tJS=" + js+"\tDIFF=" + 
			   Long.toHexString(Math.abs(Double.doubleToLongBits(js)-
						     Double.doubleToLongBits(dcrs))) +
			   " RELATIVE_DIFF=" + relative_diff);
      }
    }
  }

  public static void main (String[] args) {
    
    // testBuggyHighOrderBit();
    // testHighOrderBit();
    // testRelativeDiff();
    testSine();
  }
}