73 return Double.longBitsToDouble((transX & 0xFFFF_FFFF_0000_0000L)|low ); 74 } 75 76 /** 77 * Return the high-order 32 bits of the double argument as an int. 78 */ 79 private static int __HI(double x) { 80 long transducer = Double.doubleToRawLongBits(x); 81 return (int)(transducer >> 32); 82 } 83 84 /** 85 * Return a double with its high-order bits of the second argument 86 * and the low-order bits of the first argument.. 87 */ 88 private static double __HI(double x, int high) { 89 long transX = Double.doubleToRawLongBits(x); 90 return Double.longBitsToDouble((transX & 0x0000_0000_FFFF_FFFFL)|( ((long)high)) << 32 ); 91 } 92 93 /** 94 * Compute x**y 95 * n 96 * Method: Let x = 2 * (1+f) 97 * 1. Compute and return log2(x) in two pieces: 98 * log2(x) = w1 + w2, 99 * where w1 has 53 - 24 = 29 bit trailing zeros. 100 * 2. Perform y*log2(x) = n+y' by simulating muti-precision 101 * arithmetic, where |y'| <= 0.5. 102 * 3. Return x**y = 2**n*exp(y'*log2) 103 * 104 * Special cases: 105 * 1. (anything) ** 0 is 1 106 * 2. (anything) ** 1 is itself 107 * 3. (anything) ** NAN is NAN 108 * 4. NAN ** (anything except 0) is NAN 109 * 5. +-(|x| > 1) ** +INF is +INF 110 * 6. +-(|x| > 1) ** -INF is +0 111 * 7. +-(|x| < 1) ** +INF is +0 112 * 8. +-(|x| < 1) ** -INF is +INF | 73 return Double.longBitsToDouble((transX & 0xFFFF_FFFF_0000_0000L)|low ); 74 } 75 76 /** 77 * Return the high-order 32 bits of the double argument as an int. 78 */ 79 private static int __HI(double x) { 80 long transducer = Double.doubleToRawLongBits(x); 81 return (int)(transducer >> 32); 82 } 83 84 /** 85 * Return a double with its high-order bits of the second argument 86 * and the low-order bits of the first argument.. 87 */ 88 private static double __HI(double x, int high) { 89 long transX = Double.doubleToRawLongBits(x); 90 return Double.longBitsToDouble((transX & 0x0000_0000_FFFF_FFFFL)|( ((long)high)) << 32 ); 91 } 92 93 /** 94 * hypot(x,y) 95 * 96 * Method : 97 * If (assume round-to-nearest) z=x*x+y*y 98 * has error less than sqrt(2)/2 ulp, than 99 * sqrt(z) has error less than 1 ulp (exercise). 100 * 101 * So, compute sqrt(x*x+y*y) with some care as 102 * follows to get the error below 1 ulp: 103 * 104 * Assume x>y>0; 105 * (if possible, set rounding to round-to-nearest) 106 * 1. if x > 2y use 107 * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y 108 * where x1 = x with lower 32 bits cleared, x2 = x-x1; else 109 * 2. if x <= 2y use 110 * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y)) 111 * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, 112 * y1= y with lower 32 bits chopped, y2 = y-y1. 113 * 114 * NOTE: scaling may be necessary if some argument is too 115 * large or too tiny 116 * 117 * Special cases: 118 * hypot(x,y) is INF if x or y is +INF or -INF; else 119 * hypot(x,y) is NAN if x or y is NAN. 120 * 121 * Accuracy: 122 * hypot(x,y) returns sqrt(x^2+y^2) with error less 123 * than 1 ulps (units in the last place) 124 */ 125 public static class Hypot { 126 public static double compute(double x, double y) { 127 double a=x,b=y,t1,t2,y1,y2,w; 128 int j,k,ha,hb; 129 130 ha = __HI(x)&0x7fffffff; /* high word of x */ 131 hb = __HI(y)&0x7fffffff; /* high word of y */ 132 if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} 133 a = __HI(a, ha); /* a <- |a| */ 134 b = __HI(b, hb); /* b <- |b| */ 135 if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */ 136 k=0; 137 if(ha > 0x5f300000) { /* a>2**500 */ 138 if(ha >= 0x7ff00000) { /* Inf or NaN */ 139 w = a+b; /* for sNaN */ 140 if(((ha&0xfffff)|__LO(a))==0) w = a; 141 if(((hb^0x7ff00000)|__LO(b))==0) w = b; 142 return w; 143 } 144 /* scale a and b by 2**-600 */ 145 ha -= 0x25800000; hb -= 0x25800000; k += 600; 146 a = __HI(a, ha); 147 b = __HI(b, hb); 148 } 149 if(hb < 0x20b00000) { /* b < 2**-500 */ 150 if(hb <= 0x000fffff) { /* subnormal b or 0 */ 151 if((hb|(__LO(b)))==0) 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 return w; 189 } 190 } 191 192 /** 193 * Compute x**y 194 * n 195 * Method: Let x = 2 * (1+f) 196 * 1. Compute and return log2(x) in two pieces: 197 * log2(x) = w1 + w2, 198 * where w1 has 53 - 24 = 29 bit trailing zeros. 199 * 2. Perform y*log2(x) = n+y' by simulating muti-precision 200 * arithmetic, where |y'| <= 0.5. 201 * 3. Return x**y = 2**n*exp(y'*log2) 202 * 203 * Special cases: 204 * 1. (anything) ** 0 is 1 205 * 2. (anything) ** 1 is itself 206 * 3. (anything) ** NAN is NAN 207 * 4. NAN ** (anything except 0) is NAN 208 * 5. +-(|x| > 1) ** +INF is +INF 209 * 6. +-(|x| > 1) ** -INF is +0 210 * 7. +-(|x| < 1) ** +INF is +0 211 * 8. +-(|x| < 1) ** -INF is +INF |