microdouble.java

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

// $Id: MicroDouble.java,v 1.2 2004/08/03 04:57:42 Dave Exp $
/*
 * Double.java
 * Copyright (C) 2003, 2004 David Clausen
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 * Portions of this software are derived from FDLIBM, which contained the
 * following notice:
 *
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 *
 * For mor information on FDLIBM see:
 * http://netlib.bell-labs.com/netlib/fdlibm/index.html
 *
 */
package net.dclausen.microfloat;

import java.util.Random;

/**
 * A software implementation of IEEE-754 double precision math which does not
 * rely on the <code>double</code> data type. 
 * This class overloads the <code>long</code> data type by storing 
 * <code>double</code> data in it.
 * See the 
 * <a href="package-summary.html#package_description">package description</a> 
 * for more information.
 * <p>
 * @author David Clausen
 * @see <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Double.html">Double</a>
 * @see <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html">Math</a>
 * @see Float
 * @version $Revision: 1.2 $
 */
public class MicroDouble {
  
  /////////////////////////////////////////////////////////////////////////////
  // General-purpose constants
  /////////////////////////////////////////////////////////////////////////////

  /**
   * A constant holding the same value as <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Double.html#POSITIVE_INFINITY">Double.POSITIVE_INFINITY</a>
   */
  public  static final long POSITIVE_INFINITY = 0x7ff0000000000000L;
  
  /**
   * A constant holding the same value as <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Double.html#NEGATIVE_INFINITY">Double.NEGATIVE_INFINITY</a>
   */
  public  static final long NEGATIVE_INFINITY = 0xfff0000000000000L;

  /**
   * A constant holding the same value as <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Double.html#NaN">Double.NaN</a>
   */
  public  static final long NaN               = 0x7ff8000000000000L;
  
  /**
   * A constant holding the same value as <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Double.html#MAX_VALUE">Double.MAX_VALUE</a>
   */
  public  static final long MAX_VALUE         = 0x7fefffffffffffffL;
  
  /**
   * A constant holding the same value as <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Double.html#MIN_VALUE">Double.MIN_VALUE</a>
   */
  public  static final long MIN_VALUE         = 0x0000000000000001L; 
  
  /**
   * A constant holding the same value as <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#E">Math.E</a>
   */
  public  static final long E                 = 0x4005bf0a8b145769L;
  
  /**
   * A constant holding the same value as <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#PI">Math.PI</a>
   */
  public  static final long PI                = 0x400921fb54442d18L;

  // Other constants needed internally, and exposed as a convenience.

  /** A constant holding the value of 0.0d */
  public static final long ZERO              = 0x0000000000000000L;

  /** A constant holding the value of -0.0d */
  public static final long NEGATIVE_ZERO     = 0x8000000000000000L;
  
  /** A constant holding the value of 1.0d */
  public static final long ONE               = 0x3ff0000000000000L;
  
  /** A constant holding the value of -1.0d */
  public static final long NEGATIVE_ONE      = 0xbff0000000000000L;
  
  /** A constant holding the value of 2.0d */
  public static final long TWO               = 0x4000000000000000L;
  
  /** A constant holding the value of 3.0d */
  public static final long THREE             = 0x4008000000000000L;
  
  /** A constant holding the value of 4.0d */
  public static final long FOUR              = 0x4010000000000000L;
  
  /** A constant holding the value of 5.0d */
  public static final long FIVE              = 0x4014000000000000L;
  
  /** A constant holding the value of 6.0d */
  public static final long SIX               = 0x4018000000000000L;
  
  /** A constant holding the value of 8.0d */
  public static final long EIGHT             = 0x4020000000000000L;
  
  /** A constant holding the value of 10.0d */
  public static final long TEN               = 0x4024000000000000L;
  
  /** A constant holding the value of 100.0d */
  public static final long ONE_HUNDRED       = 0x4059000000000000L;
  
  /** A constant holding the value of 1.5d */
  public static final long THREE_HALVES      = 0x3ff8000000000000L;
  
  /** A constant holding the value of 0.5d */
  public static final long ONE_HALF          = 0x3fe0000000000000L;
  
  /** A constant holding the value of (1.0d / 3.0d) */
  public static final long ONE_THIRD         = 0x3fd5555555555555L;
  
  /** A constant holding the value of 0.25d */
  public static final long ONE_FOURTH        = 0x3fd0000000000000L;
  
  /** A constant holding the value of 0.125d */
  public static final long ONE_EIGHTH        = 0x3fc0000000000000L;
  
  /** A constant holding the natural logarithm of 2 */
  public static final long LN2               = 0x3fe62e42fefa39efL;

  
  /////////////////////////////////////////////////////////////////////////////
  // Packing and unpacking the IEEE-754 double precision format
  /////////////////////////////////////////////////////////////////////////////

  private static final long ABS_MASK          = 0x7fffffffffffffffL;
  private static final long SIGN_MASK         = 0x8000000000000000L; // 1 bit
  private static final long EXPONENT_MASK     = 0x7ff0000000000000L; // 11 bits
  private static final long FRACTION_MASK     = 0x000fffffffffffffL; // 52 bits
  private static final long IMPLIED_ONE       = 0x0010000000000000L; // 53rd bit

  /** @return true iff d is negative */
  static boolean unpackSign(long d) {
    return (d < 0L);
  }

  /** @return an integer in the range [-1075, 972] */
  static int unpackExponent(long d) {
    return (((int) (d >> 52)) & 0x7ff) - 1075;
  }

  /** @return a long in the range [0, 0x001fffffffffffffL] */
  static long unpackMantissa(long d) {
    if ((d & EXPONENT_MASK) == 0) {
      return ((d & FRACTION_MASK) << 1);
    } else {
      return ((d & FRACTION_MASK) | IMPLIED_ONE);
    }
  }

  /** 
   * @return the double which most closely represents the given base-2 mantissa
   *         and exponent
   */
  static long pack(boolean negative, int exponent, long mantissa) {
    // reduce precision of mantissa, rounding if necessary
    if (mantissa != 0) {
      // left align mantissa
      int shift = BitUtils.countLeadingZeros(mantissa);
      mantissa <<= shift;
      exponent -= shift;

      if (exponent < -1085) {
        // subnormal
        mantissa = BitUtils.roundingRightShift(mantissa, -1074 - exponent);
      } else {
        // normal
        mantissa = BitUtils.roundingRightShift(mantissa, 11);
        if (mantissa == 0x20000000000000L) {
          // oops, rounding carried into the 54th bit
          mantissa = 0x10000000000000L;
          exponent++;
        }
        // pack the exponent
        if (exponent > 960) {
          mantissa = POSITIVE_INFINITY;
        } else {
          mantissa ^= IMPLIED_ONE;
          mantissa |= ((long) (exponent + 1086)) << 52;
        }
      }
    }
    
    // pack the sign bit
    if (negative) {
      mantissa |= SIGN_MASK;
    }
    
    return mantissa;
  }

  
  /////////////////////////////////////////////////////////////////////////////
  // Simple tests 
  /////////////////////////////////////////////////////////////////////////////

  /**
   * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Double.html#isNaN(double)">Double.isNaN(double)</a>
   */
  public static boolean isNaN(long d) {
    return ((d & ABS_MASK) > POSITIVE_INFINITY);
  }

  /**
   * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Double.html#isInfinite(double)">Double.isInfinite(double)</a>
   */
  public static boolean isInfinite(long d) {
    return ((d & ABS_MASK) == POSITIVE_INFINITY);
  }
  
  /**
   * Returns <code>true</code> if the specified number has zero
   * magnitude, <code>false</code> otherwise.
   *
   * @param   d   the <code>double</code> value to be tested.
   * @return  <code>true</code> if the value of the argument is positive
   *          zero or negative zero; <code>false</code> otherwise.
   */
  public static boolean isZero(long d) {
    return ((d & ABS_MASK) == ZERO);
  }

  
  /////////////////////////////////////////////////////////////////////////////
  // Sign changes
  /////////////////////////////////////////////////////////////////////////////

  /**
   * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#abs(double)">Math.abs(double)</a>
   */
  public static long abs(long d) {
    //if (isNaN(d)) {
    //  return NaN;
    //}
    return (d & ABS_MASK);
  }

  /**
   * Returns the negation of a <code>double</code> value.
   * Special cases:
   * <ul>
   * <li>If the argument is negative zero, the result is positive zero.
   * <li>If the argument is positive zero, the result is negative zero.
   * <li>If the argument is negative infinity, the result is positive infinity.
   * <li>If the argument is positive infinity, the result is negative infinity.
   * <li>If the argument is NaN, the result is NaN.</ul>
   * <p>
   * This method takes the place of the unary <code>-</code> operator.
   *
   * @param   d   the <code>double</code> value whose negated value is to be 
   *              determined
   * @return  the negation of the argument.
   */
  public static long negate(long d) {
    if (isNaN(d)) {
      return NaN;
    }
    return (d ^ SIGN_MASK);
  }
  

  /////////////////////////////////////////////////////////////////////////////
  // Comparison
  /////////////////////////////////////////////////////////////////////////////

  /**
   * Returns <code>true</code> if the specified numbers are considered equal
   * according to 
   * <a href="http://java.sun.com/docs/books/jls/second_edition/html/expressions.doc.html#5198">section 15.21.1
   * of the JLS</a>.  Special cases:
   * <ul>
   * <li>If either operand is NaN, then the result is false
   * <li>Positive zero and negative zero are considered equal
   * </ul>
   * <p>
   * This method takes the place of the <code>==</code> operator.
   *
   * @param   d1   the first <code>double</code> value to be compared.
   * @param   d2   the second <code>double</code> value to be compared.
   * @return  <code>true</code> if the two values are considered equal;
   *          <code>false</code> otherwise.
   */
  public static boolean eq(long d1, long d2) {
    return (((d1 == d2) && (! isNaN(d1))) || (isZero(d1) && isZero(d2)));
  }

  /**
   * Returns <code>true</code> if the specified numbers are considered unequal
   * according to 
   * <a href="http://java.sun.com/docs/books/jls/second_edition/html/expressions.doc.html#5198">section
   * 15.21.1 of the JLS</a>.  Special cases:
   * <ul>
   * <li>If either operand is NaN, then the result is true
   * <li>Positive zero and negative zero are considered equal
   * </ul>
   * The value returned by <code>ne</code> is always the opposite of the value
   * returned by <code>eq</code> for the same arguments.
   * <p>
   * This method takes the place of the <code>!=</code> operator.
   *
   * @param   d1   the first <code>double</code> value to be compared.
   * @param   d2   the second <code>double</code> value to be compared.
   * @return  <code>true</code> if the two values are considered equal;
   *          <code>false</code> otherwise.
   */
  public static boolean ne(long d1, long d2) {
    return (! eq(d1, d2));
  }

  /**
   * Returns <code>true</code> if the first argument is considered less than
   * the second argument according to 
   * <a href="http://java.sun.com/docs/books/jls/second_edition/html/expressions.doc.html#153654">section
   * 15.20.1 of the JLS</a>.  Special cases:
   * <ul>
   * <li>If either operand is NaN, then the result is false
   * <li>Positive zero and negative zero are considered equal
   * <li>Negative infinity is conisdered less than all other values except NaN
   * <li>Positive infinity is conisdered greater than all other values except NaN
   * </ul>
   * <p>
   * This method takes the place of the <code>&lt;</code> operator.
   *
   * @param   d1   the first <code>double</code> value to be compared.