/*

    * Copyright 2003 Sun Microsystems, Inc.  All Rights Reserved.

    * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.

    *

    * This code is free software; you can redistribute it and/or modify it

    * under the terms of the GNU General Public License version 2 only, as

    * published by the Free Software Foundation.  Sun designates this

    * particular file as subject to the "Classpath" exception as provided

    * by Sun in the LICENSE file that accompanied this code.

   *

   * This code 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

   * version 2 for more details (a copy is included in the LICENSE file that

   * accompanied this code).

   *

   * You should have received a copy of the GNU General Public License version

   * 2 along with this work; if not, write to the Free Software Foundation,

   * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.

   *

   * Please contact Sun Microsystems, Inc., 4150 Network Circle, Santa Clara,

   * CA 95054 USA or visit www.sun.com if you need additional information or

   * have any questions.

   */

  package sun.misc;

  import sun.misc.FloatConsts;

  import sun.misc.DoubleConsts;

  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.

      */

    

     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.

     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);

     }

    

     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;

     }

    

     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);

     }

    

     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);

     }

    

     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)));

     }

    

     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)));

     }

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

    

     public static boolean isFinite(double d) {

         return Math.abs(d) <= DoubleConsts.MAX_VALUE;

     }

    

      public static boolean isFinite(float f) {

         return Math.abs(f) <= FloatConsts.MAX_VALUE;

     }

    

     public static boolean isInfinite(double d) {

         return Double.isInfinite(d);

     }

    

      public static boolean isInfinite(float f) {

          return Float.isInfinite(f);

     }

    

     public static boolean isNaN(double d) {

         return Double.isNaN(d);

     }

    

      public static boolean isNaN(float f) {

         return Float.isNaN(f);

     }

    

     public static boolean isUnordered(double arg1, double arg2) {

         return isNaN(arg1) || isNaN(arg2);

     }

    

      public static boolean isUnordered(float arg1, float arg2) {

         return isNaN(arg1) || isNaN(arg2);

     }

    

     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

         // break;

         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;

             }

         // break;

         default:

             assert( exponent >= DoubleConsts.MIN_EXPONENT &&

                     exponent <= DoubleConsts.MAX_EXPONENT);

             return exponent;

         // break;

         }

     }

    

      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

         // break;

         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;

             }

         // break;

         default:

             assert( exponent >= FloatConsts.MIN_EXPONENT &&

                     exponent <= FloatConsts.MAX_EXPONENT);

             return exponent;

         // break;

         }

     }

     /*

      * 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.

      */

    

     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;

     }

    

      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));

     }

    

     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);

         }

     }

    

      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);

         }

     }

    

     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));

         }

     }

    

      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) +