/*
 * Copyright (c) 2003, 2018 Oracle and/or its affiliates. All rights reserved.
 *
 * This program and the accompanying materials are made available under the
 * terms of the Eclipse Public License v. 2.0, which is available at
 * http://www.eclipse.org/legal/epl-2.0.
 *
 * This Source Code may also be made available under the following Secondary
 * Licenses when the conditions for such availability set forth in the
 * Eclipse Public License v. 2.0 are satisfied: GNU General Public License,
 * version 2 with the GNU Classpath Exception, which is available at
 * https://www.gnu.org/software/classpath/license.html.
 *
 * SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0
 */

package com.sun.ts.lib.util;

import com.sun.ts.lib.util.FloatConsts;
import com.sun.ts.lib.util.DoubleConsts;

/**
 * The class {@code FpUtils} contains static utility methods for manipulating
 * and inspecting {@code float} and {@code double} floating-point numbers. These
 * methods include functionality recommended or required by the IEEE 754
 * floating-point standard.
 *
 * @author Joseph D. Darcy
 */

public class FpUtils {
  /*
   * The methods in this class are reasonably implemented using direct or
   * indirect bit-level manipulation of floating-point values. However, having
   * access to the IEEE 754 recommended functions would obviate the need for
   * most programmers to engage in floating-point bit-twiddling.
   *
   * An IEEE 754 number has three fields, from most significant bit to to least
   * significant, sign, exponent, and significand.
   *
   * msb lsb [sign|exponent| fractional_significand]
   *
   * Using some encoding cleverness, explained below, the high order bit of the
   * logical significand does not need to be explicitly stored, thus
   * "fractional_significand" instead of simply "significand" in the figure
   * above.
   *
   * For finite normal numbers, the numerical value encoded is
   *
   * (-1)^sign * 2^(exponent)*(1.fractional_significand)
   *
   * Most finite floating-point numbers are normalized; the exponent value is
   * reduced until the leading significand bit is 1. Therefore, the leading 1 is
   * redundant and is not explicitly stored. If a numerical value is so small it
   * cannot be normalized, it has a subnormal representation. Subnormal numbers
   * don't have a leading 1 in their significand; subnormals are encoding using
   * a special exponent value. In other words, the high-order bit of the logical
   * significand can be elided in from the representation in either case since
   * the bit's value is implicit from the exponent value.
   *
   * The exponent field uses a biased representation; if the bits of the
   * exponent are interpreted as a unsigned integer E, the exponent represented
   * is E - E_bias where E_bias depends on the floating-point format. E can
   * range between E_min and E_max, constants which depend on the floating-point
   * format. E_min and E_max are -126 and +127 for float, -1022 and +1023 for
   * double.
   *
   * The 32-bit float format has 1 sign bit, 8 exponent bits, and 23 bits for
   * the significand (which is logically 24 bits wide because of the implicit
   * bit). The 64-bit double format has 1 sign bit, 11 exponent bits, and 52
   * bits for the significand (logically 53 bits).
   *
   * Subnormal numbers and zero have the special exponent value E_min -1; the
   * numerical value represented by a subnormal is:
   *
   * (-1)^sign * 2^(E_min)*(0.fractional_significand)
   *
   * Zero is represented by all zero bits in the exponent and all zero bits in
   * the significand; zero can have either sign.
   *
   * Infinity and NaN are encoded using the exponent value E_max + 1. Signed
   * infinities have all significand bits zero; NaNs have at least one non-zero
   * significand bit.
   *
   * The details of IEEE 754 floating-point encoding will be used in the methods
   * below without further comment. For further exposition on IEEE 754 numbers,
   * see "IEEE Standard for Binary Floating-Point Arithmetic" ANSI/IEEE Std
   * 754-1985 or William Kahan's "Lecture Notes on the Status of IEEE Standard
   * 754 for Binary Floating-Point Arithmetic",
   * http://www.cs.berkeley.edu/~wkahan/ieee754status/ieee754.ps.
   *
   * Many of this class's methods are members of the set of IEEE 754 recommended
   * functions or similar functions recommended or required by IEEE 754R.
   * Discussion of various implementation techniques for these functions have
   * occurred in:
   *
   * W.J. Cody and Jerome T. Coonen, "Algorithm 772 Functions to Support the
   * IEEE Standard for Binary Floating-Point Arithmetic," ACM Transactions on
   * Mathematical Software, vol. 19, no. 4, December 1993, pp. 443-451.
   *
   * Joseph D. Darcy, "Writing robust IEEE recommended functions in ``100% Pure
   * Java''(TM)," University of California, Berkeley technical report
   * UCB//CSD-98-1009.
   */

  /**
   * Don't let anyone instantiate this class.
   */
  private FpUtils() {
  }

  // Constants used in scalb
  static double twoToTheDoubleScaleUp = powerOfTwoD(512);

  static double twoToTheDoubleScaleDown = powerOfTwoD(-512);

  // Helper Methods

  // The following helper methods are used in the implementation of
  // the public recommended functions; they generally omit certain
  // tests for exception cases.

  /**
   * Returns unbiased exponent of a {@code double}.
   */
  public static int getExponent(double d) {
    /*
     * Bitwise convert d to long, mask out exponent bits, shift to the right and
     * then subtract out double's bias adjust to get true exponent value.
     */
    return (int) (((Double.doubleToRawLongBits(d)
        & DoubleConsts.EXP_BIT_MASK) >> (DoubleConsts.SIGNIFICAND_WIDTH - 1))
        - DoubleConsts.EXP_BIAS);
  }

  /**
   * Returns unbiased exponent of a {@code float}.
   */
  public static int getExponent(float f) {
    /*
     * Bitwise convert f to integer, mask out exponent bits, shift to the right
     * and then subtract out float's bias adjust to get true exponent value
     */
    return ((Float.floatToRawIntBits(f)
        & FloatConsts.EXP_BIT_MASK) >> (FloatConsts.SIGNIFICAND_WIDTH - 1))
        - FloatConsts.EXP_BIAS;
  }

  /**
   * Returns a floating-point power of two in the normal range.
   */
  static double powerOfTwoD(int n) {
    assert (n >= DoubleConsts.MIN_EXPONENT && n <= DoubleConsts.MAX_EXPONENT);
    return Double.longBitsToDouble((((long) n
        + (long) DoubleConsts.EXP_BIAS) << (DoubleConsts.SIGNIFICAND_WIDTH - 1))
        & DoubleConsts.EXP_BIT_MASK);
  }

  /**
   * Returns a floating-point power of two in the normal range.
   */
  static float powerOfTwoF(int n) {
    assert (n >= FloatConsts.MIN_EXPONENT && n <= FloatConsts.MAX_EXPONENT);
    return Float.intBitsToFloat(
        ((n + FloatConsts.EXP_BIAS) << (FloatConsts.SIGNIFICAND_WIDTH - 1))
            & FloatConsts.EXP_BIT_MASK);
  }

  /**
   * Returns the first floating-point argument with the sign of the second
   * floating-point argument. Note that unlike the
   * {@link FpUtils#copySign(double, double) copySign} method, this method does
   * not require NaN {@code sign} arguments to be treated as positive values;
   * implementations are permitted to treat some NaN arguments as positive and
   * other NaN arguments as negative to allow greater performance.
   *
   * @param magnitude
   *          the parameter providing the magnitude of the result
   * @param sign
   *          the parameter providing the sign of the result
   * @return a value with the magnitude of {@code magnitude} and the sign of
   *         {@code sign}.
   * @author Joseph D. Darcy
   */
  public static double rawCopySign(double magnitude, double sign) {
    return Double.longBitsToDouble(
        (Double.doubleToRawLongBits(sign) & (DoubleConsts.SIGN_BIT_MASK))
            | (Double.doubleToRawLongBits(magnitude)
                & (DoubleConsts.EXP_BIT_MASK | DoubleConsts.SIGNIF_BIT_MASK)));
  }

  /**
   * Returns the first floating-point argument with the sign of the second
   * floating-point argument. Note that unlike the
   * {@link FpUtils#copySign(float, float) copySign} method, this method does
   * not require NaN {@code sign} arguments to be treated as positive values;
   * implementations are permitted to treat some NaN arguments as positive and
   * other NaN arguments as negative to allow greater performance.
   *
   * @param magnitude
   *          the parameter providing the magnitude of the result
   * @param sign
   *          the parameter providing the sign of the result
   * @return a value with the magnitude of {@code magnitude} and the sign of
   *         {@code sign}.
   * @author Joseph D. Darcy
   */
  public static float rawCopySign(float magnitude, float sign) {
    return Float.intBitsToFloat(
        (Float.floatToRawIntBits(sign) & (FloatConsts.SIGN_BIT_MASK))
            | (Float.floatToRawIntBits(magnitude)
                & (FloatConsts.EXP_BIT_MASK | FloatConsts.SIGNIF_BIT_MASK)));
  }

  /* ***************************************************************** */

  /**
   * Returns {@code true} if the argument is a finite floating-point value;
   * returns {@code false} otherwise (for NaN and infinity arguments).
   *
   * @param d
   *          the {@code double} value to be tested
   * @return {@code true} if the argument is a finite floating-point value,
   *         {@code false} otherwise.
   */
  public static boolean isFinite(double d) {
    return Math.abs(d) <= DoubleConsts.MAX_VALUE;
  }

  /**
   * Returns {@code true} if the argument is a finite floating-point value;
   * returns {@code false} otherwise (for NaN and infinity arguments).
   *
   * @param f
   *          the {@code float} value to be tested
   * @return {@code true} if the argument is a finite floating-point value,
   *         {@code false} otherwise.
   */
  public static boolean isFinite(float f) {
    return Math.abs(f) <= FloatConsts.MAX_VALUE;
  }

  /**
   * Returns {@code true} if the specified number is infinitely large in
   * magnitude, {@code false} otherwise.
   *
   * <p>
   * Note that this method is equivalent to the {@link Double#isInfinite(double)
   * Double.isInfinite} method; the functionality is included in this class for
   * convenience.
   *
   * @param d
   *          the value to be tested.
   * @return {@code true} if the value of the argument is positive infinity or
   *         negative infinity; {@code false} otherwise.
   */
  public static boolean isInfinite(double d) {
    return Double.isInfinite(d);
  }

  /**
   * Returns {@code true} if the specified number is infinitely large in
   * magnitude, {@code false} otherwise.
   *
   * <p>
   * Note that this method is equivalent to the {@link Float#isInfinite(float)
   * Float.isInfinite} method; the functionality is included in this class for
   * convenience.
   *
   * @param f
   *          the value to be tested.
   * @return {@code true} if the argument is positive infinity or negative
   *         infinity; {@code false} otherwise.
   */
  public static boolean isInfinite(float f) {
    return Float.isInfinite(f);
  }

  /**
   * Returns {@code true} if the specified number is a Not-a-Number (NaN) value,
   * {@code false} otherwise.
   *
   * <p>
   * Note that this method is equivalent to the {@link Double#isNaN(double)
   * Double.isNaN} method; the functionality is included in this class for
   * convenience.
   *
   * @param d
   *          the value to be tested.
   * @return {@code true} if the value of the argument is NaN; {@code false}
   *         otherwise.
   */
  public static boolean isNaN(double d) {
    return Double.isNaN(d);
  }

  /**
   * Returns {@code true} if the specified number is a Not-a-Number (NaN) value,
   * {@code false} otherwise.
   *
   * <p>
   * Note that this method is equivalent to the {@link Float#isNaN(float)
   * Float.isNaN} method; the functionality is included in this class for
   * convenience.
   *
   * @param f
   *          the value to be tested.
   * @return {@code true} if the argument is NaN; {@code false} otherwise.
   */
  public static boolean isNaN(float f) {
    return Float.isNaN(f);
  }

  /**
   * Returns {@code true} if the unordered relation holds between the two
   * arguments. When two floating-point values are unordered, one value is
   * neither less than, equal to, nor greater than the other. For the unordered
   * relation to be true, at least one argument must be a {@code NaN}.
   *
   * @param arg1
   *          the first argument
   * @param arg2
   *          the second argument
   * @return {@code true} if at least one argument is a NaN, {@code false}
   *         otherwise.
   */
  public static boolean isUnordered(double arg1, double arg2) {
    return isNaN(arg1) || isNaN(arg2);
  }

  /**
   * Returns {@code true} if the unordered relation holds between the two
   * arguments. When two floating-point values are unordered, one value is
   * neither less than, equal to, nor greater than the other. For the unordered
   * relation to be true, at least one argument must be a {@code NaN}.
   *
   * @param arg1
   *          the first argument
   * @param arg2
   *          the second argument
   * @return {@code true} if at least one argument is a NaN, {@code false}
   *         otherwise.
   */
  public static boolean isUnordered(float arg1, float arg2) {
    return isNaN(arg1) || isNaN(arg2);
  }

  /**
   * Returns unbiased exponent of a {@code double}; for subnormal values, the
   * number is treated as if it were normalized. That is for all finite,
   * non-zero, positive numbers <i>x</i>,
   * <code>scalb(<i>x</i>, -ilogb(<i>x</i>))</code> is always in the range [1,
   * 2).
   * <p>
   * Special cases:
   * <ul>
   * <li>If the argument is NaN, then the result is 2<sup>30</sup>.
   * <li>If the argument is infinite, then the result is 2<sup>28</sup>.
   * <li>If the argument is zero, then the result is -(2<sup>28</sup>).
   * </ul>
   *
   * @param d
   *          floating-point number whose exponent is to be extracted
   * @return unbiased exponent of the argument.
   * @author Joseph D. Darcy
   */
  public static int ilogb(double d) {
    int exponent = getExponent(d);

    switch (exponent) {
    case DoubleConsts.MAX_EXPONENT + 1: // NaN or infinity
      if (isNaN(d))
        return (1 << 30); // 2^30
      else // infinite value
        return (1 << 28); // 2^28

    case DoubleConsts.MIN_EXPONENT - 1: // zero or subnormal
      if (d == 0.0) {
        return -(1 << 28); // -(2^28)
      } else {
        long transducer = Double.doubleToRawLongBits(d);

        /*
         * To avoid causing slow arithmetic on subnormals, the scaling to
         * determine when d's significand is normalized is done in integer
         * arithmetic. (there must be at least one "1" bit in the significand
         * since zero has been screened out.
         */

        // isolate significand bits
        transducer &= DoubleConsts.SIGNIF_BIT_MASK;
        assert (transducer != 0L);

        // This loop is simple and functional. We might be
        // able to do something more clever that was faster;
        // e.g. number of leading zero detection on
        // (transducer << (# exponent and sign bits).
        while (transducer < (1L << (DoubleConsts.SIGNIFICAND_WIDTH - 1))) {
          transducer *= 2;
          exponent--;
        }
        exponent++;
        assert (exponent >= DoubleConsts.MIN_EXPONENT
            - (DoubleConsts.SIGNIFICAND_WIDTH - 1)
            && exponent < DoubleConsts.MIN_EXPONENT);
        return exponent;
      }

    default:
      assert (exponent >= DoubleConsts.MIN_EXPONENT
          && exponent <= DoubleConsts.MAX_EXPONENT);
      return exponent;
    }
  }

  /**
   * Returns unbiased exponent of a {@code float}; for subnormal values, the
   * number is treated as if it were normalized. That is for all finite,
   * non-zero, positive numbers <i>x</i>,
   * <code>scalb(<i>x</i>, -ilogb(<i>x</i>))</code> is always in the range [1,
   * 2).
   * <p>
   * Special cases:
   * <ul>
   * <li>If the argument is NaN, then the result is 2<sup>30</sup>.
   * <li>If the argument is infinite, then the result is 2<sup>28</sup>.
   * <li>If the argument is zero, then the result is -(2<sup>28</sup>).
   * </ul>
   *
   * @param f
   *          floating-point number whose exponent is to be extracted
   * @return unbiased exponent of the argument.
   * @author Joseph D. Darcy
   */
  public static int ilogb(float f) {
    int exponent = getExponent(f);

    switch (exponent) {
    case FloatConsts.MAX_EXPONENT + 1: // NaN or infinity
      if (isNaN(f))
        return (1 << 30); // 2^30
      else // infinite value
        return (1 << 28); // 2^28

    case FloatConsts.MIN_EXPONENT - 1: // zero or subnormal
      if (f == 0.0f) {
        return -(1 << 28); // -(2^28)
      } else {
        int transducer = Float.floatToRawIntBits(f);

        /*
         * To avoid causing slow arithmetic on subnormals, the scaling to
         * determine when f's significand is normalized is done in integer
         * arithmetic. (there must be at least one "1" bit in the significand
         * since zero has been screened out.
         */

        // isolate significand bits
        transducer &= FloatConsts.SIGNIF_BIT_MASK;
        assert (transducer != 0);

        // This loop is simple and functional. We might be
        // able to do something more clever that was faster;
        // e.g. number of leading zero detection on
        // (transducer << (# exponent and sign bits).
        while (transducer < (1 << (FloatConsts.SIGNIFICAND_WIDTH - 1))) {
          transducer *= 2;
          exponent--;
        }
        exponent++;
        assert (exponent >= FloatConsts.MIN_EXPONENT
            - (FloatConsts.SIGNIFICAND_WIDTH - 1)
            && exponent < FloatConsts.MIN_EXPONENT);
        return exponent;
      }

    default:
      assert (exponent >= FloatConsts.MIN_EXPONENT
          && exponent <= FloatConsts.MAX_EXPONENT);
      return exponent;
    }
  }

  /*
   * The scalb operation should be reasonably fast; however, there are tradeoffs
   * in writing a method to minimize the worst case performance and writing a
   * method to minimize the time for expected common inputs. Some processors
   * operate very slowly on subnormal operands, taking hundreds or thousands of
   * cycles for one floating-point add or multiply as opposed to, say, four
   * cycles for normal operands. For processors with very slow subnormal
   * execution, scalb would be fastest if written entirely with integer
   * operations; in other words, scalb would need to include the logic of
   * performing correct rounding of subnormal values. This could be reasonably
   * done in at most a few hundred cycles. However, this approach may penalize
   * normal operations since at least the exponent of the floating-point
   * argument must be examined.
   *
   * The approach taken in this implementation is a compromise. Floating-point
   * multiplication is used to do most of the work; but knowingly multiplying by
   * a subnormal scaling factor is avoided. However, the floating-point argument
   * is not examined to see whether or not it is subnormal since subnormal
   * inputs are assumed to be rare. At most three multiplies are needed to scale
   * from the largest to smallest exponent ranges (scaling down, at most two
   * multiplies are needed if subnormal scaling factors are allowed). However,
   * in this implementation an expensive integer remainder operation is avoided
   * at the cost of requiring five floating-point multiplies in the worst case,
   * which should still be a performance win.
   *
   * If scaling of entire arrays is a concern, it would probably be more
   * efficient to provide a double[] scalb(double[], int) version of scalb to
   * avoid having to recompute the needed scaling factors for each
   * floating-point value.
   */

  /**
   * Return {@code d} &times; 2<sup>{@code scale_factor}</sup> rounded as if
   * performed by a single correctly rounded floating-point multiply to a member
   * of the double value set. See section 4.2.3 of <cite>The Java&trade;
   * Language Specification</cite> for a discussion of floating-point value
   * sets. If the exponent of the result is between the {@code double}'s minimum
   * exponent and maximum exponent, the answer is calculated exactly. If the
   * exponent of the result would be larger than {@code doubles}'s maximum
   * exponent, an infinity is returned. Note that if the result is subnormal,
   * precision may be lost; that is, when {@code scalb(x,
   * n)} is subnormal, {@code scalb(scalb(x, n), -n)} may not equal <i>x</i>.
   * When the result is non-NaN, the result has the same sign as {@code d}.
   *
   * <p>
   * Special cases:
   * <ul>
   * <li>If the first argument is NaN, NaN is returned.
   * <li>If the first argument is infinite, then an infinity of the same sign is
   * returned.
   * <li>If the first argument is zero, then a zero of the same sign is
   * returned.
   * </ul>
   *
   * @param d
   *          number to be scaled by a power of two.
   * @param scale_factor
   *          power of 2 used to scale {@code d}
   * @return {@code d * }2<sup>{@code scale_factor}</sup>
   * @author Joseph D. Darcy
   */
  public static double scalb(double d, int scale_factor) {
    /*
     * This method does not need to be declared strictfp to compute the same
     * correct result on all platforms. When scaling up, it does not matter what
     * order the multiply-store operations are done; the result will be finite
     * or overflow regardless of the operation ordering. However, to get the
     * correct result when scaling down, a particular ordering must be used.
     *
     * When scaling down, the multiply-store operations are sequenced so that it
     * is not possible for two consecutive multiply-stores to return subnormal
     * results. If one multiply-store result is subnormal, the next multiply
     * will round it away to zero. This is done by first multiplying by 2 ^
     * (scale_factor % n) and then multiplying several times by by 2^n as needed
     * where n is the exponent of number that is a covenient power of two. In
     * this way, at most one real rounding error occurs. If the double value set
     * is being used exclusively, the rounding will occur on a multiply. If the
     * double-extended-exponent value set is being used, the products will
     * (perhaps) be exact but the stores to d are guaranteed to round to the
     * double value set.
     *
     * It is _not_ a valid implementation to first multiply d by 2^MIN_EXPONENT
     * and then by 2 ^ (scale_factor % MIN_EXPONENT) since even in a strictfp
     * program double rounding on underflow could occur; e.g. if the
     * scale_factor argument was (MIN_EXPONENT - n) and the exponent of d was a
     * little less than -(MIN_EXPONENT - n), meaning the final result would be
     * subnormal.
     *
     * Since exact reproducibility of this method can be achieved without any
     * undue performance burden, there is no compelling reason to allow double
     * rounding on underflow in scalb.
     */

    // magnitude of a power of two so large that scaling a finite
    // nonzero value by it would be guaranteed to over or
    // underflow; due to rounding, scaling down takes takes an
    // additional power of two which is reflected here
    final int MAX_SCALE = DoubleConsts.MAX_EXPONENT + -DoubleConsts.MIN_EXPONENT
        + DoubleConsts.SIGNIFICAND_WIDTH + 1;
    int exp_adjust = 0;
    int scale_increment = 0;
    double exp_delta = Double.NaN;

    // Make sure scaling factor is in a reasonable range

    if (scale_factor < 0) {
      scale_factor = Math.max(scale_factor, -MAX_SCALE);
      scale_increment = -512;
      exp_delta = twoToTheDoubleScaleDown;
    } else {
      scale_factor = Math.min(scale_factor, MAX_SCALE);
      scale_increment = 512;
      exp_delta = twoToTheDoubleScaleUp;
    }

    // Calculate (scale_factor % +/-512), 512 = 2^9, using
    // technique from "Hacker's Delight" section 10-2.
    int t = (scale_factor >> 9 - 1) >>> 32 - 9;
    exp_adjust = ((scale_factor + t) & (512 - 1)) - t;

    d *= powerOfTwoD(exp_adjust);
    scale_factor -= exp_adjust;

    while (scale_factor != 0) {
      d *= exp_delta;
      scale_factor -= scale_increment;
    }
    return d;
  }

  /**
   * Return {@code f} &times; 2<sup>{@code scale_factor}</sup> rounded as if
   * performed by a single correctly rounded floating-point multiply to a member
   * of the float value set. See section 4.2.3 of <cite>The Java&trade; Language
   * Specification</cite> for a discussion of floating-point value sets. If the
   * exponent of the result is between the {@code float}'s minimum exponent and
   * maximum exponent, the answer is calculated exactly. If the exponent of the
   * result would be larger than {@code float}'s maximum exponent, an infinity
   * is returned. Note that if the result is subnormal, precision may be lost;
   * that is, when {@code scalb(x, n)} is subnormal,
   * {@code scalb(scalb(x, n), -n)} may not equal <i>x</i>. When the result is
   * non-NaN, the result has the same sign as {@code f}.
   *
   * <p>
   * Special cases:
   * <ul>
   * <li>If the first argument is NaN, NaN is returned.
   * <li>If the first argument is infinite, then an infinity of the same sign is
   * returned.
   * <li>If the first argument is zero, then a zero of the same sign is
   * returned.
   * </ul>
   *
   * @param f
   *          number to be scaled by a power of two.
   * @param scale_factor
   *          power of 2 used to scale {@code f}
   * @return {@code f * }2<sup>{@code scale_factor}</sup>
   * @author Joseph D. Darcy
   */
  public static float scalb(float f, int scale_factor) {
    // magnitude of a power of two so large that scaling a finite
    // nonzero value by it would be guaranteed to over or
    // underflow; due to rounding, scaling down takes takes an
    // additional power of two which is reflected here
    final int MAX_SCALE = FloatConsts.MAX_EXPONENT + -FloatConsts.MIN_EXPONENT
        + FloatConsts.SIGNIFICAND_WIDTH + 1;

    // Make sure scaling factor is in a reasonable range
    scale_factor = Math.max(Math.min(scale_factor, MAX_SCALE), -MAX_SCALE);

    /*
     * Since + MAX_SCALE for float fits well within the double exponent range
     * and + float -> double conversion is exact the multiplication below will
     * be exact. Therefore, the rounding that occurs when the double product is
     * cast to float will be the correctly rounded float result. Since all
     * operations other than the final multiply will be exact, it is not
     * necessary to declare this method strictfp.
     */
    return (float) ((double) f * powerOfTwoD(scale_factor));
  }

  /**
   * Returns the floating-point number adjacent to the first argument in the
   * direction of the second argument. If both arguments compare as equal the
   * second argument is returned.
   *
   * <p>
   * Special cases:
   * <ul>
   * <li>If either argument is a NaN, then NaN is returned.
   *
   * <li>If both arguments are signed zeros, {@code direction} is returned
   * unchanged (as implied by the requirement of returning the second argument
   * if the arguments compare as equal).
   *
   * <li>If {@code start} is &plusmn;{@code Double.MIN_VALUE} and
   * {@code direction} has a value such that the result should have a smaller
   * magnitude, then a zero with the same sign as {@code start} is returned.
   *
   * <li>If {@code start} is infinite and {@code direction} has a value such
   * that the result should have a smaller magnitude, {@code Double.MAX_VALUE}
   * with the same sign as {@code start} is returned.
   *
   * <li>If {@code start} is equal to &plusmn; {@code Double.MAX_VALUE} and
   * {@code direction} has a value such that the result should have a larger
   * magnitude, an infinity with same sign as {@code start} is returned.
   * </ul>
   *
   * @param start
   *          starting floating-point value
   * @param direction
   *          value indicating which of {@code start}'s neighbors or
   *          {@code start} should be returned
   * @return The floating-point number adjacent to {@code start} in the
   *         direction of {@code direction}.
   * @author Joseph D. Darcy
   */
  public static double nextAfter(double start, double direction) {
    /*
     * The cases:
     *
     * nextAfter(+infinity, 0) == MAX_VALUE nextAfter(+infinity, +infinity) ==
     * +infinity nextAfter(-infinity, 0) == -MAX_VALUE nextAfter(-infinity,
     * -infinity) == -infinity
     *
     * are naturally handled without any additional testing
     */

    // First check for NaN values
    if (isNaN(start) || isNaN(direction)) {
      // return a NaN derived from the input NaN(s)
      return start + direction;
    } else if (start == direction) {
      return direction;
    } else { // start > direction or start < direction
      // Add +0.0 to get rid of a -0.0 (+0.0 + -0.0 => +0.0)
      // then bitwise convert start to integer.
      long transducer = Double.doubleToRawLongBits(start + 0.0d);

      /*
       * IEEE 754 floating-point numbers are lexicographically ordered if
       * treated as signed- magnitude integers . Since Java's integers are two's
       * complement, incrementing" the two's complement representation of a
       * logically negative floating-point value *decrements* the
       * signed-magnitude representation. Therefore, when the integer
       * representation of a floating-point values is less than zero, the
       * adjustment to the representation is in the opposite direction than
       * would be expected at first .
       */
      if (direction > start) { // Calculate next greater value
        transducer = transducer + (transducer >= 0L ? 1L : -1L);
      } else { // Calculate next lesser value
        assert direction < start;
        if (transducer > 0L)
          --transducer;
        else if (transducer < 0L)
          ++transducer;
        /*
         * transducer==0, the result is -MIN_VALUE
         *
         * The transition from zero (implicitly positive) to the smallest
         * negative signed magnitude value must be done explicitly.
         */
        else
          transducer = DoubleConsts.SIGN_BIT_MASK | 1L;
      }

      return Double.longBitsToDouble(transducer);
    }
  }

  /**
   * Returns the floating-point number adjacent to the first argument in the
   * direction of the second argument. If both arguments compare as equal, the
   * second argument is returned.
   *
   * <p>
   * Special cases:
   * <ul>
   * <li>If either argument is a NaN, then NaN is returned.
   *
   * <li>If both arguments are signed zeros, a {@code float} zero with the same
   * sign as {@code direction} is returned (as implied by the requirement of
   * returning the second argument if the arguments compare as equal).
   *
   * <li>If {@code start} is &plusmn;{@code Float.MIN_VALUE} and
   * {@code direction} has a value such that the result should have a smaller
   * magnitude, then a zero with the same sign as {@code start} is returned.
   *
   * <li>If {@code start} is infinite and {@code direction} has a value such
   * that the result should have a smaller magnitude, {@code Float.MAX_VALUE}
   * with the same sign as {@code start} is returned.
   *
   * <li>If {@code start} is equal to &plusmn; {@code Float.MAX_VALUE} and
   * {@code direction} has a value such that the result should have a larger
   * magnitude, an infinity with same sign as {@code start} is returned.
   * </ul>
   *
   * @param start
   *          starting floating-point value
   * @param direction
   *          value indicating which of {@code start}'s neighbors or
   *          {@code start} should be returned
   * @return The floating-point number adjacent to {@code start} in the
   *         direction of {@code direction}.
   * @author Joseph D. Darcy
   */
  public static float nextAfter(float start, double direction) {
    /*
     * The cases:
     *
     * nextAfter(+infinity, 0) == MAX_VALUE nextAfter(+infinity, +infinity) ==
     * +infinity nextAfter(-infinity, 0) == -MAX_VALUE nextAfter(-infinity,
     * -infinity) == -infinity
     *
     * are naturally handled without any additional testing
     */

    // First check for NaN values
    if (isNaN(start) || isNaN(direction)) {
      // return a NaN derived from the input NaN(s)
      return start + (float) direction;
    } else if (start == direction) {
      return (float) direction;
    } else { // start > direction or start < direction
      // Add +0.0 to get rid of a -0.0 (+0.0 + -0.0 => +0.0)
      // then bitwise convert start to integer.
      int transducer = Float.floatToRawIntBits(start + 0.0f);

      /*
       * IEEE 754 floating-point numbers are lexicographically ordered if
       * treated as signed- magnitude integers . Since Java's integers are two's
       * complement, incrementing" the two's complement representation of a
       * logically negative floating-point value *decrements* the
       * signed-magnitude representation. Therefore, when the integer
       * representation of a floating-point values is less than zero, the
       * adjustment to the representation is in the opposite direction than
       * would be expected at first.
       */
      if (direction > start) {// Calculate next greater value
        transducer = transducer + (transducer >= 0 ? 1 : -1);
      } else { // Calculate next lesser value
        assert direction < start;
        if (transducer > 0)
          --transducer;
        else if (transducer < 0)
          ++transducer;
        /*
         * transducer==0, the result is -MIN_VALUE
         *
         * The transition from zero (implicitly positive) to the smallest
         * negative signed magnitude value must be done explicitly.
         */
        else
          transducer = FloatConsts.SIGN_BIT_MASK | 1;
      }

      return Float.intBitsToFloat(transducer);
    }
  }

  /**
   * Returns the floating-point value adjacent to {@code d} in the direction of
   * positive infinity. This method is semantically equivalent to
   * {@code nextAfter(d,
   * Double.POSITIVE_INFINITY)}; however, a {@code nextUp} implementation may
   * run faster than its equivalent {@code nextAfter} call.
   *
   * <p>
   * Special Cases:
   * <ul>
   * <li>If the argument is NaN, the result is NaN.
   *
   * <li>If the argument is positive infinity, the result is positive infinity.
   *
   * <li>If the argument is zero, the result is {@code Double.MIN_VALUE}
   *
   * </ul>
   *
   * @param d
   *          starting floating-point value
   * @return The adjacent floating-point value closer to positive infinity.
   * @author Joseph D. Darcy
   */
  public static double nextUp(double d) {
    if (isNaN(d) || d == Double.POSITIVE_INFINITY)
      return d;
    else {
      d += 0.0d;
      return Double.longBitsToDouble(
          Double.doubleToRawLongBits(d) + ((d >= 0.0d) ? +1L : -1L));
    }
  }

  /**
   * Returns the floating-point value adjacent to {@code f} in the direction of
   * positive infinity. This method is semantically equivalent to
   * {@code nextAfter(f,
   * Double.POSITIVE_INFINITY)}; however, a {@code nextUp} implementation may
   * run faster than its equivalent {@code nextAfter} call.
   *
   * <p>
   * Special Cases:
   * <ul>
   * <li>If the argument is NaN, the result is NaN.
   *
   * <li>If the argument is positive infinity, the result is positive infinity.
   *
   * <li>If the argument is zero, the result is {@code Float.MIN_VALUE}
   *
   * </ul>
   *
   * @param f
   *          starting floating-point value
   * @return The adjacent floating-point value closer to positive infinity.
   * @author Joseph D. Darcy
   */
  public static float nextUp(float f) {
    if (isNaN(f) || f == FloatConsts.POSITIVE_INFINITY)
      return f;
    else {
      f += 0.0f;
      return Float
          .intBitsToFloat(Float.floatToRawIntBits(f) + ((f >= 0.0f) ? +1 : -1));
    }
  }

  /**
   * Returns the floating-point value adjacent to {@code d} in the direction of
   * negative infinity. This method is semantically equivalent to
   * {@code nextAfter(d,
   * Double.NEGATIVE_INFINITY)}; however, a {@code nextDown} implementation may
   * run faster than its equivalent {@code nextAfter} call.
   *
   * <p>
   * Special Cases:
   * <ul>
   * <li>If the argument is NaN, the result is NaN.
   *
   * <li>If the argument is negative infinity, the result is negative infinity.
   *
   * <li>If the argument is zero, the result is {@code -Double.MIN_VALUE}
   *
   * </ul>
   *
   * @param d
   *          starting floating-point value
   * @return The adjacent floating-point value closer to negative infinity.
   * @author Joseph D. Darcy
   */
  public static double nextDown(double d) {
    if (isNaN(d) || d == Double.NEGATIVE_INFINITY)
      return d;
    else {
      if (d == 0.0)
        return -Double.MIN_VALUE;
      else
        return Double.longBitsToDouble(
            Double.doubleToRawLongBits(d) + ((d > 0.0d) ? -1L : +1L));
    }
  }

  /**
   * Returns the floating-point value adjacent to {@code f} in the direction of
   * negative infinity. This method is semantically equivalent to
   * {@code nextAfter(f,
   * Float.NEGATIVE_INFINITY)}; however, a {@code nextDown} implementation may
   * run faster than its equivalent {@code nextAfter} call.
   *
   * <p>
   * Special Cases:
   * <ul>
   * <li>If the argument is NaN, the result is NaN.
   *
   * <li>If the argument is negative infinity, the result is negative infinity.
   *
   * <li>If the argument is zero, the result is {@code -Float.MIN_VALUE}
   *
   * </ul>
   *
   * @param f
   *          starting floating-point value
   * @return The adjacent floating-point value closer to negative infinity.
   * @author Joseph D. Darcy
   */
  public static double nextDown(float f) {
    if (isNaN(f) || f == Float.NEGATIVE_INFINITY)
      return f;
    else {
      if (f == 0.0f)
        return -Float.MIN_VALUE;
      else
        return Float.intBitsToFloat(
            Float.floatToRawIntBits(f) + ((f > 0.0f) ? -1 : +1));
    }
  }

  /**
   * Returns the first floating-point argument with the sign of the second
   * floating-point argument. For this method, a NaN {@code sign} argument is
   * always treated as if it were positive.
   *
   * @param magnitude
   *          the parameter providing the magnitude of the result
   * @param sign
   *          the parameter providing the sign of the result
   * @return a value with the magnitude of {@code magnitude} and the sign of
   *         {@code sign}.
   * @author Joseph D. Darcy
   * @since 1.5
   */
  public static double copySign(double magnitude, double sign) {
    return rawCopySign(magnitude, (isNaN(sign) ? 1.0d : sign));
  }

  /**
   * Returns the first floating-point argument with the sign of the second
   * floating-point argument. For this method, a NaN {@code sign} argument is
   * always treated as if it were positive.
   *
   * @param magnitude
   *          the parameter providing the magnitude of the result
   * @param sign
   *          the parameter providing the sign of the result
   * @return a value with the magnitude of {@code magnitude} and the sign of
   *         {@code sign}.
   * @author Joseph D. Darcy
   */
  public static float copySign(float magnitude, float sign) {
    return rawCopySign(magnitude, (isNaN(sign) ? 1.0f : sign));
  }

  /**
   * Returns the size of an ulp of the argument. An ulp of a {@code double}
   * value is the positive distance between this floating-point value and the
   * {@code double} value next larger in magnitude. Note that for non-NaN
   * <i>x</i>, <code>ulp(-<i>x</i>) == ulp(<i>x</i>)</code>.
   *
   * <p>
   * Special Cases:
   * <ul>
   * <li>If the argument is NaN, then the result is NaN.
   * <li>If the argument is positive or negative infinity, then the result is
   * positive infinity.
   * <li>If the argument is positive or negative zero, then the result is
   * {@code Double.MIN_VALUE}.
   * <li>If the argument is &plusmn;{@code Double.MAX_VALUE}, then the result is
   * equal to 2<sup>971</sup>.
   * </ul>
   *
   * @param d
   *          the floating-point value whose ulp is to be returned
   * @return the size of an ulp of the argument
   * @author Joseph D. Darcy
   * @since 1.5
   */
  public static double ulp(double d) {
    int exp = getExponent(d);

    switch (exp) {
    case DoubleConsts.MAX_EXPONENT + 1: // NaN or infinity
      return Math.abs(d);

    case DoubleConsts.MIN_EXPONENT - 1: // zero or subnormal
      return Double.MIN_VALUE;

    default:
      assert exp <= DoubleConsts.MAX_EXPONENT
          && exp >= DoubleConsts.MIN_EXPONENT;

      // ulp(x) is usually 2^(SIGNIFICAND_WIDTH-1)*(2^ilogb(x))
      exp = exp - (DoubleConsts.SIGNIFICAND_WIDTH - 1);
      if (exp >= DoubleConsts.MIN_EXPONENT) {
        return powerOfTwoD(exp);
      } else {
        // return a subnormal result; left shift integer
        // representation of Double.MIN_VALUE appropriate
        // number of positions
        return Double.longBitsToDouble(1L << (exp - (DoubleConsts.MIN_EXPONENT
            - (DoubleConsts.SIGNIFICAND_WIDTH - 1))));
      }
    }
  }

  /**
   * Returns the size of an ulp of the argument. An ulp of a {@code float} value
   * is the positive distance between this floating-point value and the
   * {@code float} value next larger in magnitude. Note that for non-NaN
   * <i>x</i>, <code>ulp(-<i>x</i>) == ulp(<i>x</i>)</code>.
   *
   * <p>
   * Special Cases:
   * <ul>
   * <li>If the argument is NaN, then the result is NaN.
   * <li>If the argument is positive or negative infinity, then the result is
   * positive infinity.
   * <li>If the argument is positive or negative zero, then the result is
   * {@code Float.MIN_VALUE}.
   * <li>If the argument is &plusmn;{@code Float.MAX_VALUE}, then the result is
   * equal to 2<sup>104</sup>.
   * </ul>
   *
   * @param f
   *          the floating-point value whose ulp is to be returned
   * @return the size of an ulp of the argument
   * @author Joseph D. Darcy
   * @since 1.5
   */
  public static float ulp(float f) {
    int exp = getExponent(f);

    switch (exp) {
    case FloatConsts.MAX_EXPONENT + 1: // NaN or infinity
      return Math.abs(f);

    case FloatConsts.MIN_EXPONENT - 1: // zero or subnormal
      return FloatConsts.MIN_VALUE;

    default:
      assert exp <= FloatConsts.MAX_EXPONENT && exp >= FloatConsts.MIN_EXPONENT;

      // ulp(x) is usually 2^(SIGNIFICAND_WIDTH-1)*(2^ilogb(x))
      exp = exp - (FloatConsts.SIGNIFICAND_WIDTH - 1);
      if (exp >= FloatConsts.MIN_EXPONENT) {
        return powerOfTwoF(exp);
      } else {
        // return a subnormal result; left shift integer
        // representation of FloatConsts.MIN_VALUE appropriate
        // number of positions
        return Float.intBitsToFloat(1 << (exp - (FloatConsts.MIN_EXPONENT
            - (FloatConsts.SIGNIFICAND_WIDTH - 1))));
      }
    }
  }

  /**
   * Returns the signum function of the argument; zero if the argument is zero,
   * 1.0 if the argument is greater than zero, -1.0 if the argument is less than
   * zero.
   *
   * <p>
   * Special Cases:
   * <ul>
   * <li>If the argument is NaN, then the result is NaN.
   * <li>If the argument is positive zero or negative zero, then the result is
   * the same as the argument.
   * </ul>
   *
   * @param d
   *          the floating-point value whose signum is to be returned
   * @return the signum function of the argument
   * @author Joseph D. Darcy
   * @since 1.5
   */
  public static double signum(double d) {
    return (d == 0.0 || isNaN(d)) ? d : copySign(1.0, d);
  }

  /**
   * Returns the signum function of the argument; zero if the argument is zero,
   * 1.0f if the argument is greater than zero, -1.0f if the argument is less
   * than zero.
   *
   * <p>
   * Special Cases:
   * <ul>
   * <li>If the argument is NaN, then the result is NaN.
   * <li>If the argument is positive zero or negative zero, then the result is
   * the same as the argument.
   * </ul>
   *
   * @param f
   *          the floating-point value whose signum is to be returned
   * @return the signum function of the argument
   * @author Joseph D. Darcy
   * @since 1.5
   */
  public static float signum(float f) {
    return (f == 0.0f || isNaN(f)) ? f : copySign(1.0f, f);
  }

}
