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 }