Deleted Added
sdiff udiff text old ( 97413 ) new ( 141296 )
full compact
e_hypot.c (97413) e_hypot.c (141296)
1/* @(#)e_hypot.c 5.1 93/09/24 */
1
2/* @(#)e_hypot.c 1.3 95/01/18 */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
3/*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Developed at SunSoft, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13#ifndef lint
10 * is preserved.
11 * ====================================================
12 */
13
14#ifndef lint
14static char rcsid[] = "$FreeBSD: head/lib/msun/src/e_hypot.c 97413 2002-05-28 18:15:04Z alfred $";
15static char rcsid[] = "$FreeBSD: head/lib/msun/src/e_hypot.c 141296 2005-02-04 18:26:06Z das $";
15#endif
16
17/* __ieee754_hypot(x,y)
18 *
16#endif
17
18/* __ieee754_hypot(x,y)
19 *
19 * Method :
20 * If (assume round-to-nearest) z=x*x+y*y
21 * has error less than sqrt(2)/2 ulp, than
20 * Method :
21 * If (assume round-to-nearest) z=x*x+y*y
22 * has error less than sqrt(2)/2 ulp, than
22 * sqrt(z) has error less than 1 ulp (exercise).
23 *
23 * sqrt(z) has error less than 1 ulp (exercise).
24 *
24 * So, compute sqrt(x*x+y*y) with some care as
25 * So, compute sqrt(x*x+y*y) with some care as
25 * follows to get the error below 1 ulp:
26 *
27 * Assume x>y>0;
28 * (if possible, set rounding to round-to-nearest)
29 * 1. if x > 2y use
30 * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
31 * where x1 = x with lower 32 bits cleared, x2 = x-x1; else
32 * 2. if x <= 2y use
33 * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
26 * follows to get the error below 1 ulp:
27 *
28 * Assume x>y>0;
29 * (if possible, set rounding to round-to-nearest)
30 * 1. if x > 2y use
31 * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
32 * where x1 = x with lower 32 bits cleared, x2 = x-x1; else
33 * 2. if x <= 2y use
34 * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
34 * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
35 * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
35 * y1= y with lower 32 bits chopped, y2 = y-y1.
36 * y1= y with lower 32 bits chopped, y2 = y-y1.
36 *
37 * NOTE: scaling may be necessary if some argument is too
37 *
38 * NOTE: scaling may be necessary if some argument is too
38 * large or too tiny
39 *
40 * Special cases:
41 * hypot(x,y) is INF if x or y is +INF or -INF; else
42 * hypot(x,y) is NAN if x or y is NAN.
43 *
44 * Accuracy:
39 * large or too tiny
40 *
41 * Special cases:
42 * hypot(x,y) is INF if x or y is +INF or -INF; else
43 * hypot(x,y) is NAN if x or y is NAN.
44 *
45 * Accuracy:
45 * hypot(x,y) returns sqrt(x^2+y^2) with error less
46 * than 1 ulps (units in the last place)
46 * hypot(x,y) returns sqrt(x^2+y^2) with error less
47 * than 1 ulps (units in the last place)
47 */
48
49#include "math.h"
50#include "math_private.h"
51
52double
53__ieee754_hypot(double x, double y)
54{

--- 43 unchanged lines hidden (view full) ---

98 }
99 }
100 /* medium size a and b */
101 w = a-b;
102 if (w>b) {
103 t1 = 0;
104 SET_HIGH_WORD(t1,ha);
105 t2 = a-t1;
48 */
49
50#include "math.h"
51#include "math_private.h"
52
53double
54__ieee754_hypot(double x, double y)
55{

--- 43 unchanged lines hidden (view full) ---

99 }
100 }
101 /* medium size a and b */
102 w = a-b;
103 if (w>b) {
104 t1 = 0;
105 SET_HIGH_WORD(t1,ha);
106 t2 = a-t1;
106 w = __ieee754_sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
107 w = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
107 } else {
108 a = a+a;
109 y1 = 0;
110 SET_HIGH_WORD(y1,hb);
111 y2 = b - y1;
112 t1 = 0;
113 SET_HIGH_WORD(t1,ha+0x00100000);
114 t2 = a - t1;
108 } else {
109 a = a+a;
110 y1 = 0;
111 SET_HIGH_WORD(y1,hb);
112 y2 = b - y1;
113 t1 = 0;
114 SET_HIGH_WORD(t1,ha+0x00100000);
115 t2 = a - t1;
115 w = __ieee754_sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
116 w = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
116 }
117 if(k!=0) {
118 u_int32_t high;
119 t1 = 1.0;
120 GET_HIGH_WORD(high,t1);
121 SET_HIGH_WORD(t1,high+(k<<20));
122 return t1*w;
123 } else return w;
124}
117 }
118 if(k!=0) {
119 u_int32_t high;
120 t1 = 1.0;
121 GET_HIGH_WORD(high,t1);
122 SET_HIGH_WORD(t1,high+(k<<20));
123 return t1*w;
124 } else return w;
125}