jdk/test/java/lang/StrictMath/FdlibmTranslit.java
changeset 32928 a3f03999ed62
child 32992 19ed8781c2ba
equal deleted inserted replaced
32917:8392405ab038 32928:a3f03999ed62
       
     1 /*
       
     2  * Copyright (c) 1998, 2015, Oracle and/or its affiliates. All rights reserved.
       
     3  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
       
     4  *
       
     5  * This code is free software; you can redistribute it and/or modify it
       
     6  * under the terms of the GNU General Public License version 2 only, as
       
     7  * published by the Free Software Foundation.  Oracle designates this
       
     8  * particular file as subject to the "Classpath" exception as provided
       
     9  * by Oracle in the LICENSE file that accompanied this code.
       
    10  *
       
    11  * This code is distributed in the hope that it will be useful, but WITHOUT
       
    12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
       
    13  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
       
    14  * version 2 for more details (a copy is included in the LICENSE file that
       
    15  * accompanied this code).
       
    16  *
       
    17  * You should have received a copy of the GNU General Public License version
       
    18  * 2 along with this work; if not, write to the Free Software Foundation,
       
    19  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
       
    20  *
       
    21  * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
       
    22  * or visit www.oracle.com if you need additional information or have any
       
    23  * questions.
       
    24  */
       
    25 
       
    26 /**
       
    27  * A transliteration of the "Freely Distributable Math Library"
       
    28  * algorithms from C into Java. That is, this port of the algorithms
       
    29  * is as close to the C originals as possible while still being
       
    30  * readable legal Java.
       
    31  */
       
    32 public class FdlibmTranslit {
       
    33     private FdlibmTranslit() {
       
    34         throw new UnsupportedOperationException("No FdLibmTranslit instances for you.");
       
    35     }
       
    36 
       
    37     /**
       
    38      * Return the low-order 32 bits of the double argument as an int.
       
    39      */
       
    40     private static int __LO(double x) {
       
    41         long transducer = Double.doubleToRawLongBits(x);
       
    42         return (int)transducer;
       
    43     }
       
    44 
       
    45     /**
       
    46      * Return a double with its low-order bits of the second argument
       
    47      * and the high-order bits of the first argument..
       
    48      */
       
    49     private static double __LO(double x, int low) {
       
    50         long transX = Double.doubleToRawLongBits(x);
       
    51         return Double.longBitsToDouble((transX & 0xFFFF_FFFF_0000_0000L)|low );
       
    52     }
       
    53 
       
    54     /**
       
    55      * Return the high-order 32 bits of the double argument as an int.
       
    56      */
       
    57     private static int __HI(double x) {
       
    58         long transducer = Double.doubleToRawLongBits(x);
       
    59         return (int)(transducer >> 32);
       
    60     }
       
    61 
       
    62     /**
       
    63      * Return a double with its high-order bits of the second argument
       
    64      * and the low-order bits of the first argument..
       
    65      */
       
    66     private static double __HI(double x, int high) {
       
    67         long transX = Double.doubleToRawLongBits(x);
       
    68         return Double.longBitsToDouble((transX & 0x0000_0000_FFFF_FFFFL)|( ((long)high)) << 32 );
       
    69     }
       
    70 
       
    71     public static double hypot(double x, double y) {
       
    72         return Hypot.compute(x, y);
       
    73     }
       
    74 
       
    75     /**
       
    76      * hypot(x,y)
       
    77      *
       
    78      * Method :
       
    79      *      If (assume round-to-nearest) z = x*x + y*y
       
    80      *      has error less than sqrt(2)/2 ulp, than
       
    81      *      sqrt(z) has error less than 1 ulp (exercise).
       
    82      *
       
    83      *      So, compute sqrt(x*x + y*y) with some care as
       
    84      *      follows to get the error below 1 ulp:
       
    85      *
       
    86      *      Assume x > y > 0;
       
    87      *      (if possible, set rounding to round-to-nearest)
       
    88      *      1. if x > 2y  use
       
    89      *              x1*x1 + (y*y + (x2*(x + x1))) for x*x + y*y
       
    90      *      where x1 = x with lower 32 bits cleared, x2 = x - x1; else
       
    91      *      2. if x <= 2y use
       
    92      *              t1*y1 + ((x-y) * (x-y) + (t1*y2 + t2*y))
       
    93      *      where t1 = 2x with lower 32 bits cleared, t2 = 2x - t1,
       
    94      *      y1= y with lower 32 bits chopped, y2 = y - y1.
       
    95      *
       
    96      *      NOTE: scaling may be necessary if some argument is too
       
    97      *            large or too tiny
       
    98      *
       
    99      * Special cases:
       
   100      *      hypot(x,y) is INF if x or y is +INF or -INF; else
       
   101      *      hypot(x,y) is NAN if x or y is NAN.
       
   102      *
       
   103      * Accuracy:
       
   104      *      hypot(x,y) returns sqrt(x^2 + y^2) with error less
       
   105      *      than 1 ulps (units in the last place)
       
   106      */
       
   107     static class Hypot {
       
   108         public static double compute(double x, double y) {
       
   109             double a = x;
       
   110             double b = y;
       
   111             double t1, t2, y1, y2, w;
       
   112             int j, k, ha, hb;
       
   113 
       
   114             ha = __HI(x) & 0x7fffffff;        // high word of  x
       
   115             hb = __HI(y) & 0x7fffffff;        // high word of  y
       
   116             if(hb > ha) {
       
   117                 a = y;
       
   118                 b = x;
       
   119                 j = ha;
       
   120                 ha = hb;
       
   121                 hb = j;
       
   122             } else {
       
   123                 a = x;
       
   124                 b = y;
       
   125             }
       
   126             a = __HI(a, ha);   // a <- |a|
       
   127             b = __HI(b, hb);   // b <- |b|
       
   128             if ((ha - hb) > 0x3c00000) {
       
   129                 return a + b;  // x / y > 2**60
       
   130             }
       
   131             k=0;
       
   132             if (ha > 0x5f300000) {   // a>2**500
       
   133                 if (ha >= 0x7ff00000) {       // Inf or NaN
       
   134                     w = a + b;                // for sNaN
       
   135                     if (((ha & 0xfffff) | __LO(a)) == 0)
       
   136                         w = a;
       
   137                     if (((hb ^ 0x7ff00000) | __LO(b)) == 0)
       
   138                         w = b;
       
   139                     return w;
       
   140                 }
       
   141                 // scale a and b by 2**-600
       
   142                 ha -= 0x25800000;
       
   143                 hb -= 0x25800000;
       
   144                 k += 600;
       
   145                 a = __HI(a, ha);
       
   146                 b = __HI(b, hb);
       
   147             }
       
   148             if (hb < 0x20b00000) {   // b < 2**-500
       
   149                 if (hb <= 0x000fffff) {      // subnormal b or 0 */
       
   150                     if ((hb | (__LO(b))) == 0)
       
   151                         return a;
       
   152                     t1 = 0;
       
   153                     t1 = __HI(t1, 0x7fd00000);  // t1=2^1022
       
   154                     b *= t1;
       
   155                     a *= t1;
       
   156                     k -= 1022;
       
   157                 } else {            // scale a and b by 2^600
       
   158                     ha += 0x25800000;       // a *= 2^600
       
   159                     hb += 0x25800000;       // b *= 2^600
       
   160                     k -= 600;
       
   161                     a = __HI(a, ha);
       
   162                     b = __HI(b, hb);
       
   163                 }
       
   164             }
       
   165             // medium size a and b
       
   166             w = a - b;
       
   167             if (w > b) {
       
   168                 t1 = 0;
       
   169                 t1 = __HI(t1, ha);
       
   170                 t2 = a - t1;
       
   171                 w  = Math.sqrt(t1*t1 - (b*(-b) - t2 * (a + t1)));
       
   172             } else {
       
   173                 a  = a + a;
       
   174                 y1 = 0;
       
   175                 y1 = __HI(y1, hb);
       
   176                 y2 = b - y1;
       
   177                 t1 = 0;
       
   178                 t1 = __HI(t1, ha + 0x00100000);
       
   179                 t2 = a - t1;
       
   180                 w  = Math.sqrt(t1*y1 - (w*(-w) - (t1*y2 + t2*b)));
       
   181             }
       
   182             if (k != 0) {
       
   183                 t1 = 1.0;
       
   184                 int t1_hi = __HI(t1);
       
   185                 t1_hi += (k << 20);
       
   186                 t1 = __HI(t1, t1_hi);
       
   187                 return t1 * w;
       
   188             } else
       
   189                 return w;
       
   190         }
       
   191     }
       
   192 }