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 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 #pragma weak expm1 = __expm1
31
32 /* INDENT OFF */
33 /*
34 * expm1(x)
35 * Returns exp(x)-1, the exponential of x minus 1.
36 *
37 * Method
38 * 1. Arugment reduction:
39 * Given x, find r and integer k such that
40 *
41 * x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658
42 *
43 * Here a correction term c will be computed to compensate
44 * the error in r when rounded to a floating-point number.
45 *
46 * 2. Approximating expm1(r) by a special rational function on
47 * the interval [0,0.34658]:
48 * Since
49 * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
50 * we define R1(r*r) by
110 * expm1(INF) is INF, expm1(NaN) is NaN;
111 * expm1(-INF) is -1, and
112 * for finite argument, only expm1(0)=0 is exact.
113 *
114 * Accuracy:
115 * according to an error analysis, the error is always less than
116 * 1 ulp (unit in the last place).
117 *
118 * Misc. info.
119 * For IEEE double
120 * if x > 7.09782712893383973096e+02 then expm1(x) overflow
121 *
122 * Constants:
123 * The hexadecimal values are the intended ones for the following
124 * constants. The decimal values may be used, provided that the
125 * compiler will convert from decimal to binary accurately enough
126 * to produce the hexadecimal values shown.
127 */
128 /* INDENT ON */
129
130 #include "libm_synonyms.h" /* __expm1 */
131 #include "libm_macros.h"
132 #include <math.h>
133
134 static const double xxx[] = {
135 /* one */ 1.0,
136 /* huge */ 1.0e+300,
137 /* tiny */ 1.0e-300,
138 /* o_threshold */ 7.09782712893383973096e+02, /* 40862E42 FEFA39EF */
139 /* ln2_hi */ 6.93147180369123816490e-01, /* 3FE62E42 FEE00000 */
140 /* ln2_lo */ 1.90821492927058770002e-10, /* 3DEA39EF 35793C76 */
141 /* invln2 */ 1.44269504088896338700e+00, /* 3FF71547 652B82FE */
142 /* scaled coefficients related to expm1 */
143 /* Q1 */ -3.33333333333331316428e-02, /* BFA11111 111110F4 */
144 /* Q2 */ 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */
145 /* Q3 */ -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
146 /* Q4 */ 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
147 /* Q5 */ -2.01099218183624371326e-07 /* BE8AFDB7 6E09C32D */
148 };
149 #define one xxx[0]
150 #define huge xxx[1]
|
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 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 #pragma weak __expm1 = expm1
31
32 /* INDENT OFF */
33 /*
34 * expm1(x)
35 * Returns exp(x)-1, the exponential of x minus 1.
36 *
37 * Method
38 * 1. Arugment reduction:
39 * Given x, find r and integer k such that
40 *
41 * x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658
42 *
43 * Here a correction term c will be computed to compensate
44 * the error in r when rounded to a floating-point number.
45 *
46 * 2. Approximating expm1(r) by a special rational function on
47 * the interval [0,0.34658]:
48 * Since
49 * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
50 * we define R1(r*r) by
110 * expm1(INF) is INF, expm1(NaN) is NaN;
111 * expm1(-INF) is -1, and
112 * for finite argument, only expm1(0)=0 is exact.
113 *
114 * Accuracy:
115 * according to an error analysis, the error is always less than
116 * 1 ulp (unit in the last place).
117 *
118 * Misc. info.
119 * For IEEE double
120 * if x > 7.09782712893383973096e+02 then expm1(x) overflow
121 *
122 * Constants:
123 * The hexadecimal values are the intended ones for the following
124 * constants. The decimal values may be used, provided that the
125 * compiler will convert from decimal to binary accurately enough
126 * to produce the hexadecimal values shown.
127 */
128 /* INDENT ON */
129
130 #include "libm_macros.h"
131 #include <math.h>
132
133 static const double xxx[] = {
134 /* one */ 1.0,
135 /* huge */ 1.0e+300,
136 /* tiny */ 1.0e-300,
137 /* o_threshold */ 7.09782712893383973096e+02, /* 40862E42 FEFA39EF */
138 /* ln2_hi */ 6.93147180369123816490e-01, /* 3FE62E42 FEE00000 */
139 /* ln2_lo */ 1.90821492927058770002e-10, /* 3DEA39EF 35793C76 */
140 /* invln2 */ 1.44269504088896338700e+00, /* 3FF71547 652B82FE */
141 /* scaled coefficients related to expm1 */
142 /* Q1 */ -3.33333333333331316428e-02, /* BFA11111 111110F4 */
143 /* Q2 */ 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */
144 /* Q3 */ -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
145 /* Q4 */ 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
146 /* Q5 */ -2.01099218183624371326e-07 /* BE8AFDB7 6E09C32D */
147 };
148 #define one xxx[0]
149 #define huge xxx[1]
|