1 /* origin: FreeBSD /usr/src/lib/msun/src/e_hypot.c */
2 /*
3  * ====================================================
4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5  *
6  * Developed at SunSoft, a Sun Microsystems, Inc. business.
7  * Permission to use, copy, modify, and distribute this
8  * software is freely granted, provided that this notice
9  * is preserved.
10  * ====================================================
11  */
12 /* hypot(x,y)
13  *
14  * Method :
15  *      If (assume round-to-nearest) z=x*x+y*y
16  *      has error less than sqrt(2)/2 ulp, then
17  *      sqrt(z) has error less than 1 ulp (exercise).
18  *
19  *      So, compute sqrt(x*x+y*y) with some care as
20  *      follows to get the error below 1 ulp:
21  *
22  *      Assume x>y>0;
23  *      (if possible, set rounding to round-to-nearest)
24  *      1. if x > 2y  use
25  *              x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
26  *      where x1 = x with lower 32 bits cleared, x2 = x-x1; else
27  *      2. if x <= 2y use
28  *              t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
29  *      where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
30  *      y1= y with lower 32 bits chopped, y2 = y-y1.
31  *
32  *      NOTE: scaling may be necessary if some argument is too
33  *            large or too tiny
34  *
35  * Special cases:
36  *      hypot(x,y) is INF if x or y is +INF or -INF; else
37  *      hypot(x,y) is NAN if x or y is NAN.
38  *
39  * Accuracy:
40  *      hypot(x,y) returns sqrt(x^2+y^2) with error less
41  *      than 1 ulps (units in the last place)
42  */
44 #include "libm.h"
46 double hypot(double x, double y)
47 {
48         double a,b,t1,t2,y1,y2,w;
49         int32_t j,k,ha,hb;
51         GET_HIGH_WORD(ha, x);
52         ha &= 0x7fffffff;
53         GET_HIGH_WORD(hb, y);
54         hb &= 0x7fffffff;
55         if (hb > ha) {
56                 a = y;
57                 b = x;
58                 j=ha; ha=hb; hb=j;
59         } else {
60                 a = x;
61                 b = y;
62         }
63         a = fabs(a);
64         b = fabs(b);
65         if (ha - hb > 0x3c00000)  /* x/y > 2**60 */
66                 return a+b;
67         k = 0;
68         if (ha > 0x5f300000) {    /* a > 2**500 */
69                 if(ha >= 0x7ff00000) {  /* Inf or NaN */
70                         uint32_t low;
71                         /* Use original arg order iff result is NaN; quieten sNaNs. */
72                         w = fabs(x+0.0) - fabs(y+0.0);
73                         GET_LOW_WORD(low, a);
74                         if (((ha&0xfffff)|low) == 0) w = a;
75                         GET_LOW_WORD(low, b);
76                         if (((hb^0x7ff00000)|low) == 0) w = b;
77                         return w;
78                 }
79                 /* scale a and b by 2**-600 */
80                 ha -= 0x25800000; hb -= 0x25800000;  k += 600;
81                 SET_HIGH_WORD(a, ha);
82                 SET_HIGH_WORD(b, hb);
83         }
84         if (hb < 0x20b00000) {    /* b < 2**-500 */
85                 if (hb <= 0x000fffff) {  /* subnormal b or 0 */
86                         uint32_t low;
87                         GET_LOW_WORD(low, b);
88                         if ((hb|low) == 0)
89                                 return a;
90                         t1 = 0;
91                         SET_HIGH_WORD(t1, 0x7fd00000);  /* t1 = 2^1022 */
92                         b *= t1;
93                         a *= t1;
94                         k -= 1022;
95                 } else {            /* scale a and b by 2^600 */
96                         ha += 0x25800000;  /* a *= 2^600 */
97                         hb += 0x25800000;  /* b *= 2^600 */
98                         k -= 600;
99                         SET_HIGH_WORD(a, ha);
100                         SET_HIGH_WORD(b, hb);
101                 }
102         }
103         /* medium size a and b */
104         w = a - b;
105         if (w > b) {
106                 t1 = 0;
107                 SET_HIGH_WORD(t1, ha);
108                 t2 = a-t1;
109                 w  = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
110         } else {
111                 a  = a + a;
112                 y1 = 0;
113                 SET_HIGH_WORD(y1, hb);
114                 y2 = b - y1;
115                 t1 = 0;
116                 SET_HIGH_WORD(t1, ha+0x00100000);
117                 t2 = a - t1;
118                 w  = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
119         }
120         if (k != 0) {
121                 uint32_t high;
122                 t1 = 1.0;
123                 GET_HIGH_WORD(high, t1);
124                 SET_HIGH_WORD(t1, high+(k<<20));
125                 return t1*w;
126         }
127         return w;
128 }