1 1.1 christos /*- 2 1.1 christos * SPDX-License-Identifier: BSD-2-Clause 3 1.1 christos * 4 1.1 christos * Copyright (c) 2007-2013 Bruce D. Evans 5 1.1 christos * All rights reserved. 6 1.1 christos * 7 1.1 christos * Redistribution and use in source and binary forms, with or without 8 1.1 christos * modification, are permitted provided that the following conditions 9 1.1 christos * are met: 10 1.1 christos * 1. Redistributions of source code must retain the above copyright 11 1.1 christos * notice unmodified, this list of conditions, and the following 12 1.1 christos * disclaimer. 13 1.1 christos * 2. Redistributions in binary form must reproduce the above copyright 14 1.1 christos * notice, this list of conditions and the following disclaimer in the 15 1.1 christos * documentation and/or other materials provided with the distribution. 16 1.1 christos * 17 1.1 christos * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 18 1.1 christos * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 19 1.1 christos * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 20 1.1 christos * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 21 1.1 christos * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 22 1.1 christos * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 23 1.1 christos * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 24 1.1 christos * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 25 1.1 christos * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 26 1.1 christos * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 27 1.1 christos */ 28 1.1 christos 29 1.1 christos #include <sys/cdefs.h> 30 1.1 christos /** 31 1.1 christos * Implementation of the natural logarithm of x for Intel 80-bit format. 32 1.1 christos * 33 1.1 christos * First decompose x into its base 2 representation: 34 1.1 christos * 35 1.1 christos * log(x) = log(X * 2**k), where X is in [1, 2) 36 1.1 christos * = log(X) + k * log(2). 37 1.1 christos * 38 1.1 christos * Let X = X_i + e, where X_i is the center of one of the intervals 39 1.1 christos * [-1.0/256, 1.0/256), [1.0/256, 3.0/256), .... [2.0-1.0/256, 2.0+1.0/256) 40 1.1 christos * and X is in this interval. Then 41 1.1 christos * 42 1.1 christos * log(X) = log(X_i + e) 43 1.1 christos * = log(X_i * (1 + e / X_i)) 44 1.1 christos * = log(X_i) + log(1 + e / X_i). 45 1.1 christos * 46 1.1 christos * The values log(X_i) are tabulated below. Let d = e / X_i and use 47 1.1 christos * 48 1.1 christos * log(1 + d) = p(d) 49 1.1 christos * 50 1.1 christos * where p(d) = d - 0.5*d*d + ... is a special minimax polynomial of 51 1.1 christos * suitably high degree. 52 1.1 christos * 53 1.1 christos * To get sufficiently small roundoff errors, k * log(2), log(X_i), and 54 1.1 christos * sometimes (if |k| is not large) the first term in p(d) must be evaluated 55 1.1 christos * and added up in extra precision. Extra precision is not needed for the 56 1.1 christos * rest of p(d). In the worst case when k = 0 and log(X_i) is 0, the final 57 1.1 christos * error is controlled mainly by the error in the second term in p(d). The 58 1.1 christos * error in this term itself is at most 0.5 ulps from the d*d operation in 59 1.1 christos * it. The error in this term relative to the first term is thus at most 60 1.1 christos * 0.5 * |-0.5| * |d| < 1.0/1024 ulps. We aim for an accumulated error of 61 1.1 christos * at most twice this at the point of the final rounding step. Thus the 62 1.1 christos * final error should be at most 0.5 + 1.0/512 = 0.5020 ulps. Exhaustive 63 1.1 christos * testing of a float variant of this function showed a maximum final error 64 1.1 christos * of 0.5008 ulps. Non-exhaustive testing of a double variant of this 65 1.1 christos * function showed a maximum final error of 0.5078 ulps (near 1+1.0/256). 66 1.1 christos * 67 1.1 christos * We made the maximum of |d| (and thus the total relative error and the 68 1.1 christos * degree of p(d)) small by using a large number of intervals. Using 69 1.1 christos * centers of intervals instead of endpoints reduces this maximum by a 70 1.1 christos * factor of 2 for a given number of intervals. p(d) is special only 71 1.1 christos * in beginning with the Taylor coefficients 0 + 1*d, which tends to happen 72 1.1 christos * naturally. The most accurate minimax polynomial of a given degree might 73 1.1 christos * be different, but then we wouldn't want it since we would have to do 74 1.1 christos * extra work to avoid roundoff error (especially for P0*d instead of d). 75 1.1 christos */ 76 1.1 christos 77 1.1 christos #ifdef DEBUG 78 1.1 christos #include <fenv.h> 79 1.1 christos #endif 80 1.1 christos 81 1.1 christos #ifdef __FreeBSD__ 82 1.1 christos #include "fpmath.h" 83 1.1 christos #endif 84 1.1 christos #include "math.h" 85 1.1 christos #define i386_SSE_GOOD 86 1.1 christos #ifndef NO_STRUCT_RETURN 87 1.1 christos #define STRUCT_RETURN 88 1.1 christos #endif 89 1.1 christos #include "math_private.h" 90 1.1 christos 91 1.1 christos #if !defined(NO_UTAB) && !defined(NO_UTABL) 92 1.1 christos #define USE_UTAB 93 1.1 christos #endif 94 1.1 christos 95 1.1 christos /* 96 1.1 christos * Domain [-0.005280, 0.004838], range ~[-5.1736e-22, 5.1738e-22]: 97 1.1 christos * |log(1 + d)/d - p(d)| < 2**-70.7 98 1.1 christos */ 99 1.1 christos static const double 100 1.1 christos P2 = -0.5, 101 1.1 christos P3 = 3.3333333333333359e-1, /* 0x1555555555555a.0p-54 */ 102 1.1 christos P4 = -2.5000000000004424e-1, /* -0x1000000000031d.0p-54 */ 103 1.1 christos P5 = 1.9999999992970016e-1, /* 0x1999999972f3c7.0p-55 */ 104 1.1 christos P6 = -1.6666666072191585e-1, /* -0x15555548912c09.0p-55 */ 105 1.1 christos P7 = 1.4286227413310518e-1, /* 0x12494f9d9def91.0p-55 */ 106 1.1 christos P8 = -1.2518388626763144e-1; /* -0x1006068cc0b97c.0p-55 */ 107 1.1 christos 108 1.1 christos static volatile const double zero = 0; 109 1.1 christos 110 1.1 christos #define INTERVALS 128 111 1.1 christos #define LOG2_INTERVALS 7 112 1.1 christos #define TSIZE (INTERVALS + 1) 113 1.1 christos #define G(i) (T[(i)].G) 114 1.1 christos #define F_hi(i) (T[(i)].F_hi) 115 1.1 christos #define F_lo(i) (T[(i)].F_lo) 116 1.1 christos #define ln2_hi F_hi(TSIZE - 1) 117 1.1 christos #define ln2_lo F_lo(TSIZE - 1) 118 1.1 christos #define E(i) (U[(i)].E) 119 1.1 christos #define H(i) (U[(i)].H) 120 1.1 christos 121 1.1 christos static const struct { 122 1.1 christos float G; /* 1/(1 + i/128) rounded to 8/9 bits */ 123 1.1 christos float F_hi; /* log(1 / G_i) rounded (see below) */ 124 1.1 christos double F_lo; /* next 53 bits for log(1 / G_i) */ 125 1.1 christos } T[TSIZE] = { 126 1.1 christos /* 127 1.1 christos * ln2_hi and each F_hi(i) are rounded to a number of bits that 128 1.1 christos * makes F_hi(i) + dk*ln2_hi exact for all i and all dk. 129 1.1 christos * 130 1.1 christos * The last entry (for X just below 2) is used to define ln2_hi 131 1.1 christos * and ln2_lo, to ensure that F_hi(i) and F_lo(i) cancel exactly 132 1.1 christos * with dk*ln2_hi and dk*ln2_lo, respectively, when dk = -1. 133 1.1 christos * This is needed for accuracy when x is just below 1. (To avoid 134 1.1 christos * special cases, such x are "reduced" strangely to X just below 135 1.1 christos * 2 and dk = -1, and then the exact cancellation is needed 136 1.1 christos * because any the error from any non-exactness would be too 137 1.1 christos * large). 138 1.1 christos * 139 1.1 christos * We want to share this table between double precision and ld80, 140 1.1 christos * so the relevant range of dk is the larger one of ld80 141 1.1 christos * ([-16445, 16383]) and the relevant exactness requirement is 142 1.1 christos * the stricter one of double precision. The maximum number of 143 1.1 christos * bits in F_hi(i) that works is very dependent on i but has 144 1.1 christos * a minimum of 33. We only need about 12 bits in F_hi(i) for 145 1.1 christos * it to provide enough extra precision in double precision (11 146 1.1 christos * more than that are required for ld80). 147 1.1 christos * 148 1.1 christos * We round F_hi(i) to 24 bits so that it can have type float, 149 1.1 christos * mainly to minimize the size of the table. Using all 24 bits 150 1.1 christos * in a float for it automatically satisfies the above constraints. 151 1.1 christos */ 152 1.1 christos { 0x800000.0p-23, 0, 0 }, 153 1.1 christos { 0xfe0000.0p-24, 0x8080ac.0p-30, -0x14ee431dae6675.0p-84 }, 154 1.1 christos { 0xfc0000.0p-24, 0x8102b3.0p-29, -0x1db29ee2d83718.0p-84 }, 155 1.1 christos { 0xfa0000.0p-24, 0xc24929.0p-29, 0x1191957d173698.0p-83 }, 156 1.1 christos { 0xf80000.0p-24, 0x820aec.0p-28, 0x13ce8888e02e79.0p-82 }, 157 1.1 christos { 0xf60000.0p-24, 0xa33577.0p-28, -0x17a4382ce6eb7c.0p-82 }, 158 1.1 christos { 0xf48000.0p-24, 0xbc42cb.0p-28, -0x172a21161a1076.0p-83 }, 159 1.1 christos { 0xf30000.0p-24, 0xd57797.0p-28, -0x1e09de07cb9589.0p-82 }, 160 1.1 christos { 0xf10000.0p-24, 0xf7518e.0p-28, 0x1ae1eec1b036c5.0p-91 }, 161 1.1 christos { 0xef0000.0p-24, 0x8cb9df.0p-27, -0x1d7355325d560e.0p-81 }, 162 1.1 christos { 0xed8000.0p-24, 0x999ec0.0p-27, -0x1f9f02d256d503.0p-82 }, 163 1.1 christos { 0xec0000.0p-24, 0xa6988b.0p-27, -0x16fc0a9d12c17a.0p-83 }, 164 1.1 christos { 0xea0000.0p-24, 0xb80698.0p-27, 0x15d581c1e8da9a.0p-81 }, 165 1.1 christos { 0xe80000.0p-24, 0xc99af3.0p-27, -0x1535b3ba8f150b.0p-83 }, 166 1.1 christos { 0xe70000.0p-24, 0xd273b2.0p-27, 0x163786f5251af0.0p-85 }, 167 1.1 christos { 0xe50000.0p-24, 0xe442c0.0p-27, 0x1bc4b2368e32d5.0p-84 }, 168 1.1 christos { 0xe38000.0p-24, 0xf1b83f.0p-27, 0x1c6090f684e676.0p-81 }, 169 1.1 christos { 0xe20000.0p-24, 0xff448a.0p-27, -0x1890aa69ac9f42.0p-82 }, 170 1.1 christos { 0xe08000.0p-24, 0x8673f6.0p-26, 0x1b9985194b6b00.0p-80 }, 171 1.1 christos { 0xdf0000.0p-24, 0x8d515c.0p-26, -0x1dc08d61c6ef1e.0p-83 }, 172 1.1 christos { 0xdd8000.0p-24, 0x943a9e.0p-26, -0x1f72a2dac729b4.0p-82 }, 173 1.1 christos { 0xdc0000.0p-24, 0x9b2fe6.0p-26, -0x1fd4dfd3a0afb9.0p-80 }, 174 1.1 christos { 0xda8000.0p-24, 0xa2315d.0p-26, -0x11b26121629c47.0p-82 }, 175 1.1 christos { 0xd90000.0p-24, 0xa93f2f.0p-26, 0x1286d633e8e569.0p-81 }, 176 1.1 christos { 0xd78000.0p-24, 0xb05988.0p-26, 0x16128eba936770.0p-84 }, 177 1.1 christos { 0xd60000.0p-24, 0xb78094.0p-26, 0x16ead577390d32.0p-80 }, 178 1.1 christos { 0xd50000.0p-24, 0xbc4c6c.0p-26, 0x151131ccf7c7b7.0p-81 }, 179 1.1 christos { 0xd38000.0p-24, 0xc3890a.0p-26, -0x115e2cd714bd06.0p-80 }, 180 1.1 christos { 0xd20000.0p-24, 0xcad2d7.0p-26, -0x1847f406ebd3b0.0p-82 }, 181 1.1 christos { 0xd10000.0p-24, 0xcfb620.0p-26, 0x1c2259904d6866.0p-81 }, 182 1.1 christos { 0xcf8000.0p-24, 0xd71653.0p-26, 0x1ece57a8d5ae55.0p-80 }, 183 1.1 christos { 0xce0000.0p-24, 0xde843a.0p-26, -0x1f109d4bc45954.0p-81 }, 184 1.1 christos { 0xcd0000.0p-24, 0xe37fde.0p-26, 0x1bc03dc271a74d.0p-81 }, 185 1.1 christos { 0xcb8000.0p-24, 0xeb050c.0p-26, -0x1bf2badc0df842.0p-85 }, 186 1.1 christos { 0xca0000.0p-24, 0xf29878.0p-26, -0x18efededd89fbe.0p-87 }, 187 1.1 christos { 0xc90000.0p-24, 0xf7ad6f.0p-26, 0x1373ff977baa69.0p-81 }, 188 1.1 christos { 0xc80000.0p-24, 0xfcc8e3.0p-26, 0x196766f2fb3283.0p-80 }, 189 1.1 christos { 0xc68000.0p-24, 0x823f30.0p-25, 0x19bd076f7c434e.0p-79 }, 190 1.1 christos { 0xc58000.0p-24, 0x84d52c.0p-25, -0x1a327257af0f46.0p-79 }, 191 1.1 christos { 0xc40000.0p-24, 0x88bc74.0p-25, 0x113f23def19c5a.0p-81 }, 192 1.1 christos { 0xc30000.0p-24, 0x8b5ae6.0p-25, 0x1759f6e6b37de9.0p-79 }, 193 1.1 christos { 0xc20000.0p-24, 0x8dfccb.0p-25, 0x1ad35ca6ed5148.0p-81 }, 194 1.1 christos { 0xc10000.0p-24, 0x90a22b.0p-25, 0x1a1d71a87deba4.0p-79 }, 195 1.1 christos { 0xbf8000.0p-24, 0x94a0d8.0p-25, -0x139e5210c2b731.0p-80 }, 196 1.1 christos { 0xbe8000.0p-24, 0x974f16.0p-25, -0x18f6ebcff3ed73.0p-81 }, 197 1.1 christos { 0xbd8000.0p-24, 0x9a00f1.0p-25, -0x1aa268be39aab7.0p-79 }, 198 1.1 christos { 0xbc8000.0p-24, 0x9cb672.0p-25, -0x14c8815839c566.0p-79 }, 199 1.1 christos { 0xbb0000.0p-24, 0xa0cda1.0p-25, 0x1eaf46390dbb24.0p-81 }, 200 1.1 christos { 0xba0000.0p-24, 0xa38c6e.0p-25, 0x138e20d831f698.0p-81 }, 201 1.1 christos { 0xb90000.0p-24, 0xa64f05.0p-25, -0x1e8d3c41123616.0p-82 }, 202 1.1 christos { 0xb80000.0p-24, 0xa91570.0p-25, 0x1ce28f5f3840b2.0p-80 }, 203 1.1 christos { 0xb70000.0p-24, 0xabdfbb.0p-25, -0x186e5c0a424234.0p-79 }, 204 1.1 christos { 0xb60000.0p-24, 0xaeadef.0p-25, -0x14d41a0b2a08a4.0p-83 }, 205 1.1 christos { 0xb50000.0p-24, 0xb18018.0p-25, 0x16755892770634.0p-79 }, 206 1.1 christos { 0xb40000.0p-24, 0xb45642.0p-25, -0x16395ebe59b152.0p-82 }, 207 1.1 christos { 0xb30000.0p-24, 0xb73077.0p-25, 0x1abc65c8595f09.0p-80 }, 208 1.1 christos { 0xb20000.0p-24, 0xba0ec4.0p-25, -0x1273089d3dad89.0p-79 }, 209 1.1 christos { 0xb10000.0p-24, 0xbcf133.0p-25, 0x10f9f67b1f4bbf.0p-79 }, 210 1.1 christos { 0xb00000.0p-24, 0xbfd7d2.0p-25, -0x109fab90486409.0p-80 }, 211 1.1 christos { 0xaf0000.0p-24, 0xc2c2ac.0p-25, -0x1124680aa43333.0p-79 }, 212 1.1 christos { 0xae8000.0p-24, 0xc439b3.0p-25, -0x1f360cc4710fc0.0p-80 }, 213 1.1 christos { 0xad8000.0p-24, 0xc72afd.0p-25, -0x132d91f21d89c9.0p-80 }, 214 1.1 christos { 0xac8000.0p-24, 0xca20a2.0p-25, -0x16bf9b4d1f8da8.0p-79 }, 215 1.1 christos { 0xab8000.0p-24, 0xcd1aae.0p-25, 0x19deb5ce6a6a87.0p-81 }, 216 1.1 christos { 0xaa8000.0p-24, 0xd0192f.0p-25, 0x1a29fb48f7d3cb.0p-79 }, 217 1.1 christos { 0xaa0000.0p-24, 0xd19a20.0p-25, 0x1127d3c6457f9d.0p-81 }, 218 1.1 christos { 0xa90000.0p-24, 0xd49f6a.0p-25, -0x1ba930e486a0ac.0p-81 }, 219 1.1 christos { 0xa80000.0p-24, 0xd7a94b.0p-25, -0x1b6e645f31549e.0p-79 }, 220 1.1 christos { 0xa70000.0p-24, 0xdab7d0.0p-25, 0x1118a425494b61.0p-80 }, 221 1.1 christos { 0xa68000.0p-24, 0xdc40d5.0p-25, 0x1966f24d29d3a3.0p-80 }, 222 1.1 christos { 0xa58000.0p-24, 0xdf566d.0p-25, -0x1d8e52eb2248f1.0p-82 }, 223 1.1 christos { 0xa48000.0p-24, 0xe270ce.0p-25, -0x1ee370f96e6b68.0p-80 }, 224 1.1 christos { 0xa40000.0p-24, 0xe3ffce.0p-25, 0x1d155324911f57.0p-80 }, 225 1.1 christos { 0xa30000.0p-24, 0xe72179.0p-25, -0x1fe6e2f2f867d9.0p-80 }, 226 1.1 christos { 0xa20000.0p-24, 0xea4812.0p-25, 0x1b7be9add7f4d4.0p-80 }, 227 1.1 christos { 0xa18000.0p-24, 0xebdd3d.0p-25, 0x1b3cfb3f7511dd.0p-79 }, 228 1.1 christos { 0xa08000.0p-24, 0xef0b5b.0p-25, -0x1220de1f730190.0p-79 }, 229 1.1 christos { 0xa00000.0p-24, 0xf0a451.0p-25, -0x176364c9ac81cd.0p-80 }, 230 1.1 christos { 0x9f0000.0p-24, 0xf3da16.0p-25, 0x1eed6b9aafac8d.0p-81 }, 231 1.1 christos { 0x9e8000.0p-24, 0xf576e9.0p-25, 0x1d593218675af2.0p-79 }, 232 1.1 christos { 0x9d8000.0p-24, 0xf8b47c.0p-25, -0x13e8eb7da053e0.0p-84 }, 233 1.1 christos { 0x9d0000.0p-24, 0xfa553f.0p-25, 0x1c063259bcade0.0p-79 }, 234 1.1 christos { 0x9c0000.0p-24, 0xfd9ac5.0p-25, 0x1ef491085fa3c1.0p-79 }, 235 1.1 christos { 0x9b8000.0p-24, 0xff3f8c.0p-25, 0x1d607a7c2b8c53.0p-79 }, 236 1.1 christos { 0x9a8000.0p-24, 0x814697.0p-24, -0x12ad3817004f3f.0p-78 }, 237 1.1 christos { 0x9a0000.0p-24, 0x821b06.0p-24, -0x189fc53117f9e5.0p-81 }, 238 1.1 christos { 0x990000.0p-24, 0x83c5f8.0p-24, 0x14cf15a048907b.0p-79 }, 239 1.1 christos { 0x988000.0p-24, 0x849c7d.0p-24, 0x1cbb1d35fb8287.0p-78 }, 240 1.1 christos { 0x978000.0p-24, 0x864ba6.0p-24, 0x1128639b814f9c.0p-78 }, 241 1.1 christos { 0x970000.0p-24, 0x87244c.0p-24, 0x184733853300f0.0p-79 }, 242 1.1 christos { 0x968000.0p-24, 0x87fdaa.0p-24, 0x109d23aef77dd6.0p-80 }, 243 1.1 christos { 0x958000.0p-24, 0x89b293.0p-24, -0x1a81ef367a59de.0p-78 }, 244 1.1 christos { 0x950000.0p-24, 0x8a8e20.0p-24, -0x121ad3dbb2f452.0p-78 }, 245 1.1 christos { 0x948000.0p-24, 0x8b6a6a.0p-24, -0x1cfb981628af72.0p-79 }, 246 1.1 christos { 0x938000.0p-24, 0x8d253a.0p-24, -0x1d21730ea76cfe.0p-79 }, 247 1.1 christos { 0x930000.0p-24, 0x8e03c2.0p-24, 0x135cc00e566f77.0p-78 }, 248 1.1 christos { 0x928000.0p-24, 0x8ee30d.0p-24, -0x10fcb5df257a26.0p-80 }, 249 1.1 christos { 0x918000.0p-24, 0x90a3ee.0p-24, -0x16e171b15433d7.0p-79 }, 250 1.1 christos { 0x910000.0p-24, 0x918587.0p-24, -0x1d050da07f3237.0p-79 }, 251 1.1 christos { 0x908000.0p-24, 0x9267e7.0p-24, 0x1be03669a5268d.0p-79 }, 252 1.1 christos { 0x8f8000.0p-24, 0x942f04.0p-24, 0x10b28e0e26c337.0p-79 }, 253 1.1 christos { 0x8f0000.0p-24, 0x9513c3.0p-24, 0x1a1d820da57cf3.0p-78 }, 254 1.1 christos { 0x8e8000.0p-24, 0x95f950.0p-24, -0x19ef8f13ae3cf1.0p-79 }, 255 1.1 christos { 0x8e0000.0p-24, 0x96dfab.0p-24, -0x109e417a6e507c.0p-78 }, 256 1.1 christos { 0x8d0000.0p-24, 0x98aed2.0p-24, 0x10d01a2c5b0e98.0p-79 }, 257 1.1 christos { 0x8c8000.0p-24, 0x9997a2.0p-24, -0x1d6a50d4b61ea7.0p-78 }, 258 1.1 christos { 0x8c0000.0p-24, 0x9a8145.0p-24, 0x1b3b190b83f952.0p-78 }, 259 1.1 christos { 0x8b8000.0p-24, 0x9b6bbf.0p-24, 0x13a69fad7e7abe.0p-78 }, 260 1.1 christos { 0x8b0000.0p-24, 0x9c5711.0p-24, -0x11cd12316f576b.0p-78 }, 261 1.1 christos { 0x8a8000.0p-24, 0x9d433b.0p-24, 0x1c95c444b807a2.0p-79 }, 262 1.1 christos { 0x898000.0p-24, 0x9f1e22.0p-24, -0x1b9c224ea698c3.0p-79 }, 263 1.1 christos { 0x890000.0p-24, 0xa00ce1.0p-24, 0x125ca93186cf0f.0p-81 }, 264 1.1 christos { 0x888000.0p-24, 0xa0fc80.0p-24, -0x1ee38a7bc228b3.0p-79 }, 265 1.1 christos { 0x880000.0p-24, 0xa1ed00.0p-24, -0x1a0db876613d20.0p-78 }, 266 1.1 christos { 0x878000.0p-24, 0xa2de62.0p-24, 0x193224e8516c01.0p-79 }, 267 1.1 christos { 0x870000.0p-24, 0xa3d0a9.0p-24, 0x1fa28b4d2541ad.0p-79 }, 268 1.1 christos { 0x868000.0p-24, 0xa4c3d6.0p-24, 0x1c1b5760fb4572.0p-78 }, 269 1.1 christos { 0x858000.0p-24, 0xa6acea.0p-24, 0x1fed5d0f65949c.0p-80 }, 270 1.1 christos { 0x850000.0p-24, 0xa7a2d4.0p-24, 0x1ad270c9d74936.0p-80 }, 271 1.1 christos { 0x848000.0p-24, 0xa899ab.0p-24, 0x199ff15ce53266.0p-79 }, 272 1.1 christos { 0x840000.0p-24, 0xa99171.0p-24, 0x1a19e15ccc45d2.0p-79 }, 273 1.1 christos { 0x838000.0p-24, 0xaa8a28.0p-24, -0x121a14ec532b36.0p-80 }, 274 1.1 christos { 0x830000.0p-24, 0xab83d1.0p-24, 0x1aee319980bff3.0p-79 }, 275 1.1 christos { 0x828000.0p-24, 0xac7e6f.0p-24, -0x18ffd9e3900346.0p-80 }, 276 1.1 christos { 0x820000.0p-24, 0xad7a03.0p-24, -0x1e4db102ce29f8.0p-80 }, 277 1.1 christos { 0x818000.0p-24, 0xae768f.0p-24, 0x17c35c55a04a83.0p-81 }, 278 1.1 christos { 0x810000.0p-24, 0xaf7415.0p-24, 0x1448324047019b.0p-78 }, 279 1.1 christos { 0x808000.0p-24, 0xb07298.0p-24, -0x1750ee3915a198.0p-78 }, 280 1.1 christos { 0x800000.0p-24, 0xb17218.0p-24, -0x105c610ca86c39.0p-81 }, 281 1.1 christos }; 282 1.1 christos 283 1.1 christos #ifdef USE_UTAB 284 1.1 christos static const struct { 285 1.1 christos float H; /* 1 + i/INTERVALS (exact) */ 286 1.1 christos float E; /* H(i) * G(i) - 1 (exact) */ 287 1.1 christos } U[TSIZE] = { 288 1.1 christos { 0x800000.0p-23, 0 }, 289 1.1 christos { 0x810000.0p-23, -0x800000.0p-37 }, 290 1.1 christos { 0x820000.0p-23, -0x800000.0p-35 }, 291 1.1 christos { 0x830000.0p-23, -0x900000.0p-34 }, 292 1.1 christos { 0x840000.0p-23, -0x800000.0p-33 }, 293 1.1 christos { 0x850000.0p-23, -0xc80000.0p-33 }, 294 1.1 christos { 0x860000.0p-23, -0xa00000.0p-36 }, 295 1.1 christos { 0x870000.0p-23, 0x940000.0p-33 }, 296 1.1 christos { 0x880000.0p-23, 0x800000.0p-35 }, 297 1.1 christos { 0x890000.0p-23, -0xc80000.0p-34 }, 298 1.1 christos { 0x8a0000.0p-23, 0xe00000.0p-36 }, 299 1.1 christos { 0x8b0000.0p-23, 0x900000.0p-33 }, 300 1.1 christos { 0x8c0000.0p-23, -0x800000.0p-35 }, 301 1.1 christos { 0x8d0000.0p-23, -0xe00000.0p-33 }, 302 1.1 christos { 0x8e0000.0p-23, 0x880000.0p-33 }, 303 1.1 christos { 0x8f0000.0p-23, -0xa80000.0p-34 }, 304 1.1 christos { 0x900000.0p-23, -0x800000.0p-35 }, 305 1.1 christos { 0x910000.0p-23, 0x800000.0p-37 }, 306 1.1 christos { 0x920000.0p-23, 0x900000.0p-35 }, 307 1.1 christos { 0x930000.0p-23, 0xd00000.0p-35 }, 308 1.1 christos { 0x940000.0p-23, 0xe00000.0p-35 }, 309 1.1 christos { 0x950000.0p-23, 0xc00000.0p-35 }, 310 1.1 christos { 0x960000.0p-23, 0xe00000.0p-36 }, 311 1.1 christos { 0x970000.0p-23, -0x800000.0p-38 }, 312 1.1 christos { 0x980000.0p-23, -0xc00000.0p-35 }, 313 1.1 christos { 0x990000.0p-23, -0xd00000.0p-34 }, 314 1.1 christos { 0x9a0000.0p-23, 0x880000.0p-33 }, 315 1.1 christos { 0x9b0000.0p-23, 0xe80000.0p-35 }, 316 1.1 christos { 0x9c0000.0p-23, -0x800000.0p-35 }, 317 1.1 christos { 0x9d0000.0p-23, 0xb40000.0p-33 }, 318 1.1 christos { 0x9e0000.0p-23, 0x880000.0p-34 }, 319 1.1 christos { 0x9f0000.0p-23, -0xe00000.0p-35 }, 320 1.1 christos { 0xa00000.0p-23, 0x800000.0p-33 }, 321 1.1 christos { 0xa10000.0p-23, -0x900000.0p-36 }, 322 1.1 christos { 0xa20000.0p-23, -0xb00000.0p-33 }, 323 1.1 christos { 0xa30000.0p-23, -0xa00000.0p-36 }, 324 1.1 christos { 0xa40000.0p-23, 0x800000.0p-33 }, 325 1.1 christos { 0xa50000.0p-23, -0xf80000.0p-35 }, 326 1.1 christos { 0xa60000.0p-23, 0x880000.0p-34 }, 327 1.1 christos { 0xa70000.0p-23, -0x900000.0p-33 }, 328 1.1 christos { 0xa80000.0p-23, -0x800000.0p-35 }, 329 1.1 christos { 0xa90000.0p-23, 0x900000.0p-34 }, 330 1.1 christos { 0xaa0000.0p-23, 0xa80000.0p-33 }, 331 1.1 christos { 0xab0000.0p-23, -0xac0000.0p-34 }, 332 1.1 christos { 0xac0000.0p-23, -0x800000.0p-37 }, 333 1.1 christos { 0xad0000.0p-23, 0xf80000.0p-35 }, 334 1.1 christos { 0xae0000.0p-23, 0xf80000.0p-34 }, 335 1.1 christos { 0xaf0000.0p-23, -0xac0000.0p-33 }, 336 1.1 christos { 0xb00000.0p-23, -0x800000.0p-33 }, 337 1.1 christos { 0xb10000.0p-23, -0xb80000.0p-34 }, 338 1.1 christos { 0xb20000.0p-23, -0x800000.0p-34 }, 339 1.1 christos { 0xb30000.0p-23, -0xb00000.0p-35 }, 340 1.1 christos { 0xb40000.0p-23, -0x800000.0p-35 }, 341 1.1 christos { 0xb50000.0p-23, -0xe00000.0p-36 }, 342 1.1 christos { 0xb60000.0p-23, -0x800000.0p-35 }, 343 1.1 christos { 0xb70000.0p-23, -0xb00000.0p-35 }, 344 1.1 christos { 0xb80000.0p-23, -0x800000.0p-34 }, 345 1.1 christos { 0xb90000.0p-23, -0xb80000.0p-34 }, 346 1.1 christos { 0xba0000.0p-23, -0x800000.0p-33 }, 347 1.1 christos { 0xbb0000.0p-23, -0xac0000.0p-33 }, 348 1.1 christos { 0xbc0000.0p-23, 0x980000.0p-33 }, 349 1.1 christos { 0xbd0000.0p-23, 0xbc0000.0p-34 }, 350 1.1 christos { 0xbe0000.0p-23, 0xe00000.0p-36 }, 351 1.1 christos { 0xbf0000.0p-23, -0xb80000.0p-35 }, 352 1.1 christos { 0xc00000.0p-23, -0x800000.0p-33 }, 353 1.1 christos { 0xc10000.0p-23, 0xa80000.0p-33 }, 354 1.1 christos { 0xc20000.0p-23, 0x900000.0p-34 }, 355 1.1 christos { 0xc30000.0p-23, -0x800000.0p-35 }, 356 1.1 christos { 0xc40000.0p-23, -0x900000.0p-33 }, 357 1.1 christos { 0xc50000.0p-23, 0x820000.0p-33 }, 358 1.1 christos { 0xc60000.0p-23, 0x800000.0p-38 }, 359 1.1 christos { 0xc70000.0p-23, -0x820000.0p-33 }, 360 1.1 christos { 0xc80000.0p-23, 0x800000.0p-33 }, 361 1.1 christos { 0xc90000.0p-23, -0xa00000.0p-36 }, 362 1.1 christos { 0xca0000.0p-23, -0xb00000.0p-33 }, 363 1.1 christos { 0xcb0000.0p-23, 0x840000.0p-34 }, 364 1.1 christos { 0xcc0000.0p-23, -0xd00000.0p-34 }, 365 1.1 christos { 0xcd0000.0p-23, 0x800000.0p-33 }, 366 1.1 christos { 0xce0000.0p-23, -0xe00000.0p-35 }, 367 1.1 christos { 0xcf0000.0p-23, 0xa60000.0p-33 }, 368 1.1 christos { 0xd00000.0p-23, -0x800000.0p-35 }, 369 1.1 christos { 0xd10000.0p-23, 0xb40000.0p-33 }, 370 1.1 christos { 0xd20000.0p-23, -0x800000.0p-35 }, 371 1.1 christos { 0xd30000.0p-23, 0xaa0000.0p-33 }, 372 1.1 christos { 0xd40000.0p-23, -0xe00000.0p-35 }, 373 1.1 christos { 0xd50000.0p-23, 0x880000.0p-33 }, 374 1.1 christos { 0xd60000.0p-23, -0xd00000.0p-34 }, 375 1.1 christos { 0xd70000.0p-23, 0x9c0000.0p-34 }, 376 1.1 christos { 0xd80000.0p-23, -0xb00000.0p-33 }, 377 1.1 christos { 0xd90000.0p-23, -0x800000.0p-38 }, 378 1.1 christos { 0xda0000.0p-23, 0xa40000.0p-33 }, 379 1.1 christos { 0xdb0000.0p-23, -0xdc0000.0p-34 }, 380 1.1 christos { 0xdc0000.0p-23, 0xc00000.0p-35 }, 381 1.1 christos { 0xdd0000.0p-23, 0xca0000.0p-33 }, 382 1.1 christos { 0xde0000.0p-23, -0xb80000.0p-34 }, 383 1.1 christos { 0xdf0000.0p-23, 0xd00000.0p-35 }, 384 1.1 christos { 0xe00000.0p-23, 0xc00000.0p-33 }, 385 1.1 christos { 0xe10000.0p-23, -0xf40000.0p-34 }, 386 1.1 christos { 0xe20000.0p-23, 0x800000.0p-37 }, 387 1.1 christos { 0xe30000.0p-23, 0x860000.0p-33 }, 388 1.1 christos { 0xe40000.0p-23, -0xc80000.0p-33 }, 389 1.1 christos { 0xe50000.0p-23, -0xa80000.0p-34 }, 390 1.1 christos { 0xe60000.0p-23, 0xe00000.0p-36 }, 391 1.1 christos { 0xe70000.0p-23, 0x880000.0p-33 }, 392 1.1 christos { 0xe80000.0p-23, -0xe00000.0p-33 }, 393 1.1 christos { 0xe90000.0p-23, -0xfc0000.0p-34 }, 394 1.1 christos { 0xea0000.0p-23, -0x800000.0p-35 }, 395 1.1 christos { 0xeb0000.0p-23, 0xe80000.0p-35 }, 396 1.1 christos { 0xec0000.0p-23, 0x900000.0p-33 }, 397 1.1 christos { 0xed0000.0p-23, 0xe20000.0p-33 }, 398 1.1 christos { 0xee0000.0p-23, -0xac0000.0p-33 }, 399 1.1 christos { 0xef0000.0p-23, -0xc80000.0p-34 }, 400 1.1 christos { 0xf00000.0p-23, -0x800000.0p-35 }, 401 1.1 christos { 0xf10000.0p-23, 0x800000.0p-35 }, 402 1.1 christos { 0xf20000.0p-23, 0xb80000.0p-34 }, 403 1.1 christos { 0xf30000.0p-23, 0x940000.0p-33 }, 404 1.1 christos { 0xf40000.0p-23, 0xc80000.0p-33 }, 405 1.1 christos { 0xf50000.0p-23, -0xf20000.0p-33 }, 406 1.1 christos { 0xf60000.0p-23, -0xc80000.0p-33 }, 407 1.1 christos { 0xf70000.0p-23, -0xa20000.0p-33 }, 408 1.1 christos { 0xf80000.0p-23, -0x800000.0p-33 }, 409 1.1 christos { 0xf90000.0p-23, -0xc40000.0p-34 }, 410 1.1 christos { 0xfa0000.0p-23, -0x900000.0p-34 }, 411 1.1 christos { 0xfb0000.0p-23, -0xc80000.0p-35 }, 412 1.1 christos { 0xfc0000.0p-23, -0x800000.0p-35 }, 413 1.1 christos { 0xfd0000.0p-23, -0x900000.0p-36 }, 414 1.1 christos { 0xfe0000.0p-23, -0x800000.0p-37 }, 415 1.1 christos { 0xff0000.0p-23, -0x800000.0p-39 }, 416 1.1 christos { 0x800000.0p-22, 0 }, 417 1.1 christos }; 418 1.1 christos #endif /* USE_UTAB */ 419 1.1 christos 420 1.1 christos #ifdef STRUCT_RETURN 421 1.1 christos #define RETURN1(rp, v) do { \ 422 1.1 christos (rp)->hi = (v); \ 423 1.1 christos (rp)->lo_set = 0; \ 424 1.1 christos return; \ 425 1.1 christos } while (0) 426 1.1 christos 427 1.1 christos #define RETURN2(rp, h, l) do { \ 428 1.1 christos (rp)->hi = (h); \ 429 1.1 christos (rp)->lo = (l); \ 430 1.1 christos (rp)->lo_set = 1; \ 431 1.1 christos return; \ 432 1.1 christos } while (0) 433 1.1 christos 434 1.1 christos struct ld { 435 1.1 christos long double hi; 436 1.1 christos long double lo; 437 1.1 christos int lo_set; 438 1.1 christos }; 439 1.1 christos #else 440 1.1 christos #define RETURN1(rp, v) RETURNF(v) 441 1.1 christos #define RETURN2(rp, h, l) RETURNI((h) + (l)) 442 1.1 christos #endif 443 1.1 christos 444 1.1 christos #ifdef STRUCT_RETURN 445 1.1 christos static inline __always_inline void 446 1.1 christos k_logl(long double x, struct ld *rp) 447 1.1 christos #else 448 1.1 christos long double 449 1.1 christos logl(long double x) 450 1.1 christos #endif 451 1.1 christos { 452 1.1 christos long double d, dk, val_hi, val_lo, z; 453 1.1 christos uint64_t ix, lx; 454 1.1 christos int i, k; 455 1.1 christos uint16_t hx; 456 1.1 christos 457 1.1 christos EXTRACT_LDBL80_WORDS(hx, lx, x); 458 1.1 christos k = -16383; 459 1.1 christos #if 0 /* Hard to do efficiently. Don't do it until we support all modes. */ 460 1.1 christos if (x == 1) 461 1.1 christos RETURN1(rp, 0); /* log(1) = +0 in all rounding modes */ 462 1.1 christos #endif 463 1.1 christos if (hx == 0 || hx >= 0x8000) { /* zero, negative or subnormal? */ 464 1.1 christos if (((hx & 0x7fff) | lx) == 0) 465 1.1 christos RETURN1(rp, -1 / zero); /* log(+-0) = -Inf */ 466 1.1 christos if (hx != 0) 467 1.1 christos /* log(neg or [pseudo-]NaN) = qNaN: */ 468 1.1 christos RETURN1(rp, (x - x) / zero); 469 1.1 christos x *= 0x1.0p65; /* subnormal; scale up x */ 470 1.1 christos /* including pseudo-subnormals */ 471 1.1 christos EXTRACT_LDBL80_WORDS(hx, lx, x); 472 1.1 christos k = -16383 - 65; 473 1.1 christos } else if (hx >= 0x7fff || (lx & 0x8000000000000000ULL) == 0) 474 1.1 christos RETURN1(rp, x + x); /* log(Inf or NaN) = Inf or qNaN */ 475 1.1 christos /* log(pseudo-Inf) = qNaN */ 476 1.1 christos /* log(pseudo-NaN) = qNaN */ 477 1.1 christos /* log(unnormal) = qNaN */ 478 1.1 christos #ifndef STRUCT_RETURN 479 1.1 christos ENTERI(); 480 1.1 christos #endif 481 1.1 christos k += hx; 482 1.1 christos ix = lx & 0x7fffffffffffffffULL; 483 1.1 christos dk = k; 484 1.1 christos 485 1.1 christos /* Scale x to be in [1, 2). */ 486 1.1 christos SET_LDBL_EXPSIGN(x, 0x3fff); 487 1.1 christos 488 1.1 christos /* 0 <= i <= INTERVALS: */ 489 1.1 christos #define L2I (64 - LOG2_INTERVALS) 490 1.1 christos i = (ix + (1LL << (L2I - 2))) >> (L2I - 1); 491 1.1 christos 492 1.1 christos /* 493 1.1 christos * -0.005280 < d < 0.004838. In particular, the infinite- 494 1.1 christos * precision |d| is <= 2**-7. Rounding of G(i) to 8 bits 495 1.1 christos * ensures that d is representable without extra precision for 496 1.1 christos * this bound on |d| (since when this calculation is expressed 497 1.1 christos * as x*G(i)-1, the multiplication needs as many extra bits as 498 1.1 christos * G(i) has and the subtraction cancels 8 bits). But for 499 1.1 christos * most i (107 cases out of 129), the infinite-precision |d| 500 1.1 christos * is <= 2**-8. G(i) is rounded to 9 bits for such i to give 501 1.1 christos * better accuracy (this works by improving the bound on |d|, 502 1.1 christos * which in turn allows rounding to 9 bits in more cases). 503 1.1 christos * This is only important when the original x is near 1 -- it 504 1.1 christos * lets us avoid using a special method to give the desired 505 1.1 christos * accuracy for such x. 506 1.1 christos */ 507 1.1 christos if (/*CONSTCOND*/0) 508 1.1 christos /*NOTREACHED*/ 509 1.1 christos d = x * G(i) - 1; 510 1.1 christos else { 511 1.1 christos #ifdef USE_UTAB 512 1.1 christos d = (x - H(i)) * G(i) + E(i); 513 1.1 christos #else 514 1.1 christos long double x_hi, x_lo; 515 1.1 christos float fx_hi; 516 1.1 christos 517 1.1 christos /* 518 1.1 christos * Split x into x_hi + x_lo to calculate x*G(i)-1 exactly. 519 1.1 christos * G(i) has at most 9 bits, so the splitting point is not 520 1.1 christos * critical. 521 1.1 christos */ 522 1.1 christos SET_FLOAT_WORD(fx_hi, (lx >> 40) | 0x3f800000); 523 1.1 christos x_hi = fx_hi; 524 1.1 christos x_lo = x - x_hi; 525 1.1 christos d = x_hi * G(i) - 1 + x_lo * G(i); 526 1.1 christos #endif 527 1.1 christos } 528 1.1 christos 529 1.1 christos /* 530 1.1 christos * Our algorithm depends on exact cancellation of F_lo(i) and 531 1.1 christos * F_hi(i) with dk*ln_2_lo and dk*ln2_hi when k is -1 and i is 532 1.1 christos * at the end of the table. This and other technical complications 533 1.1 christos * make it difficult to avoid the double scaling in (dk*ln2) * 534 1.1 christos * log(base) for base != e without losing more accuracy and/or 535 1.1 christos * efficiency than is gained. 536 1.1 christos */ 537 1.1 christos z = d * d; 538 1.1 christos val_lo = z * d * z * (z * (d * P8 + P7) + (d * P6 + P5)) + 539 1.1 christos (F_lo(i) + dk * ln2_lo + z * d * (d * P4 + P3)) + z * P2; 540 1.1 christos val_hi = d; 541 1.1 christos #ifdef DEBUG 542 1.1 christos if (fetestexcept(FE_UNDERFLOW)) 543 1.1 christos breakpoint(); 544 1.1 christos #endif 545 1.1 christos 546 1.1 christos _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi); 547 1.1 christos RETURN2(rp, val_hi, val_lo); 548 1.1 christos } 549 1.1 christos 550 1.1 christos long double 551 1.1 christos log1pl(long double x) 552 1.1 christos { 553 1.1 christos long double d, d_hi, d_lo, dk, f_lo, val_hi, val_lo, z; 554 1.1 christos long double f_hi, twopminusk; 555 1.1 christos uint64_t ix, lx; 556 1.1 christos int i, k; 557 1.1 christos int16_t ax, hx; 558 1.1 christos 559 1.1 christos EXTRACT_LDBL80_WORDS(hx, lx, x); 560 1.1 christos if (hx < 0x3fff) { /* x < 1, or x neg NaN */ 561 1.1 christos ax = hx & 0x7fff; 562 1.1 christos if (ax >= 0x3fff) { /* x <= -1, or x neg NaN */ 563 1.1 christos if (ax == 0x3fff && lx == 0x8000000000000000ULL) 564 1.1 christos RETURNF(-1 / zero); /* log1p(-1) = -Inf */ 565 1.1 christos /* log1p(x < 1, or x [pseudo-]NaN) = qNaN: */ 566 1.1 christos RETURNF((x - x) / (x - x)); 567 1.1 christos } 568 1.1 christos if (ax <= 0x3fbe) { /* |x| < 2**-64 */ 569 1.1 christos if ((int)x == 0) 570 1.1 christos RETURNF(x); /* x with inexact if x != 0 */ 571 1.1 christos } 572 1.1 christos f_hi = 1; 573 1.1 christos f_lo = x; 574 1.1 christos } else if (hx >= 0x7fff) { /* x +Inf or non-neg NaN */ 575 1.1 christos RETURNF(x + x); /* log1p(Inf or NaN) = Inf or qNaN */ 576 1.1 christos /* log1p(pseudo-Inf) = qNaN */ 577 1.1 christos /* log1p(pseudo-NaN) = qNaN */ 578 1.1 christos /* log1p(unnormal) = qNaN */ 579 1.1 christos } else if (hx < 0x407f) { /* 1 <= x < 2**128 */ 580 1.1 christos f_hi = x; 581 1.1 christos f_lo = 1; 582 1.1 christos } else { /* 2**128 <= x < +Inf */ 583 1.1 christos f_hi = x; 584 1.1 christos f_lo = 0; /* avoid underflow of the P5 term */ 585 1.1 christos } 586 1.1 christos ENTERI(); 587 1.1 christos x = f_hi + f_lo; 588 1.1 christos f_lo = (f_hi - x) + f_lo; 589 1.1 christos 590 1.1 christos EXTRACT_LDBL80_WORDS(hx, lx, x); 591 1.1 christos k = -16383; 592 1.1 christos 593 1.1 christos k += hx; 594 1.1 christos ix = lx & 0x7fffffffffffffffULL; 595 1.1 christos dk = k; 596 1.1 christos 597 1.1 christos SET_LDBL_EXPSIGN(x, 0x3fff); 598 1.1 christos twopminusk = 1; 599 1.1 christos SET_LDBL_EXPSIGN(twopminusk, 0x7ffe - (hx & 0x7fff)); 600 1.1 christos f_lo *= twopminusk; 601 1.1 christos 602 1.1 christos i = (ix + (1LL << (L2I - 2))) >> (L2I - 1); 603 1.1 christos 604 1.1 christos /* 605 1.1 christos * x*G(i)-1 (with a reduced x) can be represented exactly, as 606 1.1 christos * above, but now we need to evaluate the polynomial on d = 607 1.1 christos * (x+f_lo)*G(i)-1 and extra precision is needed for that. 608 1.1 christos * Since x+x_lo is a hi+lo decomposition and subtracting 1 609 1.1 christos * doesn't lose too many bits, an inexact calculation for 610 1.1 christos * f_lo*G(i) is good enough. 611 1.1 christos */ 612 1.1 christos if (/*CONSTCOND*/0) 613 1.1 christos /*NOTREACHED*/ 614 1.1 christos d_hi = x * G(i) - 1; 615 1.1 christos else { 616 1.1 christos #ifdef USE_UTAB 617 1.1 christos d_hi = (x - H(i)) * G(i) + E(i); 618 1.1 christos #else 619 1.1 christos long double x_hi, x_lo; 620 1.1 christos float fx_hi; 621 1.1 christos 622 1.1 christos SET_FLOAT_WORD(fx_hi, (lx >> 40) | 0x3f800000); 623 1.1 christos x_hi = fx_hi; 624 1.1 christos x_lo = x - x_hi; 625 1.1 christos d_hi = x_hi * G(i) - 1 + x_lo * G(i); 626 1.1 christos #endif 627 1.1 christos } 628 1.1 christos d_lo = f_lo * G(i); 629 1.1 christos 630 1.1 christos /* 631 1.1 christos * This is _2sumF(d_hi, d_lo) inlined. The condition 632 1.1 christos * (d_hi == 0 || |d_hi| >= |d_lo|) for using _2sumF() is not 633 1.1 christos * always satisifed, so it is not clear that this works, but 634 1.1 christos * it works in practice. It works even if it gives a wrong 635 1.1 christos * normalized d_lo, since |d_lo| > |d_hi| implies that i is 636 1.1 christos * nonzero and d is tiny, so the F(i) term dominates d_lo. 637 1.1 christos * In float precision: 638 1.1 christos * (By exhaustive testing, the worst case is d_hi = 0x1.bp-25. 639 1.1 christos * And if d is only a little tinier than that, we would have 640 1.1 christos * another underflow problem for the P3 term; this is also ruled 641 1.1 christos * out by exhaustive testing.) 642 1.1 christos */ 643 1.1 christos d = d_hi + d_lo; 644 1.1 christos d_lo = d_hi - d + d_lo; 645 1.1 christos d_hi = d; 646 1.1 christos 647 1.1 christos z = d * d; 648 1.1 christos val_lo = z * d * z * (z * (d * P8 + P7) + (d * P6 + P5)) + 649 1.1 christos (F_lo(i) + dk * ln2_lo + d_lo + z * d * (d * P4 + P3)) + z * P2; 650 1.1 christos val_hi = d_hi; 651 1.1 christos #ifdef DEBUG 652 1.1 christos if (fetestexcept(FE_UNDERFLOW)) 653 1.1 christos breakpoint(); 654 1.1 christos #endif 655 1.1 christos 656 1.1 christos _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi); 657 1.1 christos RETURNI(val_hi + val_lo); 658 1.1 christos } 659 1.1 christos 660 1.1 christos #ifdef STRUCT_RETURN 661 1.1 christos 662 1.1 christos long double 663 1.1 christos logl(long double x) 664 1.1 christos { 665 1.1 christos struct ld r; 666 1.1 christos 667 1.1 christos ENTERI(); 668 1.1 christos k_logl(x, &r); 669 1.1 christos RETURNSPI(&r); 670 1.1 christos } 671 1.1 christos 672 1.1 christos /* Use macros since GCC < 8 rejects static const expressions in initializers. */ 673 1.1 christos #define invln10_hi 4.3429448190317999e-1 /* 0x1bcb7b1526e000.0p-54 */ 674 1.1 christos #define invln10_lo 7.1842412889749798e-14 /* 0x1438ca9aadd558.0p-96 */ 675 1.1 christos #define invln2_hi 1.4426950408887933e0 /* 0x171547652b8000.0p-52 */ 676 1.1 christos #define invln2_lo 1.7010652264631490e-13 /* 0x17f0bbbe87fed0.0p-95 */ 677 1.1 christos /* Let the compiler pre-calculate this sum to avoid FE_INEXACT at run time. */ 678 1.1 christos static const double invln10_lo_plus_hi = invln10_lo + invln10_hi; 679 1.1 christos static const double invln2_lo_plus_hi = invln2_lo + invln2_hi; 680 1.1 christos 681 1.1 christos long double 682 1.1 christos log10l(long double x) 683 1.1 christos { 684 1.1 christos struct ld r; 685 1.1 christos long double hi, lo; 686 1.1 christos 687 1.1 christos ENTERI(); 688 1.1 christos k_logl(x, &r); 689 1.1 christos if (!r.lo_set) 690 1.1 christos RETURNI(r.hi); 691 1.1 christos _2sumF(r.hi, r.lo); 692 1.1 christos hi = (float)r.hi; 693 1.1 christos lo = r.lo + (r.hi - hi); 694 1.1 christos RETURNI(invln10_hi * hi + 695 1.1 christos (invln10_lo_plus_hi * lo + invln10_lo * hi)); 696 1.1 christos } 697 1.1 christos 698 1.1 christos long double 699 1.1 christos log2l(long double x) 700 1.1 christos { 701 1.1 christos struct ld r; 702 1.1 christos long double hi, lo; 703 1.1 christos 704 1.1 christos ENTERI(); 705 1.1 christos k_logl(x, &r); 706 1.1 christos if (!r.lo_set) 707 1.1 christos RETURNI(r.hi); 708 1.1 christos _2sumF(r.hi, r.lo); 709 1.1 christos hi = (float)r.hi; 710 1.1 christos lo = r.lo + (r.hi - hi); 711 1.1 christos RETURNI(invln2_hi * hi + 712 1.1 christos (invln2_lo_plus_hi * lo + invln2_lo * hi)); 713 1.1 christos } 714 1.1 christos 715 1.1 christos #endif /* STRUCT_RETURN */ 716