1 | /* |
---|
2 | (C) Copyright IBM Corp. 2009 |
---|
3 | |
---|
4 | All rights reserved. |
---|
5 | |
---|
6 | Redistribution and use in source and binary forms, with or without |
---|
7 | modification, are permitted provided that the following conditions are met: |
---|
8 | |
---|
9 | * Redistributions of source code must retain the above copyright notice, |
---|
10 | this list of conditions and the following disclaimer. |
---|
11 | * Redistributions in binary form must reproduce the above copyright |
---|
12 | notice, this list of conditions and the following disclaimer in the |
---|
13 | documentation and/or other materials provided with the distribution. |
---|
14 | * Neither the name of IBM nor the names of its contributors may be |
---|
15 | used to endorse or promote products derived from this software without |
---|
16 | specific prior written permission. |
---|
17 | |
---|
18 | THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" |
---|
19 | AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
---|
20 | IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
---|
21 | ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE |
---|
22 | LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR |
---|
23 | CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF |
---|
24 | SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS |
---|
25 | INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN |
---|
26 | CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) |
---|
27 | ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
---|
28 | POSSIBILITY OF SUCH DAMAGE. |
---|
29 | */ |
---|
30 | |
---|
31 | #include <math.h> |
---|
32 | #include "local.h" |
---|
33 | |
---|
34 | #ifdef _LDBL_EQ_DBL |
---|
35 | /* On platforms where long double is as wide as double. */ |
---|
36 | long double |
---|
37 | sqrtl (long double x) |
---|
38 | { |
---|
39 | return sqrt(x); |
---|
40 | } |
---|
41 | |
---|
42 | #else |
---|
43 | |
---|
44 | /* This code is based upon the version in the BSD math's library. |
---|
45 | That code is... |
---|
46 | * |
---|
47 | * Copyright (c) 2007 Steven G. Kargl |
---|
48 | * All rights reserved. |
---|
49 | * |
---|
50 | * Redistribution and use in source and binary forms, with or without |
---|
51 | * modification, are permitted provided that the following conditions |
---|
52 | * are met: |
---|
53 | * 1. Redistributions of source code must retain the above copyright |
---|
54 | * notice unmodified, this list of conditions, and the following |
---|
55 | * disclaimer. |
---|
56 | * 2. Redistributions in binary form must reproduce the above copyright |
---|
57 | * notice, this list of conditions and the following disclaimer in the |
---|
58 | * documentation and/or other materials provided with the distribution. |
---|
59 | * |
---|
60 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR |
---|
61 | * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES |
---|
62 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. |
---|
63 | * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, |
---|
64 | * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT |
---|
65 | * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, |
---|
66 | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY |
---|
67 | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |
---|
68 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF |
---|
69 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
---|
70 | */ |
---|
71 | #include <float.h> |
---|
72 | #include "ieeefp.h" |
---|
73 | |
---|
74 | #ifndef LDBL_NBIT |
---|
75 | #define LDBL_NBIT 0 |
---|
76 | #endif |
---|
77 | |
---|
78 | #ifndef LDBL_MAX_EXP |
---|
79 | #define LDBL_MAX_EXP DBL_MAX_EXP |
---|
80 | #endif |
---|
81 | |
---|
82 | /* Return (x + ulp) for normal positive x. Assumes no overflow. */ |
---|
83 | |
---|
84 | static inline long double |
---|
85 | inc (long double x) |
---|
86 | { |
---|
87 | union ieee_ext_u ux = { .extu_ld = x, }; |
---|
88 | |
---|
89 | if (++ux.extu_ext.ext_fracl == 0) |
---|
90 | { |
---|
91 | if (++ux.extu_ext.ext_frach == 0) |
---|
92 | { |
---|
93 | ux.extu_ext.ext_exp++; |
---|
94 | ux.extu_ext.ext_frach |= LDBL_NBIT; |
---|
95 | } |
---|
96 | } |
---|
97 | |
---|
98 | return ux.extu_ld; |
---|
99 | } |
---|
100 | |
---|
101 | /* Return (x - ulp) for normal positive x. Assumes no underflow. */ |
---|
102 | |
---|
103 | static inline long double |
---|
104 | dec (long double x) |
---|
105 | { |
---|
106 | union ieee_ext_u ux = { .extu_ld = x, }; |
---|
107 | |
---|
108 | if (ux.extu_ext.ext_fracl-- == 0) |
---|
109 | { |
---|
110 | if (ux.extu_ext.ext_frach-- == LDBL_NBIT) |
---|
111 | { |
---|
112 | ux.extu_ext.ext_exp--; |
---|
113 | ux.extu_ext.ext_frach |= LDBL_NBIT; |
---|
114 | } |
---|
115 | } |
---|
116 | |
---|
117 | return ux.extu_ld; |
---|
118 | } |
---|
119 | |
---|
120 | /* This is slow, but simple and portable. */ |
---|
121 | |
---|
122 | long double |
---|
123 | sqrtl (long double x) |
---|
124 | { |
---|
125 | union ieee_ext_u ux = { .extu_ld = x, }; |
---|
126 | int k, r; |
---|
127 | long double lo, xn; |
---|
128 | |
---|
129 | /* If x = NaN, then sqrt(x) = NaN. */ |
---|
130 | /* If x = Inf, then sqrt(x) = Inf. */ |
---|
131 | /* If x = -Inf, then sqrt(x) = NaN. */ |
---|
132 | if (ux.extu_ext.ext_exp == LDBL_MAX_EXP * 2 - 1) |
---|
133 | return (x * x + x); |
---|
134 | |
---|
135 | /* If x = +-0, then sqrt(x) = +-0. */ |
---|
136 | if (x == 0.0L || x == -0.0L) |
---|
137 | return x; |
---|
138 | |
---|
139 | /* If x < 0, then raise invalid and return NaN. */ |
---|
140 | if (ux.extu_ext.ext_sign) |
---|
141 | return ((x - x) / (x - x)); |
---|
142 | |
---|
143 | if (ux.extu_ext.ext_exp == 0) |
---|
144 | { |
---|
145 | /* Adjust subnormal numbers. */ |
---|
146 | ux.extu_ld *= 0x1.0p514; |
---|
147 | k = -514; |
---|
148 | } |
---|
149 | else |
---|
150 | k = 0; |
---|
151 | |
---|
152 | /* ux.extu_ld is a normal number, so break it into ux.extu_ld = e*2^n where |
---|
153 | ux.extu_ld = (2*e)*2^2k for odd n and ux.extu_ld = (4*e)*2^2k for even n. */ |
---|
154 | |
---|
155 | if ((ux.extu_ext.ext_exp - EXT_EXP_BIAS) & 1) |
---|
156 | { |
---|
157 | /* n is even. */ |
---|
158 | k += ux.extu_ext.ext_exp - EXT_EXP_BIAS - 1; /* 2k = n - 2. */ |
---|
159 | ux.extu_ext.ext_exp = EXT_EXP_BIAS + 1; /* ux.extu_ld in [2,4). */ |
---|
160 | } |
---|
161 | else |
---|
162 | { |
---|
163 | k += ux.extu_ext.ext_exp - EXT_EXP_BIAS; /* 2k = n - 1. */ |
---|
164 | ux.extu_ext.ext_exp = EXT_EXP_BIAS; /* ux.extu_ld in [1,2). */ |
---|
165 | } |
---|
166 | |
---|
167 | /* Newton's iteration. |
---|
168 | Split ux.extu_ld into a high and low part to achieve additional precision. */ |
---|
169 | |
---|
170 | xn = sqrt ((double) ux.extu_ld); /* 53-bit estimate of sqrtl(x). */ |
---|
171 | |
---|
172 | #if LDBL_MANT_DIG > 100 |
---|
173 | xn = (xn + (ux.extu_ld / xn)) * 0.5; /* 106-bit estimate. */ |
---|
174 | #endif |
---|
175 | |
---|
176 | lo = ux.extu_ld; |
---|
177 | ux.extu_ext.ext_fracl = 0; /* Zero out lower bits. */ |
---|
178 | lo = (lo - ux.extu_ld) / xn; /* Low bits divided by xn. */ |
---|
179 | xn = xn + (ux.extu_ld / xn); /* High portion of estimate. */ |
---|
180 | ux.extu_ld = xn + lo; /* Combine everything. */ |
---|
181 | ux.extu_ext.ext_exp += (k >> 1) - 1; |
---|
182 | |
---|
183 | xn = x / ux.extu_ld; /* Chopped quotient (inexact?). */ |
---|
184 | |
---|
185 | /* For simplicity we round to nearest. */ |
---|
186 | xn = inc (xn); /* xn = xn + ulp. */ |
---|
187 | |
---|
188 | ux.extu_ld = ux.extu_ld + xn; /* Chopped sum. */ |
---|
189 | ux.extu_ext.ext_exp--; |
---|
190 | |
---|
191 | return ux.extu_ld; |
---|
192 | } |
---|
193 | #endif /* ! _LDBL_EQ_DBL */ |
---|
194 | |
---|
195 | |
---|