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 2005 Sun Microsystems, Inc.  All rights reserved.
  26  * Use is subject to license terms.
  27  */
  28 
  29 #pragma weak __log = log
  30 
  31 /* INDENT OFF */
  32 /*
  33  * log(x)
  34  * Table look-up algorithm with product polynomial approximation.
  35  * By K.C. Ng, Oct 23, 2004. Updated Oct 18, 2005.
  36  *
  37  * (a). For x in [1-0.125, 1+0.1328125], using a special approximation:
  38  *      Let f = x - 1 and z = f*f.
  39  *      return f + ((a1*z) *
  40  *                 ((a2 + (a3*f)*(a4+f)) + (f*z)*(a5+f))) *
  41  *                 (((a6 + f*(a7+f)) + (f*z)*(a8+f)) *
  42  *                 ((a9 + (a10*f)*(a11+f)) + (f*z)*(a12+f)))
  43  * a1   -6.88821452420390473170286327331268694251775741577e-0002,
  44  * a2    1.97493380704769294631262255279580131173133850098e+0000,
  45  * a3    2.24963218866067560242072431719861924648284912109e+0000,
  46  * a4   -9.02975906958474405783476868236903101205825805664e-0001,
  47  * a5   -1.47391630715542865104339398385491222143173217773e+0000,
  48  * a6    1.86846544648220058704168877738993614912033081055e+0000,
  49  * a7    1.82277370459347465292410106485476717352867126465e+0000,
  50  * a8    1.25295479915214102994980294170090928673744201660e+0000,
  51  * a9    1.96709676945198275177517643896862864494323730469e+0000,


  78  *      with remez error |log(1+s) - P(s)| <= 2**-63.5
  79  *
  80  * (c). Otherwise, get "n", the exponent of x, and then normalize x to
  81  *      z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
  82  *      significant bits. Then
  83  *          log(x) = n*ln2 + log(Y[i]) + log(z/Y[i]).
  84  *
  85  * Special cases:
  86  *      log(x) is NaN with signal if x < 0 (including -INF) ;
  87  *      log(+INF) is +INF; log(0) is -INF with signal;
  88  *      log(NaN) is that NaN with no signal.
  89  *
  90  * Maximum error observed: less than 0.90 ulp
  91  *
  92  * Constants:
  93  * The hexadecimal values are the intended ones for the following constants.
  94  * The decimal values may be used, provided that the compiler will convert
  95  * from decimal to binary accurately enough to produce the hexadecimal values
  96  * shown.
  97  */
  98 /* INDENT ON */
  99 
 100 #include "libm.h"
 101 
 102 extern const double _TBL_log[];
 103 
 104 static const double P[] = {
 105 /* ONE   */  1.0,

 106 /* TWO52 */  4503599627370496.0,
 107 /* LN2HI */  6.93147180369123816490e-01,        /* 3fe62e42, fee00000 */
 108 /* LN2LO */  1.90821492927058770002e-10,        /* 3dea39ef, 35793c76 */
 109 /* A1    */ -6.88821452420390473170286327331268694251775741577e-0002,
 110 /* A2    */  1.97493380704769294631262255279580131173133850098e+0000,
 111 /* A3    */  2.24963218866067560242072431719861924648284912109e+0000,
 112 /* A4    */ -9.02975906958474405783476868236903101205825805664e-0001,
 113 /* A5    */ -1.47391630715542865104339398385491222143173217773e+0000,
 114 /* A6    */  1.86846544648220058704168877738993614912033081055e+0000,
 115 /* A7    */  1.82277370459347465292410106485476717352867126465e+0000,
 116 /* A8    */  1.25295479915214102994980294170090928673744201660e+0000,
 117 /* A9    */  1.96709676945198275177517643896862864494323730469e+0000,
 118 /* A10   */ -4.00127989749189894030934055990655906498432159424e-0001,
 119 /* A11   */  3.01675528558798333733648178167641162872314453125e+0000,
 120 /* A12   */ -9.52325445049240770778453679668018594384193420410e-0001,
 121 /* B1    */ -1.25041641589283658575482149899471551179885864258e-0001,
 122 /* B2    */  1.87161713283355151891381127914642725337613123482e+0000,
 123 /* B3    */ -1.89082956295731507978530316904652863740921020508e+0000,
 124 /* B4    */ -2.50562891673640253387134180229622870683670043945e+0000,
 125 /* B5    */  1.64822828085258366037635369139024987816810607910e+0000,


 137 #define A3    P[6]
 138 #define A4    P[7]
 139 #define A5    P[8]
 140 #define A6    P[9]
 141 #define A7    P[10]
 142 #define A8    P[11]
 143 #define A9    P[12]
 144 #define A10   P[13]
 145 #define A11   P[14]
 146 #define A12   P[15]
 147 #define B1    P[16]
 148 #define B2    P[17]
 149 #define B3    P[18]
 150 #define B4    P[19]
 151 #define B5    P[20]
 152 #define B6    P[21]
 153 #define B7    P[22]
 154 #define B8    P[23]
 155 
 156 double
 157 log(double x) {

 158         double  *tb, dn, dn1, s, z, r, w;
 159         int     i, hx, ix, n, lx;
 160 
 161         n = 0;
 162         hx = ((int *)&x)[HIWORD];
 163         ix = hx & 0x7fffffff;
 164         lx = ((int *)&x)[LOWORD];
 165 
 166         /* subnormal,0,negative,inf,nan */
 167         if ((hx + 0x100000) < 0x200000) {
 168                 if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0)) /* nan */
 169                         return (x * x);

 170                 if (((hx << 1) | lx) == 0)                /* zero */
 171                         return (_SVID_libm_err(x, x, 16));

 172                 if (hx < 0)                          /* negative */
 173                         return (_SVID_libm_err(x, x, 17));

 174                 if (((hx - 0x7ff00000) | lx) == 0)      /* +inf */
 175                         return (x);
 176 
 177                 /* x must be positive and subnormal */
 178                 x *= TWO52;
 179                 n = -52;
 180                 ix = ((int *)&x)[HIWORD];
 181                 lx = ((int *)&x)[LOWORD];
 182         }
 183 
 184         i = ix >> 19;

 185         if (i >= 0x7f7 && i <= 0x806) {
 186                 /* 0.09375 (0x3fb80000) <= x < 24 (0x40380000) */
 187                 if (ix >= 0x3fec0000 && ix < 0x3ff22000) {
 188                         /* 0.875 <= x < 1.125 */
 189                         s = x - ONE;
 190                         z = s * s;

 191                         if (((ix - 0x3ff00000) | lx) == 0) /* x = 1 */
 192                                 return (z);

 193                         r = (A10 * s) * (A11 + s);
 194                         w = z * s;
 195                         return (s + ((A1 * z) *
 196                                 (A2 + ((A3 * s) * (A4 + s) + w * (A5 + s)))) *
 197                                 ((A6 + (s * (A7 + s) + w * (A8 + s))) *
 198                                 (A9 + (r + w * (A12 + s)))));
 199                 } else {
 200                         i = (ix - 0x3fb80000) >> 15;
 201                         tb = (double *)_TBL_log + (i + i + i);
 202                         s = (x - tb[0]) * tb[1];
 203                         return (tb[2] +  ((B1 * s) * (B2 + s * (B3 + s))) *
 204                                 (((B4 + s * B5) + (s * s) * (B6 + s)) *
 205                                 (B7 + s * (B8 + s))));
 206                 }
 207         } else {
 208                 dn = (double)(n + ((ix >> 20) - 0x3ff));
 209                 dn1 = dn * LN2HI;
 210                 i = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
 211                 ((int *)&x)[HIWORD] = i;
 212                 i = (i - 0x3fb80000) >> 15;
 213                 tb = (double *)_TBL_log + (i + i + i);
 214                 s = (x - tb[0]) * tb[1];
 215                 dn = dn * LN2LO + tb[2];
 216                 return (dn1 + (dn + ((B1 * s) * (B2 + s * (B3 + s))) *
 217                         (((B4 + s * B5) + (s * s) * (B6 + s)) *
 218                         (B7 + s * (B8 + s)))));
   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 2005 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __log = log
  32 
  33 
  34 /*
  35  * log(x)
  36  * Table look-up algorithm with product polynomial approximation.
  37  * By K.C. Ng, Oct 23, 2004. Updated Oct 18, 2005.
  38  *
  39  * (a). For x in [1-0.125, 1+0.1328125], using a special approximation:
  40  *      Let f = x - 1 and z = f*f.
  41  *      return f + ((a1*z) *
  42  *                 ((a2 + (a3*f)*(a4+f)) + (f*z)*(a5+f))) *
  43  *                 (((a6 + f*(a7+f)) + (f*z)*(a8+f)) *
  44  *                 ((a9 + (a10*f)*(a11+f)) + (f*z)*(a12+f)))
  45  * a1   -6.88821452420390473170286327331268694251775741577e-0002,
  46  * a2    1.97493380704769294631262255279580131173133850098e+0000,
  47  * a3    2.24963218866067560242072431719861924648284912109e+0000,
  48  * a4   -9.02975906958474405783476868236903101205825805664e-0001,
  49  * a5   -1.47391630715542865104339398385491222143173217773e+0000,
  50  * a6    1.86846544648220058704168877738993614912033081055e+0000,
  51  * a7    1.82277370459347465292410106485476717352867126465e+0000,
  52  * a8    1.25295479915214102994980294170090928673744201660e+0000,
  53  * a9    1.96709676945198275177517643896862864494323730469e+0000,


  80  *      with remez error |log(1+s) - P(s)| <= 2**-63.5
  81  *
  82  * (c). Otherwise, get "n", the exponent of x, and then normalize x to
  83  *      z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
  84  *      significant bits. Then
  85  *          log(x) = n*ln2 + log(Y[i]) + log(z/Y[i]).
  86  *
  87  * Special cases:
  88  *      log(x) is NaN with signal if x < 0 (including -INF) ;
  89  *      log(+INF) is +INF; log(0) is -INF with signal;
  90  *      log(NaN) is that NaN with no signal.
  91  *
  92  * Maximum error observed: less than 0.90 ulp
  93  *
  94  * Constants:
  95  * The hexadecimal values are the intended ones for the following constants.
  96  * The decimal values may be used, provided that the compiler will convert
  97  * from decimal to binary accurately enough to produce the hexadecimal values
  98  * shown.
  99  */

 100 
 101 #include "libm.h"
 102 
 103 extern const double _TBL_log[];
 104 
 105 static const double P[] = {
 106 /* ONE   */
 107         1.0,
 108 /* TWO52 */ 4503599627370496.0,
 109 /* LN2HI */ 6.93147180369123816490e-01, /* 3fe62e42, fee00000 */
 110 /* LN2LO */ 1.90821492927058770002e-10, /* 3dea39ef, 35793c76 */
 111 /* A1    */ -6.88821452420390473170286327331268694251775741577e-0002,
 112 /* A2    */ 1.97493380704769294631262255279580131173133850098e+0000,
 113 /* A3    */ 2.24963218866067560242072431719861924648284912109e+0000,
 114 /* A4    */ -9.02975906958474405783476868236903101205825805664e-0001,
 115 /* A5    */ -1.47391630715542865104339398385491222143173217773e+0000,
 116 /* A6    */ 1.86846544648220058704168877738993614912033081055e+0000,
 117 /* A7    */ 1.82277370459347465292410106485476717352867126465e+0000,
 118 /* A8    */ 1.25295479915214102994980294170090928673744201660e+0000,
 119 /* A9    */ 1.96709676945198275177517643896862864494323730469e+0000,
 120 /* A10   */ -4.00127989749189894030934055990655906498432159424e-0001,
 121 /* A11   */ 3.01675528558798333733648178167641162872314453125e+0000,
 122 /* A12   */ -9.52325445049240770778453679668018594384193420410e-0001,
 123 /* B1    */ -1.25041641589283658575482149899471551179885864258e-0001,
 124 /* B2    */ 1.87161713283355151891381127914642725337613123482e+0000,
 125 /* B3    */ -1.89082956295731507978530316904652863740921020508e+0000,
 126 /* B4    */ -2.50562891673640253387134180229622870683670043945e+0000,
 127 /* B5    */ 1.64822828085258366037635369139024987816810607910e+0000,


 139 #define A3              P[6]
 140 #define A4              P[7]
 141 #define A5              P[8]
 142 #define A6              P[9]
 143 #define A7              P[10]
 144 #define A8              P[11]
 145 #define A9              P[12]
 146 #define A10             P[13]
 147 #define A11             P[14]
 148 #define A12             P[15]
 149 #define B1              P[16]
 150 #define B2              P[17]
 151 #define B3              P[18]
 152 #define B4              P[19]
 153 #define B5              P[20]
 154 #define B6              P[21]
 155 #define B7              P[22]
 156 #define B8              P[23]
 157 
 158 double
 159 log(double x)
 160 {
 161         double *tb, dn, dn1, s, z, r, w;
 162         int i, hx, ix, n, lx;
 163 
 164         n = 0;
 165         hx = ((int *)&x)[HIWORD];
 166         ix = hx & 0x7fffffff;
 167         lx = ((int *)&x)[LOWORD];
 168 
 169         /* subnormal,0,negative,inf,nan */
 170         if ((hx + 0x100000) < 0x200000) {
 171                 if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0)) /* nan */
 172                         return (x * x);
 173 
 174                 if (((hx << 1) | lx) == 0) /* zero */
 175                         return (_SVID_libm_err(x, x, 16));
 176 
 177                 if (hx < 0)  /* negative */
 178                         return (_SVID_libm_err(x, x, 17));
 179 
 180                 if (((hx - 0x7ff00000) | lx) == 0) /* +inf */
 181                         return (x);
 182 
 183                 /* x must be positive and subnormal */
 184                 x *= TWO52;
 185                 n = -52;
 186                 ix = ((int *)&x)[HIWORD];
 187                 lx = ((int *)&x)[LOWORD];
 188         }
 189 
 190         i = ix >> 19;
 191 
 192         if (i >= 0x7f7 && i <= 0x806) {
 193                 /* 0.09375 (0x3fb80000) <= x < 24 (0x40380000) */
 194                 if (ix >= 0x3fec0000 && ix < 0x3ff22000) {
 195                         /* 0.875 <= x < 1.125 */
 196                         s = x - ONE;
 197                         z = s * s;
 198 
 199                         if (((ix - 0x3ff00000) | lx) == 0)      /* x = 1 */
 200                                 return (z);
 201 
 202                         r = (A10 * s) * (A11 + s);
 203                         w = z * s;
 204                         return (s + ((A1 * z) * (A2 + ((A3 * s) * (A4 + s) + w *
 205                             (A5 + s)))) * ((A6 + (s * (A7 + s) + w *
 206                             (A8 + s))) * (A9 + (r + w * (A12 + s)))));

 207                 } else {
 208                         i = (ix - 0x3fb80000) >> 15;
 209                         tb = (double *)_TBL_log + (i + i + i);
 210                         s = (x - tb[0]) * tb[1];
 211                         return (tb[2] + ((B1 * s) * (B2 + s * (B3 + s))) *
 212                             (((B4 + s * B5) + (s * s) * (B6 + s)) *
 213                             (B7 + s * (B8 + s))));
 214                 }
 215         } else {
 216                 dn = (double)(n + ((ix >> 20) - 0x3ff));
 217                 dn1 = dn * LN2HI;
 218                 i = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
 219                 ((int *)&x)[HIWORD] = i;
 220                 i = (i - 0x3fb80000) >> 15;
 221                 tb = (double *)_TBL_log + (i + i + i);
 222                 s = (x - tb[0]) * tb[1];
 223                 dn = dn * LN2LO + tb[2];
 224                 return (dn1 + (dn + ((B1 * s) * (B2 + s * (B3 + s))) *
 225                     (((B4 + s * B5) + (s * s) * (B6 + s)) *
 226                     (B7 + s * (B8 + s)))));