Print this page
11210 libm should be cstyle(1ONBLD) clean
   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 }