1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
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 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
23 */
24 /*
25 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
26 * Use is subject to license terms.
27 */
28
29 #pragma weak __exp = exp
30
31 /*
32 * exp(x)
33 * Hybrid algorithm of Peter Tang's Table driven method (for large
34 * arguments) and an accurate table (for small arguments).
35 * Written by K.C. Ng, November 1988.
36 * Method (large arguments):
37 * 1. Argument Reduction: given the input x, find r and integer k
38 * and j such that
39 * x = (k+j/32)*(ln2) + r, |r| <= (1/64)*ln2
40 *
41 * 2. exp(x) = 2^k * (2^(j/32) + 2^(j/32)*expm1(r))
42 * a. expm1(r) is approximated by a polynomial:
43 * expm1(r) ~ r + t1*r^2 + t2*r^3 + ... + t5*r^6
266 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
267 7.45133219101941108420e+02, /* 0x40874910, 0xD52D3051 */
268 5.55111512312578270212e-17, /* 0x3c900000, 0x00000000 */
269 };
270
271 #define half C[0]
272 #define invln2_32 C[1]
273 #define ln2_32hi C[2]
274 #define ln2_32lo C[3]
275 #define t2 C[4]
276 #define t3 C[5]
277 #define t4 C[6]
278 #define t5 C[7]
279 #define one C[8]
280 #define zero C[9]
281 #define threshold1 C[10]
282 #define threshold2 C[11]
283 #define twom54 C[12]
284
285 double
286 exp(double x) {
287 double y, z, t;
288 int hx, ix, k, j, m;
289
290 ix = ((int *)&x)[HIWORD];
291 hx = ix & ~0x80000000;
292
293 if (hx < 0x3ff0a2b2) { /* |x| < 3/2 ln 2 */
294 if (hx < 0x3f862e42) { /* |x| < 1/64 ln 2 */
295 if (hx < 0x3ed00000) { /* |x| < 2^-18 */
296 volatile int dummy __unused;
297
298 dummy = (int)x; /* raise inexact if x != 0 */
299 #ifdef lint
300 dummy = dummy;
301 #endif
302 if (hx < 0x3e300000)
303 return (one + x);
304 return (one + x * (one + half * x));
305 }
306 t = x * x;
307 y = x + (t * (half + x * t2) +
308 (t * t) * (t3 + x * t4 + t * t5));
309 return (one + y);
310 }
311
312 /* find the multiple of 2^-6 nearest x */
313 k = hx >> 20;
314 j = (0x00100000 | (hx & 0x000fffff)) >> (0x40c - k);
315 j = (j - 1) & ~1;
316 if (ix < 0)
317 j += 134;
318 z = x - TBL2[j];
319 t = z * z;
320 y = z + (t * (half + z * t2) +
321 (t * t) * (t3 + z * t4 + t * t5));
322 return (TBL2[j+1] + TBL2[j+1] * y);
323 }
324
325 if (hx >= 0x40862e42) { /* x is large, infinite, or nan */
326 if (hx >= 0x7ff00000) {
327 if (ix == 0xfff00000 && ((int *)&x)[LOWORD] == 0)
328 return (zero);
329 return (x * x);
330 }
331 if (x > threshold1)
332 return (_SVID_libm_err(x, x, 6));
333 if (-x > threshold2)
334 return (_SVID_libm_err(x, x, 7));
335 }
336
337 t = invln2_32 * x;
338 if (ix < 0)
339 t -= half;
340 else
341 t += half;
342 k = (int)t;
343 j = (k & 0x1f) << 1;
344 m = k >> 5;
345 z = (x - k * ln2_32hi) - k * ln2_32lo;
346
347 /* z is now in primary range */
348 t = z * z;
349 y = z + (t * (half + z * t2) + (t * t) * (t3 + z * t4 + t * t5));
350 y = TBL[j] + (TBL[j+1] + TBL[j] * y);
351 if (m < -1021) {
352 ((int *)&y)[HIWORD] += (m + 54) << 20;
353 return (twom54 * y);
354 }
355 ((int *)&y)[HIWORD] += m << 20;
356 return (y);
357 }
|
1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
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 #pragma weak __exp = exp
32
33 /*
34 * exp(x)
35 * Hybrid algorithm of Peter Tang's Table driven method (for large
36 * arguments) and an accurate table (for small arguments).
37 * Written by K.C. Ng, November 1988.
38 * Method (large arguments):
39 * 1. Argument Reduction: given the input x, find r and integer k
40 * and j such that
41 * x = (k+j/32)*(ln2) + r, |r| <= (1/64)*ln2
42 *
43 * 2. exp(x) = 2^k * (2^(j/32) + 2^(j/32)*expm1(r))
44 * a. expm1(r) is approximated by a polynomial:
45 * expm1(r) ~ r + t1*r^2 + t2*r^3 + ... + t5*r^6
268 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
269 7.45133219101941108420e+02, /* 0x40874910, 0xD52D3051 */
270 5.55111512312578270212e-17, /* 0x3c900000, 0x00000000 */
271 };
272
273 #define half C[0]
274 #define invln2_32 C[1]
275 #define ln2_32hi C[2]
276 #define ln2_32lo C[3]
277 #define t2 C[4]
278 #define t3 C[5]
279 #define t4 C[6]
280 #define t5 C[7]
281 #define one C[8]
282 #define zero C[9]
283 #define threshold1 C[10]
284 #define threshold2 C[11]
285 #define twom54 C[12]
286
287 double
288 exp(double x)
289 {
290 double y, z, t;
291 int hx, ix, k, j, m;
292
293 ix = ((int *)&x)[HIWORD];
294 hx = ix & ~0x80000000;
295
296 if (hx < 0x3ff0a2b2) { /* |x| < 3/2 ln 2 */
297 if (hx < 0x3f862e42) { /* |x| < 1/64 ln 2 */
298 if (hx < 0x3ed00000) { /* |x| < 2^-18 */
299 volatile int dummy __unused;
300
301 dummy = (int)x; /* raise inexact if x != 0 */
302 #ifdef lint
303 dummy = dummy;
304 #endif
305
306 if (hx < 0x3e300000)
307 return (one + x);
308
309 return (one + x * (one + half * x));
310 }
311
312 t = x * x;
313 y = x + (t * (half + x * t2) + (t * t) * (t3 + x * t4 +
314 t * t5));
315 return (one + y);
316 }
317
318 /* find the multiple of 2^-6 nearest x */
319 k = hx >> 20;
320 j = (0x00100000 | (hx & 0x000fffff)) >> (0x40c - k);
321 j = (j - 1) & ~1;
322
323 if (ix < 0)
324 j += 134;
325
326 z = x - TBL2[j];
327 t = z * z;
328 y = z + (t * (half + z * t2) + (t * t) * (t3 + z * t4 + t *
329 t5));
330 return (TBL2[j + 1] + TBL2[j + 1] * y);
331 }
332
333 if (hx >= 0x40862e42) { /* x is large, infinite, or nan */
334 if (hx >= 0x7ff00000) {
335 if (ix == 0xfff00000 && ((int *)&x)[LOWORD] == 0)
336 return (zero);
337
338 return (x * x);
339 }
340
341 if (x > threshold1)
342 return (_SVID_libm_err(x, x, 6));
343
344 if (-x > threshold2)
345 return (_SVID_libm_err(x, x, 7));
346 }
347
348 t = invln2_32 * x;
349
350 if (ix < 0)
351 t -= half;
352 else
353 t += half;
354
355 k = (int)t;
356 j = (k & 0x1f) << 1;
357 m = k >> 5;
358 z = (x - k * ln2_32hi) - k * ln2_32lo;
359
360 /* z is now in primary range */
361 t = z * z;
362 y = z + (t * (half + z * t2) + (t * t) * (t3 + z * t4 + t * t5));
363 y = TBL[j] + (TBL[j + 1] + TBL[j] * y);
364
365 if (m < -1021) {
366 ((int *)&y)[HIWORD] += (m + 54) << 20;
367 return (twom54 * y);
368 }
369
370 ((int *)&y)[HIWORD] += m << 20;
371 return (y);
372 }
|