1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2021 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/*********************************************************************/
50double
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
148libm_alias_finite (__ieee754_sqrt, __sqrt)
149#endif
150