microdouble.java

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

      return NaN;
    }

    boolean n1 = unpackSign(d1);
    boolean n2 = unpackSign(d2);
    
    // special handling of infinity
    boolean i1 = isInfinite(d1);
    boolean i2 = isInfinite(d2);
    if (i1 || i2) {
      if (i1 && i2) {
        if (n1 != n2) {
          // infinites of opposite sign -> NaN
          return NaN;
        } else {
          // infinites of same sign -> infinity the same sign
          return d1;
        }
      } else if (i1) {
        return d1; // infinite + finite = infinite
      } else {
        return d2; // finite + infinite = infinite
      }
    }
    
    // special handling of zero
    boolean z1 = isZero(d1);
    boolean z2 = isZero(d2);
    if (z1 || z2) {
      if (z1 && z2) {
        if (n1 != n2) {
          // zeros of opposite sign -> positive zero
          return ZERO;
        } else {
          return d1; // zeros of same sign -> zero of the same sign
        }
      } else if (z1) {
        return d2; // zero + nonzero = nonzero
      } else {
        return d1; // nonzero + zero = nonzero
      }
    }
    
    // unpack, and add 3 guard digits
    long m1 = unpackMantissa(d1) << 3;
    int x1 = unpackExponent(d1) - 3;
    long m2 = unpackMantissa(d2) << 3;
    int x2 = unpackExponent(d2) - 3;
    
    // make exponents equal
    int dx = x1 - x2;
    if (dx > 0) {
      m2 = BitUtils.stickyRightShift(m2, dx);
      x2 = x1;
    } else if (dx < 0) {
      m1 = BitUtils.stickyRightShift(m1, -dx);
      x1 = x2;
    }

    // if the signs are different, negate the smaller mantissa and choose
    // the sign of the larger
    if (n1 ^ n2) { 
      if (m1 > m2) {
        m2 = -m2;
      } else {
        m1 = -m1;
        n1 = n2;
      }
    }
    
    // add (or subtract) mantissas
    m1 += m2;

    // pack result, and handle special case of zero (which always returns +0.0) 
    long d = pack(n1, x1, m1);
    if (d == NEGATIVE_ZERO) {
      return ZERO;
    }
    return d;
  }
  
  /**
   * Returns the difference of the two <code>double</code> arguments according to 
   * <a href="http://java.sun.com/docs/books/jls/second_edition/html/expressions.doc.html#13510">section
   * 15.18.2 of the JLS</a>.
   * <p>
   * This method takes the place of the binary <code>-</code> operator.
   *
   * @param   d1   the first <code>double</code> value 
   * @param   d2   the second <code>double</code> value
   * @return  the difference of the two arguments
   */
  public static long sub(long d1, long d2) {
    return add(d1, negate(d2));
  }

  /**
   * Returns the product of the two <code>double</code> arguments according to 
   * <a href="http://java.sun.com/docs/books/jls/second_edition/html/expressions.doc.html#5036">section
   * 15.17.1 of the JLS</a>.
   * <p>
   * This method takes the place of the <code>*</code> operator.
   *
   * @param   d1   the first <code>double</code> value
   * @param   d2   the second <code>double</code> value
   * @return  the product of the two arguments
   */
  public static long mul(long d1, long d2) {
    if (isNaN(d1) || isNaN(d2)) {
      return NaN;
    }

    boolean negative = unpackSign(d1) ^ unpackSign(d2);
    
    // special handling of infinity
    if (isInfinite(d1) || isInfinite(d2)) {
      if (isZero(d1) || isZero(d2)) {
        return NaN;
      } else {
        return (negative ? NEGATIVE_INFINITY : POSITIVE_INFINITY);
      }
    }
    
    // unpack
    long m1 = unpackMantissa(d1);
    int x1 = unpackExponent(d1);
    long m2 = unpackMantissa(d2);
    int x2 = unpackExponent(d2);
    
    // compute the resultant exponent
    x1 += x2;
    
    // compute the resultant mantissa using double-precision integer 
    // multiplication with 28 bit words
    long m11 = m1 & 0x0fffffff;
    long m12 = m1 >> 28;
    long m21 = m2 & 0x0fffffff;
    long m22 = m2 >> 28;
    
    long t1 = m11 * m21;
    long t2 = (m11 * m22) + (m12 * m21);
    long t3 = m12 * m22;
    
    t1 += (t2 & 0x0fffffff) << 28;
    t3 += t2 >>> 28;
    t3 += t1 >>> 56;
    t1 <<= 8;
    // the 128 bit result is now in t3t1

    if (t3 == 0) {
      // the high 64 bits are zero and can be ignored.
      return pack(negative, x1, t1);
    }
    
    // shift the result left into t3 and discard excess precision
    int s = BitUtils.countLeadingZeros(t3);
    x1 += 56 - s;
    t3 <<= s;
    t3 |= t1 >>> (64 - s);
    if ((t1 << s) != 0) {
      // discarded low bits go into the sticky bit
      t3 |= 1;
    }
    
    // round and pack the result
    return pack(negative, x1, t3);
  }
  
  /**
   * Returns the quotient of the two <code>double</code> arguments according to 
   * <a href="http://java.sun.com/docs/books/jls/second_edition/html/expressions.doc.html#5047">section
   * 15.17.2 of the JLS</a>.
   * <p>
   * This method takes the place of the <code>/</code> operator.
   *
   * @param   d1   the <code>double</code> dividend 
   * @param   d2   the <code>double</code> divisor
   * @return  the quotient of the two arguments
   */
  public static long div(long d1, long d2) {
    if (isNaN(d1) || isNaN(d2)) {
      return NaN;
    }

    boolean negative = unpackSign(d1) ^ unpackSign(d2);
    
    // special handling of infinity
    boolean n1 = isInfinite(d1);
    boolean n2 = isInfinite(d2);
    if (n1 || n2) {
      if (n1 && n2) {
        return NaN;
      } else if (n1) {
        return (negative ? NEGATIVE_INFINITY : POSITIVE_INFINITY);
      } else {
        return (negative ? NEGATIVE_ZERO : ZERO);
      }
    }
    // neither value is infinite
    
    // special handling of zero
    n1 = isZero(d1);
    n2 = isZero(d2);
    if (n1 || n2) {
      if (n1 && n2) {
        return NaN;
      } else if (n1) {
        return (negative ? NEGATIVE_ZERO : ZERO);
      } else {
        return (negative ? NEGATIVE_INFINITY : POSITIVE_INFINITY);
      }
    }
    // neither value is zero
    
    // unpack
    long m1 = unpackMantissa(d1);
    int x1 = unpackExponent(d1);
    long m2 = unpackMantissa(d2);
    int x2 = unpackExponent(d2);

    // shift, divide, mod, repeat
    long m = 0;
    x1 -= x2;
    while (true) {
      int s = Math.min(BitUtils.countLeadingZeros(m1) - 1, 
              BitUtils.countLeadingZeros(m));
      if (s <= 8) {
        if (m1 != 0) {
          m |= 1;
        }
        break;
      }
      m1 <<= s;
      m <<= s;
      x1 -= s;
      m |= m1 / m2;
      m1 %= m2;
    }
    return pack(negative, x1, m);
  }
  
  /**
   * Returns the remainder of the two <code>double</code> arguments according to 
   * <a href="http://java.sun.com/docs/books/jls/second_edition/html/expressions.doc.html#24956">section
   * 15.17.3 of the JLS</a>.
   * <p>
   * This method takes the place of the <code>%</code> operator.
   *
   * @param   d1   the <code>double</code> dividend 
   * @param   d2   the <code>double</code> divisor
   * @return  the remainder of the two arguments
   * @see #IEEEremainder(long, long)
   */
  public static long mod(long d1, long d2) {
    if (isNaN(d1) || isNaN(d2) || isInfinite(d1) || isZero(d2)) {
      return NaN;
    } else if (isZero(d1) || isInfinite(d2)) {
      return d1;
    }
    
    // unpack
    int x1 = unpackExponent(d1);
    int x2 = unpackExponent(d2);
    if (x1 < x2) {
      return d1;
    }
    boolean n = unpackSign(d1);
    long m1 = unpackMantissa(d1);
    long m2 = unpackMantissa(d2);
    if (x1 == x2) {
      m1 %= m2;
    } else {
      // reduce m1 by left shifting and modding until the exponents x1 and x2 are 
      // equal
      while (x1 != x2) {
        int s = Math.min(BitUtils.countLeadingZeros(m1) - 1, x1 - x2);
        x1 -= s;
        m1 = (m1 << s) % m2;
      }
    }
    return pack(n, x1, m1);
  }

  
  /////////////////////////////////////////////////////////////////////////////
  // Rounding
  /////////////////////////////////////////////////////////////////////////////

  
  /**
   * Returns the <code>double</code> of greatest magnitude (furthest from zero)
   * that is equal to a mathematical integer and which has a mignitude not
   * greater than the argument's magnitude.  Special cases:
   * <ul><li>If the argument value is already equal to a mathematical 
   * integer, then the result is the same as the argument. 
   * <li>If the argument is NaN or an infinity or positive zero or 
   * negative zero, then the result is the same as the argument.</ul>
   *
   * @param   d   a <code>double</code> value.
   * @return the <code>double</code> of greatest magnitude (furthest from zero)
   *         whose magnitude is not greater than the argument's and which 
   *         is equal to a mathematical integer.
   */
  public static long truncate(long d) {
    return round(d, false, unpackSign(d));
  }

  /**
   * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#rint(double)">Math.rint(double)</a>.
   */
  public static long rint(long d) {
    return round(d, true, false);
  }
  
  /**
   * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#floor(double)">Math.floor(double)</a>.
   */
  public static long floor(long d) {
    return round(d, false, false);
  }
  
  /**
   * Mimcs <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#ceil(double)">Math.ceil(double)</a>.
   */
  public static long ceil(long d) {
    return round(d, false, true);
  }

  /**
   * Mimcs <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#round(double)">Math.round(double)</a>.
   */
  public static long round(long d) {
    return longValue(floor(add(d, ONE_HALF)));
  }

  private static long round(long d, boolean round, boolean ceil) {
    if (isNaN(d)) {
      return NaN;
    } else if (isZero(d) || isInfinite(d)) {
      return d;
    }
    int x = unpackExponent(d);
    if (x >= 0) {
      return d;
    }
    boolean n = unpackSign(d);
    long m = unpackMantissa(d);
    if (round) {
      m = BitUtils.roundingRightShift(m, -x);
    } else {
      long r;
      if (x <= -64) {
        r = m;
        m = 0;
      } else {