[337] | 1 | /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_expl.c */
|
---|
| 2 | /*
|
---|
| 3 | * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
|
---|
| 4 | *
|
---|
| 5 | * Permission to use, copy, modify, and distribute this software for any
|
---|
| 6 | * purpose with or without fee is hereby granted, provided that the above
|
---|
| 7 | * copyright notice and this permission notice appear in all copies.
|
---|
| 8 | *
|
---|
| 9 | * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
|
---|
| 10 | * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
|
---|
| 11 | * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
|
---|
| 12 | * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
|
---|
| 13 | * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
|
---|
| 14 | * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
|
---|
| 15 | * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
|
---|
| 16 | */
|
---|
| 17 | /*
|
---|
| 18 | * Exponential function, long double precision
|
---|
| 19 | *
|
---|
| 20 | *
|
---|
| 21 | * SYNOPSIS:
|
---|
| 22 | *
|
---|
| 23 | * long double x, y, expl();
|
---|
| 24 | *
|
---|
| 25 | * y = expl( x );
|
---|
| 26 | *
|
---|
| 27 | *
|
---|
| 28 | * DESCRIPTION:
|
---|
| 29 | *
|
---|
| 30 | * Returns e (2.71828...) raised to the x power.
|
---|
| 31 | *
|
---|
| 32 | * Range reduction is accomplished by separating the argument
|
---|
| 33 | * into an integer k and fraction f such that
|
---|
| 34 | *
|
---|
| 35 | * x k f
|
---|
| 36 | * e = 2 e.
|
---|
| 37 | *
|
---|
| 38 | * A Pade' form of degree 5/6 is used to approximate exp(f) - 1
|
---|
| 39 | * in the basic range [-0.5 ln 2, 0.5 ln 2].
|
---|
| 40 | *
|
---|
| 41 | *
|
---|
| 42 | * ACCURACY:
|
---|
| 43 | *
|
---|
| 44 | * Relative error:
|
---|
| 45 | * arithmetic domain # trials peak rms
|
---|
| 46 | * IEEE +-10000 50000 1.12e-19 2.81e-20
|
---|
| 47 | *
|
---|
| 48 | *
|
---|
| 49 | * Error amplification in the exponential function can be
|
---|
| 50 | * a serious matter. The error propagation involves
|
---|
| 51 | * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
|
---|
| 52 | * which shows that a 1 lsb error in representing X produces
|
---|
| 53 | * a relative error of X times 1 lsb in the function.
|
---|
| 54 | * While the routine gives an accurate result for arguments
|
---|
| 55 | * that are exactly represented by a long double precision
|
---|
| 56 | * computer number, the result contains amplified roundoff
|
---|
| 57 | * error for large arguments not exactly represented.
|
---|
| 58 | *
|
---|
| 59 | *
|
---|
| 60 | * ERROR MESSAGES:
|
---|
| 61 | *
|
---|
| 62 | * message condition value returned
|
---|
| 63 | * exp underflow x < MINLOG 0.0
|
---|
| 64 | * exp overflow x > MAXLOG MAXNUM
|
---|
| 65 | *
|
---|
| 66 | */
|
---|
| 67 |
|
---|
| 68 | #include "libm.h"
|
---|
| 69 |
|
---|
| 70 | #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
|
---|
| 71 | long double expl(long double x)
|
---|
| 72 | {
|
---|
| 73 | return exp(x);
|
---|
| 74 | }
|
---|
| 75 | #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
|
---|
| 76 |
|
---|
| 77 | static const long double P[3] = {
|
---|
| 78 | 1.2617719307481059087798E-4L,
|
---|
| 79 | 3.0299440770744196129956E-2L,
|
---|
| 80 | 9.9999999999999999991025E-1L,
|
---|
| 81 | };
|
---|
| 82 | static const long double Q[4] = {
|
---|
| 83 | 3.0019850513866445504159E-6L,
|
---|
| 84 | 2.5244834034968410419224E-3L,
|
---|
| 85 | 2.2726554820815502876593E-1L,
|
---|
| 86 | 2.0000000000000000000897E0L,
|
---|
| 87 | };
|
---|
| 88 | static const long double
|
---|
| 89 | LN2HI = 6.9314575195312500000000E-1L,
|
---|
| 90 | LN2LO = 1.4286068203094172321215E-6L,
|
---|
| 91 | LOG2E = 1.4426950408889634073599E0L;
|
---|
| 92 |
|
---|
| 93 | long double expl(long double x)
|
---|
| 94 | {
|
---|
| 95 | long double px, xx;
|
---|
| 96 | int k;
|
---|
| 97 |
|
---|
| 98 | if (isnan(x))
|
---|
| 99 | return x;
|
---|
| 100 | if (x > 11356.5234062941439488L) /* x > ln(2^16384 - 0.5) */
|
---|
| 101 | return x * 0x1p16383L;
|
---|
| 102 | if (x < -11399.4985314888605581L) /* x < ln(2^-16446) */
|
---|
| 103 | return -0x1p-16445L/x;
|
---|
| 104 |
|
---|
| 105 | /* Express e**x = e**f 2**k
|
---|
| 106 | * = e**(f + k ln(2))
|
---|
| 107 | */
|
---|
| 108 | px = floorl(LOG2E * x + 0.5);
|
---|
| 109 | k = px;
|
---|
| 110 | x -= px * LN2HI;
|
---|
| 111 | x -= px * LN2LO;
|
---|
| 112 |
|
---|
| 113 | /* rational approximation of the fractional part:
|
---|
| 114 | * e**x = 1 + 2x P(x**2)/(Q(x**2) - x P(x**2))
|
---|
| 115 | */
|
---|
| 116 | xx = x * x;
|
---|
| 117 | px = x * __polevll(xx, P, 2);
|
---|
| 118 | x = px/(__polevll(xx, Q, 3) - px);
|
---|
| 119 | x = 1.0 + 2.0 * x;
|
---|
| 120 | return scalbnl(x, k);
|
---|
| 121 | }
|
---|
| 122 | #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
|
---|
| 123 | // TODO: broken implementation to make things compile
|
---|
| 124 | long double expl(long double x)
|
---|
| 125 | {
|
---|
| 126 | return exp(x);
|
---|
| 127 | }
|
---|
| 128 | #endif
|
---|