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
|
---|