newlib/libm/ld80/b_expl.c

   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 "../ld/fpmath.h"
  39 #include "../ld/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 }