2362N/A * Copyright (c) 1998, 2001, Oracle and/or its affiliates. All rights reserved. 0N/A * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. 0N/A * This code is free software; you can redistribute it and/or modify it 0N/A * under the terms of the GNU General Public License version 2 only, as 2362N/A * published by the Free Software Foundation. Oracle designates this 0N/A * particular file as subject to the "Classpath" exception as provided 2362N/A * by Oracle in the LICENSE file that accompanied this code. 0N/A * This code is distributed in the hope that it will be useful, but WITHOUT 0N/A * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 0N/A * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 0N/A * version 2 for more details (a copy is included in the LICENSE file that 0N/A * accompanied this code). 0N/A * You should have received a copy of the GNU General Public License version 0N/A * 2 along with this work; if not, write to the Free Software Foundation, 0N/A * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 2362N/A * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA 0N/A/* __ieee754_hypot(x,y) 0N/A * If (assume round-to-nearest) z=x*x+y*y 0N/A * has error less than sqrt(2)/2 ulp, than 0N/A * sqrt(z) has error less than 1 ulp (exercise). 0N/A * So, compute sqrt(x*x+y*y) with some care as 0N/A * follows to get the error below 1 ulp: 0N/A * (if possible, set rounding to round-to-nearest) 0N/A * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y 0N/A * where x1 = x with lower 32 bits cleared, x2 = x-x1; else 0N/A * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y)) 0N/A * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, 0N/A * y1= y with lower 32 bits chopped, y2 = y-y1. 0N/A * NOTE: scaling may be necessary if some argument is too 0N/A * hypot(x,y) is INF if x or y is +INF or -INF; else 0N/A * hypot(x,y) is NAN if x or y is NAN. 0N/A * hypot(x,y) returns sqrt(x^2+y^2) with error less 0N/A * than 1 ulps (units in the last place) 0N/A ha =
__HI(x)&
0x7fffffff;
/* high word of x */ 0N/A hb =
__HI(y)&
0x7fffffff;
/* high word of y */ 0N/A if((
ha-
hb)>
0x3c00000) {
return a+b;}
/* x/y > 2**60 */ 0N/A if(
ha >
0x5f300000) {
/* a>2**500 */ 0N/A if(
ha >=
0x7ff00000) {
/* Inf or NaN */ 0N/A w = a+b;
/* for sNaN */ 0N/A /* scale a and b by 2**-600 */ 0N/A ha -=
0x25800000;
hb -=
0x25800000; k +=
600;
0N/A if(
hb <
0x20b00000) {
/* b < 2**-500 */ 0N/A if(
hb <=
0x000fffff) {
/* subnormal b or 0 */ 0N/A }
else {
/* scale a and b by 2^600 */ 0N/A ha +=
0x25800000;
/* a *= 2^600 */ 0N/A hb +=
0x25800000;
/* b *= 2^600 */ 0N/A /* medium size a and b */