Print this page
11210 libm should be cstyle(1ONBLD) clean


   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 }