1175460Sdas/*- 2175460Sdas * Copyright (c) 2005-2008 David Schultz <das@FreeBSD.ORG> 3175460Sdas * All rights reserved. 4175460Sdas * 5175460Sdas * Redistribution and use in source and binary forms, with or without 6175460Sdas * modification, are permitted provided that the following conditions 7175460Sdas * are met: 8175460Sdas * 1. Redistributions of source code must retain the above copyright 9175460Sdas * notice, this list of conditions and the following disclaimer. 10175460Sdas * 2. Redistributions in binary form must reproduce the above copyright 11175460Sdas * notice, this list of conditions and the following disclaimer in the 12175460Sdas * documentation and/or other materials provided with the distribution. 13175460Sdas * 14175460Sdas * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15175460Sdas * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16175460Sdas * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17175460Sdas * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18175460Sdas * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19175460Sdas * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20175460Sdas * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21175460Sdas * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22175460Sdas * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23175460Sdas * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24175460Sdas * SUCH DAMAGE. 25175460Sdas */ 26175460Sdas 27175460Sdas#include <sys/cdefs.h> 28175460Sdas__FBSDID("$FreeBSD$"); 29175460Sdas 30175460Sdas#include <float.h> 31175460Sdas#include <stdint.h> 32175460Sdas 33175460Sdas#ifdef __i386__ 34175460Sdas#include <ieeefp.h> 35175460Sdas#endif 36175460Sdas 37175460Sdas#include "fpmath.h" 38175460Sdas#include "math.h" 39175460Sdas 40175460Sdas#define TBLBITS 7 41175460Sdas#define TBLSIZE (1 << TBLBITS) 42175460Sdas 43175460Sdas#define BIAS (LDBL_MAX_EXP - 1) 44175460Sdas#define EXPMASK (BIAS + LDBL_MAX_EXP) 45175460Sdas 46175460Sdasstatic const long double huge = 0x1p10000L; 47175460Sdas#if 0 /* XXX Prevent gcc from erroneously constant folding this. */ 48175460Sdasstatic const long double twom10000 = 0x1p-10000L; 49175460Sdas#else 50175460Sdasstatic volatile long double twom10000 = 0x1p-10000L; 51175460Sdas#endif 52175460Sdas 53175460Sdasstatic const double 54175460Sdas redux = 0x1.8p63 / TBLSIZE, 55175460Sdas P1 = 0x1.62e42fefa39efp-1, 56175460Sdas P2 = 0x1.ebfbdff82c58fp-3, 57175460Sdas P3 = 0x1.c6b08d7049fap-5, 58175460Sdas P4 = 0x1.3b2ab6fba4da5p-7, 59175460Sdas P5 = 0x1.5d8804780a736p-10, 60175460Sdas P6 = 0x1.430918835e33dp-13; 61175460Sdas 62175460Sdasstatic const double tbl[TBLSIZE * 2] = { 63175460Sdas 0x1.6a09e667f3bcdp-1, -0x1.bdd3413b2648p-55, 64175460Sdas 0x1.6c012750bdabfp-1, -0x1.2895667ff0cp-57, 65175460Sdas 0x1.6dfb23c651a2fp-1, -0x1.bbe3a683c88p-58, 66175460Sdas 0x1.6ff7df9519484p-1, -0x1.83c0f25860fp-56, 67175460Sdas 0x1.71f75e8ec5f74p-1, -0x1.16e4786887bp-56, 68175460Sdas 0x1.73f9a48a58174p-1, -0x1.0a8d96c65d5p-55, 69175460Sdas 0x1.75feb564267c9p-1, -0x1.0245957316ep-55, 70175460Sdas 0x1.780694fde5d3fp-1, 0x1.866b80a0216p-55, 71175460Sdas 0x1.7a11473eb0187p-1, -0x1.41577ee0499p-56, 72175460Sdas 0x1.7c1ed0130c132p-1, 0x1.f124cd1164ep-55, 73175460Sdas 0x1.7e2f336cf4e62p-1, 0x1.05d02ba157ap-57, 74175460Sdas 0x1.80427543e1a12p-1, -0x1.27c86626d97p-55, 75175460Sdas 0x1.82589994cce13p-1, -0x1.d4c1dd41533p-55, 76175460Sdas 0x1.8471a4623c7adp-1, -0x1.8d684a341cep-56, 77175460Sdas 0x1.868d99b4492edp-1, -0x1.fc6f89bd4f68p-55, 78175460Sdas 0x1.88ac7d98a6699p-1, 0x1.994c2f37cb5p-55, 79175460Sdas 0x1.8ace5422aa0dbp-1, 0x1.6e9f156864bp-55, 80175460Sdas 0x1.8cf3216b5448cp-1, -0x1.0d55e32e9e4p-57, 81175460Sdas 0x1.8f1ae99157736p-1, 0x1.5cc13a2e397p-56, 82175460Sdas 0x1.9145b0b91ffc6p-1, -0x1.dd6792e5825p-55, 83175460Sdas 0x1.93737b0cdc5e5p-1, -0x1.75fc781b58p-58, 84175460Sdas 0x1.95a44cbc8520fp-1, -0x1.64b7c96a5fp-57, 85175460Sdas 0x1.97d829fde4e5p-1, -0x1.d185b7c1b86p-55, 86175460Sdas 0x1.9a0f170ca07bap-1, -0x1.173bd91cee6p-55, 87175460Sdas 0x1.9c49182a3f09p-1, 0x1.c7c46b071f2p-57, 88175460Sdas 0x1.9e86319e32323p-1, 0x1.824ca78e64cp-57, 89175460Sdas 0x1.a0c667b5de565p-1, -0x1.359495d1cd5p-55, 90175460Sdas 0x1.a309bec4a2d33p-1, 0x1.6305c7ddc368p-55, 91175460Sdas 0x1.a5503b23e255dp-1, -0x1.d2f6edb8d42p-55, 92175460Sdas 0x1.a799e1330b358p-1, 0x1.bcb7ecac564p-55, 93175460Sdas 0x1.a9e6b5579fdbfp-1, 0x1.0fac90ef7fdp-55, 94175460Sdas 0x1.ac36bbfd3f37ap-1, -0x1.f9234cae76dp-56, 95175460Sdas 0x1.ae89f995ad3adp-1, 0x1.7a1cd345dcc8p-55, 96175460Sdas 0x1.b0e07298db666p-1, -0x1.bdef54c80e4p-55, 97175460Sdas 0x1.b33a2b84f15fbp-1, -0x1.2805e3084d8p-58, 98175460Sdas 0x1.b59728de5593ap-1, -0x1.c71dfbbba6ep-55, 99175460Sdas 0x1.b7f76f2fb5e47p-1, -0x1.5584f7e54acp-57, 100175460Sdas 0x1.ba5b030a1064ap-1, -0x1.efcd30e5429p-55, 101175460Sdas 0x1.bcc1e904bc1d2p-1, 0x1.23dd07a2d9fp-56, 102175460Sdas 0x1.bf2c25bd71e09p-1, -0x1.efdca3f6b9c8p-55, 103175460Sdas 0x1.c199bdd85529cp-1, 0x1.11065895049p-56, 104175460Sdas 0x1.c40ab5fffd07ap-1, 0x1.b4537e083c6p-55, 105175460Sdas 0x1.c67f12e57d14bp-1, 0x1.2884dff483c8p-55, 106175460Sdas 0x1.c8f6d9406e7b5p-1, 0x1.1acbc48805cp-57, 107175460Sdas 0x1.cb720dcef9069p-1, 0x1.503cbd1e94ap-57, 108175460Sdas 0x1.cdf0b555dc3fap-1, -0x1.dd83b53829dp-56, 109175460Sdas 0x1.d072d4a07897cp-1, -0x1.cbc3743797a8p-55, 110175460Sdas 0x1.d2f87080d89f2p-1, -0x1.d487b719d858p-55, 111175460Sdas 0x1.d5818dcfba487p-1, 0x1.2ed02d75b37p-56, 112175460Sdas 0x1.d80e316c98398p-1, -0x1.11ec18bedep-55, 113175460Sdas 0x1.da9e603db3285p-1, 0x1.c2300696db5p-55, 114175460Sdas 0x1.dd321f301b46p-1, 0x1.2da5778f019p-55, 115175460Sdas 0x1.dfc97337b9b5fp-1, -0x1.1a5cd4f184b8p-55, 116175460Sdas 0x1.e264614f5a129p-1, -0x1.7b627817a148p-55, 117175460Sdas 0x1.e502ee78b3ff6p-1, 0x1.39e8980a9cdp-56, 118175460Sdas 0x1.e7a51fbc74c83p-1, 0x1.2d522ca0c8ep-55, 119175460Sdas 0x1.ea4afa2a490dap-1, -0x1.e9c23179c288p-55, 120175460Sdas 0x1.ecf482d8e67f1p-1, -0x1.c93f3b411ad8p-55, 121175460Sdas 0x1.efa1bee615a27p-1, 0x1.dc7f486a4b68p-55, 122175460Sdas 0x1.f252b376bba97p-1, 0x1.3a1a5bf0d8e8p-55, 123175460Sdas 0x1.f50765b6e454p-1, 0x1.9d3e12dd8a18p-55, 124175460Sdas 0x1.f7bfdad9cbe14p-1, -0x1.dbb12d00635p-55, 125175460Sdas 0x1.fa7c1819e90d8p-1, 0x1.74853f3a593p-56, 126175460Sdas 0x1.fd3c22b8f71f1p-1, 0x1.2eb74966578p-58, 127175460Sdas 0x1p+0, 0x0p+0, 128175460Sdas 0x1.0163da9fb3335p+0, 0x1.b61299ab8cd8p-54, 129175460Sdas 0x1.02c9a3e778061p+0, -0x1.19083535b08p-56, 130175460Sdas 0x1.04315e86e7f85p+0, -0x1.0a31c1977c98p-54, 131175460Sdas 0x1.059b0d3158574p+0, 0x1.d73e2a475b4p-55, 132175460Sdas 0x1.0706b29ddf6dep+0, -0x1.c91dfe2b13cp-55, 133175460Sdas 0x1.0874518759bc8p+0, 0x1.186be4bb284p-57, 134175460Sdas 0x1.09e3ecac6f383p+0, 0x1.14878183161p-54, 135175460Sdas 0x1.0b5586cf9890fp+0, 0x1.8a62e4adc61p-54, 136175460Sdas 0x1.0cc922b7247f7p+0, 0x1.01edc16e24f8p-54, 137175460Sdas 0x1.0e3ec32d3d1a2p+0, 0x1.03a1727c58p-59, 138175460Sdas 0x1.0fb66affed31bp+0, -0x1.b9bedc44ebcp-57, 139175460Sdas 0x1.11301d0125b51p+0, -0x1.6c51039449bp-54, 140175460Sdas 0x1.12abdc06c31ccp+0, -0x1.1b514b36ca8p-58, 141175460Sdas 0x1.1429aaea92dep+0, -0x1.32fbf9af1368p-54, 142175460Sdas 0x1.15a98c8a58e51p+0, 0x1.2406ab9eeabp-55, 143175460Sdas 0x1.172b83c7d517bp+0, -0x1.19041b9d78ap-55, 144175460Sdas 0x1.18af9388c8deap+0, -0x1.11023d1970f8p-54, 145175460Sdas 0x1.1a35beb6fcb75p+0, 0x1.e5b4c7b4969p-55, 146175460Sdas 0x1.1bbe084045cd4p+0, -0x1.95386352ef6p-54, 147175460Sdas 0x1.1d4873168b9aap+0, 0x1.e016e00a264p-54, 148175460Sdas 0x1.1ed5022fcd91dp+0, -0x1.1df98027bb78p-54, 149175460Sdas 0x1.2063b88628cd6p+0, 0x1.dc775814a85p-55, 150175460Sdas 0x1.21f49917ddc96p+0, 0x1.2a97e9494a6p-55, 151175460Sdas 0x1.2387a6e756238p+0, 0x1.9b07eb6c7058p-54, 152175460Sdas 0x1.251ce4fb2a63fp+0, 0x1.ac155bef4f5p-55, 153175460Sdas 0x1.26b4565e27cddp+0, 0x1.2bd339940eap-55, 154175460Sdas 0x1.284dfe1f56381p+0, -0x1.a4c3a8c3f0d8p-54, 155175460Sdas 0x1.29e9df51fdee1p+0, 0x1.612e8afad12p-55, 156175460Sdas 0x1.2b87fd0dad99p+0, -0x1.10adcd6382p-59, 157175460Sdas 0x1.2d285a6e4030bp+0, 0x1.0024754db42p-54, 158175460Sdas 0x1.2ecafa93e2f56p+0, 0x1.1ca0f45d524p-56, 159175460Sdas 0x1.306fe0a31b715p+0, 0x1.6f46ad23183p-55, 160175460Sdas 0x1.32170fc4cd831p+0, 0x1.a9ce78e1804p-55, 161175460Sdas 0x1.33c08b26416ffp+0, 0x1.327218436598p-54, 162175460Sdas 0x1.356c55f929ff1p+0, -0x1.b5cee5c4e46p-55, 163175460Sdas 0x1.371a7373aa9cbp+0, -0x1.63aeabf42ebp-54, 164175460Sdas 0x1.38cae6d05d866p+0, -0x1.e958d3c99048p-54, 165175460Sdas 0x1.3a7db34e59ff7p+0, -0x1.5e436d661f6p-56, 166175460Sdas 0x1.3c32dc313a8e5p+0, -0x1.efff8375d2ap-54, 167175460Sdas 0x1.3dea64c123422p+0, 0x1.ada0911f09fp-55, 168175460Sdas 0x1.3fa4504ac801cp+0, -0x1.7d023f956fap-54, 169175460Sdas 0x1.4160a21f72e2ap+0, -0x1.ef3691c309p-58, 170175460Sdas 0x1.431f5d950a897p+0, -0x1.1c7dde35f7ap-55, 171175460Sdas 0x1.44e086061892dp+0, 0x1.89b7a04ef8p-59, 172175460Sdas 0x1.46a41ed1d0057p+0, 0x1.c944bd1648a8p-54, 173175460Sdas 0x1.486a2b5c13cdp+0, 0x1.3c1a3b69062p-56, 174175460Sdas 0x1.4a32af0d7d3dep+0, 0x1.9cb62f3d1be8p-54, 175175460Sdas 0x1.4bfdad5362a27p+0, 0x1.d4397afec42p-56, 176175460Sdas 0x1.4dcb299fddd0dp+0, 0x1.8ecdbbc6a78p-54, 177175460Sdas 0x1.4f9b2769d2ca7p+0, -0x1.4b309d25958p-54, 178175460Sdas 0x1.516daa2cf6642p+0, -0x1.f768569bd94p-55, 179175460Sdas 0x1.5342b569d4f82p+0, -0x1.07abe1db13dp-55, 180175460Sdas 0x1.551a4ca5d920fp+0, -0x1.d689cefede6p-55, 181175460Sdas 0x1.56f4736b527dap+0, 0x1.9bb2c011d938p-54, 182175460Sdas 0x1.58d12d497c7fdp+0, 0x1.295e15b9a1ep-55, 183175460Sdas 0x1.5ab07dd485429p+0, 0x1.6324c0546478p-54, 184175460Sdas 0x1.5c9268a5946b7p+0, 0x1.c4b1b81698p-60, 185175460Sdas 0x1.5e76f15ad2148p+0, 0x1.ba6f93080e68p-54, 186175460Sdas 0x1.605e1b976dc09p+0, -0x1.3e2429b56de8p-54, 187175460Sdas 0x1.6247eb03a5585p+0, -0x1.383c17e40b48p-54, 188175460Sdas 0x1.6434634ccc32p+0, -0x1.c483c759d89p-55, 189175460Sdas 0x1.6623882552225p+0, -0x1.bb60987591cp-54, 190175460Sdas 0x1.68155d44ca973p+0, 0x1.038ae44f74p-57, 191175460Sdas}; 192175460Sdas 193175460Sdas/* 194175460Sdas * exp2l(x): compute the base 2 exponential of x 195175460Sdas * 196175460Sdas * Accuracy: Peak error < 0.511 ulp. 197175460Sdas * 198175460Sdas * Method: (equally-spaced tables) 199175460Sdas * 200175460Sdas * Reduce x: 201175460Sdas * x = 2**k + y, for integer k and |y| <= 1/2. 202175460Sdas * Thus we have exp2l(x) = 2**k * exp2(y). 203175460Sdas * 204175460Sdas * Reduce y: 205175460Sdas * y = i/TBLSIZE + z for integer i near y * TBLSIZE. 206175460Sdas * Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z), 207175460Sdas * with |z| <= 2**-(TBLBITS+1). 208175460Sdas * 209175460Sdas * We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a 210175460Sdas * degree-6 minimax polynomial with maximum error under 2**-69. 211175460Sdas * The table entries each have 104 bits of accuracy, encoded as 212175460Sdas * a pair of double precision values. 213175460Sdas */ 214175460Sdaslong double 215175460Sdasexp2l(long double x) 216175460Sdas{ 217176074Sbde union IEEEl2bits u, v; 218176074Sbde long double r, twopk, twopkp10000, z; 219175460Sdas uint32_t hx, ix, i0; 220175460Sdas int k; 221175460Sdas 222175460Sdas /* Filter out exceptional cases. */ 223175460Sdas u.e = x; 224175460Sdas hx = u.xbits.expsign; 225175460Sdas ix = hx & EXPMASK; 226175460Sdas if (ix >= BIAS + 14) { /* |x| >= 16384 or x is NaN */ 227175460Sdas if (ix == BIAS + LDBL_MAX_EXP) { 228175460Sdas if (u.xbits.man != 1ULL << 63 || (hx & 0x8000) == 0) 229176231Sbde return (x + x); /* x is +Inf or NaN */ 230175460Sdas else 231175460Sdas return (0.0); /* x is -Inf */ 232175460Sdas } 233175460Sdas if (x >= 16384) 234175460Sdas return (huge * huge); /* overflow */ 235175460Sdas if (x <= -16446) 236175460Sdas return (twom10000 * twom10000); /* underflow */ 237175460Sdas } else if (ix <= BIAS - 66) { /* |x| < 0x1p-66 */ 238175460Sdas return (1.0 + x); 239175460Sdas } 240175460Sdas 241175460Sdas#ifdef __i386__ 242175460Sdas /* 243175460Sdas * The default precision on i386 is 53 bits, so long doubles are 244175460Sdas * broken. Call exp2() to get an accurate (double precision) result. 245175460Sdas */ 246175460Sdas if (fpgetprec() != FP_PE) 247175460Sdas return (exp2(x)); 248175460Sdas#endif 249175460Sdas 250175460Sdas /* 251175460Sdas * Reduce x, computing z, i0, and k. The low bits of x + redux 252175460Sdas * contain the 16-bit integer part of the exponent (k) followed by 253175460Sdas * TBLBITS fractional bits (i0). We use bit tricks to extract these 254175460Sdas * as integers, then set z to the remainder. 255175460Sdas * 256175460Sdas * Example: Suppose x is 0xabc.123456p0 and TBLBITS is 8. 257175460Sdas * Then the low-order word of x + redux is 0x000abc12, 258175460Sdas * We split this into k = 0xabc and i0 = 0x12 (adjusted to 259175460Sdas * index into the table), then we compute z = 0x0.003456p0. 260175460Sdas * 261175460Sdas * XXX If the exponent is negative, the computation of k depends on 262175460Sdas * '>>' doing sign extension. 263175460Sdas */ 264175460Sdas u.e = x + redux; 265175460Sdas i0 = u.bits.manl + TBLSIZE / 2; 266175460Sdas k = (int)i0 >> TBLBITS; 267175460Sdas i0 = (i0 & (TBLSIZE - 1)) << 1; 268175460Sdas u.e -= redux; 269175460Sdas z = x - u.e; 270176074Sbde v.xbits.man = 1ULL << 63; 271176074Sbde if (k >= LDBL_MIN_EXP) { 272176074Sbde v.xbits.expsign = LDBL_MAX_EXP - 1 + k; 273176074Sbde twopk = v.e; 274176074Sbde } else { 275176074Sbde v.xbits.expsign = LDBL_MAX_EXP - 1 + k + 10000; 276176074Sbde twopkp10000 = v.e; 277176074Sbde } 278175460Sdas 279175460Sdas /* Compute r = exp2l(y) = exp2lt[i0] * p(z). */ 280175460Sdas long double t_hi = tbl[i0]; 281175460Sdas long double t_lo = tbl[i0 + 1]; 282175460Sdas /* XXX This gives > 1 ulp errors outside of FE_TONEAREST mode */ 283175460Sdas r = t_lo + (t_hi + t_lo) * z * (P1 + z * (P2 + z * (P3 + z * (P4 284175460Sdas + z * (P5 + z * P6))))) + t_hi; 285175460Sdas 286175460Sdas /* Scale by 2**k. */ 287175460Sdas if (k >= LDBL_MIN_EXP) { 288176074Sbde if (k == LDBL_MAX_EXP) 289176074Sbde return (r * 2.0 * 0x1p16383L); 290176074Sbde return (r * twopk); 291175460Sdas } else { 292176074Sbde return (r * twopkp10000 * twom10000); 293175460Sdas } 294175460Sdas} 295