]>
Commit | Line | Data |
---|---|---|
f7eac6eb | 1 | /* |
e4d82761 | 2 | * IBM Accurate Mathematical Library |
aeb25823 | 3 | * written by International Business Machines Corp. |
04277e02 | 4 | * Copyright (C) 2001-2019 Free Software Foundation, Inc. |
f7eac6eb | 5 | * |
e4d82761 UD |
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 | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
6d52618b | 10 | * |
e4d82761 UD |
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 | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
6d52618b | 15 | * |
e4d82761 | 16 | * You should have received a copy of the GNU Lesser General Public License |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
f7eac6eb | 18 | */ |
e4d82761 UD |
19 | /*********************************************************************/ |
20 | /* MODULE_NAME: uroot.c */ | |
21 | /* */ | |
22 | /* FUNCTION: usqrt */ | |
23 | /* */ | |
60f435bb | 24 | /* FILES NEEDED: dla.h endian.h mydefs.h */ |
e4d82761 UD |
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" | |
c8b3296b | 37 | #include <dla.h> |
e4d82761 UD |
38 | #include "MathLib.h" |
39 | #include "root.tbl" | |
b4d5b8b0 | 40 | #include <math-barriers.h> |
1ed0291c | 41 | #include <math_private.h> |
70e2ba33 | 42 | #include <fenv_private.h> |
e4d82761 UD |
43 | |
44 | /*********************************************************************/ | |
7c370086 | 45 | /* An ultimate sqrt routine. Given an IEEE double machine number x */ |
e4d82761 UD |
46 | /* it computes the correctly rounded (to nearest) value of square */ |
47 | /* root of x. */ | |
48 | /*********************************************************************/ | |
c5d5d574 OB |
49 | double |
50 | __ieee754_sqrt (double x) | |
51 | { | |
e4d82761 UD |
52 | static const double |
53 | rt0 = 9.99999999859990725855365213134618E-01, | |
54 | rt1 = 4.99999999495955425917856814202739E-01, | |
55 | rt2 = 3.75017500867345182581453026130850E-01, | |
56 | rt3 = 3.12523626554518656309172508769531E-01; | |
c5d5d574 OB |
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 } }; | |
7ea88fb4 | 60 | int4 k; |
e4d82761 | 61 | |
c5d5d574 OB |
62 | a.x = x; |
63 | k = a.i[HIGH_HALF]; | |
64 | a.i[HIGH_HALF] = (k & 0x001fffff) | 0x3fe00000; | |
65 | t = inroot[(k & 0x001fffff) >> 14]; | |
66 | s = a.x; | |
e4d82761 | 67 | /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/ |
c5d5d574 OB |
68 | if (k > 0x000fffff && k < 0x7ff00000) |
69 | { | |
b93c2205 | 70 | int rm = __fegetround (); |
c5d5d574 | 71 | fenv_t env; |
3c1c46a6 | 72 | libc_feholdexcept_setround (&env, FE_TONEAREST); |
c5d5d574 OB |
73 | double ret; |
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); | |
77 | y = t * s; | |
78 | hy = (y + big) - big; | |
79 | del = 0.5 * t * ((s - hy * hy) - (y - hy) * (y + hy)); | |
80 | res = y + del; | |
81 | if (res == (res + 1.002 * ((y - res) + del))) | |
82 | ret = res * c.x; | |
83 | else | |
84 | { | |
85 | res1 = res + 1.5 * ((y - res) + del); | |
86 | EMULV (res, res1, z, zz, p, hx, tx, hy, ty); /* (z+zz)=res*res1 */ | |
3c1c46a6 JM |
87 | res = ((((z - s) + zz) < 0) ? max (res, res1) : |
88 | min (res, res1)); | |
89 | ret = res * c.x; | |
c5d5d574 OB |
90 | } |
91 | math_force_eval (ret); | |
92 | libc_fesetenv (&env); | |
3c1c46a6 JM |
93 | double dret = x / ret; |
94 | if (dret != ret) | |
c5d5d574 OB |
95 | { |
96 | double force_inexact = 1.0 / 3.0; | |
97 | math_force_eval (force_inexact); | |
3c1c46a6 JM |
98 | /* The square root is inexact, ret is the round-to-nearest |
99 | value which may need adjusting for other rounding | |
100 | modes. */ | |
101 | switch (rm) | |
102 | { | |
103 | #ifdef FE_UPWARD | |
104 | case FE_UPWARD: | |
105 | if (dret > ret) | |
106 | ret = (res + 0x1p-1022) * c.x; | |
107 | break; | |
108 | #endif | |
109 | ||
110 | #ifdef FE_DOWNWARD | |
111 | case FE_DOWNWARD: | |
112 | #endif | |
113 | #ifdef FE_TOWARDZERO | |
114 | case FE_TOWARDZERO: | |
115 | #endif | |
116 | #if defined FE_DOWNWARD || defined FE_TOWARDZERO | |
117 | if (dret < ret) | |
118 | ret = (res - 0x1p-1022) * c.x; | |
119 | break; | |
120 | #endif | |
121 | ||
122 | default: | |
123 | break; | |
124 | } | |
c5d5d574 OB |
125 | } |
126 | /* Otherwise (x / ret == ret), either the square root was exact or | |
127 | the division was inexact. */ | |
128 | return ret; | |
129 | } | |
130 | else | |
131 | { | |
132 | if ((k & 0x7ff00000) == 0x7ff00000) | |
133 | return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */ | |
134 | if (x == 0) | |
135 | return x; /* sqrt(+0)=+0, sqrt(-0)=-0 */ | |
136 | if (k < 0) | |
137 | return (x - x) / (x - x); /* sqrt(-ve)=sNaN */ | |
60f435bb | 138 | return 0x1p-256 * __ieee754_sqrt (x * 0x1p512); |
e4d82761 | 139 | } |
f7eac6eb | 140 | } |
0ac5ae23 | 141 | strong_alias (__ieee754_sqrt, __sqrt_finite) |