microdouble.java

This is a Java library for performing floating-point calculations on small devices such as mobile p

          -85, -81, -78, -75, -71, -68, -65, -62, 
          -58, -55, -52, -48, -45, -42, -38, -35, 
          -32, -28, -25, -22, -18, -15, -12, -9, 
          -5, -2, 1, 5, 8, 11, 15, 18, 
          21, 25, 28, 31, 35, 38, 41, 44, 
          48, 51, 54, 58, 61, 64, 68, 71, 
          74, 78, 81, 84, 87, 91, 94, 97, 
          101, 104, 107, 111, 114, 117, 121, 124, 
          127, 131, 134, 137, 140, 144, 147, 150, 
          154, 157, 160, 164, 167, 170, 174, 177, 
          180, 184, 187, 190, 193, 197, 200, 203, 
          207, 210, 213, 217, 220, 223, 227, 230, 
          233, 237, 240, 243, 246, 250, 253, 256, 
          260, 263, 266, 270, 273, 276, 280, 283, 
          286, 289, 293 
  } ;

  /**
   * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Double.html#toString(double)">Double.toString(double)</a>.
   * <p>
   * String conversion is a bit of a gray area.  The J2SE implementation of
   * this function (<code>Double.toString(double)</code> has some problems.  
   * Often times it does not return the shortest valid String, even though it 
   * claims to do so, and it has a few
   * corner cases where it behaves oddly (e.g. 0.001 gets converted to 
   * the String "0.0010").
   * <p>
   * The implementation in MicroDouble uses a much simpler table-based 
   * algorithm.  It frequently returns slightly different results than 
   * <code>Double.toString(double)</code>.  Sometimes the results are better,
   * and sometimes worse.  Ususally the difference is confined to the last
   * character, which may be different or missing in one of the results.
   */
  public static String toString(long d) {
    return toString(d, 100);
  }
  
  /**
   * Returns a string representation of the double argument, rounded so that
   * the returned <code>String</code> is no longer than
   * <code>maxStringLength</code> characters (or 9 characters, if 
   * <code>maxStringLength</code> is less than 9).  
   *
   * @param      d   the <code>double</code> to be converted.
   * @param      maxStringLength the maximum length of the returned string 
   * @return     a string representation of the argument.
   *
   * @see #toString(long)
   */
  public static String toString(long d, int maxStringLength) {
    if (isNaN(d)) {
      return "NaN";
    }
    boolean n = unpackSign(d);
    if (isZero(d)) {
      return (n ? "-0.0" : "0.0");
    } else if (isInfinite(d)) {
      return (n ? "-Infinity" : "Infinity");
    }
    if (maxStringLength < 9) {
      maxStringLength = 9;
    }
    // convert from base 2 to base 10
    int base2x = unpackExponent(d);
    long base2m = unpackMantissa(d);
    int idx = base2x + 1075;
    int dx = idx % 11;
    base2m <<= dx;
    idx /= 11;
    int base10x = pow2x[idx];
    while (base2m <= 0xcccccccccccccccL) {
      base2m = (base2m << 3) + (base2m << 1); // base2m *= 10;
      base10x--;
    }
    long base10m = dpMul(base2m, pow2m[idx]);
    boolean roundedUp = false;
    while (true) {
      int r = (int) (base10m % 10);
      long mt = base10m / 10;
      int xt = base10x + 1;
      if (r != 0) {
        boolean rut;
        if ((r > 5) || ((r == 5) && (! roundedUp))) {
          rut = true;
          mt++;
        } else {
          rut = false;
        }
        long dt = decToDouble(n, xt, mt);
        if (dt != d) {
          if (rut) {
            mt--;
          } else {
            mt++;
          }
          rut ^= true;
          dt = decToDouble(n, xt, mt);
          if (dt != d) {
            break;
          }
        }
        roundedUp = rut;
      }
      base10m = mt;
      base10x = xt;
    }
    
    while (true) {
      String s = toString(n, base10x, base10m);
      if (s.length() <= maxStringLength) {
        return s;
      }
      int r = (int) (base10m % 10);
      base10m /= 10;
      base10x++;
      if ((r > 5) || ((r == 5) && (! roundedUp))) {
        roundedUp = true;
        base10m++;
      } else {
        roundedUp = false;
      }      
      while ((base10m % 10) == 0) {
        base10m /= 10;
        base10x++;
      }
    }
  }
  
  private static String toString(boolean negative, int base10x, long base10m) {
    StringBuffer sb = new StringBuffer(26);
    if (negative) {
      sb.append('-');
    }
    String s = Long.toString(base10m);
    base10x += s.length() - 1;
    boolean scientific = ((base10x < -3) || (base10x >= 7));
    int dp; // index of decimal point in final string
    if (scientific) {
      dp = 1;
    } else {
      dp = base10x + 1;
      if (dp < 1) {
        sb.append('0');
      }
    }
    for (int i=0; i<dp; i++) {
      if (i < s.length()) {
        sb.append(s.charAt(i));
      } else {
        sb.append('0');
      }
    }
    sb.append('.');
    if (dp >= s.length()) {
      sb.append('0');
    } else {
      for (int i=dp; i<s.length(); i++) {
        if (i < 0) {
          sb.append('0');
        } else {
          sb.append(s.charAt(i));
        }
      }
    }
    if (scientific) {
      sb.append('E');
      sb.append(Integer.toString(base10x));
    }
    return sb.toString();
  }

  private static final long ONE_EIGHTY =           0x4066800000000000L;
  private static final long TWO_HUNDRED =          0x4069000000000000L;
  
  /**
   * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#toDegrees(double)">Math.toDegrees(double)</a>.
   */
  public static long toDegrees(long angrad) {
    return div(mul(angrad, ONE_EIGHTY), PI);
  }
  
  /**
   * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#toRadians(double)">Math.toRadians(double)</a>.
   */
  public static long toRadians(long angdeg) {
    return mul(div(angdeg, ONE_EIGHTY), PI);
  }

  public static long toGradians(long angrad) {
    return div(mul(angrad, TWO_HUNDRED), PI);
  }

  public static long gradiansToRadians(long anggrad) {
    return mul(div(anggrad, TWO_HUNDRED), PI);
  }
  
  /////////////////////////////////////////////////////////////////////////////
  // Elementary functions.  Most are ported directly from fdlibm.
  /////////////////////////////////////////////////////////////////////////////

  private static long set(int newHiPart, int newLowPart) {
    return ((((long) newHiPart) << 32) | newLowPart);
  }

  private static long setLO(long d, int newLowPart) {
    return ((d & 0xFFFFFFFF00000000L) | newLowPart);
  }

  private static long setHI(long d, int newHiPart) {
    return ((d & 0x00000000FFFFFFFFL) | (((long) newHiPart) << 32));
  }

  private static int getHI(long d) {
    return ((int) (d >> 32));
  }
  
  private static int getLO(long d) {
    return ((int) d);
  }

  private static int ilogb(long d) {
    if (isZero(d)) {
      return 0x80000001;
    } else if (isNaN(d) || (isInfinite(d))) {
      return 0x7fffffff;
    }
    int x = (((int) (d >> 52)) & 0x7ff);
    if (x == 0) {
      long m = (d & FRACTION_MASK);
      while (m < IMPLIED_ONE) {
        m <<= 1;
        x--;
      }
    }
    return x - 1023;
  }
  
  /**
   * @return the magnitude of x with the sign of y
   */
  private static long copySign(long x, long y) {
    return (x & 0x7fffffffffffffffL) | (y & 0x8000000000000000L);
  }

  /** 
   * Returns the value of the first argument, multiplied by 2 raised to the
   * power of the second argument.  Note that the second argument is really
   * an <code>int</code>, not a <code>float</code> or <code>double</code>.
   *
   * @param d   a <code>double</code> value.
   * @param n   an <code>int</code> value.
   * @return  the value <code>d * 2<sup>n</sup></code>.
   */
  public static long scalbn(long d, int n) {
    if (isNaN(d)) {
      return NaN;
    } else if ((n == 0) || isInfinite(d) || isZero(d)) {
      return d;
    } else if (n >= 2098) {
      return copySign(POSITIVE_INFINITY, d);
    } else if (n <= -2099) {
      return copySign(ZERO, d);
    }
    int x = ((int) (d >> 52) & 0x7ff);
    int x2 = x + n;
    if ((x == 0) || (x2 <= 0)) { // argument and/or return value are subnormal 
      return pack(unpackSign(d), x2 - 1075, unpackMantissa(d));
    } else if (x2 >= 0x7ff) { // overflow
      return copySign(POSITIVE_INFINITY, d);
    }
    return ((d & 0x800fffffffffffffL) | (((long) x2) << 52));
  }

  /**
   * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#IEEEremainder(double, double)">Math.IEEEremainder(double, double)</a>.
   */
  public static long IEEEremainder(long d1, long d2) {
    if (isNaN(d1) || isNaN(d2) || isInfinite(d1) || isZero(d2)) {
      return NaN;
    } else if (isZero(d1) || isInfinite(d2)) {
      return d1;
    }
    int hx = getHI(d1); // high word of x 
    int lx = getLO(d1); // low  word of x 
    int hp = getHI(d2); // high word of p 
    int lp = getLO(d2); // low  word of p 
    boolean negative = unpackSign(d1);
    hp &= 0x7fffffff;
    hx &= 0x7fffffff;

    if (hp<=0x7fdfffff) d1 = mod(d1,scalbn(d2, 1)); // now x < 2p 
    if (((hx-hp)|(lx-lp))==0) return ZERO; //zero*x;
    d1  = abs(d1);
    d2  = abs(d2);
    if (hp<0x00200000) {
      if(gt(scalbn(d1, 1), d2)) {
        d1 = sub(d1, d2);
        if (ge(scalbn(d1, 1), d2)) d1 = sub(d1, d2);
      }
    } else {
      long p_half = scalbn(d2, -1);
      if (gt(d1, p_half)) {
        d1 = sub(d1, d2);
        if (ge(d1, p_half)) d1 = sub(d1, d2);
      }
    }
    if (negative) {
      return negate(d1);
    }
    return d1;
  }

  /**
   * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#sqrt(double)">Math.sqrt(double)</a>.
   */
  public static long sqrt(long d) {
    if (isZero(d)) {
      return d;
    } else if (unpackSign(d) || isNaN(d)) {
      return NaN;
    } else if (d == POSITIVE_INFINITY) {
      return d;
    }
    // f is positive, nonzero, and finite

    // unpack
    int x = unpackExponent(d);
    long m = unpackMantissa(d);
    // normalize 
    while (m < IMPLIED_ONE) {
      m <<= 1;
      x--;
    }
    // make exponent even
    if ((x & 1) != 0) {
      m <<= 1;
    }
    // compute final exponent
    x = (x >> 1) - 26;
    
    // generate sqrt(x) bit by bit
    m <<= 1;
    long q = 0L; // q = sqrt(x)
    long s = 0L;
    long r = 0x0020000000000000L;
    while (r != 0) {
      long t = s + r;
      if (t < m) {
        s = t + r;
        m -= t;
        q |= r;
      }