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 __expf = expf
30
31 /* INDENT OFF */
32 /*
33 * float expf(float x);
34 * Code by K.C. Ng for SUN 5.0 libmopt
35 * 11/5/99
36 * Method :
37 * 1. For |x| >= 2^7, either underflow/overflow.
38 * More precisely:
39 * x > 88.722839355...(0x42B17218) => overflow;
40 * x < -103.97207642..(0xc2CFF1B4) => underflow.
41 * 2. For |x| < 2^-6, use polynomail
42 * exp(x) = 1 + x + p1*x^2 + p2*x^3
43 * 3. Otherwise, write |x|=(1+r)*2^n, where 0<=r<1.
44 * Let t = 2^n * (1+r) .... x > 0;
45 * t = 2^n * (1-r) .... x < 0. (x= -2**(n+1)+t)
46 * Since -6 <= n <= 6, we may break t into
47 * six 6-bits chunks:
48 * -5 -11 -17 -23 -29
49 * t=j *2+j *2 +j *2 +j *2 +j *2 +j *2
50 * 1 2 3 4 5 6
51 *
83 * 5. If x < 0, return exp(-2 )* exp(t). Note that
84 * -6 <= n <= 6. Let k = n - 6, then we can
85 * precompute
86 * k-5 n+1
87 * EN[k] = exp(-2 ) = exp(-2 ) for k=0,1,...,12.
88 *
89 *
90 * Special cases:
91 * exp(INF) is INF, exp(NaN) is NaN;
92 * exp(-INF) = 0;
93 * for finite argument, only exp(0) = 1 is exact.
94 *
95 * Accuracy:
96 * All calculations are done in double precision except for
97 * the case |x| < 2^-6. When |x| < 2^-6, the error is less
98 * than 0.55 ulp. When |x| >= 2^-6 and the result is normal,
99 * the error is less than 0.51 ulp. When FDTOS_TRAPS_... is
100 * defined and the result is subnormal, the error can be as
101 * large as 0.75 ulp.
102 */
103 /* INDENT ON */
104
105 #include "libm.h"
106
107 /*
108 * ET[k] = exp(j*2^(7-6i)) , where j = k mod 64, i = k/64
109 */
110 static const double ET[] = {
111 1.00000000000000000000e+00, 1.00000000186264514923e+00,
112 1.00000000372529029846e+00, 1.00000000558793544769e+00,
113 1.00000000745058059692e+00, 1.00000000931322574615e+00,
114 1.00000001117587089539e+00, 1.00000001303851604462e+00,
115 1.00000001490116119385e+00, 1.00000001676380656512e+00,
116 1.00000001862645171435e+00, 1.00000002048909686359e+00,
117 1.00000002235174201282e+00, 1.00000002421438716205e+00,
118 1.00000002607703253332e+00, 1.00000002793967768255e+00,
119 1.00000002980232283178e+00, 1.00000003166496798102e+00,
120 1.00000003352761335229e+00, 1.00000003539025850152e+00,
121 1.00000003725290365075e+00, 1.00000003911554879998e+00,
122 1.00000004097819417126e+00, 1.00000004284083932049e+00,
123 1.00000004470348446972e+00, 1.00000004656612984100e+00,
325 5.0000000951292138e-01F,
326 1.6666518897347284e-01F,
327 3.4028234663852885981170E+38F,
328 1.1754943508222875079688E-38F,
329 #if defined(FDTOS_TRAPS_INCOMPLETE_IN_FNS_MODE)
330 8.67361737988403547205962240695953369140625e-19F
331 #endif
332 };
333
334 #define zero F[0]
335 #define one F[1]
336 #define p1 F[2]
337 #define p2 F[3]
338 #define big F[4]
339 #define tiny F[5]
340 #if defined(FDTOS_TRAPS_INCOMPLETE_IN_FNS_MODE)
341 #define twom60 F[6]
342 #endif
343
344 float
345 expf(float xf) {
346 double w, p, q;
347 int hx, ix, n;
348
349 hx = *(int *)&xf;
350 ix = hx & ~0x80000000;
351
352 if (ix < 0x3c800000) { /* |x| < 2**-6 */
353 if (ix < 0x38800000) /* |x| < 2**-14 */
354 return (one + xf);
355 return (one + (xf + (xf * xf) * (p1 + xf * p2)));
356 }
357
358 n = ix >> 23; /* biased exponent */
359
360 if (n >= 0x86) { /* |x| >= 2^7 */
361 if (n >= 0xff) { /* x is nan of +-inf */
362 if (hx == 0xff800000)
363 return (zero); /* exp(-inf)=0 */
364 return (xf * xf); /* exp(nan/inf) is nan or inf */
365 }
366 if (hx > 0)
367 return (big * big); /* overflow */
368 else
369 return (tiny * tiny); /* underflow */
370 }
371
372 ix -= n << 23;
373 if (hx > 0)
374 ix += 0x800000;
375 else
376 ix = 0x800000 - ix;
377 if (n >= 0x7f) { /* n >= 0 */
378 ix <<= n - 0x7f;
379 w = ET[(ix & 0x3f) + 64] * ET[((ix >> 6) & 0x3f) + 128];
380 p = ET[((ix >> 12) & 0x3f) + 192] *
381 ET[((ix >> 18) & 0x3f) + 256];
382 q = ET[((ix >> 24) & 0x3f) + 320];
383 } else {
384 ix <<= n - 0x79;
385 w = ET[ix & 0x3f] * ET[((ix >> 6) & 0x3f) + 64];
386 p = ET[((ix >> 12) & 0x3f) + 128] *
387 ET[((ix >> 18) & 0x3f) + 192];
388 q = ET[((ix >> 24) & 0x3f) + 256];
389 }
390 xf = (float)((w * p) * (hx < 0 ? q * EN[n - 0x79] : q));
391 #if defined(FDTOS_TRAPS_INCOMPLETE_IN_FNS_MODE)
392 if ((unsigned)hx >= 0xc2800000u) {
393 if ((unsigned)hx >= 0xc2aeac50) { /* force underflow */
394 volatile float t = tiny;
395 t *= t;
396 }
397 return (xf * twom60);
398 }
399 #endif
400 return (xf);
401 }
|
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 __expf = expf
32
33
34 /*
35 * float expf(float x);
36 * Code by K.C. Ng for SUN 5.0 libmopt
37 * 11/5/99
38 * Method :
39 * 1. For |x| >= 2^7, either underflow/overflow.
40 * More precisely:
41 * x > 88.722839355...(0x42B17218) => overflow;
42 * x < -103.97207642..(0xc2CFF1B4) => underflow.
43 * 2. For |x| < 2^-6, use polynomail
44 * exp(x) = 1 + x + p1*x^2 + p2*x^3
45 * 3. Otherwise, write |x|=(1+r)*2^n, where 0<=r<1.
46 * Let t = 2^n * (1+r) .... x > 0;
47 * t = 2^n * (1-r) .... x < 0. (x= -2**(n+1)+t)
48 * Since -6 <= n <= 6, we may break t into
49 * six 6-bits chunks:
50 * -5 -11 -17 -23 -29
51 * t=j *2+j *2 +j *2 +j *2 +j *2 +j *2
52 * 1 2 3 4 5 6
53 *
85 * 5. If x < 0, return exp(-2 )* exp(t). Note that
86 * -6 <= n <= 6. Let k = n - 6, then we can
87 * precompute
88 * k-5 n+1
89 * EN[k] = exp(-2 ) = exp(-2 ) for k=0,1,...,12.
90 *
91 *
92 * Special cases:
93 * exp(INF) is INF, exp(NaN) is NaN;
94 * exp(-INF) = 0;
95 * for finite argument, only exp(0) = 1 is exact.
96 *
97 * Accuracy:
98 * All calculations are done in double precision except for
99 * the case |x| < 2^-6. When |x| < 2^-6, the error is less
100 * than 0.55 ulp. When |x| >= 2^-6 and the result is normal,
101 * the error is less than 0.51 ulp. When FDTOS_TRAPS_... is
102 * defined and the result is subnormal, the error can be as
103 * large as 0.75 ulp.
104 */
105
106 #include "libm.h"
107
108 /*
109 * ET[k] = exp(j*2^(7-6i)) , where j = k mod 64, i = k/64
110 */
111 static const double ET[] = {
112 1.00000000000000000000e+00, 1.00000000186264514923e+00,
113 1.00000000372529029846e+00, 1.00000000558793544769e+00,
114 1.00000000745058059692e+00, 1.00000000931322574615e+00,
115 1.00000001117587089539e+00, 1.00000001303851604462e+00,
116 1.00000001490116119385e+00, 1.00000001676380656512e+00,
117 1.00000001862645171435e+00, 1.00000002048909686359e+00,
118 1.00000002235174201282e+00, 1.00000002421438716205e+00,
119 1.00000002607703253332e+00, 1.00000002793967768255e+00,
120 1.00000002980232283178e+00, 1.00000003166496798102e+00,
121 1.00000003352761335229e+00, 1.00000003539025850152e+00,
122 1.00000003725290365075e+00, 1.00000003911554879998e+00,
123 1.00000004097819417126e+00, 1.00000004284083932049e+00,
124 1.00000004470348446972e+00, 1.00000004656612984100e+00,
326 5.0000000951292138e-01F,
327 1.6666518897347284e-01F,
328 3.4028234663852885981170E+38F,
329 1.1754943508222875079688E-38F,
330 #if defined(FDTOS_TRAPS_INCOMPLETE_IN_FNS_MODE)
331 8.67361737988403547205962240695953369140625e-19F
332 #endif
333 };
334
335 #define zero F[0]
336 #define one F[1]
337 #define p1 F[2]
338 #define p2 F[3]
339 #define big F[4]
340 #define tiny F[5]
341 #if defined(FDTOS_TRAPS_INCOMPLETE_IN_FNS_MODE)
342 #define twom60 F[6]
343 #endif
344
345 float
346 expf(float xf)
347 {
348 double w, p, q;
349 int hx, ix, n;
350
351 hx = *(int *)&xf;
352 ix = hx & ~0x80000000;
353
354 if (ix < 0x3c800000) { /* |x| < 2**-6 */
355 if (ix < 0x38800000) /* |x| < 2**-14 */
356 return (one + xf);
357
358 return (one + (xf + (xf * xf) * (p1 + xf * p2)));
359 }
360
361 n = ix >> 23; /* biased exponent */
362
363 if (n >= 0x86) { /* |x| >= 2^7 */
364 if (n >= 0xff) { /* x is nan of +-inf */
365 if (hx == 0xff800000)
366 return (zero); /* exp(-inf)=0 */
367
368 return (xf * xf); /* exp(nan/inf) is nan or inf */
369 }
370
371 if (hx > 0)
372 return (big * big); /* overflow */
373 else
374 return (tiny * tiny); /* underflow */
375 }
376
377 ix -= n << 23;
378
379 if (hx > 0)
380 ix += 0x800000;
381 else
382 ix = 0x800000 - ix;
383
384 if (n >= 0x7f) { /* n >= 0 */
385 ix <<= n - 0x7f;
386 w = ET[(ix & 0x3f) + 64] * ET[((ix >> 6) & 0x3f) + 128];
387 p = ET[((ix >> 12) & 0x3f) + 192] * ET[((ix >> 18) & 0x3f) +
388 256];
389 q = ET[((ix >> 24) & 0x3f) + 320];
390 } else {
391 ix <<= n - 0x79;
392 w = ET[ix & 0x3f] * ET[((ix >> 6) & 0x3f) + 64];
393 p = ET[((ix >> 12) & 0x3f) + 128] * ET[((ix >> 18) & 0x3f) +
394 192];
395 q = ET[((ix >> 24) & 0x3f) + 256];
396 }
397
398 xf = (float)((w * p) * (hx < 0 ? q * EN[n - 0x79] : q));
399 #if defined(FDTOS_TRAPS_INCOMPLETE_IN_FNS_MODE)
400 if ((unsigned)hx >= 0xc2800000u) {
401 if ((unsigned)hx >= 0xc2aeac50) { /* force underflow */
402 volatile float t = tiny;
403
404 t *= t;
405 }
406
407 return (xf * twom60);
408 }
409 #endif
410 return (xf);
411 }
|