1/*- 2 * SPDX-License-Identifier: BSD-3-Clause 3 * 4 * Copyright (c) 1985, 1993 5 * The Regents of the University of California. 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, this list of conditions and the following disclaimer. 12 * 2. Redistributions in binary form must reproduce the above copyright 13 * notice, this list of conditions and the following disclaimer in the 14 * documentation and/or other materials provided with the distribution. 15 * 3. Neither the name of the University nor the names of its contributors 16 * may be used to endorse or promote products derived from this software 17 * without specific prior written permission. 18 * 19 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 20 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 22 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 23 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 24 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 25 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 26 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 27 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 28 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 29 * SUCH DAMAGE. 30 */ 31 32/* 33 * See bsdsrc/b_exp.c for implementation details. 34 * 35 * bsdrc/b_exp.c converted to long double by Steven G. Kargl. 36 */ 37 38#include "math_private.h" 39 40static const union ieee_ext_u 41 p0u = LD80C(0xaaaaaaaaaaaaaaab, -3, 1.66666666666666666671e-01L), 42 p1u = LD80C(0xb60b60b60b60b59a, -9, -2.77777777777777775377e-03L), 43 p2u = LD80C(0x8ab355e008a3cfce, -14, 6.61375661375629297465e-05L), 44 p3u = LD80C(0xddebbc994b0c1376, -20, -1.65343915327882529784e-06L), 45 p4u = LD80C(0xb354784cb4ef4c41, -25, 4.17535101591534118469e-08L), 46 p5u = LD80C(0x913e8a718382ce75, -30, -1.05679137034774806475e-09L), 47 p6u = LD80C(0xe8f0042aa134502e, -36, 2.64819349895429516863e-11L); 48#define p1 (p0u.extu_ld) 49#define p2 (p1u.extu_ld) 50#define p3 (p2u.extu_ld) 51#define p4 (p3u.extu_ld) 52#define p5 (p4u.extu_ld) 53#define p6 (p5u.extu_ld) 54#define p7 (p6u.extu_ld) 55 56/* 57 * lnhuge = (LDBL_MAX_EXP + 9) * log(2.) 58 * lntiny = (LDBL_MIN_EXP - 64 - 10) * log(2.) 59 * invln2 = 1 / log(2.) 60 */ 61static const union ieee_ext_u 62ln2hiu = LD80C(0xb17217f700000000, -1, 6.93147180369123816490e-01L), 63ln2lou = LD80C(0xd1cf79abc9e3b398, -33, 1.90821492927058781614e-10L), 64lnhugeu = LD80C(0xb18b0c0330a8fad9, 13, 1.13627617309191834574e+04L), 65lntinyu = LD80C(0xb236f28a68bc3bd7, 13, -1.14057368561139000667e+04L), 66invln2u = LD80C(0xb8aa3b295c17f0bc, 0, 1.44269504088896340739e+00L); 67#define ln2hi (ln2hiu.extu_ld) 68#define ln2lo (ln2lou.extu_ld) 69#define lnhuge (lnhugeu.extu_ld) 70#define lntiny (lntinyu.extu_ld) 71#define invln2 (invln2u.extu_ld) 72 73/* returns exp(r = x + c) for |c| < |x| with no overlap. */ 74 75static long double 76__exp__D(long double x, long double c) 77{ 78 long double hi, lo, z; 79 int k; 80 81 if (x != x) /* x is NaN. */ 82 return(x); 83 84 if (x <= lnhuge) { 85 if (x >= lntiny) { 86 /* argument reduction: x --> x - k*ln2 */ 87 z = invln2 * x; 88 k = z + copysignl(0.5L, x); 89 90 /* 91 * Express (x + c) - k * ln2 as hi - lo. 92 * Let x = hi - lo rounded. 93 */ 94 hi = x - k * ln2hi; /* Exact. */ 95 lo = k * ln2lo - c; 96 x = hi - lo; 97 98 /* Return 2^k*[1+x+x*c/(2+c)] */ 99 z = x * x; 100 c = x - z * (p1 + z * (p2 + z * (p3 + z * (p4 + 101 z * (p5 + z * (p6 + z * p7)))))); 102 c = (x * c) / (2 - c); 103 104 return (ldexpl(1 + (hi - (lo - c)), k)); 105 } else { 106 /* exp(-INF) is 0. exp(-big) underflows to 0. */ 107 return (isfinite(x) ? ldexpl(1., -5000) : 0); 108 } 109 } else 110 /* exp(INF) is INF, exp(+big#) overflows to INF */ 111 return (isfinite(x) ? ldexpl(1., 5000) : x); 112} 113