1 /*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2023 Free Software Foundation, Inc.
5 *
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.
10 *
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.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <https://www.gnu.org/licenses/>.
18 */
19 /*********************************************************************/
20 /* MODULE_NAME: uroot.c */
21 /* */
22 /* FUNCTION: usqrt */
23 /* */
24 /* FILES NEEDED: dla.h endian.h mydefs.h */
25 /* uroot.tbl */
26 /* */
27 /* An ultimate sqrt routine. Given an IEEE double machine number x */
28 /* it computes the correctly rounded (to nearest) value of square */
29 /* root of x. */
30 /* Assumption: Machine arithmetic operations are performed in */
31 /* round to nearest mode of IEEE 754 standard. */
32 /* */
33 /*********************************************************************/
34
35 #include "endian.h"
36 #include "mydefs.h"
37 #include <dla.h>
38 #include "root.tbl"
39 #include <math-barriers.h>
40 #include <math_private.h>
41 #include <fenv_private.h>
42 #include <libm-alias-finite.h>
43 #include <math-use-builtins.h>
44
45 /*********************************************************************/
46 /* An ultimate sqrt routine. Given an IEEE double machine number x */
47 /* it computes the correctly rounded (to nearest) value of square */
48 /* root of x. */
49 /*********************************************************************/
50 double
51 __ieee754_sqrt (double x)
52 {
53 #if USE_SQRT_BUILTIN
54 return __builtin_sqrt (x);
55 #else
56 /* Use generic implementation. */
57 static const double
58 rt0 = 9.99999999859990725855365213134618E-01,
59 rt1 = 4.99999999495955425917856814202739E-01,
60 rt2 = 3.75017500867345182581453026130850E-01,
61 rt3 = 3.12523626554518656309172508769531E-01;
62 static const double big = 134217728.0;
63 double y, t, del, res, res1, hy, z, zz, s;
64 mynumber a, c = { { 0, 0 } };
65 int4 k;
66
67 a.x = x;
68 k = a.i[HIGH_HALF];
69 a.i[HIGH_HALF] = (k & 0x001fffff) | 0x3fe00000;
70 t = inroot[(k & 0x001fffff) >> 14];
71 s = a.x;
72 /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/
73 if (k > 0x000fffff && k < 0x7ff00000)
74 {
75 int rm = __fegetround ();
76 fenv_t env;
77 libc_feholdexcept_setround (&env, FE_TONEAREST);
78 double ret;
79 y = 1.0 - t * (t * s);
80 t = t * (rt0 + y * (rt1 + y * (rt2 + y * rt3)));
81 c.i[HIGH_HALF] = 0x20000000 + ((k & 0x7fe00000) >> 1);
82 y = t * s;
83 hy = (y + big) - big;
84 del = 0.5 * t * ((s - hy * hy) - (y - hy) * (y + hy));
85 res = y + del;
86 if (res == (res + 1.002 * ((y - res) + del)))
87 ret = res * c.x;
88 else
89 {
90 res1 = res + 1.5 * ((y - res) + del);
91 EMULV (res, res1, z, zz); /* (z+zz)=res*res1 */
92 res = ((((z - s) + zz) < 0) ? max (res, res1) :
93 min (res, res1));
94 ret = res * c.x;
95 }
96 math_force_eval (ret);
97 libc_fesetenv (&env);
98 double dret = x / ret;
99 if (dret != ret)
100 {
101 double force_inexact = 1.0 / 3.0;
102 math_force_eval (force_inexact);
103 /* The square root is inexact, ret is the round-to-nearest
104 value which may need adjusting for other rounding
105 modes. */
106 switch (rm)
107 {
108 #ifdef FE_UPWARD
109 case FE_UPWARD:
110 if (dret > ret)
111 ret = (res + 0x1p-1022) * c.x;
112 break;
113 #endif
114
115 #ifdef FE_DOWNWARD
116 case FE_DOWNWARD:
117 #endif
118 #ifdef FE_TOWARDZERO
119 case FE_TOWARDZERO:
120 #endif
121 #if defined FE_DOWNWARD || defined FE_TOWARDZERO
122 if (dret < ret)
123 ret = (res - 0x1p-1022) * c.x;
124 break;
125 #endif
126
127 default:
128 break;
129 }
130 }
131 /* Otherwise (x / ret == ret), either the square root was exact or
132 the division was inexact. */
133 return ret;
134 }
135 else
136 {
137 if ((k & 0x7ff00000) == 0x7ff00000)
138 return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
139 if (x == 0)
140 return x; /* sqrt(+0)=+0, sqrt(-0)=-0 */
141 if (k < 0)
142 return (x - x) / (x - x); /* sqrt(-ve)=sNaN */
143 return 0x1p-256 * __ieee754_sqrt (x * 0x1p512);
144 }
145 #endif /* ! USE_SQRT_BUILTIN */
146 }
147 #ifndef __ieee754_sqrt
148 libm_alias_finite (__ieee754_sqrt, __sqrt)
149 #endif