1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2020 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 "MathLib.h"
39#include "root.tbl"
40#include <math-barriers.h>
41#include <math_private.h>
42#include <fenv_private.h>
43#include <libm-alias-finite.h>
44#include <math-use-builtins.h>
45
46/*********************************************************************/
47/* An ultimate sqrt routine. Given an IEEE double machine number x */
48/* it computes the correctly rounded (to nearest) value of square */
49/* root of x. */
50/*********************************************************************/
51double
52__ieee754_sqrt (double x)
53{
54#if USE_SQRT_BUILTIN
55 return __builtin_sqrt (x);
56#else
57 /* Use generic implementation. */
58 static const double
59 rt0 = 9.99999999859990725855365213134618E-01,
60 rt1 = 4.99999999495955425917856814202739E-01,
61 rt2 = 3.75017500867345182581453026130850E-01,
62 rt3 = 3.12523626554518656309172508769531E-01;
63 static const double big = 134217728.0;
64 double y, t, del, res, res1, hy, z, zz, s;
65 mynumber a, c = { { 0, 0 } };
66 int4 k;
67
68 a.x = x;
69 k = a.i[HIGH_HALF];
70 a.i[HIGH_HALF] = (k & 0x001fffff) | 0x3fe00000;
71 t = inroot[(k & 0x001fffff) >> 14];
72 s = a.x;
73 /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/
74 if (k > 0x000fffff && k < 0x7ff00000)
75 {
76 int rm = __fegetround ();
77 fenv_t env;
78 libc_feholdexcept_setround (&env, FE_TONEAREST);
79 double ret;
80 y = 1.0 - t * (t * s);
81 t = t * (rt0 + y * (rt1 + y * (rt2 + y * rt3)));
82 c.i[HIGH_HALF] = 0x20000000 + ((k & 0x7fe00000) >> 1);
83 y = t * s;
84 hy = (y + big) - big;
85 del = 0.5 * t * ((s - hy * hy) - (y - hy) * (y + hy));
86 res = y + del;
87 if (res == (res + 1.002 * ((y - res) + del)))
88 ret = res * c.x;
89 else
90 {
91 res1 = res + 1.5 * ((y - res) + del);
92 EMULV (res, res1, z, zz); /* (z+zz)=res*res1 */
93 res = ((((z - s) + zz) < 0) ? max (res, res1) :
94 min (res, res1));
95 ret = res * c.x;
96 }
97 math_force_eval (ret);
98 libc_fesetenv (&env);
99 double dret = x / ret;
100 if (dret != ret)
101 {
102 double force_inexact = 1.0 / 3.0;
103 math_force_eval (force_inexact);
104 /* The square root is inexact, ret is the round-to-nearest
105 value which may need adjusting for other rounding
106 modes. */
107 switch (rm)
108 {
109#ifdef FE_UPWARD
110 case FE_UPWARD:
111 if (dret > ret)
112 ret = (res + 0x1p-1022) * c.x;
113 break;
114#endif
115
116#ifdef FE_DOWNWARD
117 case FE_DOWNWARD:
118#endif
119#ifdef FE_TOWARDZERO
120 case FE_TOWARDZERO:
121#endif
122#if defined FE_DOWNWARD || defined FE_TOWARDZERO
123 if (dret < ret)
124 ret = (res - 0x1p-1022) * c.x;
125 break;
126#endif
127
128 default:
129 break;
130 }
131 }
132 /* Otherwise (x / ret == ret), either the square root was exact or
133 the division was inexact. */
134 return ret;
135 }
136 else
137 {
138 if ((k & 0x7ff00000) == 0x7ff00000)
139 return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
140 if (x == 0)
141 return x; /* sqrt(+0)=+0, sqrt(-0)=-0 */
142 if (k < 0)
143 return (x - x) / (x - x); /* sqrt(-ve)=sNaN */
144 return 0x1p-256 * __ieee754_sqrt (x * 0x1p512);
145 }
146#endif /* ! USE_SQRT_BUILTIN */
147}
148#ifndef __ieee754_sqrt
149libm_alias_finite (__ieee754_sqrt, __sqrt)
150#endif
151