[352] | 1 | /* origin: FreeBSD /usr/src/lib/msun/ld80/s_exp2l.c and /usr/src/lib/msun/ld128/s_exp2l.c */
|
---|
| 2 | /*-
|
---|
| 3 | * Copyright (c) 2005-2008 David Schultz <das@FreeBSD.ORG>
|
---|
| 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
|
---|
| 8 | * are met:
|
---|
| 9 | * 1. Redistributions of source code must retain the above copyright
|
---|
| 10 | * notice, this list of conditions and the following disclaimer.
|
---|
| 11 | * 2. 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 | *
|
---|
| 15 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
|
---|
| 16 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
---|
| 17 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
---|
| 18 | * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
|
---|
| 19 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
|
---|
| 20 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
|
---|
| 21 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
|
---|
| 22 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
|
---|
| 23 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
|
---|
| 24 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
|
---|
| 25 | * SUCH DAMAGE.
|
---|
| 26 | */
|
---|
| 27 |
|
---|
| 28 | #include "libm.h"
|
---|
| 29 |
|
---|
| 30 | #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
|
---|
| 31 | long double exp2l(long double x)
|
---|
| 32 | {
|
---|
| 33 | return exp2(x);
|
---|
| 34 | }
|
---|
| 35 | #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
|
---|
| 36 | #define TBLBITS 7
|
---|
| 37 | #define TBLSIZE (1 << TBLBITS)
|
---|
| 38 |
|
---|
| 39 | static const double
|
---|
| 40 | redux = 0x1.8p63 / TBLSIZE,
|
---|
| 41 | P1 = 0x1.62e42fefa39efp-1,
|
---|
| 42 | P2 = 0x1.ebfbdff82c58fp-3,
|
---|
| 43 | P3 = 0x1.c6b08d7049fap-5,
|
---|
| 44 | P4 = 0x1.3b2ab6fba4da5p-7,
|
---|
| 45 | P5 = 0x1.5d8804780a736p-10,
|
---|
| 46 | P6 = 0x1.430918835e33dp-13;
|
---|
| 47 |
|
---|
| 48 | static const double tbl[TBLSIZE * 2] = {
|
---|
| 49 | 0x1.6a09e667f3bcdp-1, -0x1.bdd3413b2648p-55,
|
---|
| 50 | 0x1.6c012750bdabfp-1, -0x1.2895667ff0cp-57,
|
---|
| 51 | 0x1.6dfb23c651a2fp-1, -0x1.bbe3a683c88p-58,
|
---|
| 52 | 0x1.6ff7df9519484p-1, -0x1.83c0f25860fp-56,
|
---|
| 53 | 0x1.71f75e8ec5f74p-1, -0x1.16e4786887bp-56,
|
---|
| 54 | 0x1.73f9a48a58174p-1, -0x1.0a8d96c65d5p-55,
|
---|
| 55 | 0x1.75feb564267c9p-1, -0x1.0245957316ep-55,
|
---|
| 56 | 0x1.780694fde5d3fp-1, 0x1.866b80a0216p-55,
|
---|
| 57 | 0x1.7a11473eb0187p-1, -0x1.41577ee0499p-56,
|
---|
| 58 | 0x1.7c1ed0130c132p-1, 0x1.f124cd1164ep-55,
|
---|
| 59 | 0x1.7e2f336cf4e62p-1, 0x1.05d02ba157ap-57,
|
---|
| 60 | 0x1.80427543e1a12p-1, -0x1.27c86626d97p-55,
|
---|
| 61 | 0x1.82589994cce13p-1, -0x1.d4c1dd41533p-55,
|
---|
| 62 | 0x1.8471a4623c7adp-1, -0x1.8d684a341cep-56,
|
---|
| 63 | 0x1.868d99b4492edp-1, -0x1.fc6f89bd4f68p-55,
|
---|
| 64 | 0x1.88ac7d98a6699p-1, 0x1.994c2f37cb5p-55,
|
---|
| 65 | 0x1.8ace5422aa0dbp-1, 0x1.6e9f156864bp-55,
|
---|
| 66 | 0x1.8cf3216b5448cp-1, -0x1.0d55e32e9e4p-57,
|
---|
| 67 | 0x1.8f1ae99157736p-1, 0x1.5cc13a2e397p-56,
|
---|
| 68 | 0x1.9145b0b91ffc6p-1, -0x1.dd6792e5825p-55,
|
---|
| 69 | 0x1.93737b0cdc5e5p-1, -0x1.75fc781b58p-58,
|
---|
| 70 | 0x1.95a44cbc8520fp-1, -0x1.64b7c96a5fp-57,
|
---|
| 71 | 0x1.97d829fde4e5p-1, -0x1.d185b7c1b86p-55,
|
---|
| 72 | 0x1.9a0f170ca07bap-1, -0x1.173bd91cee6p-55,
|
---|
| 73 | 0x1.9c49182a3f09p-1, 0x1.c7c46b071f2p-57,
|
---|
| 74 | 0x1.9e86319e32323p-1, 0x1.824ca78e64cp-57,
|
---|
| 75 | 0x1.a0c667b5de565p-1, -0x1.359495d1cd5p-55,
|
---|
| 76 | 0x1.a309bec4a2d33p-1, 0x1.6305c7ddc368p-55,
|
---|
| 77 | 0x1.a5503b23e255dp-1, -0x1.d2f6edb8d42p-55,
|
---|
| 78 | 0x1.a799e1330b358p-1, 0x1.bcb7ecac564p-55,
|
---|
| 79 | 0x1.a9e6b5579fdbfp-1, 0x1.0fac90ef7fdp-55,
|
---|
| 80 | 0x1.ac36bbfd3f37ap-1, -0x1.f9234cae76dp-56,
|
---|
| 81 | 0x1.ae89f995ad3adp-1, 0x1.7a1cd345dcc8p-55,
|
---|
| 82 | 0x1.b0e07298db666p-1, -0x1.bdef54c80e4p-55,
|
---|
| 83 | 0x1.b33a2b84f15fbp-1, -0x1.2805e3084d8p-58,
|
---|
| 84 | 0x1.b59728de5593ap-1, -0x1.c71dfbbba6ep-55,
|
---|
| 85 | 0x1.b7f76f2fb5e47p-1, -0x1.5584f7e54acp-57,
|
---|
| 86 | 0x1.ba5b030a1064ap-1, -0x1.efcd30e5429p-55,
|
---|
| 87 | 0x1.bcc1e904bc1d2p-1, 0x1.23dd07a2d9fp-56,
|
---|
| 88 | 0x1.bf2c25bd71e09p-1, -0x1.efdca3f6b9c8p-55,
|
---|
| 89 | 0x1.c199bdd85529cp-1, 0x1.11065895049p-56,
|
---|
| 90 | 0x1.c40ab5fffd07ap-1, 0x1.b4537e083c6p-55,
|
---|
| 91 | 0x1.c67f12e57d14bp-1, 0x1.2884dff483c8p-55,
|
---|
| 92 | 0x1.c8f6d9406e7b5p-1, 0x1.1acbc48805cp-57,
|
---|
| 93 | 0x1.cb720dcef9069p-1, 0x1.503cbd1e94ap-57,
|
---|
| 94 | 0x1.cdf0b555dc3fap-1, -0x1.dd83b53829dp-56,
|
---|
| 95 | 0x1.d072d4a07897cp-1, -0x1.cbc3743797a8p-55,
|
---|
| 96 | 0x1.d2f87080d89f2p-1, -0x1.d487b719d858p-55,
|
---|
| 97 | 0x1.d5818dcfba487p-1, 0x1.2ed02d75b37p-56,
|
---|
| 98 | 0x1.d80e316c98398p-1, -0x1.11ec18bedep-55,
|
---|
| 99 | 0x1.da9e603db3285p-1, 0x1.c2300696db5p-55,
|
---|
| 100 | 0x1.dd321f301b46p-1, 0x1.2da5778f019p-55,
|
---|
| 101 | 0x1.dfc97337b9b5fp-1, -0x1.1a5cd4f184b8p-55,
|
---|
| 102 | 0x1.e264614f5a129p-1, -0x1.7b627817a148p-55,
|
---|
| 103 | 0x1.e502ee78b3ff6p-1, 0x1.39e8980a9cdp-56,
|
---|
| 104 | 0x1.e7a51fbc74c83p-1, 0x1.2d522ca0c8ep-55,
|
---|
| 105 | 0x1.ea4afa2a490dap-1, -0x1.e9c23179c288p-55,
|
---|
| 106 | 0x1.ecf482d8e67f1p-1, -0x1.c93f3b411ad8p-55,
|
---|
| 107 | 0x1.efa1bee615a27p-1, 0x1.dc7f486a4b68p-55,
|
---|
| 108 | 0x1.f252b376bba97p-1, 0x1.3a1a5bf0d8e8p-55,
|
---|
| 109 | 0x1.f50765b6e454p-1, 0x1.9d3e12dd8a18p-55,
|
---|
| 110 | 0x1.f7bfdad9cbe14p-1, -0x1.dbb12d00635p-55,
|
---|
| 111 | 0x1.fa7c1819e90d8p-1, 0x1.74853f3a593p-56,
|
---|
| 112 | 0x1.fd3c22b8f71f1p-1, 0x1.2eb74966578p-58,
|
---|
| 113 | 0x1p+0, 0x0p+0,
|
---|
| 114 | 0x1.0163da9fb3335p+0, 0x1.b61299ab8cd8p-54,
|
---|
| 115 | 0x1.02c9a3e778061p+0, -0x1.19083535b08p-56,
|
---|
| 116 | 0x1.04315e86e7f85p+0, -0x1.0a31c1977c98p-54,
|
---|
| 117 | 0x1.059b0d3158574p+0, 0x1.d73e2a475b4p-55,
|
---|
| 118 | 0x1.0706b29ddf6dep+0, -0x1.c91dfe2b13cp-55,
|
---|
| 119 | 0x1.0874518759bc8p+0, 0x1.186be4bb284p-57,
|
---|
| 120 | 0x1.09e3ecac6f383p+0, 0x1.14878183161p-54,
|
---|
| 121 | 0x1.0b5586cf9890fp+0, 0x1.8a62e4adc61p-54,
|
---|
| 122 | 0x1.0cc922b7247f7p+0, 0x1.01edc16e24f8p-54,
|
---|
| 123 | 0x1.0e3ec32d3d1a2p+0, 0x1.03a1727c58p-59,
|
---|
| 124 | 0x1.0fb66affed31bp+0, -0x1.b9bedc44ebcp-57,
|
---|
| 125 | 0x1.11301d0125b51p+0, -0x1.6c51039449bp-54,
|
---|
| 126 | 0x1.12abdc06c31ccp+0, -0x1.1b514b36ca8p-58,
|
---|
| 127 | 0x1.1429aaea92dep+0, -0x1.32fbf9af1368p-54,
|
---|
| 128 | 0x1.15a98c8a58e51p+0, 0x1.2406ab9eeabp-55,
|
---|
| 129 | 0x1.172b83c7d517bp+0, -0x1.19041b9d78ap-55,
|
---|
| 130 | 0x1.18af9388c8deap+0, -0x1.11023d1970f8p-54,
|
---|
| 131 | 0x1.1a35beb6fcb75p+0, 0x1.e5b4c7b4969p-55,
|
---|
| 132 | 0x1.1bbe084045cd4p+0, -0x1.95386352ef6p-54,
|
---|
| 133 | 0x1.1d4873168b9aap+0, 0x1.e016e00a264p-54,
|
---|
| 134 | 0x1.1ed5022fcd91dp+0, -0x1.1df98027bb78p-54,
|
---|
| 135 | 0x1.2063b88628cd6p+0, 0x1.dc775814a85p-55,
|
---|
| 136 | 0x1.21f49917ddc96p+0, 0x1.2a97e9494a6p-55,
|
---|
| 137 | 0x1.2387a6e756238p+0, 0x1.9b07eb6c7058p-54,
|
---|
| 138 | 0x1.251ce4fb2a63fp+0, 0x1.ac155bef4f5p-55,
|
---|
| 139 | 0x1.26b4565e27cddp+0, 0x1.2bd339940eap-55,
|
---|
| 140 | 0x1.284dfe1f56381p+0, -0x1.a4c3a8c3f0d8p-54,
|
---|
| 141 | 0x1.29e9df51fdee1p+0, 0x1.612e8afad12p-55,
|
---|
| 142 | 0x1.2b87fd0dad99p+0, -0x1.10adcd6382p-59,
|
---|
| 143 | 0x1.2d285a6e4030bp+0, 0x1.0024754db42p-54,
|
---|
| 144 | 0x1.2ecafa93e2f56p+0, 0x1.1ca0f45d524p-56,
|
---|
| 145 | 0x1.306fe0a31b715p+0, 0x1.6f46ad23183p-55,
|
---|
| 146 | 0x1.32170fc4cd831p+0, 0x1.a9ce78e1804p-55,
|
---|
| 147 | 0x1.33c08b26416ffp+0, 0x1.327218436598p-54,
|
---|
| 148 | 0x1.356c55f929ff1p+0, -0x1.b5cee5c4e46p-55,
|
---|
| 149 | 0x1.371a7373aa9cbp+0, -0x1.63aeabf42ebp-54,
|
---|
| 150 | 0x1.38cae6d05d866p+0, -0x1.e958d3c99048p-54,
|
---|
| 151 | 0x1.3a7db34e59ff7p+0, -0x1.5e436d661f6p-56,
|
---|
| 152 | 0x1.3c32dc313a8e5p+0, -0x1.efff8375d2ap-54,
|
---|
| 153 | 0x1.3dea64c123422p+0, 0x1.ada0911f09fp-55,
|
---|
| 154 | 0x1.3fa4504ac801cp+0, -0x1.7d023f956fap-54,
|
---|
| 155 | 0x1.4160a21f72e2ap+0, -0x1.ef3691c309p-58,
|
---|
| 156 | 0x1.431f5d950a897p+0, -0x1.1c7dde35f7ap-55,
|
---|
| 157 | 0x1.44e086061892dp+0, 0x1.89b7a04ef8p-59,
|
---|
| 158 | 0x1.46a41ed1d0057p+0, 0x1.c944bd1648a8p-54,
|
---|
| 159 | 0x1.486a2b5c13cdp+0, 0x1.3c1a3b69062p-56,
|
---|
| 160 | 0x1.4a32af0d7d3dep+0, 0x1.9cb62f3d1be8p-54,
|
---|
| 161 | 0x1.4bfdad5362a27p+0, 0x1.d4397afec42p-56,
|
---|
| 162 | 0x1.4dcb299fddd0dp+0, 0x1.8ecdbbc6a78p-54,
|
---|
| 163 | 0x1.4f9b2769d2ca7p+0, -0x1.4b309d25958p-54,
|
---|
| 164 | 0x1.516daa2cf6642p+0, -0x1.f768569bd94p-55,
|
---|
| 165 | 0x1.5342b569d4f82p+0, -0x1.07abe1db13dp-55,
|
---|
| 166 | 0x1.551a4ca5d920fp+0, -0x1.d689cefede6p-55,
|
---|
| 167 | 0x1.56f4736b527dap+0, 0x1.9bb2c011d938p-54,
|
---|
| 168 | 0x1.58d12d497c7fdp+0, 0x1.295e15b9a1ep-55,
|
---|
| 169 | 0x1.5ab07dd485429p+0, 0x1.6324c0546478p-54,
|
---|
| 170 | 0x1.5c9268a5946b7p+0, 0x1.c4b1b81698p-60,
|
---|
| 171 | 0x1.5e76f15ad2148p+0, 0x1.ba6f93080e68p-54,
|
---|
| 172 | 0x1.605e1b976dc09p+0, -0x1.3e2429b56de8p-54,
|
---|
| 173 | 0x1.6247eb03a5585p+0, -0x1.383c17e40b48p-54,
|
---|
| 174 | 0x1.6434634ccc32p+0, -0x1.c483c759d89p-55,
|
---|
| 175 | 0x1.6623882552225p+0, -0x1.bb60987591cp-54,
|
---|
| 176 | 0x1.68155d44ca973p+0, 0x1.038ae44f74p-57,
|
---|
| 177 | };
|
---|
| 178 |
|
---|
| 179 | /*
|
---|
| 180 | * exp2l(x): compute the base 2 exponential of x
|
---|
| 181 | *
|
---|
| 182 | * Accuracy: Peak error < 0.511 ulp.
|
---|
| 183 | *
|
---|
| 184 | * Method: (equally-spaced tables)
|
---|
| 185 | *
|
---|
| 186 | * Reduce x:
|
---|
| 187 | * x = 2**k + y, for integer k and |y| <= 1/2.
|
---|
| 188 | * Thus we have exp2l(x) = 2**k * exp2(y).
|
---|
| 189 | *
|
---|
| 190 | * Reduce y:
|
---|
| 191 | * y = i/TBLSIZE + z for integer i near y * TBLSIZE.
|
---|
| 192 | * Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z),
|
---|
| 193 | * with |z| <= 2**-(TBLBITS+1).
|
---|
| 194 | *
|
---|
| 195 | * We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a
|
---|
| 196 | * degree-6 minimax polynomial with maximum error under 2**-69.
|
---|
| 197 | * The table entries each have 104 bits of accuracy, encoded as
|
---|
| 198 | * a pair of double precision values.
|
---|
| 199 | */
|
---|
| 200 | long double exp2l(long double x)
|
---|
| 201 | {
|
---|
| 202 | union ldshape u = {x};
|
---|
| 203 | int e = u.i.se & 0x7fff;
|
---|
| 204 | long double r, z;
|
---|
| 205 | uint32_t i0;
|
---|
| 206 | union {uint32_t u; int32_t i;} k;
|
---|
| 207 |
|
---|
| 208 | /* Filter out exceptional cases. */
|
---|
| 209 | if (e >= 0x3fff + 13) { /* |x| >= 8192 or x is NaN */
|
---|
| 210 | if (u.i.se >= 0x3fff + 14 && u.i.se >> 15 == 0)
|
---|
| 211 | /* overflow */
|
---|
| 212 | return x * 0x1p16383L;
|
---|
| 213 | if (e == 0x7fff) /* -inf or -nan */
|
---|
| 214 | return -1/x;
|
---|
| 215 | if (x < -16382) {
|
---|
| 216 | if (x <= -16446 || x - 0x1p63 + 0x1p63 != x)
|
---|
| 217 | /* underflow */
|
---|
| 218 | FORCE_EVAL((float)(-0x1p-149/x));
|
---|
| 219 | if (x <= -16446)
|
---|
| 220 | return 0;
|
---|
| 221 | }
|
---|
| 222 | } else if (e < 0x3fff - 64) {
|
---|
| 223 | return 1 + x;
|
---|
| 224 | }
|
---|
| 225 |
|
---|
| 226 | /*
|
---|
| 227 | * Reduce x, computing z, i0, and k. The low bits of x + redux
|
---|
| 228 | * contain the 16-bit integer part of the exponent (k) followed by
|
---|
| 229 | * TBLBITS fractional bits (i0). We use bit tricks to extract these
|
---|
| 230 | * as integers, then set z to the remainder.
|
---|
| 231 | *
|
---|
| 232 | * Example: Suppose x is 0xabc.123456p0 and TBLBITS is 8.
|
---|
| 233 | * Then the low-order word of x + redux is 0x000abc12,
|
---|
| 234 | * We split this into k = 0xabc and i0 = 0x12 (adjusted to
|
---|
| 235 | * index into the table), then we compute z = 0x0.003456p0.
|
---|
| 236 | */
|
---|
| 237 | u.f = x + redux;
|
---|
| 238 | i0 = u.i.m + TBLSIZE / 2;
|
---|
| 239 | k.u = i0 / TBLSIZE * TBLSIZE;
|
---|
| 240 | k.i /= TBLSIZE;
|
---|
| 241 | i0 %= TBLSIZE;
|
---|
| 242 | u.f -= redux;
|
---|
| 243 | z = x - u.f;
|
---|
| 244 |
|
---|
| 245 | /* Compute r = exp2l(y) = exp2lt[i0] * p(z). */
|
---|
| 246 | long double t_hi = tbl[2*i0];
|
---|
| 247 | long double t_lo = tbl[2*i0 + 1];
|
---|
| 248 | /* XXX This gives > 1 ulp errors outside of FE_TONEAREST mode */
|
---|
| 249 | r = t_lo + (t_hi + t_lo) * z * (P1 + z * (P2 + z * (P3 + z * (P4
|
---|
| 250 | + z * (P5 + z * P6))))) + t_hi;
|
---|
| 251 |
|
---|
| 252 | return scalbnl(r, k.i);
|
---|
| 253 | }
|
---|
| 254 | #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
|
---|
| 255 | #define TBLBITS 7
|
---|
| 256 | #define TBLSIZE (1 << TBLBITS)
|
---|
| 257 |
|
---|
| 258 | static const long double
|
---|
| 259 | P1 = 0x1.62e42fefa39ef35793c7673007e6p-1L,
|
---|
| 260 | P2 = 0x1.ebfbdff82c58ea86f16b06ec9736p-3L,
|
---|
| 261 | P3 = 0x1.c6b08d704a0bf8b33a762bad3459p-5L,
|
---|
| 262 | P4 = 0x1.3b2ab6fba4e7729ccbbe0b4f3fc2p-7L,
|
---|
| 263 | P5 = 0x1.5d87fe78a67311071dee13fd11d9p-10L,
|
---|
| 264 | P6 = 0x1.430912f86c7876f4b663b23c5fe5p-13L;
|
---|
| 265 |
|
---|
| 266 | static const double
|
---|
| 267 | P7 = 0x1.ffcbfc588b041p-17,
|
---|
| 268 | P8 = 0x1.62c0223a5c7c7p-20,
|
---|
| 269 | P9 = 0x1.b52541ff59713p-24,
|
---|
| 270 | P10 = 0x1.e4cf56a391e22p-28,
|
---|
| 271 | redux = 0x1.8p112 / TBLSIZE;
|
---|
| 272 |
|
---|
| 273 | static const long double tbl[TBLSIZE] = {
|
---|
| 274 | 0x1.6a09e667f3bcc908b2fb1366dfeap-1L,
|
---|
| 275 | 0x1.6c012750bdabeed76a99800f4edep-1L,
|
---|
| 276 | 0x1.6dfb23c651a2ef220e2cbe1bc0d4p-1L,
|
---|
| 277 | 0x1.6ff7df9519483cf87e1b4f3e1e98p-1L,
|
---|
| 278 | 0x1.71f75e8ec5f73dd2370f2ef0b148p-1L,
|
---|
| 279 | 0x1.73f9a48a58173bd5c9a4e68ab074p-1L,
|
---|
| 280 | 0x1.75feb564267c8bf6e9aa33a489a8p-1L,
|
---|
| 281 | 0x1.780694fde5d3f619ae02808592a4p-1L,
|
---|
| 282 | 0x1.7a11473eb0186d7d51023f6ccb1ap-1L,
|
---|
| 283 | 0x1.7c1ed0130c1327c49334459378dep-1L,
|
---|
| 284 | 0x1.7e2f336cf4e62105d02ba1579756p-1L,
|
---|
| 285 | 0x1.80427543e1a11b60de67649a3842p-1L,
|
---|
| 286 | 0x1.82589994cce128acf88afab34928p-1L,
|
---|
| 287 | 0x1.8471a4623c7acce52f6b97c6444cp-1L,
|
---|
| 288 | 0x1.868d99b4492ec80e41d90ac2556ap-1L,
|
---|
| 289 | 0x1.88ac7d98a669966530bcdf2d4cc0p-1L,
|
---|
| 290 | 0x1.8ace5422aa0db5ba7c55a192c648p-1L,
|
---|
| 291 | 0x1.8cf3216b5448bef2aa1cd161c57ap-1L,
|
---|
| 292 | 0x1.8f1ae991577362b982745c72eddap-1L,
|
---|
| 293 | 0x1.9145b0b91ffc588a61b469f6b6a0p-1L,
|
---|
| 294 | 0x1.93737b0cdc5e4f4501c3f2540ae8p-1L,
|
---|
| 295 | 0x1.95a44cbc8520ee9b483695a0e7fep-1L,
|
---|
| 296 | 0x1.97d829fde4e4f8b9e920f91e8eb6p-1L,
|
---|
| 297 | 0x1.9a0f170ca07b9ba3109b8c467844p-1L,
|
---|
| 298 | 0x1.9c49182a3f0901c7c46b071f28dep-1L,
|
---|
| 299 | 0x1.9e86319e323231824ca78e64c462p-1L,
|
---|
| 300 | 0x1.a0c667b5de564b29ada8b8cabbacp-1L,
|
---|
| 301 | 0x1.a309bec4a2d3358c171f770db1f4p-1L,
|
---|
| 302 | 0x1.a5503b23e255c8b424491caf88ccp-1L,
|
---|
| 303 | 0x1.a799e1330b3586f2dfb2b158f31ep-1L,
|
---|
| 304 | 0x1.a9e6b5579fdbf43eb243bdff53a2p-1L,
|
---|
| 305 | 0x1.ac36bbfd3f379c0db966a3126988p-1L,
|
---|
| 306 | 0x1.ae89f995ad3ad5e8734d17731c80p-1L,
|
---|
| 307 | 0x1.b0e07298db66590842acdfc6fb4ep-1L,
|
---|
| 308 | 0x1.b33a2b84f15faf6bfd0e7bd941b0p-1L,
|
---|
| 309 | 0x1.b59728de559398e3881111648738p-1L,
|
---|
| 310 | 0x1.b7f76f2fb5e46eaa7b081ab53ff6p-1L,
|
---|
| 311 | 0x1.ba5b030a10649840cb3c6af5b74cp-1L,
|
---|
| 312 | 0x1.bcc1e904bc1d2247ba0f45b3d06cp-1L,
|
---|
| 313 | 0x1.bf2c25bd71e088408d7025190cd0p-1L,
|
---|
| 314 | 0x1.c199bdd85529c2220cb12a0916bap-1L,
|
---|
| 315 | 0x1.c40ab5fffd07a6d14df820f17deap-1L,
|
---|
| 316 | 0x1.c67f12e57d14b4a2137fd20f2a26p-1L,
|
---|
| 317 | 0x1.c8f6d9406e7b511acbc48805c3f6p-1L,
|
---|
| 318 | 0x1.cb720dcef90691503cbd1e949d0ap-1L,
|
---|
| 319 | 0x1.cdf0b555dc3f9c44f8958fac4f12p-1L,
|
---|
| 320 | 0x1.d072d4a07897b8d0f22f21a13792p-1L,
|
---|
| 321 | 0x1.d2f87080d89f18ade123989ea50ep-1L,
|
---|
| 322 | 0x1.d5818dcfba48725da05aeb66dff8p-1L,
|
---|
| 323 | 0x1.d80e316c98397bb84f9d048807a0p-1L,
|
---|
| 324 | 0x1.da9e603db3285708c01a5b6d480cp-1L,
|
---|
| 325 | 0x1.dd321f301b4604b695de3c0630c0p-1L,
|
---|
| 326 | 0x1.dfc97337b9b5eb968cac39ed284cp-1L,
|
---|
| 327 | 0x1.e264614f5a128a12761fa17adc74p-1L,
|
---|
| 328 | 0x1.e502ee78b3ff6273d130153992d0p-1L,
|
---|
| 329 | 0x1.e7a51fbc74c834b548b2832378a4p-1L,
|
---|
| 330 | 0x1.ea4afa2a490d9858f73a18f5dab4p-1L,
|
---|
| 331 | 0x1.ecf482d8e67f08db0312fb949d50p-1L,
|
---|
| 332 | 0x1.efa1bee615a27771fd21a92dabb6p-1L,
|
---|
| 333 | 0x1.f252b376bba974e8696fc3638f24p-1L,
|
---|
| 334 | 0x1.f50765b6e4540674f84b762861a6p-1L,
|
---|
| 335 | 0x1.f7bfdad9cbe138913b4bfe72bd78p-1L,
|
---|
| 336 | 0x1.fa7c1819e90d82e90a7e74b26360p-1L,
|
---|
| 337 | 0x1.fd3c22b8f71f10975ba4b32bd006p-1L,
|
---|
| 338 | 0x1.0000000000000000000000000000p+0L,
|
---|
| 339 | 0x1.0163da9fb33356d84a66ae336e98p+0L,
|
---|
| 340 | 0x1.02c9a3e778060ee6f7caca4f7a18p+0L,
|
---|
| 341 | 0x1.04315e86e7f84bd738f9a20da442p+0L,
|
---|
| 342 | 0x1.059b0d31585743ae7c548eb68c6ap+0L,
|
---|
| 343 | 0x1.0706b29ddf6ddc6dc403a9d87b1ep+0L,
|
---|
| 344 | 0x1.0874518759bc808c35f25d942856p+0L,
|
---|
| 345 | 0x1.09e3ecac6f3834521e060c584d5cp+0L,
|
---|
| 346 | 0x1.0b5586cf9890f6298b92b7184200p+0L,
|
---|
| 347 | 0x1.0cc922b7247f7407b705b893dbdep+0L,
|
---|
| 348 | 0x1.0e3ec32d3d1a2020742e4f8af794p+0L,
|
---|
| 349 | 0x1.0fb66affed31af232091dd8a169ep+0L,
|
---|
| 350 | 0x1.11301d0125b50a4ebbf1aed9321cp+0L,
|
---|
| 351 | 0x1.12abdc06c31cbfb92bad324d6f84p+0L,
|
---|
| 352 | 0x1.1429aaea92ddfb34101943b2588ep+0L,
|
---|
| 353 | 0x1.15a98c8a58e512480d573dd562aep+0L,
|
---|
| 354 | 0x1.172b83c7d517adcdf7c8c50eb162p+0L,
|
---|
| 355 | 0x1.18af9388c8de9bbbf70b9a3c269cp+0L,
|
---|
| 356 | 0x1.1a35beb6fcb753cb698f692d2038p+0L,
|
---|
| 357 | 0x1.1bbe084045cd39ab1e72b442810ep+0L,
|
---|
| 358 | 0x1.1d4873168b9aa7805b8028990be8p+0L,
|
---|
| 359 | 0x1.1ed5022fcd91cb8819ff61121fbep+0L,
|
---|
| 360 | 0x1.2063b88628cd63b8eeb0295093f6p+0L,
|
---|
| 361 | 0x1.21f49917ddc962552fd29294bc20p+0L,
|
---|
| 362 | 0x1.2387a6e75623866c1fadb1c159c0p+0L,
|
---|
| 363 | 0x1.251ce4fb2a63f3582ab7de9e9562p+0L,
|
---|
| 364 | 0x1.26b4565e27cdd257a673281d3068p+0L,
|
---|
| 365 | 0x1.284dfe1f5638096cf15cf03c9fa0p+0L,
|
---|
| 366 | 0x1.29e9df51fdee12c25d15f5a25022p+0L,
|
---|
| 367 | 0x1.2b87fd0dad98ffddea46538fca24p+0L,
|
---|
| 368 | 0x1.2d285a6e4030b40091d536d0733ep+0L,
|
---|
| 369 | 0x1.2ecafa93e2f5611ca0f45d5239a4p+0L,
|
---|
| 370 | 0x1.306fe0a31b7152de8d5a463063bep+0L,
|
---|
| 371 | 0x1.32170fc4cd8313539cf1c3009330p+0L,
|
---|
| 372 | 0x1.33c08b26416ff4c9c8610d96680ep+0L,
|
---|
| 373 | 0x1.356c55f929ff0c94623476373be4p+0L,
|
---|
| 374 | 0x1.371a7373aa9caa7145502f45452ap+0L,
|
---|
| 375 | 0x1.38cae6d05d86585a9cb0d9bed530p+0L,
|
---|
| 376 | 0x1.3a7db34e59ff6ea1bc9299e0a1fep+0L,
|
---|
| 377 | 0x1.3c32dc313a8e484001f228b58cf0p+0L,
|
---|
| 378 | 0x1.3dea64c12342235b41223e13d7eep+0L,
|
---|
| 379 | 0x1.3fa4504ac801ba0bf701aa417b9cp+0L,
|
---|
| 380 | 0x1.4160a21f72e29f84325b8f3dbacap+0L,
|
---|
| 381 | 0x1.431f5d950a896dc704439410b628p+0L,
|
---|
| 382 | 0x1.44e086061892d03136f409df0724p+0L,
|
---|
| 383 | 0x1.46a41ed1d005772512f459229f0ap+0L,
|
---|
| 384 | 0x1.486a2b5c13cd013c1a3b69062f26p+0L,
|
---|
| 385 | 0x1.4a32af0d7d3de672d8bcf46f99b4p+0L,
|
---|
| 386 | 0x1.4bfdad5362a271d4397afec42e36p+0L,
|
---|
| 387 | 0x1.4dcb299fddd0d63b36ef1a9e19dep+0L,
|
---|
| 388 | 0x1.4f9b2769d2ca6ad33d8b69aa0b8cp+0L,
|
---|
| 389 | 0x1.516daa2cf6641c112f52c84d6066p+0L,
|
---|
| 390 | 0x1.5342b569d4f81df0a83c49d86bf4p+0L,
|
---|
| 391 | 0x1.551a4ca5d920ec52ec620243540cp+0L,
|
---|
| 392 | 0x1.56f4736b527da66ecb004764e61ep+0L,
|
---|
| 393 | 0x1.58d12d497c7fd252bc2b7343d554p+0L,
|
---|
| 394 | 0x1.5ab07dd48542958c93015191e9a8p+0L,
|
---|
| 395 | 0x1.5c9268a5946b701c4b1b81697ed4p+0L,
|
---|
| 396 | 0x1.5e76f15ad21486e9be4c20399d12p+0L,
|
---|
| 397 | 0x1.605e1b976dc08b076f592a487066p+0L,
|
---|
| 398 | 0x1.6247eb03a5584b1f0fa06fd2d9eap+0L,
|
---|
| 399 | 0x1.6434634ccc31fc76f8714c4ee122p+0L,
|
---|
| 400 | 0x1.66238825522249127d9e29b92ea2p+0L,
|
---|
| 401 | 0x1.68155d44ca973081c57227b9f69ep+0L,
|
---|
| 402 | };
|
---|
| 403 |
|
---|
| 404 | static const float eps[TBLSIZE] = {
|
---|
| 405 | -0x1.5c50p-101,
|
---|
| 406 | -0x1.5d00p-106,
|
---|
| 407 | 0x1.8e90p-102,
|
---|
| 408 | -0x1.5340p-103,
|
---|
| 409 | 0x1.1bd0p-102,
|
---|
| 410 | -0x1.4600p-105,
|
---|
| 411 | -0x1.7a40p-104,
|
---|
| 412 | 0x1.d590p-102,
|
---|
| 413 | -0x1.d590p-101,
|
---|
| 414 | 0x1.b100p-103,
|
---|
| 415 | -0x1.0d80p-105,
|
---|
| 416 | 0x1.6b00p-103,
|
---|
| 417 | -0x1.9f00p-105,
|
---|
| 418 | 0x1.c400p-103,
|
---|
| 419 | 0x1.e120p-103,
|
---|
| 420 | -0x1.c100p-104,
|
---|
| 421 | -0x1.9d20p-103,
|
---|
| 422 | 0x1.a800p-108,
|
---|
| 423 | 0x1.4c00p-106,
|
---|
| 424 | -0x1.9500p-106,
|
---|
| 425 | 0x1.6900p-105,
|
---|
| 426 | -0x1.29d0p-100,
|
---|
| 427 | 0x1.4c60p-103,
|
---|
| 428 | 0x1.13a0p-102,
|
---|
| 429 | -0x1.5b60p-103,
|
---|
| 430 | -0x1.1c40p-103,
|
---|
| 431 | 0x1.db80p-102,
|
---|
| 432 | 0x1.91a0p-102,
|
---|
| 433 | 0x1.dc00p-105,
|
---|
| 434 | 0x1.44c0p-104,
|
---|
| 435 | 0x1.9710p-102,
|
---|
| 436 | 0x1.8760p-103,
|
---|
| 437 | -0x1.a720p-103,
|
---|
| 438 | 0x1.ed20p-103,
|
---|
| 439 | -0x1.49c0p-102,
|
---|
| 440 | -0x1.e000p-111,
|
---|
| 441 | 0x1.86a0p-103,
|
---|
| 442 | 0x1.2b40p-103,
|
---|
| 443 | -0x1.b400p-108,
|
---|
| 444 | 0x1.1280p-99,
|
---|
| 445 | -0x1.02d8p-102,
|
---|
| 446 | -0x1.e3d0p-103,
|
---|
| 447 | -0x1.b080p-105,
|
---|
| 448 | -0x1.f100p-107,
|
---|
| 449 | -0x1.16c0p-105,
|
---|
| 450 | -0x1.1190p-103,
|
---|
| 451 | -0x1.a7d2p-100,
|
---|
| 452 | 0x1.3450p-103,
|
---|
| 453 | -0x1.67c0p-105,
|
---|
| 454 | 0x1.4b80p-104,
|
---|
| 455 | -0x1.c4e0p-103,
|
---|
| 456 | 0x1.6000p-108,
|
---|
| 457 | -0x1.3f60p-105,
|
---|
| 458 | 0x1.93f0p-104,
|
---|
| 459 | 0x1.5fe0p-105,
|
---|
| 460 | 0x1.6f80p-107,
|
---|
| 461 | -0x1.7600p-106,
|
---|
| 462 | 0x1.21e0p-106,
|
---|
| 463 | -0x1.3a40p-106,
|
---|
| 464 | -0x1.40c0p-104,
|
---|
| 465 | -0x1.9860p-105,
|
---|
| 466 | -0x1.5d40p-108,
|
---|
| 467 | -0x1.1d70p-106,
|
---|
| 468 | 0x1.2760p-105,
|
---|
| 469 | 0x0.0000p+0,
|
---|
| 470 | 0x1.21e2p-104,
|
---|
| 471 | -0x1.9520p-108,
|
---|
| 472 | -0x1.5720p-106,
|
---|
| 473 | -0x1.4810p-106,
|
---|
| 474 | -0x1.be00p-109,
|
---|
| 475 | 0x1.0080p-105,
|
---|
| 476 | -0x1.5780p-108,
|
---|
| 477 | -0x1.d460p-105,
|
---|
| 478 | -0x1.6140p-105,
|
---|
| 479 | 0x1.4630p-104,
|
---|
| 480 | 0x1.ad50p-103,
|
---|
| 481 | 0x1.82e0p-105,
|
---|
| 482 | 0x1.1d3cp-101,
|
---|
| 483 | 0x1.6100p-107,
|
---|
| 484 | 0x1.ec30p-104,
|
---|
| 485 | 0x1.f200p-108,
|
---|
| 486 | 0x1.0b40p-103,
|
---|
| 487 | 0x1.3660p-102,
|
---|
| 488 | 0x1.d9d0p-103,
|
---|
| 489 | -0x1.02d0p-102,
|
---|
| 490 | 0x1.b070p-103,
|
---|
| 491 | 0x1.b9c0p-104,
|
---|
| 492 | -0x1.01c0p-103,
|
---|
| 493 | -0x1.dfe0p-103,
|
---|
| 494 | 0x1.1b60p-104,
|
---|
| 495 | -0x1.ae94p-101,
|
---|
| 496 | -0x1.3340p-104,
|
---|
| 497 | 0x1.b3d8p-102,
|
---|
| 498 | -0x1.6e40p-105,
|
---|
| 499 | -0x1.3670p-103,
|
---|
| 500 | 0x1.c140p-104,
|
---|
| 501 | 0x1.1840p-101,
|
---|
| 502 | 0x1.1ab0p-102,
|
---|
| 503 | -0x1.a400p-104,
|
---|
| 504 | 0x1.1f00p-104,
|
---|
| 505 | -0x1.7180p-103,
|
---|
| 506 | 0x1.4ce0p-102,
|
---|
| 507 | 0x1.9200p-107,
|
---|
| 508 | -0x1.54c0p-103,
|
---|
| 509 | 0x1.1b80p-105,
|
---|
| 510 | -0x1.1828p-101,
|
---|
| 511 | 0x1.5720p-102,
|
---|
| 512 | -0x1.a060p-100,
|
---|
| 513 | 0x1.9160p-102,
|
---|
| 514 | 0x1.a280p-104,
|
---|
| 515 | 0x1.3400p-107,
|
---|
| 516 | 0x1.2b20p-102,
|
---|
| 517 | 0x1.7800p-108,
|
---|
| 518 | 0x1.cfd0p-101,
|
---|
| 519 | 0x1.2ef0p-102,
|
---|
| 520 | -0x1.2760p-99,
|
---|
| 521 | 0x1.b380p-104,
|
---|
| 522 | 0x1.0048p-101,
|
---|
| 523 | -0x1.60b0p-102,
|
---|
| 524 | 0x1.a1ccp-100,
|
---|
| 525 | -0x1.a640p-104,
|
---|
| 526 | -0x1.08a0p-101,
|
---|
| 527 | 0x1.7e60p-102,
|
---|
| 528 | 0x1.22c0p-103,
|
---|
| 529 | -0x1.7200p-106,
|
---|
| 530 | 0x1.f0f0p-102,
|
---|
| 531 | 0x1.eb4ep-99,
|
---|
| 532 | 0x1.c6e0p-103,
|
---|
| 533 | };
|
---|
| 534 |
|
---|
| 535 | /*
|
---|
| 536 | * exp2l(x): compute the base 2 exponential of x
|
---|
| 537 | *
|
---|
| 538 | * Accuracy: Peak error < 0.502 ulp.
|
---|
| 539 | *
|
---|
| 540 | * Method: (accurate tables)
|
---|
| 541 | *
|
---|
| 542 | * Reduce x:
|
---|
| 543 | * x = 2**k + y, for integer k and |y| <= 1/2.
|
---|
| 544 | * Thus we have exp2(x) = 2**k * exp2(y).
|
---|
| 545 | *
|
---|
| 546 | * Reduce y:
|
---|
| 547 | * y = i/TBLSIZE + z - eps[i] for integer i near y * TBLSIZE.
|
---|
| 548 | * Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z - eps[i]),
|
---|
| 549 | * with |z - eps[i]| <= 2**-8 + 2**-98 for the table used.
|
---|
| 550 | *
|
---|
| 551 | * We compute exp2(i/TBLSIZE) via table lookup and exp2(z - eps[i]) via
|
---|
| 552 | * a degree-10 minimax polynomial with maximum error under 2**-120.
|
---|
| 553 | * The values in exp2t[] and eps[] are chosen such that
|
---|
| 554 | * exp2t[i] = exp2(i/TBLSIZE + eps[i]), and eps[i] is a small offset such
|
---|
| 555 | * that exp2t[i] is accurate to 2**-122.
|
---|
| 556 | *
|
---|
| 557 | * Note that the range of i is +-TBLSIZE/2, so we actually index the tables
|
---|
| 558 | * by i0 = i + TBLSIZE/2.
|
---|
| 559 | *
|
---|
| 560 | * This method is due to Gal, with many details due to Gal and Bachelis:
|
---|
| 561 | *
|
---|
| 562 | * Gal, S. and Bachelis, B. An Accurate Elementary Mathematical Library
|
---|
| 563 | * for the IEEE Floating Point Standard. TOMS 17(1), 26-46 (1991).
|
---|
| 564 | */
|
---|
| 565 | long double
|
---|
| 566 | exp2l(long double x)
|
---|
| 567 | {
|
---|
| 568 | union ldshape u = {x};
|
---|
| 569 | int e = u.i.se & 0x7fff;
|
---|
| 570 | long double r, z, t;
|
---|
| 571 | uint32_t i0;
|
---|
| 572 | union {uint32_t u; int32_t i;} k;
|
---|
| 573 |
|
---|
| 574 | /* Filter out exceptional cases. */
|
---|
| 575 | if (e >= 0x3fff + 14) { /* |x| >= 16384 or x is NaN */
|
---|
| 576 | if (u.i.se >= 0x3fff + 15 && u.i.se >> 15 == 0)
|
---|
| 577 | /* overflow */
|
---|
| 578 | return x * 0x1p16383L;
|
---|
| 579 | if (e == 0x7fff) /* -inf or -nan */
|
---|
| 580 | return -1/x;
|
---|
| 581 | if (x < -16382) {
|
---|
| 582 | if (x <= -16495 || x - 0x1p112 + 0x1p112 != x)
|
---|
| 583 | /* underflow */
|
---|
| 584 | FORCE_EVAL((float)(-0x1p-149/x));
|
---|
| 585 | if (x <= -16446)
|
---|
| 586 | return 0;
|
---|
| 587 | }
|
---|
| 588 | } else if (e < 0x3fff - 114) {
|
---|
| 589 | return 1 + x;
|
---|
| 590 | }
|
---|
| 591 |
|
---|
| 592 | /*
|
---|
| 593 | * Reduce x, computing z, i0, and k. The low bits of x + redux
|
---|
| 594 | * contain the 16-bit integer part of the exponent (k) followed by
|
---|
| 595 | * TBLBITS fractional bits (i0). We use bit tricks to extract these
|
---|
| 596 | * as integers, then set z to the remainder.
|
---|
| 597 | *
|
---|
| 598 | * Example: Suppose x is 0xabc.123456p0 and TBLBITS is 8.
|
---|
| 599 | * Then the low-order word of x + redux is 0x000abc12,
|
---|
| 600 | * We split this into k = 0xabc and i0 = 0x12 (adjusted to
|
---|
| 601 | * index into the table), then we compute z = 0x0.003456p0.
|
---|
| 602 | */
|
---|
| 603 | u.f = x + redux;
|
---|
| 604 | i0 = u.i2.lo + TBLSIZE / 2;
|
---|
| 605 | k.u = i0 / TBLSIZE * TBLSIZE;
|
---|
| 606 | k.i /= TBLSIZE;
|
---|
| 607 | i0 %= TBLSIZE;
|
---|
| 608 | u.f -= redux;
|
---|
| 609 | z = x - u.f;
|
---|
| 610 |
|
---|
| 611 | /* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
|
---|
| 612 | t = tbl[i0];
|
---|
| 613 | z -= eps[i0];
|
---|
| 614 | r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * (P5 + z * (P6
|
---|
| 615 | + z * (P7 + z * (P8 + z * (P9 + z * P10)))))))));
|
---|
| 616 |
|
---|
| 617 | return scalbnl(r, k.i);
|
---|
| 618 | }
|
---|
| 619 | #endif
|
---|