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 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 /* INDENT OFF */
31 /*
32 * double __k_cexp(double x, int *n);
33 * Returns the exponential of x in the form of 2**n * y, y=__k_cexp(x,&n).
34 *
35 * Method
36 * 1. Argument reduction:
37 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
38 * Given x, find r and integer k such that
39 *
40 * x = k*ln2 + r, |r| <= 0.5*ln2.
41 *
42 * Here r will be represented as r = hi-lo for better
43 * accuracy.
44 *
45 * 2. Approximation of exp(r) by a special rational function on
46 * the interval [0,0.34658]:
47 * Write
48 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
49 * We use a special Remez algorithm on [0,0.34658] to generate
50 * a polynomial of degree 5 to approximate R. The maximum error
70 * 3. Return n = k and __k_cexp = exp(r).
71 *
72 * Special cases:
73 * exp(INF) is INF, exp(NaN) is NaN;
74 * exp(-INF) is 0, and
75 * for finite argument, only exp(0)=1 is exact.
76 *
77 * Range and Accuracy:
78 * When |x| is really big, say |x| > 50000, the accuracy
79 * is not important because the ultimate result will over or under
80 * flow. So we will simply replace n = 50000 and r = 0.0. For
81 * moderate size x, according to an error analysis, the error is
82 * always less than 1 ulp (unit in the last place).
83 *
84 * Constants:
85 * The hexadecimal values are the intended ones for the following
86 * constants. The decimal values may be used, provided that the
87 * compiler will convert from decimal to binary accurately enough
88 * to produce the hexadecimal values shown.
89 */
90 /* INDENT ON */
91
92 #include "libm.h" /* __k_cexp */
93 #include "complex_wrapper.h" /* HI_WORD/LO_WORD */
94
95 /* INDENT OFF */
96 static const double
97 one = 1.0,
98 two128 = 3.40282366920938463463e+38,
99 halF[2] = {
100 0.5, -0.5,
101 },
102 ln2HI[2] = {
103 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
104 -6.93147180369123816490e-01, /* 0xbfe62e42, 0xfee00000 */
105 },
106 ln2LO[2] = {
107 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
108 -1.90821492927058770002e-10, /* 0xbdea39ef, 0x35793c76 */
109 },
110 invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
111 P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
112 P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
113 P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
114 P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
115 P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
116 /* INDENT ON */
117
118 double
119 __k_cexp(double x, int *n) {
120 double hi = 0.0L, lo = 0.0L, c, t;
121 int k, xsb;
122 unsigned hx, lx;
123
124 hx = HI_WORD(x); /* high word of x */
125 lx = LO_WORD(x); /* low word of x */
126 xsb = (hx >> 31) & 1; /* sign bit of x */
127 hx &= 0x7fffffff; /* high word of |x| */
128
129 /* filter out non-finite argument */
130 if (hx >= 0x40e86a00) { /* if |x| > 50000 */
131 if (hx >= 0x7ff00000) {
132 *n = 1;
133 if (((hx & 0xfffff) | lx) != 0)
134 return (x + x); /* NaN */
135 else
136 return ((xsb == 0) ? x : 0.0);
137 /* exp(+-inf)={inf,0} */
138 }
139 *n = (xsb == 0) ? 50000 : -50000;
140 return (one + ln2LO[1] * ln2LO[1]); /* generate inexact */
141 }
142
143 *n = 0;
144 /* argument reduction */
145 if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
146 if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
147 hi = x - ln2HI[xsb];
148 lo = ln2LO[xsb];
149 k = 1 - xsb - xsb;
150 } else {
151 k = (int) (invln2 * x + halF[xsb]);
152 t = k;
153 hi = x - t * ln2HI[0];
154 /* t*ln2HI is exact for t<2**20 */
155 lo = t * ln2LO[0];
156 }
157 x = hi - lo;
158 *n = k;
159 } else if (hx < 0x3e300000) { /* when |x|<2**-28 */
160 return (one + x);
161 } else
162 k = 0;
163
164 /* x is now in primary range */
165 t = x * x;
166 c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
167 if (k == 0)
168 return (one - ((x * c) / (c - 2.0) - x));
169 else {
170 t = one - ((lo - (x * c) / (2.0 - c)) - hi);
171 if (k > 128) {
172 t *= two128;
173 *n = k - 128;
174 } else if (k > 0) {
175 HI_WORD(t) += (k << 20);
176 *n = 0;
177 }
178 return (t);
179 }
180 }
|
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
32 /*
33 * double __k_cexp(double x, int *n);
34 * Returns the exponential of x in the form of 2**n * y, y=__k_cexp(x,&n).
35 *
36 * Method
37 * 1. Argument reduction:
38 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
39 * Given x, find r and integer k such that
40 *
41 * x = k*ln2 + r, |r| <= 0.5*ln2.
42 *
43 * Here r will be represented as r = hi-lo for better
44 * accuracy.
45 *
46 * 2. Approximation of exp(r) by a special rational function on
47 * the interval [0,0.34658]:
48 * Write
49 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
50 * We use a special Remez algorithm on [0,0.34658] to generate
51 * a polynomial of degree 5 to approximate R. The maximum error
71 * 3. Return n = k and __k_cexp = exp(r).
72 *
73 * Special cases:
74 * exp(INF) is INF, exp(NaN) is NaN;
75 * exp(-INF) is 0, and
76 * for finite argument, only exp(0)=1 is exact.
77 *
78 * Range and Accuracy:
79 * When |x| is really big, say |x| > 50000, the accuracy
80 * is not important because the ultimate result will over or under
81 * flow. So we will simply replace n = 50000 and r = 0.0. For
82 * moderate size x, according to an error analysis, the error is
83 * always less than 1 ulp (unit in the last place).
84 *
85 * Constants:
86 * The hexadecimal values are the intended ones for the following
87 * constants. The decimal values may be used, provided that the
88 * compiler will convert from decimal to binary accurately enough
89 * to produce the hexadecimal values shown.
90 */
91
92 #include "libm.h" /* __k_cexp */
93 #include "complex_wrapper.h" /* HI_WORD/LO_WORD */
94
95 static const double one = 1.0,
96 two128 = 3.40282366920938463463e+38;
97
98 static const double halF[2] = { 0.5, -0.5, };
99
100 static const double ln2HI[2] = {
101 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
102 -6.93147180369123816490e-01, /* 0xbfe62e42, 0xfee00000 */
103 };
104
105 static const double ln2LO[2] = {
106 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
107 -1.90821492927058770002e-10, /* 0xbdea39ef, 0x35793c76 */
108 };
109
110 static const double
111 invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
112 P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
113 P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
114 P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
115 P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
116 P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
117
118
119 double
120 __k_cexp(double x, int *n)
121 {
122 double hi = 0.0L, lo = 0.0L, c, t;
123 int k, xsb;
124 unsigned hx, lx;
125
126 hx = HI_WORD(x); /* high word of x */
127 lx = LO_WORD(x); /* low word of x */
128 xsb = (hx >> 31) & 1; /* sign bit of x */
129 hx &= 0x7fffffff; /* high word of |x| */
130
131 /* filter out non-finite argument */
132 if (hx >= 0x40e86a00) { /* if |x| > 50000 */
133 if (hx >= 0x7ff00000) {
134 *n = 1;
135
136 if (((hx & 0xfffff) | lx) != 0)
137 return (x + x); /* NaN */
138 else
139 return ((xsb == 0) ? x : 0.0);
140
141 /* exp(+-inf)={inf,0} */
142 }
143
144 *n = (xsb == 0) ? 50000 : -50000;
145 return (one + ln2LO[1] * ln2LO[1]); /* generate inexact */
146 }
147
148 *n = 0;
149
150 /* argument reduction */
151 if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
152 if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
153 hi = x - ln2HI[xsb];
154 lo = ln2LO[xsb];
155 k = 1 - xsb - xsb;
156 } else {
157 k = (int)(invln2 * x + halF[xsb]);
158 t = k;
159 hi = x - t * ln2HI[0];
160 /* t*ln2HI is exact for t<2**20 */
161 lo = t * ln2LO[0];
162 }
163
164 x = hi - lo;
165 *n = k;
166 } else if (hx < 0x3e300000) { /* when |x|<2**-28 */
167 return (one + x);
168 } else {
169 k = 0;
170 }
171
172 /* x is now in primary range */
173 t = x * x;
174 c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
175
176 if (k == 0) {
177 return (one - ((x * c) / (c - 2.0) - x));
178 } else {
179 t = one - ((lo - (x * c) / (2.0 - c)) - hi);
180
181 if (k > 128) {
182 t *= two128;
183 *n = k - 128;
184 } else if (k > 0) {
185 HI_WORD(t) += (k << 20);
186 *n = 0;
187 }
188
189 return (t);
190 }
191 }
|