// File Numbers.java for lecture on computer numbers
// Contents: 
// * Binary representation of integers and floating-point numbers
// * Experiments with IEEE 754-2008 floating point numbers: 
//   infinities, NaNs, denormal numbers, etc.
// * Examples of accurate and inaccurate numeric algorithms
// * LaTeX code generator for tables of IEEE floating-point operations
// Peter Sestoft * sestoft@itu.dk * 2009-02-15, 2011-10-28, 2012-02-08

class Numbers {
  public static void main(String[] args) {
    byteExamples();
    intExamples();
    floatExamples();
    doubleExamples();
    // Works only because IEEE handles infinities correctly:
    System.out.println("R(1) = " + R(1));
    System.out.println("R(2) = " + R(2));
    System.out.println("R(3) = " + R(3));
    System.out.println("R(4) = " + R(4));
    // Negative and positive zero
    double zp = 0.0, zn = -0.0;
    System.out.println(zp == zn && 1/zp != 1/zn);  // true
    System.out.println((new Double(zp)).hashCode() + " " + (new Double(zn)).hashCode());
    System.out.println(-0.0 < 0.0);		   // false
    allDoubleOperations();
    Point[] ps1 = new Point[] { 
      new Point(2.1, 5.2), new Point(2.2, 5.4), new Point(2.4, 5.8) };
    linear1(ps1);
    linear2(ps1);
    Point[] ps2 = Point.moveAll(ps1, 1E7, 1E7);
    linear1(ps2);
    linear2(ps2);
    Point[] ps3 = Point.moveAll(ps1, 5E7, 5E7);
    linear1(ps3);
    linear2(ps3);
    // Testing sum
    final int N = 10000;
    final double[] xs = new double[2*N];
    for (int i=0; i<N; i++) {
      xs[2*i] = 1E12;
      xs[2*i+1] = -1;
    }
    final double result = N*(-1+1E12);
    System.out.printf("Exact result = %f%n", result);
    System.out.printf("Kahan error = %f; naive error = %f%n",
		      sumKahan(xs)-result, sum(xs) - result);
    quadratic1(1, 1E9, 1);
    quadratic2(1, 1E9, 1);
    // Examples of various anomalies:
    // oneTenth();
    // wrapAround();    
    // denormalFloat();
    // idempotence();
    // stewartSequence();
    cancellation();
    testnans();
    showExactDecimal();
    relativeComparison();
  }

  // ------------------------------------------------------------

  // Quadratic equation ax^2 + bx + c = 0 solved two ways 
 
  // Poor solution suffering from cancellation:
  public static void quadratic1(double a, double b, double c) {
    System.out.println("Quadratic1"); 
    double d = b * b - 4 * a * c;
    if (d < 0) 
      System.out.println("No solutions"); 
    else if (d == 0)
      System.out.printf("One solution: %g%n", -b / (2 * a)); 
    else { 
      double y = Math.sqrt(d);
      double x1 = (-b - y)/(2 * a);
      double x2 = (-b + y)/(2 * a);
      System.out.printf("Two solutions: %g and %g%n", x1, x2);
    }
  }

  // Better solution:
  public static void quadratic2(double a, double b, double c) {
    System.out.println("Quadratic2"); 
    double d = b * b - 4 * a * c;
    if (d < 0) 
      System.out.println("No solutions"); 
    else if (d == 0)
      System.out.printf("One solution: %g%n", -b / (2 * a)); 
    else {
      double y = Math.sqrt(d);
      double x1 = b > 0 ? (-b - y)/(2 * a) : (-b + y)/(2 * a);
      double x2 = c / (x1 * a);
      System.out.printf("Two solutions: %g and %g%n", x1, x2);
    }
  }

  // ------------------------------------------------------------

  public static void denormalFloat() {
    float y = (float)Math.pow(2,-126);
    float d = y/4.0f;
    float x = y+d;
    displayFloat(x);
    displayFloat(y);
    displayFloat(x-y);
  }

  // ------------------------------------------------------------

  public static void idempotence() {
    double z = Math.pow(2, 53);
    displayDouble(z);
    double zz = z+1;
    displayDouble(zz);
    System.out.println(z==zz);
  }

  // ------------------------------------------------------------
  // Cancellation, loss of significant digits
  
  public static void cancellation() {
    double v = 9876543210.2, w = 9876543210.1;
    double r = v-w;
    // The actual representable doubles are
    // v~ = 9876543210.20000076293945312500
    // w~ = 9876543210.10000038146972656250
    // The difference between these is:
    //      0000000000.10000038146972656250
    // So no surprise that the double result is
    // v-w =         0.10000038146972656000 
    System.out.printf("v = %.20f%n", v);
    displayDouble(v);
    System.out.printf("w = %.20f%n", w);
    displayDouble(w);
    System.out.printf("v-w = %.20f%n", r);
    displayDouble(r);
    displayExact("v   = ", v);
    displayExact("w   = ", w);
    displayExact("v-w = ", r);
    System.out.println(v+w-w != v);		  // true
  }

  // ------------------------------------------------------------
  // Dangerous loops

  // This doesn't terminate!
  public static void oneTenth() {
    double d = 0.0;
    while (d != 1.0) {
      d += 0.1;
      System.out.printf("%.20f%n", d);
    }
  }

  // This does terminate!
  public static void wrapAround() {
    int i = 1;
    while (i>0) 
      i++;
    System.out.print(i);
  }

  // ------------------------------------------------------------
  // Linear regression computed two ways

  public static class Point {
    public final double x, y;
    
    public Point(double x, double y) {
      this.x = x;
      this.y = y;
    }

    public Point move(double dx, double dy) {
      return new Point(x + dx, y + dy);
    }

    public static Point[] moveAll(Point[] ps, double dx, double dy) {
      Point[] res = new Point[ps.length];
      for (int i=0; i<ps.length; i++)
	res[i] = ps[i].move(dx, dy);
      return res;
    }
  }

  public static void linear1(Point[] ps) {
    final int n = ps.length;
    double SX = 0.0, SY = 0.0, SSX = 0.0, SXY = 0.0;
    for (int i=0; i<n; i++) {
      Point p = ps[i];
      SX += p.x;
      SY += p.y;
      SXY += p.x * p.y;
      SSX += p.x * p.x;
    }
    double beta = (SXY - SX*SY/n)/(SSX - SX*SX/n);
    double alpha = SY/n - SX/n * beta;
    System.out.printf("LINREG1: y = %f + %f * x%n", alpha, beta);
  }

  public static void linear2(Point[] ps) {
    final int n = ps.length;
    double SX = 0.0, SY = 0.0;
    for (int i=0; i<n; i++) {
      Point p = ps[i];
      SX += p.x;
      SY += p.y;
    }
    double EX = SX/n, EY = SY/n;
    double SDXDY = 0.0, SSDX = 0.0;
    for (int i=0; i<n; i++) {
      Point p = ps[i];
      double dx = p.x - EX, dy = p.y - EY;
      SDXDY += dx * dy;
      SSDX += dx * dx;
    }
    double beta = SDXDY/SSDX;
    double alpha = SY/n - SX/n * beta;
    System.out.printf("LINREG2: y = %f + %f * x%n", alpha, beta);
  }

  // ------------------------------------------------------------
  // From Goldberg 1991, due to Kahan; relative error only 2 eps!
  // Optimizer should not reduce C to (T-S)-Y to ((S+Y)-S)-Y to 0.0.

  public static double sumKahan(double[] xs) {
    double S = 0.0, C = 0.0;
    for (int i=0; i<xs.length; i++) {
      double Y = xs[i] - C, T = S + Y;
      C = (T - S) - Y;
      S = T;
    }
    return S;
  }

  // Naive summation of array of doubles:
  
  public static double sum(double[] xs) {
    double S = 0.0;
    for (int i=0; i<xs.length; i++) 
      S += xs[i];
    return S;
  }

  // ------------------------------------------------------------
  // Example from Goldberg 1991 of infinities; 
  // compute x/(x^2+1) as 1/(x + 1/x) to avoid overflows; 
  // this works for x==0 too:

  public static double G(double x) {
    return 1/(x + 1/x);
  }

  // ------------------------------------------------------------
  // Example from William Kahan 1981 "Why do we need..."

  public static double R(double z) {
    return 7-3/(z-2-1/(z-7+10/(z-2-2/(z-3))));
  }

  // ------------------------------------------------------------
  // Example from GW Stewart's "Afternotes..." 8.9-8.10
  // Should decrease steadily, but grows after k=20 approx.

  public static void stewartSequence() {
    double x1 = 1.0/3.0, x2 = 1.0/12.0;
    for (int k = 3; k < 40; k++) {
      double x3 = 2.25 * x2 - 0.5 * x1;
      System.out.println("x" + k + " = " + Math.log(x3));
      x1 = x2; 
      x2 = x3;
    }
  }

  // ------------------------------------------------------------
  // Fast relative comparison of doubles, modified from Bruce Dawson's
  // http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm

  // Return true if there are at most maxUlps representable doubles
  // between x and y.  Works correctly for positive and negative
  // zeros and subnormals.  Because size(int) = 32 < 52, it
  // considers all numbers unequal to infinities, and it will not
  // deem the first 2^32 NaNs, including the default NaN (produced by
  // mathematical functions), equal to anything else.  It will
  // consider any NaN equal to itself, unlike (==).

  // If d==MinValue, the largest negative long value, Math.abs(d)
  // returns d itself, so Math.abs(d) <= maxUlps holds for any
  // maxUlps, but this is not the desired result.  So this case must
  // be handled explicitly.  It arises eg. for x=2.0 and y=-2.0.  It
  // does not seem to be handled correctly by Bruce Dawson's code.
  
  public static boolean almostEquals(double x, double y, int maxUlps) {
    assert maxUlps >= 0;
    long xBits = Double.doubleToRawLongBits(x),
         yBits = Double.doubleToRawLongBits(y),
         MinValue = 1L << 63;
    // Make xyBits lexicographically ordered as two's-complement longs:
    if (xBits < 0)
      xBits = MinValue - xBits;
    if (yBits < 0)
      yBits = MinValue - yBits;
    long d = xBits - yBits;
    return d != MinValue && Math.abs(d) <= maxUlps;
  }

  // ------------------------------------------------------------
  // Test whether java.lang.Math methods preserve NaN payload

  private static void testnans() {
    for (int i=-100; i<100; i++) {
      double d = MakeNaN(i);
      if (ErrorCode(d) != i)
	System.out.printf("ErrorCode at %d%n", i);
      if (ErrorCode(Math.abs(d)) != i)
	System.out.printf("abs at %d%n", i);
      if (ErrorCode(Math.acos(d)) != i)
	System.out.printf("acos at %d%n", i);
      if (ErrorCode(Math.asin(d)) != i)
	System.out.printf("asin at %d%n", i);
      if (ErrorCode(Math.atan(d)) != i)
	System.out.printf("atan at %d%n", i);
      if (ErrorCode(Math.cbrt(d)) != i)
	System.out.printf("cbrt at %d%n", i);
      if (ErrorCode(Math.ceil(d)) != i)
	System.out.printf("ceil at %d%n", i);
      if (ErrorCode(Math.cos(d)) != i)
	System.out.printf("cos at %d%n", i);
      if (ErrorCode(Math.exp(d)) != i)
	System.out.printf("exp at %d%n", i);
      if (ErrorCode(Math.floor(d)) != i)
	System.out.printf("floor at %d%n", i);
      if (ErrorCode(Math.log(d)) != i)
	System.out.printf("log at %d%n", i);
      if (ErrorCode(Math.log10(d)) != i)
	System.out.printf("log10 at %d%n", i);
      if (ErrorCode(Math.rint(d)) != i)
	System.out.printf("rint at %d%n", i);
      if (ErrorCode(Math.sin(d)) != i)
	System.out.printf("sin at %d%n", i);
      if (ErrorCode(Math.sqrt(d)) != i)
	System.out.printf("sqrt at %d%n", i);
      if (ErrorCode(Math.tan(d)) != i)
	System.out.printf("tan at %d%n", i);
      // Arithmetic operations
      if (ErrorCode(d + 1) != i)
	System.out.printf("NaN+1 at %d%n", i);
      if (ErrorCode(1 + d) != i)
	System.out.printf("1+NaN at %d%n", i);
      if (ErrorCode(d - 1) != i)
	System.out.printf("NaN-1 at %d%n", i);
      if (ErrorCode(1 - d) != i)
	System.out.printf("1-NaN at %d%n", i);
      if (ErrorCode(d * 1) != i)
	System.out.printf("NaN*1 at %d%n", i);
      if (ErrorCode(1 * d) != i)
	System.out.printf("1*NaN at %d%n", i);
      if (ErrorCode(d / 1) != i)
	System.out.printf("NaN/1 at %d%n", i);
      if (ErrorCode(1 / d) != i)
	System.out.printf("1/NaN at %d%n", i);
      if (ErrorCode(d % 1) != i)
	System.out.printf("NaN%%1 at %d%n", i);
      if (ErrorCode(1 % d) != i)
	System.out.printf("1%%NaN at %d%n", i);
      // Two-argument Math functions
      if (ErrorCode(Math.atan2(d, 1)) != i)
	System.out.printf("atan2(NaN, 1) at %d%n", i);
      if (ErrorCode(Math.atan2(1, d)) != i)
	System.out.printf("atan2(1, NaN) at %d%n", i);
      if (ErrorCode(Math.IEEEremainder(d, 1)) != i)
  	System.out.printf("IEEEremainder(NaN, 1) at %d%n", i);
      if (ErrorCode(Math.IEEEremainder(1, d)) != i)
  	System.out.printf("IEEEremainder(1, NaN) at %d%n", i);
      if (ErrorCode(Math.max(d, 1)) != i)
	System.out.printf("max(NaN, 1) at %d%n", i);
      if (ErrorCode(Math.max(1, d)) != i)
	System.out.printf("max(1, NaN) at %d%n", i);
      if (ErrorCode(Math.min(d, 1)) != i)
	System.out.printf("min(NaN, 1) at %d%n", i);
      if (ErrorCode(Math.min(1, d)) != i)
	System.out.printf("min(1, NaN) at %d%n", i);
      if (ErrorCode(Math.pow(d, 1)) != i)
  	System.out.printf("pow(NaN, 1) at %d%n", i); 
      if (ErrorCode(Math.pow(1, d)) != i)
  	System.out.printf("pow(1, NaN) at %d%n", i);
    }
  }

  private static void showExactDecimal() {
    displayExact("double 0.125", 0.125);
    displayExact("float 0.125f", 0.125f);
    displayExact("double 0.1", 0.1);
    displayExact("float 0.1f", 0.1f);
    displayExact("double 0.01", 0.01);
    displayExact("float 0.01f", 0.01f);
  }

  private static void relativeComparison() {
    System.out.println("1.0 == 0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1 is " + 
		       (1.0 == 0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1));
    System.out.println("almostEquals(1.0, 0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1, 16) is " +
		       almostEquals(1.0, 0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1, 16));
  }

  // ------------------------------------------------------------
  // Manipulating NaNs 

  public static int ErrorCode(double d)
  {
    return (int)Double.doubleToRawLongBits(d);
  }

  public static double MakeNaN(int errorNumber)
  {
    long nanbits = Double.doubleToRawLongBits(Double.NaN);
    return Double.longBitsToDouble(nanbits | errorNumber);
  }
  
  // ------------------------------------------------------------
  // Generate LaTeX code for tables of arithmetic operations
  // on double = IEEE binary64

  public static void allDoubleOperations() {
    System.out.println("\\section{Results}");
    System.out.println("Results are as expected, except perhaps for ");
    System.out.println("\\texttt{Math.pow(NaN, +/-0.0)}, \\texttt{Math.log(-0.0)}, \\texttt{Math.log10(-0.0)}, and \\texttt{Math.sqrt(-0.0)}. \\\\");
    System.out.println("Source code in java/Numbers.java.");
    System.out.println("\\subsection{Arithmetic operators}");
    System.out.println("\\tt");
    displayAll2("+", new IFunction() { 
	public String invoke(double x, double y) { return (x+y) + ""; }});
    displayAll2("-", new IFunction() { 
	public String invoke(double x, double y) { return (x-y) + ""; }});
    displayAll2("*", new IFunction() { 
	public String invoke(double x, double y) { return (x*y) + ""; }});
    displayAll2("/", new IFunction() { 
	public String invoke(double x, double y) { return (x/y) + ""; }});
    displayAll2("\\%", new IFunction() { 
	public String invoke(double x, double y) { return (x%y) + ""; }});
    System.out.println("\\subsection{Comparison operators}");
    displayAll2("==", new IFunction() { 
	public String invoke(double x, double y) { return (x==y) + ""; }});
    displayAll2("!=", new IFunction() { 
	public String invoke(double x, double y) { return (x!=y) + ""; }});
    displayAll2("<", new IFunction() { 
	public String invoke(double x, double y) { return (x<y) + ""; }});
    displayAll2("<=", new IFunction() { 
	public String invoke(double x, double y) { return (x<=y) + ""; }});
    displayAll2(">", new IFunction() { 
	public String invoke(double x, double y) { return (x>y) + ""; }});
    displayAll2(">=", new IFunction() { 
	public String invoke(double x, double y) { return (x>=y) + ""; }});
    System.out.println("\\subsection{Two-argument mathematical functions}");
    displayAll2("Math.atan2", new IFunction() { 
	public String invoke(double x, double y) { return String.format("%.3f", Math.atan2(x,y)); }});
    displayAll2("Math.IEEEremainder", new IFunction() { 
	public String invoke(double x, double y) { return Math.IEEEremainder(x,y)+""; }});
    displayAll2("Math.max", new IFunction() { 
	public String invoke(double x, double y) { return Math.max(x,y)+""; }});
    displayAll2("Math.min", new IFunction() { 
	public String invoke(double x, double y) { return Math.min(x,y)+""; }});
    displayAll2("Math.pow", new IFunction() { 
	public String invoke(double x, double y) { return Math.pow(x,y)+""; }});
    System.out.println("\\subsection{One-argument mathematical functions}");
    displayAll1(
       new Function1("Math.abs") { 
	   public double fun(double x) { return Math.abs(x); } },
       new Function1("Math.acos") { 
	   public double fun(double x) { return Math.acos(x); } },
       new Function1("Math.asin") { 
	   public double fun(double x) { return Math.asin(x); } },
       new Function1("Math.atan") { 
	   public double fun(double x) { return Math.atan(x); } },
       new Function1("Math.ceil") { 
	   public double fun(double x) { return Math.ceil(x); } },
       new Function1("Math.cbrt") { 
	   public double fun(double x) { return Math.cbrt(x); } },
       new Function1("Math.cos") { 
	   public double fun(double x) { return Math.cos(x); } },
       new Function1("Math.exp") { 
	   public double fun(double x) { return Math.exp(x); } },
       new Function1("Math.floor") { 
	   public double fun(double x) { return Math.floor(x); } },
       new Function1("Math.log") { 
	   public double fun(double x) { return Math.log(x); } },
       new Function1("Math.log10") { 
	   public double fun(double x) { return Math.log10(x); } },
       new Function1("Math.rint") { 
	   public double fun(double x) { return Math.rint(x); } },
       new Function1("Math.sin") { 
	   public double fun(double x) { return Math.sin(x); } },
       new Function1("Math.signum") { 
	   public double fun(double x) { return Math.signum(x); } },
       new Function1("Math.sqrt") { 
	   public double fun(double x) { return Math.sqrt(x); } },
       new Function1("Math.tan") { 
	   public double fun(double x) { return Math.tan(x); } }
       );
  }

  interface IFunction {
    public String invoke(double x, double y);
  }

  abstract static class Function1 {
    public final String name;
    public Function1(String name) {
      this.name = name;
    }
    public String invoke(double x) {
      return String.format("%.3f", fun(x));
    }
    abstract public double fun(double x);
  }

  private final static double[] values 
    = { Double.NEGATIVE_INFINITY, -2.0, -0.0, 
	0.0, +2.0, Double.POSITIVE_INFINITY, Double.NaN };
  
  private static void displayAll2(String name, IFunction function) {
    System.out.println("\\noindent");
    System.out.println("\\begin{tabular}{r|rrrrrrr}");
    System.out.print("\\multicolumn{1}{c|}{" + name + "}");
    for (int i=0; i<values.length; i++) 
      System.out.print(" & " + fix("" + values[i]));
    System.out.println("\\\\\\hline");
    for (int i=0; i<values.length; i++) { 
      System.out.print(fix("" + values[i]));
      for (int j=0; j<values.length; j++) { 
	double x = values[i], y = values[j];
	System.out.print(" & " + fix(function.invoke(x, y)));
      }
      System.out.println("\\\\");
    }
    System.out.println("\\end{tabular}\n\\vspace{0.5cm}\n");
  }

  private static void displayAll1(Function1... functions) {
    System.out.println("\\noindent");
    System.out.println("\\begin{tabular}{l|rrrrrrr}");
    for (int i=0; i<values.length; i++) 
      System.out.print(" & " + fix("" + String.format("%.3f", values[i])));
    System.out.println("\\\\\\hline");
    for (int i=0; i<functions.length; i++) { 
      Function1 function = functions[i];
      System.out.print(function.name);
      for (int j=0; j<values.length; j++) { 
	double x = values[j];
	System.out.print(" & " + fix(function.invoke(x)));
      }
      System.out.println("\\\\");
    }
    System.out.println("\\end{tabular}\n\\vspace{0.5cm}\n");
  }

  private static String fix(String s) {
    if (s.equals("Infinity"))
      return "+Inf";
    else if (s.equals("-Infinity"))
      return "-Inf";
    else
      return s;
  }

  private static String repeat(String s, int n) {
    StringBuilder sb = new StringBuilder();
    for (int i=0; i<n; i++)
      sb.append(s);
    return sb.toString();
  }

  // ------------------------------------------------------------
  // Conversions to bit patterns

  public static void byteExamples() {
    System.out.println("\n8-bit two's complement integers, byte:"); 
    displayByte(0);
    displayByte(1);
    displayByte(-1);
    displayByte(2);
    displayByte(-2);
    displayByte(6);
    displayByte(15);
    displayByte(127);
    displayByte(-127);
    displayByte(-128);
    displayByte(127+1);
  }

  public static void intExamples() {
    System.out.println("\n32-bit two's complement integers, int:"); 
    displayInt(0);
    displayInt(1);
    displayInt(-1);
    displayInt(2);
    displayInt(-2);
    displayInt(6);
    displayInt(15);
    displayInt(Integer.MAX_VALUE);
    displayInt(-Integer.MAX_VALUE);
    displayInt(Integer.MIN_VALUE);
    displayInt(Integer.MAX_VALUE+1);
  }

  public static void floatExamples() {
    System.out.println("\nIEEE 754 32-bit, float:"); 
    displayFloat(0.0f);
    displayFloat(-0.0f);
    displayFloat(1.0f);
    displayFloat(0.5f);
    displayFloat(-118.625f);
    displayFloat(Float.POSITIVE_INFINITY);
    displayFloat(Float.NEGATIVE_INFINITY);
    displayFloat(Float.NaN);
    displayFloat(Float.MIN_NORMAL/16.0f);
    displayFloat(0.1f);
    displayFloat(0.1f+0.1f+0.1f+0.1f+0.1f+0.1f+0.1f+0.1f+0.1f+0.1f);
  }

  public static void doubleExamples() {
    System.out.println("\nIEEE 754 64-bit, double:"); 
    displayDouble(0.0);
    displayDouble(-0.0);
    displayDouble(1.0);
    displayDouble(0.5);
    displayDouble(-118.625);
    displayDouble(Double.POSITIVE_INFINITY);
    displayDouble(Double.NEGATIVE_INFINITY);
    displayDouble(Double.NaN);
    displayDouble(Double.MIN_NORMAL/16.0);
    displayDouble(Math.log(0.0));
    displayDouble(Math.log(-1));
    displayDouble(Math.sqrt(-1));
    displayDouble(Math.asin(1.5));
    displayDouble(0.1);
    displayDouble(0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1);
  }

  public static void displayByte(int n) {
    byte b = (byte)n;
    System.out.println(toBits(b) + " = " + b);
  }

  public static void displayInt(int n) {
    System.out.println(toBits(n) + " = " + n);
  }

  public static void displayFloat(float f) {
    System.out.println(toBits(f) + " = " + f);
  }

  public static void displayDouble(double d) {
    System.out.println(toBits(d) + " = " + d);
  }

  public static void displayExact(String s, double d) {
    System.out.println(s + " = " + new java.math.BigDecimal(d).toString());
  }


  public static String toBits(float f) {
    StringBuilder sb = insertBlanks(toBits(Float.floatToRawIntBits(f), 32), 
				    9, 1);
    return sb.toString();
  }
  
  public static String toBits(double d) {
    StringBuilder sb = insertBlanks(toBits(Double.doubleToRawLongBits(d), 64),
				    12, 1);
    return sb.toString();
   }
  
  public static String toBits(byte b) {
    return toBits(b, 8);
  }

  public static String toBits(int n) {
    StringBuilder sb = insertBlanks(toBits(n, 32), 24, 16, 8);
    return sb.toString();
  }

  public static String toBits(long n) {
    StringBuilder sb = insertBlanks(toBits(n, 64), 
				    56, 48, 40, 32, 24, 16, 8);
    return sb.toString();
  }

  private static StringBuilder insertBlanks(String s, int... pos) {
    StringBuilder sb = new StringBuilder(s);
    for (int i=0; i<pos.length; i++)
      sb.insert((int)pos[i], ' ');
    return sb;
  }

  public static String toBits(long n, int size) {
    char[] cs = new char[size];
    for (int i=1; i<=size; i++) {
      cs[size-i] = (n & 1) != 0 ? '1' : '0';
      n >>= 1;
    }
    return new String(cs);
  }
}