1 /*-
2 * SPDX-License-Identifier: BSD-2-Clause
3 *
4 * Copyright (c) 2007-2013 Bruce D. Evans
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
9 * are met:
10 * 1. Redistributions of source code must retain the above copyright
11 * notice unmodified, this list of conditions, and the following
12 * disclaimer.
13 * 2. Redistributions in binary form must reproduce the above copyright
14 * notice, this list of conditions and the following disclaimer in the
15 * documentation and/or other materials provided with the distribution.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
18 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
19 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
20 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
21 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
22 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
23 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
24 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
26 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 */
28
29 #include <sys/cdefs.h>
30 /**
31 * Implementation of the natural logarithm of x for 128-bit format.
32 *
33 * First decompose x into its base 2 representation:
34 *
35 * log(x) = log(X * 2**k), where X is in [1, 2)
36 * = log(X) + k * log(2).
37 *
38 * Let X = X_i + e, where X_i is the center of one of the intervals
39 * [-1.0/256, 1.0/256), [1.0/256, 3.0/256), .... [2.0-1.0/256, 2.0+1.0/256)
40 * and X is in this interval. Then
41 *
42 * log(X) = log(X_i + e)
43 * = log(X_i * (1 + e / X_i))
44 * = log(X_i) + log(1 + e / X_i).
45 *
46 * The values log(X_i) are tabulated below. Let d = e / X_i and use
47 *
48 * log(1 + d) = p(d)
49 *
50 * where p(d) = d - 0.5*d*d + ... is a special minimax polynomial of
51 * suitably high degree.
52 *
53 * To get sufficiently small roundoff errors, k * log(2), log(X_i), and
54 * sometimes (if |k| is not large) the first term in p(d) must be evaluated
55 * and added up in extra precision. Extra precision is not needed for the
56 * rest of p(d). In the worst case when k = 0 and log(X_i) is 0, the final
57 * error is controlled mainly by the error in the second term in p(d). The
58 * error in this term itself is at most 0.5 ulps from the d*d operation in
59 * it. The error in this term relative to the first term is thus at most
60 * 0.5 * |-0.5| * |d| < 1.0/1024 ulps. We aim for an accumulated error of
61 * at most twice this at the point of the final rounding step. Thus the
62 * final error should be at most 0.5 + 1.0/512 = 0.5020 ulps. Exhaustive
63 * testing of a float variant of this function showed a maximum final error
64 * of 0.5008 ulps. Non-exhaustive testing of a double variant of this
65 * function showed a maximum final error of 0.5078 ulps (near 1+1.0/256).
66 *
67 * We made the maximum of |d| (and thus the total relative error and the
68 * degree of p(d)) small by using a large number of intervals. Using
69 * centers of intervals instead of endpoints reduces this maximum by a
70 * factor of 2 for a given number of intervals. p(d) is special only
71 * in beginning with the Taylor coefficients 0 + 1*d, which tends to happen
72 * naturally. The most accurate minimax polynomial of a given degree might
73 * be different, but then we wouldn't want it since we would have to do
74 * extra work to avoid roundoff error (especially for P0*d instead of d).
75 */
76
77 #ifdef DEBUG
78 #include <assert.h>
79 #include <fenv.h>
80 #endif
81
82 #include "fpmath.h"
83 #include "math.h"
84 #ifndef NO_STRUCT_RETURN
85 #define STRUCT_RETURN
86 #endif
87 #include "math_private.h"
88
89 #if !defined(NO_UTAB) && !defined(NO_UTABL)
90 #define USE_UTAB
91 #endif
92
93 /*
94 * Domain [-0.005280, 0.004838], range ~[-1.1577e-37, 1.1582e-37]:
95 * |log(1 + d)/d - p(d)| < 2**-122.7
96 */
97 static const long double
98 P2 = -0.5L,
99 P3 = 3.33333333333333333333333333333233795e-1L, /* 0x15555555555555555555555554d42.0p-114L */
100 P4 = -2.49999999999999999999999999941139296e-1L, /* -0x1ffffffffffffffffffffffdab14e.0p-115L */
101 P5 = 2.00000000000000000000000085468039943e-1L, /* 0x19999999999999999999a6d3567f4.0p-115L */
102 P6 = -1.66666666666666666666696142372698408e-1L, /* -0x15555555555555555567267a58e13.0p-115L */
103 P7 = 1.42857142857142857119522943477166120e-1L, /* 0x1249249249249248ed79a0ae434de.0p-115L */
104 P8 = -1.24999999999999994863289015033581301e-1L; /* -0x1fffffffffffffa13e91765e46140.0p-116L */
105 /* Double precision gives ~ 53 + log2(P9 * max(|d|)**8) ~= 120 bits. */
106 static const double
107 P9 = 1.1111111111111401e-1, /* 0x1c71c71c71c7ed.0p-56 */
108 P10 = -1.0000000000040135e-1, /* -0x199999999a0a92.0p-56 */
109 P11 = 9.0909090728136258e-2, /* 0x1745d173962111.0p-56 */
110 P12 = -8.3333318851855284e-2, /* -0x1555551722c7a3.0p-56 */
111 P13 = 7.6928634666404178e-2, /* 0x13b1985204a4ae.0p-56 */
112 P14 = -7.1626810078462499e-2; /* -0x12562276cdc5d0.0p-56 */
113
114 static volatile const double zero = 0;
115
116 #define INTERVALS 128
117 #define LOG2_INTERVALS 7
118 #define TSIZE (INTERVALS + 1)
119 #define G(i) (T[(i)].G)
120 #define F_hi(i) (T[(i)].F_hi)
121 #define F_lo(i) (T[(i)].F_lo)
122 #define ln2_hi F_hi(TSIZE - 1)
123 #define ln2_lo F_lo(TSIZE - 1)
124 #define E(i) (U[(i)].E)
125 #define H(i) (U[(i)].H)
126
127 static const struct {
128 float G; /* 1/(1 + i/128) rounded to 8/9 bits */
129 float F_hi; /* log(1 / G_i) rounded (see below) */
130 /* The compiler will insert 8 bytes of padding here. */
131 long double F_lo; /* next 113 bits for log(1 / G_i) */
132 } T[TSIZE] = {
133 /*
134 * ln2_hi and each F_hi(i) are rounded to a number of bits that
135 * makes F_hi(i) + dk*ln2_hi exact for all i and all dk.
136 *
137 * The last entry (for X just below 2) is used to define ln2_hi
138 * and ln2_lo, to ensure that F_hi(i) and F_lo(i) cancel exactly
139 * with dk*ln2_hi and dk*ln2_lo, respectively, when dk = -1.
140 * This is needed for accuracy when x is just below 1. (To avoid
141 * special cases, such x are "reduced" strangely to X just below
142 * 2 and dk = -1, and then the exact cancellation is needed
143 * because any the error from any non-exactness would be too
144 * large).
145 *
146 * The relevant range of dk is [-16445, 16383]. The maximum number
147 * of bits in F_hi(i) that works is very dependent on i but has
148 * a minimum of 93. We only need about 12 bits in F_hi(i) for
149 * it to provide enough extra precision.
150 *
151 * We round F_hi(i) to 24 bits so that it can have type float,
152 * mainly to minimize the size of the table. Using all 24 bits
153 * in a float for it automatically satisfies the above constraints.
154 */
155 {0x800000.0p-23, 0, 0},
156 {0xfe0000.0p-24, 0x8080ac.0p-30, -0x14ee431dae6674afa0c4bfe16e8fd.0p-144L},
157 {0xfc0000.0p-24, 0x8102b3.0p-29, -0x1db29ee2d83717be918e1119642ab.0p-144L},
158 {0xfa0000.0p-24, 0xc24929.0p-29, 0x1191957d173697cf302cc9476f561.0p-143L},
159 {0xf80000.0p-24, 0x820aec.0p-28, 0x13ce8888e02e78eba9b1113bc1c18.0p-142L},
160 {0xf60000.0p-24, 0xa33577.0p-28, -0x17a4382ce6eb7bfa509bec8da5f22.0p-142L},
161 {0xf48000.0p-24, 0xbc42cb.0p-28, -0x172a21161a107674986dcdca6709c.0p-143L},
162 {0xf30000.0p-24, 0xd57797.0p-28, -0x1e09de07cb958897a3ea46e84abb3.0p-142L},
163 {0xf10000.0p-24, 0xf7518e.0p-28, 0x1ae1eec1b036c484993c549c4bf40.0p-151L},
164 {0xef0000.0p-24, 0x8cb9df.0p-27, -0x1d7355325d560d9e9ab3d6ebab580.0p-141L},
165 {0xed8000.0p-24, 0x999ec0.0p-27, -0x1f9f02d256d5037108f4ec21e48cd.0p-142L},
166 {0xec0000.0p-24, 0xa6988b.0p-27, -0x16fc0a9d12c17a70f7a684c596b12.0p-143L},
167 {0xea0000.0p-24, 0xb80698.0p-27, 0x15d581c1e8da99ded322fb08b8462.0p-141L},
168 {0xe80000.0p-24, 0xc99af3.0p-27, -0x1535b3ba8f150ae09996d7bb4653e.0p-143L},
169 {0xe70000.0p-24, 0xd273b2.0p-27, 0x163786f5251aefe0ded34c8318f52.0p-145L},
170 {0xe50000.0p-24, 0xe442c0.0p-27, 0x1bc4b2368e32d56699c1799a244d4.0p-144L},
171 {0xe38000.0p-24, 0xf1b83f.0p-27, 0x1c6090f684e6766abceccab1d7174.0p-141L},
172 {0xe20000.0p-24, 0xff448a.0p-27, -0x1890aa69ac9f4215f93936b709efb.0p-142L},
173 {0xe08000.0p-24, 0x8673f6.0p-26, 0x1b9985194b6affd511b534b72a28e.0p-140L},
174 {0xdf0000.0p-24, 0x8d515c.0p-26, -0x1dc08d61c6ef1d9b2ef7e68680598.0p-143L},
175 {0xdd8000.0p-24, 0x943a9e.0p-26, -0x1f72a2dac729b3f46662238a9425a.0p-142L},
176 {0xdc0000.0p-24, 0x9b2fe6.0p-26, -0x1fd4dfd3a0afb9691aed4d5e3df94.0p-140L},
177 {0xda8000.0p-24, 0xa2315d.0p-26, -0x11b26121629c46c186384993e1c93.0p-142L},
178 {0xd90000.0p-24, 0xa93f2f.0p-26, 0x1286d633e8e5697dc6a402a56fce1.0p-141L},
179 {0xd78000.0p-24, 0xb05988.0p-26, 0x16128eba9367707ebfa540e45350c.0p-144L},
180 {0xd60000.0p-24, 0xb78094.0p-26, 0x16ead577390d31ef0f4c9d43f79b2.0p-140L},
181 {0xd50000.0p-24, 0xbc4c6c.0p-26, 0x151131ccf7c7b75e7d900b521c48d.0p-141L},
182 {0xd38000.0p-24, 0xc3890a.0p-26, -0x115e2cd714bd06508aeb00d2ae3e9.0p-140L},
183 {0xd20000.0p-24, 0xcad2d7.0p-26, -0x1847f406ebd3af80485c2f409633c.0p-142L},
184 {0xd10000.0p-24, 0xcfb620.0p-26, 0x1c2259904d686581799fbce0b5f19.0p-141L},
185 {0xcf8000.0p-24, 0xd71653.0p-26, 0x1ece57a8d5ae54f550444ecf8b995.0p-140L},
186 {0xce0000.0p-24, 0xde843a.0p-26, -0x1f109d4bc4595412b5d2517aaac13.0p-141L},
187 {0xcd0000.0p-24, 0xe37fde.0p-26, 0x1bc03dc271a74d3a85b5b43c0e727.0p-141L},
188 {0xcb8000.0p-24, 0xeb050c.0p-26, -0x1bf2badc0df841a71b79dd5645b46.0p-145L},
189 {0xca0000.0p-24, 0xf29878.0p-26, -0x18efededd89fbe0bcfbe6d6db9f66.0p-147L},
190 {0xc90000.0p-24, 0xf7ad6f.0p-26, 0x1373ff977baa6911c7bafcb4d84fb.0p-141L},
191 {0xc80000.0p-24, 0xfcc8e3.0p-26, 0x196766f2fb328337cc050c6d83b22.0p-140L},
192 {0xc68000.0p-24, 0x823f30.0p-25, 0x19bd076f7c434e5fcf1a212e2a91e.0p-139L},
193 {0xc58000.0p-24, 0x84d52c.0p-25, -0x1a327257af0f465e5ecab5f2a6f81.0p-139L},
194 {0xc40000.0p-24, 0x88bc74.0p-25, 0x113f23def19c5a0fe396f40f1dda9.0p-141L},
195 {0xc30000.0p-24, 0x8b5ae6.0p-25, 0x1759f6e6b37de945a049a962e66c6.0p-139L},
196 {0xc20000.0p-24, 0x8dfccb.0p-25, 0x1ad35ca6ed5147bdb6ddcaf59c425.0p-141L},
197 {0xc10000.0p-24, 0x90a22b.0p-25, 0x1a1d71a87deba46bae9827221dc98.0p-139L},
198 {0xbf8000.0p-24, 0x94a0d8.0p-25, -0x139e5210c2b730e28aba001a9b5e0.0p-140L},
199 {0xbe8000.0p-24, 0x974f16.0p-25, -0x18f6ebcff3ed72e23e13431adc4a5.0p-141L},
200 {0xbd8000.0p-24, 0x9a00f1.0p-25, -0x1aa268be39aab7148e8d80caa10b7.0p-139L},
201 {0xbc8000.0p-24, 0x9cb672.0p-25, -0x14c8815839c5663663d15faed7771.0p-139L},
202 {0xbb0000.0p-24, 0xa0cda1.0p-25, 0x1eaf46390dbb2438273918db7df5c.0p-141L},
203 {0xba0000.0p-24, 0xa38c6e.0p-25, 0x138e20d831f698298adddd7f32686.0p-141L},
204 {0xb90000.0p-24, 0xa64f05.0p-25, -0x1e8d3c41123615b147a5d47bc208f.0p-142L},
205 {0xb80000.0p-24, 0xa91570.0p-25, 0x1ce28f5f3840b263acb4351104631.0p-140L},
206 {0xb70000.0p-24, 0xabdfbb.0p-25, -0x186e5c0a42423457e22d8c650b355.0p-139L},
207 {0xb60000.0p-24, 0xaeadef.0p-25, -0x14d41a0b2a08a465dc513b13f567d.0p-143L},
208 {0xb50000.0p-24, 0xb18018.0p-25, 0x16755892770633947ffe651e7352f.0p-139L},
209 {0xb40000.0p-24, 0xb45642.0p-25, -0x16395ebe59b15228bfe8798d10ff0.0p-142L},
210 {0xb30000.0p-24, 0xb73077.0p-25, 0x1abc65c8595f088b61a335f5b688c.0p-140L},
211 {0xb20000.0p-24, 0xba0ec4.0p-25, -0x1273089d3dad88e7d353e9967d548.0p-139L},
212 {0xb10000.0p-24, 0xbcf133.0p-25, 0x10f9f67b1f4bbf45de06ecebfaf6d.0p-139L},
213 {0xb00000.0p-24, 0xbfd7d2.0p-25, -0x109fab904864092b34edda19a831e.0p-140L},
214 {0xaf0000.0p-24, 0xc2c2ac.0p-25, -0x1124680aa43333221d8a9b475a6ba.0p-139L},
215 {0xae8000.0p-24, 0xc439b3.0p-25, -0x1f360cc4710fbfe24b633f4e8d84d.0p-140L},
216 {0xad8000.0p-24, 0xc72afd.0p-25, -0x132d91f21d89c89c45003fc5d7807.0p-140L},
217 {0xac8000.0p-24, 0xca20a2.0p-25, -0x16bf9b4d1f8da8002f2449e174504.0p-139L},
218 {0xab8000.0p-24, 0xcd1aae.0p-25, 0x19deb5ce6a6a8717d5626e16acc7d.0p-141L},
219 {0xaa8000.0p-24, 0xd0192f.0p-25, 0x1a29fb48f7d3ca87dabf351aa41f4.0p-139L},
220 {0xaa0000.0p-24, 0xd19a20.0p-25, 0x1127d3c6457f9d79f51dcc73014c9.0p-141L},
221 {0xa90000.0p-24, 0xd49f6a.0p-25, -0x1ba930e486a0ac42d1bf9199188e7.0p-141L},
222 {0xa80000.0p-24, 0xd7a94b.0p-25, -0x1b6e645f31549dd1160bcc45c7e2c.0p-139L},
223 {0xa70000.0p-24, 0xdab7d0.0p-25, 0x1118a425494b610665377f15625b6.0p-140L},
224 {0xa68000.0p-24, 0xdc40d5.0p-25, 0x1966f24d29d3a2d1b2176010478be.0p-140L},
225 {0xa58000.0p-24, 0xdf566d.0p-25, -0x1d8e52eb2248f0c95dd83626d7333.0p-142L},
226 {0xa48000.0p-24, 0xe270ce.0p-25, -0x1ee370f96e6b67ccb006a5b9890ea.0p-140L},
227 {0xa40000.0p-24, 0xe3ffce.0p-25, 0x1d155324911f56db28da4d629d00a.0p-140L},
228 {0xa30000.0p-24, 0xe72179.0p-25, -0x1fe6e2f2f867d8f4d60c713346641.0p-140L},
229 {0xa20000.0p-24, 0xea4812.0p-25, 0x1b7be9add7f4d3b3d406b6cbf3ce5.0p-140L},
230 {0xa18000.0p-24, 0xebdd3d.0p-25, 0x1b3cfb3f7511dd73692609040ccc2.0p-139L},
231 {0xa08000.0p-24, 0xef0b5b.0p-25, -0x1220de1f7301901b8ad85c25afd09.0p-139L},
232 {0xa00000.0p-24, 0xf0a451.0p-25, -0x176364c9ac81cc8a4dfb804de6867.0p-140L},
233 {0x9f0000.0p-24, 0xf3da16.0p-25, 0x1eed6b9aafac8d42f78d3e65d3727.0p-141L},
234 {0x9e8000.0p-24, 0xf576e9.0p-25, 0x1d593218675af269647b783d88999.0p-139L},
235 {0x9d8000.0p-24, 0xf8b47c.0p-25, -0x13e8eb7da053e063714615f7cc91d.0p-144L},
236 {0x9d0000.0p-24, 0xfa553f.0p-25, 0x1c063259bcade02951686d5373aec.0p-139L},
237 {0x9c0000.0p-24, 0xfd9ac5.0p-25, 0x1ef491085fa3c1649349630531502.0p-139L},
238 {0x9b8000.0p-24, 0xff3f8c.0p-25, 0x1d607a7c2b8c5320619fb9433d841.0p-139L},
239 {0x9a8000.0p-24, 0x814697.0p-24, -0x12ad3817004f3f0bdff99f932b273.0p-138L},
240 {0x9a0000.0p-24, 0x821b06.0p-24, -0x189fc53117f9e54e78103a2bc1767.0p-141L},
241 {0x990000.0p-24, 0x83c5f8.0p-24, 0x14cf15a048907b7d7f47ddb45c5a3.0p-139L},
242 {0x988000.0p-24, 0x849c7d.0p-24, 0x1cbb1d35fb82873b04a9af1dd692c.0p-138L},
243 {0x978000.0p-24, 0x864ba6.0p-24, 0x1128639b814f9b9770d8cb6573540.0p-138L},
244 {0x970000.0p-24, 0x87244c.0p-24, 0x184733853300f002e836dfd47bd41.0p-139L},
245 {0x968000.0p-24, 0x87fdaa.0p-24, 0x109d23aef77dd5cd7cc94306fb3ff.0p-140L},
246 {0x958000.0p-24, 0x89b293.0p-24, -0x1a81ef367a59de2b41eeebd550702.0p-138L},
247 {0x950000.0p-24, 0x8a8e20.0p-24, -0x121ad3dbb2f45275c917a30df4ac9.0p-138L},
248 {0x948000.0p-24, 0x8b6a6a.0p-24, -0x1cfb981628af71a89df4e6df2e93b.0p-139L},
249 {0x938000.0p-24, 0x8d253a.0p-24, -0x1d21730ea76cfdec367828734cae5.0p-139L},
250 {0x930000.0p-24, 0x8e03c2.0p-24, 0x135cc00e566f76b87333891e0dec4.0p-138L},
251 {0x928000.0p-24, 0x8ee30d.0p-24, -0x10fcb5df257a263e3bf446c6e3f69.0p-140L},
252 {0x918000.0p-24, 0x90a3ee.0p-24, -0x16e171b15433d723a4c7380a448d8.0p-139L},
253 {0x910000.0p-24, 0x918587.0p-24, -0x1d050da07f3236f330972da2a7a87.0p-139L},
254 {0x908000.0p-24, 0x9267e7.0p-24, 0x1be03669a5268d21148c6002becd3.0p-139L},
255 {0x8f8000.0p-24, 0x942f04.0p-24, 0x10b28e0e26c336af90e00533323ba.0p-139L},
256 {0x8f0000.0p-24, 0x9513c3.0p-24, 0x1a1d820da57cf2f105a89060046aa.0p-138L},
257 {0x8e8000.0p-24, 0x95f950.0p-24, -0x19ef8f13ae3cf162409d8ea99d4c0.0p-139L},
258 {0x8e0000.0p-24, 0x96dfab.0p-24, -0x109e417a6e507b9dc10dac743ad7a.0p-138L},
259 {0x8d0000.0p-24, 0x98aed2.0p-24, 0x10d01a2c5b0e97c4990b23d9ac1f5.0p-139L},
260 {0x8c8000.0p-24, 0x9997a2.0p-24, -0x1d6a50d4b61ea74540bdd2aa99a42.0p-138L},
261 {0x8c0000.0p-24, 0x9a8145.0p-24, 0x1b3b190b83f9527e6aba8f2d783c1.0p-138L},
262 {0x8b8000.0p-24, 0x9b6bbf.0p-24, 0x13a69fad7e7abe7ba81c664c107e0.0p-138L},
263 {0x8b0000.0p-24, 0x9c5711.0p-24, -0x11cd12316f576aad348ae79867223.0p-138L},
264 {0x8a8000.0p-24, 0x9d433b.0p-24, 0x1c95c444b807a246726b304ccae56.0p-139L},
265 {0x898000.0p-24, 0x9f1e22.0p-24, -0x1b9c224ea698c2f9b47466d6123fe.0p-139L},
266 {0x890000.0p-24, 0xa00ce1.0p-24, 0x125ca93186cf0f38b4619a2483399.0p-141L},
267 {0x888000.0p-24, 0xa0fc80.0p-24, -0x1ee38a7bc228b3597043be78eaf49.0p-139L},
268 {0x880000.0p-24, 0xa1ed00.0p-24, -0x1a0db876613d204147dc69a07a649.0p-138L},
269 {0x878000.0p-24, 0xa2de62.0p-24, 0x193224e8516c008d3602a7b41c6e8.0p-139L},
270 {0x870000.0p-24, 0xa3d0a9.0p-24, 0x1fa28b4d2541aca7d5844606b2421.0p-139L},
271 {0x868000.0p-24, 0xa4c3d6.0p-24, 0x1c1b5760fb4571acbcfb03f16daf4.0p-138L},
272 {0x858000.0p-24, 0xa6acea.0p-24, 0x1fed5d0f65949c0a345ad743ae1ae.0p-140L},
273 {0x850000.0p-24, 0xa7a2d4.0p-24, 0x1ad270c9d749362382a7688479e24.0p-140L},
274 {0x848000.0p-24, 0xa899ab.0p-24, 0x199ff15ce532661ea9643a3a2d378.0p-139L},
275 {0x840000.0p-24, 0xa99171.0p-24, 0x1a19e15ccc45d257530a682b80490.0p-139L},
276 {0x838000.0p-24, 0xaa8a28.0p-24, -0x121a14ec532b35ba3e1f868fd0b5e.0p-140L},
277 {0x830000.0p-24, 0xab83d1.0p-24, 0x1aee319980bff3303dd481779df69.0p-139L},
278 {0x828000.0p-24, 0xac7e6f.0p-24, -0x18ffd9e3900345a85d2d86161742e.0p-140L},
279 {0x820000.0p-24, 0xad7a03.0p-24, -0x1e4db102ce29f79b026b64b42caa1.0p-140L},
280 {0x818000.0p-24, 0xae768f.0p-24, 0x17c35c55a04a82ab19f77652d977a.0p-141L},
281 {0x810000.0p-24, 0xaf7415.0p-24, 0x1448324047019b48d7b98c1cf7234.0p-138L},
282 {0x808000.0p-24, 0xb07298.0p-24, -0x1750ee3915a197e9c7359dd94152f.0p-138L},
283 {0x800000.0p-24, 0xb17218.0p-24, -0x105c610ca86c3898cff81a12a17e2.0p-141L},
284 };
285
286 #ifdef USE_UTAB
287 static const struct {
288 float H; /* 1 + i/INTERVALS (exact) */
289 float E; /* H(i) * G(i) - 1 (exact) */
290 } U[TSIZE] = {
291 {0x800000.0p-23, 0},
292 {0x810000.0p-23, -0x800000.0p-37},
293 {0x820000.0p-23, -0x800000.0p-35},
294 {0x830000.0p-23, -0x900000.0p-34},
295 {0x840000.0p-23, -0x800000.0p-33},
296 {0x850000.0p-23, -0xc80000.0p-33},
297 {0x860000.0p-23, -0xa00000.0p-36},
298 {0x870000.0p-23, 0x940000.0p-33},
299 {0x880000.0p-23, 0x800000.0p-35},
300 {0x890000.0p-23, -0xc80000.0p-34},
301 {0x8a0000.0p-23, 0xe00000.0p-36},
302 {0x8b0000.0p-23, 0x900000.0p-33},
303 {0x8c0000.0p-23, -0x800000.0p-35},
304 {0x8d0000.0p-23, -0xe00000.0p-33},
305 {0x8e0000.0p-23, 0x880000.0p-33},
306 {0x8f0000.0p-23, -0xa80000.0p-34},
307 {0x900000.0p-23, -0x800000.0p-35},
308 {0x910000.0p-23, 0x800000.0p-37},
309 {0x920000.0p-23, 0x900000.0p-35},
310 {0x930000.0p-23, 0xd00000.0p-35},
311 {0x940000.0p-23, 0xe00000.0p-35},
312 {0x950000.0p-23, 0xc00000.0p-35},
313 {0x960000.0p-23, 0xe00000.0p-36},
314 {0x970000.0p-23, -0x800000.0p-38},
315 {0x980000.0p-23, -0xc00000.0p-35},
316 {0x990000.0p-23, -0xd00000.0p-34},
317 {0x9a0000.0p-23, 0x880000.0p-33},
318 {0x9b0000.0p-23, 0xe80000.0p-35},
319 {0x9c0000.0p-23, -0x800000.0p-35},
320 {0x9d0000.0p-23, 0xb40000.0p-33},
321 {0x9e0000.0p-23, 0x880000.0p-34},
322 {0x9f0000.0p-23, -0xe00000.0p-35},
323 {0xa00000.0p-23, 0x800000.0p-33},
324 {0xa10000.0p-23, -0x900000.0p-36},
325 {0xa20000.0p-23, -0xb00000.0p-33},
326 {0xa30000.0p-23, -0xa00000.0p-36},
327 {0xa40000.0p-23, 0x800000.0p-33},
328 {0xa50000.0p-23, -0xf80000.0p-35},
329 {0xa60000.0p-23, 0x880000.0p-34},
330 {0xa70000.0p-23, -0x900000.0p-33},
331 {0xa80000.0p-23, -0x800000.0p-35},
332 {0xa90000.0p-23, 0x900000.0p-34},
333 {0xaa0000.0p-23, 0xa80000.0p-33},
334 {0xab0000.0p-23, -0xac0000.0p-34},
335 {0xac0000.0p-23, -0x800000.0p-37},
336 {0xad0000.0p-23, 0xf80000.0p-35},
337 {0xae0000.0p-23, 0xf80000.0p-34},
338 {0xaf0000.0p-23, -0xac0000.0p-33},
339 {0xb00000.0p-23, -0x800000.0p-33},
340 {0xb10000.0p-23, -0xb80000.0p-34},
341 {0xb20000.0p-23, -0x800000.0p-34},
342 {0xb30000.0p-23, -0xb00000.0p-35},
343 {0xb40000.0p-23, -0x800000.0p-35},
344 {0xb50000.0p-23, -0xe00000.0p-36},
345 {0xb60000.0p-23, -0x800000.0p-35},
346 {0xb70000.0p-23, -0xb00000.0p-35},
347 {0xb80000.0p-23, -0x800000.0p-34},
348 {0xb90000.0p-23, -0xb80000.0p-34},
349 {0xba0000.0p-23, -0x800000.0p-33},
350 {0xbb0000.0p-23, -0xac0000.0p-33},
351 {0xbc0000.0p-23, 0x980000.0p-33},
352 {0xbd0000.0p-23, 0xbc0000.0p-34},
353 {0xbe0000.0p-23, 0xe00000.0p-36},
354 {0xbf0000.0p-23, -0xb80000.0p-35},
355 {0xc00000.0p-23, -0x800000.0p-33},
356 {0xc10000.0p-23, 0xa80000.0p-33},
357 {0xc20000.0p-23, 0x900000.0p-34},
358 {0xc30000.0p-23, -0x800000.0p-35},
359 {0xc40000.0p-23, -0x900000.0p-33},
360 {0xc50000.0p-23, 0x820000.0p-33},
361 {0xc60000.0p-23, 0x800000.0p-38},
362 {0xc70000.0p-23, -0x820000.0p-33},
363 {0xc80000.0p-23, 0x800000.0p-33},
364 {0xc90000.0p-23, -0xa00000.0p-36},
365 {0xca0000.0p-23, -0xb00000.0p-33},
366 {0xcb0000.0p-23, 0x840000.0p-34},
367 {0xcc0000.0p-23, -0xd00000.0p-34},
368 {0xcd0000.0p-23, 0x800000.0p-33},
369 {0xce0000.0p-23, -0xe00000.0p-35},
370 {0xcf0000.0p-23, 0xa60000.0p-33},
371 {0xd00000.0p-23, -0x800000.0p-35},
372 {0xd10000.0p-23, 0xb40000.0p-33},
373 {0xd20000.0p-23, -0x800000.0p-35},
374 {0xd30000.0p-23, 0xaa0000.0p-33},
375 {0xd40000.0p-23, -0xe00000.0p-35},
376 {0xd50000.0p-23, 0x880000.0p-33},
377 {0xd60000.0p-23, -0xd00000.0p-34},
378 {0xd70000.0p-23, 0x9c0000.0p-34},
379 {0xd80000.0p-23, -0xb00000.0p-33},
380 {0xd90000.0p-23, -0x800000.0p-38},
381 {0xda0000.0p-23, 0xa40000.0p-33},
382 {0xdb0000.0p-23, -0xdc0000.0p-34},
383 {0xdc0000.0p-23, 0xc00000.0p-35},
384 {0xdd0000.0p-23, 0xca0000.0p-33},
385 {0xde0000.0p-23, -0xb80000.0p-34},
386 {0xdf0000.0p-23, 0xd00000.0p-35},
387 {0xe00000.0p-23, 0xc00000.0p-33},
388 {0xe10000.0p-23, -0xf40000.0p-34},
389 {0xe20000.0p-23, 0x800000.0p-37},
390 {0xe30000.0p-23, 0x860000.0p-33},
391 {0xe40000.0p-23, -0xc80000.0p-33},
392 {0xe50000.0p-23, -0xa80000.0p-34},
393 {0xe60000.0p-23, 0xe00000.0p-36},
394 {0xe70000.0p-23, 0x880000.0p-33},
395 {0xe80000.0p-23, -0xe00000.0p-33},
396 {0xe90000.0p-23, -0xfc0000.0p-34},
397 {0xea0000.0p-23, -0x800000.0p-35},
398 {0xeb0000.0p-23, 0xe80000.0p-35},
399 {0xec0000.0p-23, 0x900000.0p-33},
400 {0xed0000.0p-23, 0xe20000.0p-33},
401 {0xee0000.0p-23, -0xac0000.0p-33},
402 {0xef0000.0p-23, -0xc80000.0p-34},
403 {0xf00000.0p-23, -0x800000.0p-35},
404 {0xf10000.0p-23, 0x800000.0p-35},
405 {0xf20000.0p-23, 0xb80000.0p-34},
406 {0xf30000.0p-23, 0x940000.0p-33},
407 {0xf40000.0p-23, 0xc80000.0p-33},
408 {0xf50000.0p-23, -0xf20000.0p-33},
409 {0xf60000.0p-23, -0xc80000.0p-33},
410 {0xf70000.0p-23, -0xa20000.0p-33},
411 {0xf80000.0p-23, -0x800000.0p-33},
412 {0xf90000.0p-23, -0xc40000.0p-34},
413 {0xfa0000.0p-23, -0x900000.0p-34},
414 {0xfb0000.0p-23, -0xc80000.0p-35},
415 {0xfc0000.0p-23, -0x800000.0p-35},
416 {0xfd0000.0p-23, -0x900000.0p-36},
417 {0xfe0000.0p-23, -0x800000.0p-37},
418 {0xff0000.0p-23, -0x800000.0p-39},
419 {0x800000.0p-22, 0},
420 };
421 #endif /* USE_UTAB */
422
423 #ifdef STRUCT_RETURN
424 #define RETURN1(rp, v) do { \
425 (rp)->hi = (v); \
426 (rp)->lo_set = 0; \
427 return; \
428 } while (0)
429
430 #define RETURN2(rp, h, l) do { \
431 (rp)->hi = (h); \
432 (rp)->lo = (l); \
433 (rp)->lo_set = 1; \
434 return; \
435 } while (0)
436
437 struct ld {
438 long double hi;
439 long double lo;
440 int lo_set;
441 };
442 #else
443 #define RETURN1(rp, v) RETURNF(v)
444 #define RETURN2(rp, h, l) RETURNI((h) + (l))
445 #endif
446
447 #ifdef STRUCT_RETURN
448 static inline __always_inline void
k_logl(long double x,struct ld * rp)449 k_logl(long double x, struct ld *rp)
450 #else
451 long double
452 logl(long double x)
453 #endif
454 {
455 long double d, val_hi, val_lo;
456 double dd, dk;
457 uint64_t lx, llx;
458 int i, k;
459 uint16_t hx;
460
461 EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
462 k = -16383;
463 #if 0 /* Hard to do efficiently. Don't do it until we support all modes. */
464 if (x == 1)
465 RETURN1(rp, 0); /* log(1) = +0 in all rounding modes */
466 #endif
467 if (hx == 0 || hx >= 0x8000) { /* zero, negative or subnormal? */
468 if (((hx & 0x7fff) | lx | llx) == 0)
469 RETURN1(rp, -1 / zero); /* log(+-0) = -Inf */
470 if (hx != 0)
471 /* log(neg or NaN) = qNaN: */
472 RETURN1(rp, (x - x) / zero);
473 x *= 0x1.0p113; /* subnormal; scale up x */
474 EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
475 k = -16383 - 113;
476 } else if (hx >= 0x7fff)
477 RETURN1(rp, x + x); /* log(Inf or NaN) = Inf or qNaN */
478 #ifndef STRUCT_RETURN
479 ENTERI();
480 #endif
481 k += hx;
482 dk = k;
483
484 /* Scale x to be in [1, 2). */
485 SET_LDBL_EXPSIGN(x, 0x3fff);
486
487 /* 0 <= i <= INTERVALS: */
488 #define L2I (49 - LOG2_INTERVALS)
489 i = (lx + (1LL << (L2I - 2))) >> (L2I - 1);
490
491 /*
492 * -0.005280 < d < 0.004838. In particular, the infinite-
493 * precision |d| is <= 2**-7. Rounding of G(i) to 8 bits
494 * ensures that d is representable without extra precision for
495 * this bound on |d| (since when this calculation is expressed
496 * as x*G(i)-1, the multiplication needs as many extra bits as
497 * G(i) has and the subtraction cancels 8 bits). But for
498 * most i (107 cases out of 129), the infinite-precision |d|
499 * is <= 2**-8. G(i) is rounded to 9 bits for such i to give
500 * better accuracy (this works by improving the bound on |d|,
501 * which in turn allows rounding to 9 bits in more cases).
502 * This is only important when the original x is near 1 -- it
503 * lets us avoid using a special method to give the desired
504 * accuracy for such x.
505 */
506 if (0)
507 d = x * G(i) - 1;
508 else {
509 #ifdef USE_UTAB
510 d = (x - H(i)) * G(i) + E(i);
511 #else
512 long double x_hi;
513 double x_lo;
514
515 /*
516 * Split x into x_hi + x_lo to calculate x*G(i)-1 exactly.
517 * G(i) has at most 9 bits, so the splitting point is not
518 * critical.
519 */
520 INSERT_LDBL128_WORDS(x_hi, 0x3fff, lx,
521 llx & 0xffffffffff000000ULL);
522 x_lo = x - x_hi;
523 d = x_hi * G(i) - 1 + x_lo * G(i);
524 #endif
525 }
526
527 /*
528 * Our algorithm depends on exact cancellation of F_lo(i) and
529 * F_hi(i) with dk*ln_2_lo and dk*ln2_hi when k is -1 and i is
530 * at the end of the table. This and other technical complications
531 * make it difficult to avoid the double scaling in (dk*ln2) *
532 * log(base) for base != e without losing more accuracy and/or
533 * efficiency than is gained.
534 */
535 /*
536 * Use double precision operations wherever possible, since
537 * long double operations are emulated and were very slow on
538 * the old sparc64 and unknown on the newer aarch64 and riscv
539 * machines. Also, don't try to improve parallelism by
540 * increasing the number of operations, since any parallelism
541 * on such machines is needed for the emulation. Horner's
542 * method is good for this, and is also good for accuracy.
543 * Horner's method doesn't handle the `lo' term well, either
544 * for efficiency or accuracy. However, for accuracy we
545 * evaluate d * d * P2 separately to take advantage of by P2
546 * being exact, and this gives a good place to sum the 'lo'
547 * term too.
548 */
549 dd = (double)d;
550 val_lo = d * d * d * (P3 +
551 d * (P4 + d * (P5 + d * (P6 + d * (P7 + d * (P8 +
552 dd * (P9 + dd * (P10 + dd * (P11 + dd * (P12 + dd * (P13 +
553 dd * P14))))))))))) + (F_lo(i) + dk * ln2_lo) + d * d * P2;
554 val_hi = d;
555 #ifdef DEBUG
556 if (fetestexcept(FE_UNDERFLOW))
557 breakpoint();
558 #endif
559
560 _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
561 RETURN2(rp, val_hi, val_lo);
562 }
563
564 long double
log1pl(long double x)565 log1pl(long double x)
566 {
567 long double d, d_hi, f_lo, val_hi, val_lo;
568 long double f_hi, twopminusk;
569 double d_lo, dd, dk;
570 uint64_t lx, llx;
571 int i, k;
572 int16_t ax, hx;
573
574 DOPRINT_START(&x);
575 EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
576 if (hx < 0x3fff) { /* x < 1, or x neg NaN */
577 ax = hx & 0x7fff;
578 if (ax >= 0x3fff) { /* x <= -1, or x neg NaN */
579 if (ax == 0x3fff && (lx | llx) == 0)
580 RETURNP(-1 / zero); /* log1p(-1) = -Inf */
581 /* log1p(x < 1, or x NaN) = qNaN: */
582 RETURNP((x - x) / (x - x));
583 }
584 if (ax <= 0x3f8d) { /* |x| < 2**-113 */
585 if ((int)x == 0)
586 RETURNP(x); /* x with inexact if x != 0 */
587 }
588 f_hi = 1;
589 f_lo = x;
590 } else if (hx >= 0x7fff) { /* x +Inf or non-neg NaN */
591 RETURNP(x + x); /* log1p(Inf or NaN) = Inf or qNaN */
592 } else if (hx < 0x40e1) { /* 1 <= x < 2**226 */
593 f_hi = x;
594 f_lo = 1;
595 } else { /* 2**226 <= x < +Inf */
596 f_hi = x;
597 f_lo = 0; /* avoid underflow of the P3 term */
598 }
599 ENTERI();
600 x = f_hi + f_lo;
601 f_lo = (f_hi - x) + f_lo;
602
603 EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
604 k = -16383;
605
606 k += hx;
607 dk = k;
608
609 SET_LDBL_EXPSIGN(x, 0x3fff);
610 twopminusk = 1;
611 SET_LDBL_EXPSIGN(twopminusk, 0x7ffe - (hx & 0x7fff));
612 f_lo *= twopminusk;
613
614 i = (lx + (1LL << (L2I - 2))) >> (L2I - 1);
615
616 /*
617 * x*G(i)-1 (with a reduced x) can be represented exactly, as
618 * above, but now we need to evaluate the polynomial on d =
619 * (x+f_lo)*G(i)-1 and extra precision is needed for that.
620 * Since x+x_lo is a hi+lo decomposition and subtracting 1
621 * doesn't lose too many bits, an inexact calculation for
622 * f_lo*G(i) is good enough.
623 */
624 if (0)
625 d_hi = x * G(i) - 1;
626 else {
627 #ifdef USE_UTAB
628 d_hi = (x - H(i)) * G(i) + E(i);
629 #else
630 long double x_hi;
631 double x_lo;
632
633 INSERT_LDBL128_WORDS(x_hi, 0x3fff, lx,
634 llx & 0xffffffffff000000ULL);
635 x_lo = x - x_hi;
636 d_hi = x_hi * G(i) - 1 + x_lo * G(i);
637 #endif
638 }
639 d_lo = f_lo * G(i);
640
641 /*
642 * This is _2sumF(d_hi, d_lo) inlined. The condition
643 * (d_hi == 0 || |d_hi| >= |d_lo|) for using _2sumF() is not
644 * always satisifed, so it is not clear that this works, but
645 * it works in practice. It works even if it gives a wrong
646 * normalized d_lo, since |d_lo| > |d_hi| implies that i is
647 * nonzero and d is tiny, so the F(i) term dominates d_lo.
648 * In float precision:
649 * (By exhaustive testing, the worst case is d_hi = 0x1.bp-25.
650 * And if d is only a little tinier than that, we would have
651 * another underflow problem for the P3 term; this is also ruled
652 * out by exhaustive testing.)
653 */
654 d = d_hi + d_lo;
655 d_lo = d_hi - d + d_lo;
656 d_hi = d;
657
658 dd = (double)d;
659 val_lo = d * d * d * (P3 +
660 d * (P4 + d * (P5 + d * (P6 + d * (P7 + d * (P8 +
661 dd * (P9 + dd * (P10 + dd * (P11 + dd * (P12 + dd * (P13 +
662 dd * P14))))))))))) + (F_lo(i) + dk * ln2_lo + d_lo) + d * d * P2;
663 val_hi = d_hi;
664 #ifdef DEBUG
665 if (fetestexcept(FE_UNDERFLOW))
666 breakpoint();
667 #endif
668
669 _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
670 RETURN2PI(val_hi, val_lo);
671 }
672
673 #ifdef STRUCT_RETURN
674
675 long double
logl(long double x)676 logl(long double x)
677 {
678 struct ld r;
679
680 ENTERI();
681 DOPRINT_START(&x);
682 k_logl(x, &r);
683 RETURNSPI(&r);
684 }
685
686 /*
687 * 29+113 bit decompositions. The bits are distributed so that the products
688 * of the hi terms are exact in double precision. The types are chosen so
689 * that the products of the hi terms are done in at least double precision,
690 * without any explicit conversions. More natural choices would require a
691 * slow long double precision multiplication.
692 */
693 static const double
694 invln10_hi = 4.3429448176175356e-1, /* 0x1bcb7b15000000.0p-54 */
695 invln2_hi = 1.4426950402557850e0; /* 0x17154765000000.0p-52 */
696 static const long double
697 invln10_lo = 1.41498268538580090791605082294397000e-10L, /* 0x137287195355baaafad33dc323ee3.0p-145L */
698 invln2_lo = 6.33178418956604368501892137426645911e-10L, /* 0x15c17f0bbbe87fed0691d3e88eb57.0p-143L */
699 invln10_lo_plus_hi = invln10_lo + invln10_hi,
700 invln2_lo_plus_hi = invln2_lo + invln2_hi;
701
702 long double
log10l(long double x)703 log10l(long double x)
704 {
705 struct ld r;
706 long double hi, lo;
707
708 ENTERI();
709 DOPRINT_START(&x);
710 k_logl(x, &r);
711 if (!r.lo_set)
712 RETURNPI(r.hi);
713 _2sumF(r.hi, r.lo);
714 hi = (float)r.hi;
715 lo = r.lo + (r.hi - hi);
716 RETURN2PI(invln10_hi * hi,
717 invln10_lo_plus_hi * lo + invln10_lo * hi);
718 }
719
720 long double
log2l(long double x)721 log2l(long double x)
722 {
723 struct ld r;
724 long double hi, lo;
725
726 ENTERI();
727 DOPRINT_START(&x);
728 k_logl(x, &r);
729 if (!r.lo_set)
730 RETURNPI(r.hi);
731 _2sumF(r.hi, r.lo);
732 hi = (float)r.hi;
733 lo = r.lo + (r.hi - hi);
734 RETURN2PI(invln2_hi * hi,
735 invln2_lo_plus_hi * lo + invln2_lo * hi);
736 }
737
738 #endif /* STRUCT_RETURN */
739