125c28e83SPiotr Jasiukajtis /*
225c28e83SPiotr Jasiukajtis * CDDL HEADER START
325c28e83SPiotr Jasiukajtis *
425c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e83SPiotr Jasiukajtis *
825c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e83SPiotr Jasiukajtis * and limitations under the License.
1225c28e83SPiotr Jasiukajtis *
1325c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e83SPiotr Jasiukajtis *
1925c28e83SPiotr Jasiukajtis * CDDL HEADER END
2025c28e83SPiotr Jasiukajtis */
2125c28e83SPiotr Jasiukajtis
2225c28e83SPiotr Jasiukajtis /*
2325c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
2425c28e83SPiotr Jasiukajtis */
2525c28e83SPiotr Jasiukajtis /*
2625c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
2725c28e83SPiotr Jasiukajtis * Use is subject to license terms.
2825c28e83SPiotr Jasiukajtis */
2925c28e83SPiotr Jasiukajtis
30*ddc0e0b5SRichard Lowe #pragma weak __expm1l = expm1l
31*ddc0e0b5SRichard Lowe
3225c28e83SPiotr Jasiukajtis #if !defined(__sparc)
3325c28e83SPiotr Jasiukajtis #error Unsupported architecture
3425c28e83SPiotr Jasiukajtis #endif
3525c28e83SPiotr Jasiukajtis
3625c28e83SPiotr Jasiukajtis /*
3725c28e83SPiotr Jasiukajtis * expm1l(x)
3825c28e83SPiotr Jasiukajtis *
3925c28e83SPiotr Jasiukajtis * Table driven method
4025c28e83SPiotr Jasiukajtis * Written by K.C. Ng, June 1995.
4125c28e83SPiotr Jasiukajtis * Algorithm :
4225c28e83SPiotr Jasiukajtis * 1. expm1(x) = x if x<2**-114
4325c28e83SPiotr Jasiukajtis * 2. if |x| <= 0.0625 = 1/16, use approximation
4425c28e83SPiotr Jasiukajtis * expm1(x) = x + x*P/(2-P)
4525c28e83SPiotr Jasiukajtis * where
4625c28e83SPiotr Jasiukajtis * P = x - z*(P1+z*(P2+z*(P3+z*(P4+z*(P5+z*P6+z*P7))))), z = x*x;
4725c28e83SPiotr Jasiukajtis * (this formula is derived from
4825c28e83SPiotr Jasiukajtis * 2-P+x = R = x*(exp(x)+1)/(exp(x)-1) ~ 2 + x*x/6 - x^4/360 + ...)
4925c28e83SPiotr Jasiukajtis *
5025c28e83SPiotr Jasiukajtis * P1 = 1.66666666666666666666666666666638500528074603030e-0001
5125c28e83SPiotr Jasiukajtis * P2 = -2.77777777777777777777777759668391122822266551158e-0003
5225c28e83SPiotr Jasiukajtis * P3 = 6.61375661375661375657437408890138814721051293054e-0005
5325c28e83SPiotr Jasiukajtis * P4 = -1.65343915343915303310185228411892601606669528828e-0006
5425c28e83SPiotr Jasiukajtis * P5 = 4.17535139755122945763580609663414647067443411178e-0008
5525c28e83SPiotr Jasiukajtis * P6 = -1.05683795988668526689182102605260986731620026832e-0009
5625c28e83SPiotr Jasiukajtis * P7 = 2.67544168821852702827123344217198187229611470514e-0011
5725c28e83SPiotr Jasiukajtis *
5825c28e83SPiotr Jasiukajtis * Accuracy: |R-x*(exp(x)+1)/(exp(x)-1)|<=2**-119.13
5925c28e83SPiotr Jasiukajtis *
6025c28e83SPiotr Jasiukajtis * 3. For 1/16 < |x| < 1.125, choose x(+-i) ~ +-(i+4.5)/64, i=0,..,67
6125c28e83SPiotr Jasiukajtis * since
6225c28e83SPiotr Jasiukajtis * exp(x) = exp(xi+(x-xi))= exp(xi)*exp((x-xi))
6325c28e83SPiotr Jasiukajtis * we have
6425c28e83SPiotr Jasiukajtis * expm1(x) = expm1(xi)+(exp(xi))*(expm1(x-xi))
6525c28e83SPiotr Jasiukajtis * where
6625c28e83SPiotr Jasiukajtis * |s=x-xi| <= 1/128
6725c28e83SPiotr Jasiukajtis * and
6825c28e83SPiotr Jasiukajtis * expm1(s)=2s/(2-R), R= s-s^2*(T1+s^2*(T2+s^2*(T3+s^2*(T4+s^2*T5))))
6925c28e83SPiotr Jasiukajtis *
7025c28e83SPiotr Jasiukajtis * T1 = 1.666666666666666666666666666660876387437e-1L,
7125c28e83SPiotr Jasiukajtis * T2 = -2.777777777777777777777707812093173478756e-3L,
7225c28e83SPiotr Jasiukajtis * T3 = 6.613756613756613482074280932874221202424e-5L,
7325c28e83SPiotr Jasiukajtis * T4 = -1.653439153392139954169609822742235851120e-6L,
7425c28e83SPiotr Jasiukajtis * T5 = 4.175314851769539751387852116610973796053e-8L;
7525c28e83SPiotr Jasiukajtis *
7625c28e83SPiotr Jasiukajtis * 4. For |x| >= 1.125, return exp(x)-1.
7725c28e83SPiotr Jasiukajtis * (see algorithm for exp)
7825c28e83SPiotr Jasiukajtis *
7925c28e83SPiotr Jasiukajtis * Special cases:
8025c28e83SPiotr Jasiukajtis * expm1l(INF) is INF, expm1l(NaN) is NaN;
8125c28e83SPiotr Jasiukajtis * expm1l(-INF)= -1;
8225c28e83SPiotr Jasiukajtis * for finite argument, only expm1l(0)=0 is exact.
8325c28e83SPiotr Jasiukajtis *
8425c28e83SPiotr Jasiukajtis * Accuracy:
8525c28e83SPiotr Jasiukajtis * according to an error analysis, the error is always less than
8625c28e83SPiotr Jasiukajtis * 2 ulp (unit in the last place).
8725c28e83SPiotr Jasiukajtis *
8825c28e83SPiotr Jasiukajtis * Misc. info.
8925c28e83SPiotr Jasiukajtis * For 113 bit long double
9025c28e83SPiotr Jasiukajtis * if x > 1.135652340629414394949193107797076342845e+4
9125c28e83SPiotr Jasiukajtis * then expm1l(x) overflow;
9225c28e83SPiotr Jasiukajtis *
9325c28e83SPiotr Jasiukajtis * Constants:
9425c28e83SPiotr Jasiukajtis * Only decimal values are given. We assume that the compiler will convert
9525c28e83SPiotr Jasiukajtis * from decimal to binary accurately enough to produce the correct
9625c28e83SPiotr Jasiukajtis * hexadecimal values.
9725c28e83SPiotr Jasiukajtis */
9825c28e83SPiotr Jasiukajtis
9925c28e83SPiotr Jasiukajtis #include "libm.h"
10025c28e83SPiotr Jasiukajtis
10125c28e83SPiotr Jasiukajtis extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
10225c28e83SPiotr Jasiukajtis extern const long double _TBL_expm1lx[], _TBL_expm1l[];
10325c28e83SPiotr Jasiukajtis
10425c28e83SPiotr Jasiukajtis static const long double
10525c28e83SPiotr Jasiukajtis zero = +0.0L,
10625c28e83SPiotr Jasiukajtis one = +1.0L,
10725c28e83SPiotr Jasiukajtis two = +2.0L,
10825c28e83SPiotr Jasiukajtis ln2_64 = +1.083042469624914545964425189778400898568e-2L,
10925c28e83SPiotr Jasiukajtis ovflthreshold = +1.135652340629414394949193107797076342845e+4L,
11025c28e83SPiotr Jasiukajtis invln2_32 = +4.616624130844682903551758979206054839765e+1L,
11125c28e83SPiotr Jasiukajtis ln2_32hi = +2.166084939249829091928849858592451515688e-2L,
11225c28e83SPiotr Jasiukajtis ln2_32lo = +5.209643502595475652782654157501186731779e-27L,
11325c28e83SPiotr Jasiukajtis huge = +1.0e4000L,
11425c28e83SPiotr Jasiukajtis tiny = +1.0e-4000L,
11525c28e83SPiotr Jasiukajtis P1 = +1.66666666666666666666666666666638500528074603030e-0001L,
11625c28e83SPiotr Jasiukajtis P2 = -2.77777777777777777777777759668391122822266551158e-0003L,
11725c28e83SPiotr Jasiukajtis P3 = +6.61375661375661375657437408890138814721051293054e-0005L,
11825c28e83SPiotr Jasiukajtis P4 = -1.65343915343915303310185228411892601606669528828e-0006L,
11925c28e83SPiotr Jasiukajtis P5 = +4.17535139755122945763580609663414647067443411178e-0008L,
12025c28e83SPiotr Jasiukajtis P6 = -1.05683795988668526689182102605260986731620026832e-0009L,
12125c28e83SPiotr Jasiukajtis P7 = +2.67544168821852702827123344217198187229611470514e-0011L,
12225c28e83SPiotr Jasiukajtis /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
12325c28e83SPiotr Jasiukajtis T1 = +1.666666666666666666666666666660876387437e-1L,
12425c28e83SPiotr Jasiukajtis T2 = -2.777777777777777777777707812093173478756e-3L,
12525c28e83SPiotr Jasiukajtis T3 = +6.613756613756613482074280932874221202424e-5L,
12625c28e83SPiotr Jasiukajtis T4 = -1.653439153392139954169609822742235851120e-6L,
12725c28e83SPiotr Jasiukajtis T5 = +4.175314851769539751387852116610973796053e-8L;
12825c28e83SPiotr Jasiukajtis
12925c28e83SPiotr Jasiukajtis long double
expm1l(long double x)13025c28e83SPiotr Jasiukajtis expm1l(long double x) {
13125c28e83SPiotr Jasiukajtis int hx, ix, j, k, m;
13225c28e83SPiotr Jasiukajtis long double t, r, s, w;
13325c28e83SPiotr Jasiukajtis
13425c28e83SPiotr Jasiukajtis hx = ((int *) &x)[HIXWORD];
13525c28e83SPiotr Jasiukajtis ix = hx & ~0x80000000;
13625c28e83SPiotr Jasiukajtis if (ix >= 0x7fff0000) {
13725c28e83SPiotr Jasiukajtis if (x != x)
13825c28e83SPiotr Jasiukajtis return (x + x); /* NaN */
13925c28e83SPiotr Jasiukajtis if (x < zero)
14025c28e83SPiotr Jasiukajtis return (-one); /* -inf */
14125c28e83SPiotr Jasiukajtis return (x); /* +inf */
14225c28e83SPiotr Jasiukajtis }
14325c28e83SPiotr Jasiukajtis if (ix < 0x3fff4000) { /* |x| < 1.25 */
14425c28e83SPiotr Jasiukajtis if (ix < 0x3ffb0000) { /* |x| < 0.0625 */
14525c28e83SPiotr Jasiukajtis if (ix < 0x3f8d0000) {
14625c28e83SPiotr Jasiukajtis if ((int) x == 0)
14725c28e83SPiotr Jasiukajtis return (x); /* |x|<2^-114 */
14825c28e83SPiotr Jasiukajtis }
14925c28e83SPiotr Jasiukajtis t = x * x;
15025c28e83SPiotr Jasiukajtis r = (x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t *
15125c28e83SPiotr Jasiukajtis (P5 + t * (P6 + t * P7)))))));
15225c28e83SPiotr Jasiukajtis return (x + (x * r) / (two - r));
15325c28e83SPiotr Jasiukajtis }
15425c28e83SPiotr Jasiukajtis /* compute i = [64*x] */
15525c28e83SPiotr Jasiukajtis m = 0x4009 - (ix >> 16);
15625c28e83SPiotr Jasiukajtis j = ((ix & 0x0000ffff) | 0x10000) >> m; /* j=4,...,67 */
15725c28e83SPiotr Jasiukajtis if (hx < 0)
15825c28e83SPiotr Jasiukajtis j += 82; /* negative */
15925c28e83SPiotr Jasiukajtis s = x - _TBL_expm1lx[j];
16025c28e83SPiotr Jasiukajtis t = s * s;
16125c28e83SPiotr Jasiukajtis r = s - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))));
16225c28e83SPiotr Jasiukajtis r = (s + s) / (two - r);
16325c28e83SPiotr Jasiukajtis w = _TBL_expm1l[j];
16425c28e83SPiotr Jasiukajtis return (w + (w + one) * r);
16525c28e83SPiotr Jasiukajtis }
16625c28e83SPiotr Jasiukajtis if (hx > 0) {
16725c28e83SPiotr Jasiukajtis if (x > ovflthreshold)
16825c28e83SPiotr Jasiukajtis return (huge * huge);
16925c28e83SPiotr Jasiukajtis k = (int) (invln2_32 * (x + ln2_64));
17025c28e83SPiotr Jasiukajtis } else {
17125c28e83SPiotr Jasiukajtis if (x < -80.0)
17225c28e83SPiotr Jasiukajtis return (tiny - x / x);
17325c28e83SPiotr Jasiukajtis k = (int) (invln2_32 * (x - ln2_64));
17425c28e83SPiotr Jasiukajtis }
17525c28e83SPiotr Jasiukajtis j = k & 0x1f;
17625c28e83SPiotr Jasiukajtis m = k >> 5;
17725c28e83SPiotr Jasiukajtis t = (long double) k;
17825c28e83SPiotr Jasiukajtis x = (x - t * ln2_32hi) - t * ln2_32lo;
17925c28e83SPiotr Jasiukajtis t = x * x;
18025c28e83SPiotr Jasiukajtis r = (x - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))))) - two;
18125c28e83SPiotr Jasiukajtis x = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r -
18225c28e83SPiotr Jasiukajtis _TBL_expl_lo[j]);
18325c28e83SPiotr Jasiukajtis return (scalbnl(x, m) - one);
18425c28e83SPiotr Jasiukajtis }
185