1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 
  26 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __expm1l = expm1l
  32 
  33 #if !defined(__sparc)
  34 #error Unsupported architecture
  35 #endif
  36 
  37 /*
  38  * expm1l(x)
  39  *
  40  * Table driven method
  41  * Written by K.C. Ng, June 1995.
  42  * Algorithm :
  43  *      1. expm1(x) = x if x<2**-114
  44  *      2. if |x| <= 0.0625 = 1/16, use approximation
  45  *              expm1(x) = x + x*P/(2-P)
  46  * where
  47  *      P = x - z*(P1+z*(P2+z*(P3+z*(P4+z*(P5+z*P6+z*P7))))), z = x*x;
  48  * (this formula is derived from
  49  *      2-P+x = R = x*(exp(x)+1)/(exp(x)-1) ~ 2 + x*x/6 - x^4/360 + ...)
  50  *
  51  * P1 =   1.66666666666666666666666666666638500528074603030e-0001
  52  * P2 =  -2.77777777777777777777777759668391122822266551158e-0003
  53  * P3 =   6.61375661375661375657437408890138814721051293054e-0005
  54  * P4 =  -1.65343915343915303310185228411892601606669528828e-0006
  55  * P5 =   4.17535139755122945763580609663414647067443411178e-0008
  56  * P6 =  -1.05683795988668526689182102605260986731620026832e-0009
  57  * P7 =   2.67544168821852702827123344217198187229611470514e-0011
  58  *
  59  * Accuracy: |R-x*(exp(x)+1)/(exp(x)-1)|<=2**-119.13
  60  *
  61  *      3. For 1/16 < |x| < 1.125, choose x(+-i) ~ +-(i+4.5)/64, i=0,..,67
  62  *         since
  63  *              exp(x) = exp(xi+(x-xi))= exp(xi)*exp((x-xi))
  64  *         we have
  65  *              expm1(x) = expm1(xi)+(exp(xi))*(expm1(x-xi))
  66  *         where
  67  *              |s=x-xi| <= 1/128
  68  *         and
  69  *      expm1(s)=2s/(2-R), R= s-s^2*(T1+s^2*(T2+s^2*(T3+s^2*(T4+s^2*T5))))
  70  *
  71  * T1 =   1.666666666666666666666666666660876387437e-1L,
  72  * T2 =  -2.777777777777777777777707812093173478756e-3L,
  73  * T3 =   6.613756613756613482074280932874221202424e-5L,
  74  * T4 =  -1.653439153392139954169609822742235851120e-6L,
  75  * T5 =   4.175314851769539751387852116610973796053e-8L;
  76  *
  77  *      4. For |x| >= 1.125, return exp(x)-1.
  78  *          (see algorithm for exp)
  79  *
  80  * Special cases:
  81  *      expm1l(INF) is INF, expm1l(NaN) is NaN;
  82  *      expm1l(-INF)= -1;
  83  *      for finite argument, only expm1l(0)=0 is exact.
  84  *
  85  * Accuracy:
  86  *      according to an error analysis, the error is always less than
  87  *      2 ulp (unit in the last place).
  88  *
  89  * Misc. info.
  90  *      For 113 bit long double
  91  *              if x >  1.135652340629414394949193107797076342845e+4
  92  *      then expm1l(x) overflow;
  93  *
  94  * Constants:
  95  * Only decimal values are given. We assume that the compiler will convert
  96  * from decimal to binary accurately enough to produce the correct
  97  * hexadecimal values.
  98  */
  99 
 100 #include "libm.h"
 101 
 102 extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
 103 extern const long double _TBL_expm1lx[], _TBL_expm1l[];
 104 static const long double zero = +0.0L,
 105         one = +1.0L,
 106         two = +2.0L,
 107         ln2_64 = +1.083042469624914545964425189778400898568e-2L,
 108         ovflthreshold = +1.135652340629414394949193107797076342845e+4L,
 109         invln2_32 = +4.616624130844682903551758979206054839765e+1L,
 110         ln2_32hi = +2.166084939249829091928849858592451515688e-2L,
 111         ln2_32lo = +5.209643502595475652782654157501186731779e-27L,
 112         huge = +1.0e4000L,
 113         tiny = +1.0e-4000L,
 114         P1 = +1.66666666666666666666666666666638500528074603030e-0001L,
 115         P2 = -2.77777777777777777777777759668391122822266551158e-0003L,
 116         P3 = +6.61375661375661375657437408890138814721051293054e-0005L,
 117         P4 = -1.65343915343915303310185228411892601606669528828e-0006L,
 118         P5 = +4.17535139755122945763580609663414647067443411178e-0008L,
 119         P6 = -1.05683795988668526689182102605260986731620026832e-0009L,
 120         P7 = +2.67544168821852702827123344217198187229611470514e-0011L,
 121 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
 122         T1 = +1.666666666666666666666666666660876387437e-1L,
 123         T2 = -2.777777777777777777777707812093173478756e-3L,
 124         T3 = +6.613756613756613482074280932874221202424e-5L,
 125         T4 = -1.653439153392139954169609822742235851120e-6L,
 126         T5 = +4.175314851769539751387852116610973796053e-8L;
 127 
 128 long double
 129 expm1l(long double x)
 130 {
 131         int hx, ix, j, k, m;
 132         long double t, r, s, w;
 133 
 134         hx = ((int *)&x)[HIXWORD];
 135         ix = hx & ~0x80000000;
 136 
 137         if (ix >= 0x7fff0000) {
 138                 if (x != x)
 139                         return (x + x);                 /* NaN */
 140 
 141                 if (x < zero)
 142                         return (-one);                  /* -inf */
 143 
 144                 return (x);                             /* +inf */
 145         }
 146 
 147         if (ix < 0x3fff4000) {                               /* |x| < 1.25 */
 148                 if (ix < 0x3ffb0000) {                       /* |x| < 0.0625 */
 149                         if (ix < 0x3f8d0000) {
 150                                 if ((int)x == 0)
 151                                         return (x);     /* |x|<2^-114 */
 152                         }
 153 
 154                         t = x * x;
 155                         r = (x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t *
 156                             (P5 + t * (P6 + t * P7)))))));
 157                         return (x + (x * r) / (two - r));
 158                 }
 159 
 160                 /* compute i = [64*x] */
 161                 m = 0x4009 - (ix >> 16);
 162                 j = ((ix & 0x0000ffff) | 0x10000) >> m;       /* j=4,...,67 */
 163 
 164                 if (hx < 0)
 165                         j += 82;                        /* negative */
 166 
 167                 s = x - _TBL_expm1lx[j];
 168                 t = s * s;
 169                 r = s - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))));
 170                 r = (s + s) / (two - r);
 171                 w = _TBL_expm1l[j];
 172                 return (w + (w + one) * r);
 173         }
 174 
 175         if (hx > 0) {
 176                 if (x > ovflthreshold)
 177                         return (huge * huge);
 178 
 179                 k = (int)(invln2_32 * (x + ln2_64));
 180         } else {
 181                 if (x < -80.0)
 182                         return (tiny - x / x);
 183 
 184                 k = (int)(invln2_32 * (x - ln2_64));
 185         }
 186 
 187         j = k & 0x1f;
 188         m = k >> 5;
 189         t = (long double)k;
 190         x = (x - t * ln2_32hi) - t * ln2_32lo;
 191         t = x * x;
 192         r = (x - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))))) - two;
 193         x = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r -
 194             _TBL_expl_lo[j]);
 195         return (scalbnl(x, m) - one);
 196 }