]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/dbl-64/e_sqrt.c
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2018 Free Software Foundation, Inc.
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <http://www.gnu.org/licenses/>.
19 /*********************************************************************/
20 /* MODULE_NAME: uroot.c */
24 /* FILES NEEDED: dla.h endian.h mydefs.h */
27 /* An ultimate sqrt routine. Given an IEEE double machine number x */
28 /* it computes the correctly rounded (to nearest) value of square */
30 /* Assumption: Machine arithmetic operations are performed in */
31 /* round to nearest mode of IEEE 754 standard. */
33 /*********************************************************************/
40 #include <math-barriers.h>
41 #include <math_private.h>
42 #include <fenv_private.h>
44 /*********************************************************************/
45 /* An ultimate sqrt routine. Given an IEEE double machine number x */
46 /* it computes the correctly rounded (to nearest) value of square */
48 /*********************************************************************/
50 __ieee754_sqrt (double x
)
53 rt0
= 9.99999999859990725855365213134618E-01,
54 rt1
= 4.99999999495955425917856814202739E-01,
55 rt2
= 3.75017500867345182581453026130850E-01,
56 rt3
= 3.12523626554518656309172508769531E-01;
57 static const double big
= 134217728.0;
58 double y
, t
, del
, res
, res1
, hy
, z
, zz
, p
, hx
, tx
, ty
, s
;
59 mynumber a
, c
= { { 0, 0 } };
64 a
.i
[HIGH_HALF
] = (k
& 0x001fffff) | 0x3fe00000;
65 t
= inroot
[(k
& 0x001fffff) >> 14];
67 /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/
68 if (k
> 0x000fffff && k
< 0x7ff00000)
70 int rm
= __fegetround ();
72 libc_feholdexcept_setround (&env
, FE_TONEAREST
);
74 y
= 1.0 - t
* (t
* s
);
75 t
= t
* (rt0
+ y
* (rt1
+ y
* (rt2
+ y
* rt3
)));
76 c
.i
[HIGH_HALF
] = 0x20000000 + ((k
& 0x7fe00000) >> 1);
79 del
= 0.5 * t
* ((s
- hy
* hy
) - (y
- hy
) * (y
+ hy
));
81 if (res
== (res
+ 1.002 * ((y
- res
) + del
)))
85 res1
= res
+ 1.5 * ((y
- res
) + del
);
86 EMULV (res
, res1
, z
, zz
, p
, hx
, tx
, hy
, ty
); /* (z+zz)=res*res1 */
87 res
= ((((z
- s
) + zz
) < 0) ? max (res
, res1
) :
91 math_force_eval (ret
);
93 double dret
= x
/ ret
;
96 double force_inexact
= 1.0 / 3.0;
97 math_force_eval (force_inexact
);
98 /* The square root is inexact, ret is the round-to-nearest
99 value which may need adjusting for other rounding
106 ret
= (res
+ 0x1p
-1022) * c
.x
;
116 #if defined FE_DOWNWARD || defined FE_TOWARDZERO
118 ret
= (res
- 0x1p
-1022) * c
.x
;
126 /* Otherwise (x / ret == ret), either the square root was exact or
127 the division was inexact. */
132 if ((k
& 0x7ff00000) == 0x7ff00000)
133 return x
* x
+ x
; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
135 return x
; /* sqrt(+0)=+0, sqrt(-0)=-0 */
137 return (x
- x
) / (x
- x
); /* sqrt(-ve)=sNaN */
138 return 0x1p
-256 * __ieee754_sqrt (x
* 0x1p
512);
141 strong_alias (__ieee754_sqrt
, __sqrt_finite
)