4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
30 #pragma weak __hypotl = hypotl
33 * long double hypotl(long double x, long double y);
35 * If z=x*x+y*y has error less than sqrt(2)/2 ulp than sqrt(z) has
36 * error less than 1 ulp.
37 * So, compute sqrt(x*x+y*y) with some care as follows:
39 * 1. save and set rounding to round-to-nearest
41 * x1*x1+(y*y+(x2*(x+x2))) for x*x+y*y
42 * where x1 = x with lower 64 bits cleared, x2 = x-x1; else
44 * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
45 * where t1 = 2x with lower 64 bits cleared, t2 = 2x-t1, y1= y with
46 * lower 64 bits chopped, y2 = y-y1.
48 * NOTE: DO NOT remove parenthsis!
51 * hypot(x,y) is INF if x or y is +INF or -INF; else
52 * hypot(x,y) is NAN if x or y is NAN.
55 * hypot(x,y) returns sqrt(x^2+y^2) with error less than 1 ulps (units
60 #include "longdouble.h"
62 extern enum fp_direction_type
__swapRD(enum fp_direction_type
);
64 static const long double zero
= 0.0L, one
= 1.0L;
67 hypotl(long double x
, long double y
) {
69 long double t1
, t2
, y1
, y2
, w
;
70 int *px
= (int *) &x
, *py
= (int *) &y
;
71 int *pt1
= (int *) &t1
, *py1
= (int *) &y1
;
72 enum fp_direction_type rd
;
75 if ((*(int *) &one
) != 0) { /* determine word ordering */
87 px
[n0
] &= 0x7fffffff; /* clear sign bit of x and y */
90 nx
= px
[n0
] & k
; /* exponent of x and y */
100 if ((nx
- ny
) >= 0x00730000)
101 return (x
+ y
); /* x/y >= 2**116 */
102 if (nx
< 0x5ff30000 && ny
> 0x205b0000) { /* medium x,y */
103 /* save and set RD to Rounding to nearest */
104 rd
= __swapRD(fp_nearest
);
109 pt1
[n2
] = pt1
[n3
] = 0;
111 x
= sqrtl(t1
* t1
- (y
* (-y
) - t2
* (x
+ t1
)));
116 py1
[n2
] = py1
[n3
] = 0;
120 pt1
[n2
] = pt1
[n3
] = 0;
122 x
= sqrtl(t1
* y1
- (w
* (-w
) - (t2
* y1
+ y2
* x
)));
124 if (rd
!= fp_nearest
)
125 (void) __swapRD(rd
); /* restore rounding mode */
128 if (nx
== k
|| ny
== k
) { /* x or y is INF or NaN */
134 t2
= x
+ y
; /* invalid if x or y is sNaN */
138 if (y
== zero
|| x
== zero
)
140 t1
= scalbnl(one
, 16381);
143 return (scalbnl(one
, -16381) * hypotl(x
, y
));
149 pt1
[n1
] = pt1
[n2
] = pt1
[n3
] = 0;
150 return (t1
* hypotl(x
, y
));