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 __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 }