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 "fpmath.h" 39 #include "math_private.h" 40 41 static const union IEEEl2bits 42 p0u = LD80C(0xaaaaaaaaaaaaaaab, -3, 1.66666666666666666671e-01L), 43 p1u = LD80C(0xb60b60b60b60b59a, -9, -2.77777777777777775377e-03L), 44 p2u = LD80C(0x8ab355e008a3cfce, -14, 6.61375661375629297465e-05L), 45 p3u = LD80C(0xddebbc994b0c1376, -20, -1.65343915327882529784e-06L), 46 p4u = LD80C(0xb354784cb4ef4c41, -25, 4.17535101591534118469e-08L), 47 p5u = LD80C(0x913e8a718382ce75, -30, -1.05679137034774806475e-09L), 48 p6u = LD80C(0xe8f0042aa134502e, -36, 2.64819349895429516863e-11L); 49 #define p1 (p0u.e) 50 #define p2 (p1u.e) 51 #define p3 (p2u.e) 52 #define p4 (p3u.e) 53 #define p5 (p4u.e) 54 #define p6 (p5u.e) 55 #define p7 (p6u.e) 56 57 /* 58 * lnhuge = (LDBL_MAX_EXP + 9) * log(2.) 59 * lntiny = (LDBL_MIN_EXP - 64 - 10) * log(2.) 60 * invln2 = 1 / log(2.) 61 */ 62 static const union IEEEl2bits 63 ln2hiu = LD80C(0xb17217f700000000, -1, 6.93147180369123816490e-01L), 64 ln2lou = LD80C(0xd1cf79abc9e3b398, -33, 1.90821492927058781614e-10L), 65 lnhugeu = LD80C(0xb18b0c0330a8fad9, 13, 1.13627617309191834574e+04L), 66 lntinyu = LD80C(0xb236f28a68bc3bd7, 13, -1.14057368561139000667e+04L), 67 invln2u = LD80C(0xb8aa3b295c17f0bc, 0, 1.44269504088896340739e+00L); 68 #define ln2hi (ln2hiu.e) 69 #define ln2lo (ln2lou.e) 70 #define lnhuge (lnhugeu.e) 71 #define lntiny (lntinyu.e) 72 #define invln2 (invln2u.e) 73 74 /* returns exp(r = x + c) for |c| < |x| with no overlap. */ 75 76 static long double 77 __exp__D(long double x, long double c) 78 { 79 long double hi, lo, z; 80 int k; 81 82 if (x != x) /* x is NaN. */ 83 return(x); 84 85 if (x <= lnhuge) { 86 if (x >= lntiny) { 87 /* argument reduction: x --> x - k*ln2 */ 88 z = invln2 * x; 89 k = z + copysignl(0.5L, x); 90 91 /* 92 * Express (x + c) - k * ln2 as hi - lo. 93 * Let x = hi - lo rounded. 94 */ 95 hi = x - k * ln2hi; /* Exact. */ 96 lo = k * ln2lo - c; 97 x = hi - lo; 98 99 /* Return 2^k*[1+x+x*c/(2+c)] */ 100 z = x * x; 101 c = x - z * (p1 + z * (p2 + z * (p3 + z * (p4 + 102 z * (p5 + z * (p6 + z * p7)))))); 103 c = (x * c) / (2 - c); 104 105 return (ldexpl(1 + (hi - (lo - c)), k)); 106 } else { 107 /* exp(-INF) is 0. exp(-big) underflows to 0. */ 108 return (isfinite(x) ? ldexpl(1., -5000) : 0); 109 } 110 } else 111 /* exp(INF) is INF, exp(+big#) overflows to INF */ 112 return (isfinite(x) ? ldexpl(1., 5000) : x); 113 } 114