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 __log10 = log10
  30 
  31 /* INDENT OFF */
  32 /*
  33  * log10(x) = log(x)/log10
  34  *
  35  * Base on Table look-up algorithm with product polynomial
  36  * approximation for log(x).
  37  *
  38  * By K.C. Ng, Nov 29, 2004
  39  *
  40  * (a). For x in [1-0.125, 1+0.125], from log.c we have
  41  *      log(x) =  f + ((a1*f^2) *
  42  *                 ((a2 + (a3*f)*(a4+f)) + (f^3)*(a5+f))) *
  43  *                 (((a6 + f*(a7+f)) + (f^3)*(a8+f)) *
  44  *                 ((a9 + (a10*f)*(a11+f)) + (f^3)*(a12+f)))
  45  *      where f = x - 1.
  46  *      (i) modify a1 <- a1 / log10
  47  *      (ii) 1/log10 = 0.4342944819...
  48  *                   = 0.4375 - 0.003205518... (7 bit shift)
  49  *           Let lgv = 0.4375 - 1/log10, then
  50  *           lgv = 0.003205518096748172348871081083395...,
  51  *      (iii) f*0.4375 is exact because f has 3 trailing zero.


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

  96 /* TWO52 */  4503599627370496.0,
  97 /* LNAHI */  3.01029995607677847147e-01,        /* 3FD34413 50900000 */
  98 /* LNALO */  5.63033480667509769841e-11,        /* 3DCEF3FD E623E256 */
  99 /* A1    */ -2.9142521960136582507385480707044582802184e-02,
 100 /* A2    */  1.99628461483039965074226529395673424005508422852e+0000,
 101 /* A3    */  2.26812367662950720159642514772713184356689453125e+0000,
 102 /* A4    */ -9.05030639084976384900471657601883634924888610840e-0001,
 103 /* A5    */ -1.48275767132434044270894446526654064655303955078e+0000,
 104 /* A6    */  1.88158320939722756293122074566781520843505859375e+0000,
 105 /* A7    */  1.83309386046986411145098827546462416648864746094e+0000,
 106 /* A8    */  1.24847063988317086291601754055591300129890441895e+0000,
 107 /* A9    */  1.98372421445537705508854742220137268304824829102e+0000,
 108 /* A10   */ -3.94711735767898475035764249696512706577777862549e-0001,
 109 /* A11   */  3.07890395362954372160402272129431366920471191406e+0000,
 110 /* A12   */ -9.60099585275022149311041630426188930869102478027e-0001,
 111 /* B1    */ -5.4304894950350052960838096752491540286689e-02,
 112 /* B2    */  1.87161713283355151891381127914642725337613123482e+0000,
 113 /* B3    */ -1.89082956295731507978530316904652863740921020508e+0000,
 114 /* B4    */ -2.50562891673640253387134180229622870683670043945e+0000,
 115 /* B5    */  1.64822828085258366037635369139024987816810607910e+0000,


 133 #define A6    P[9]
 134 #define A7    P[10]
 135 #define A8    P[11]
 136 #define A9    P[12]
 137 #define A10   P[13]
 138 #define A11   P[14]
 139 #define A12   P[15]
 140 #define B1    P[16]
 141 #define B2    P[17]
 142 #define B3    P[18]
 143 #define B4    P[19]
 144 #define B5    P[20]
 145 #define B6    P[21]
 146 #define B7    P[22]
 147 #define B8    P[23]
 148 #define LGH   P[24]
 149 #define LGL   P[25]
 150 #define LG10V P[26]
 151 
 152 double
 153 log10(double x) {

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

 166                 if (((hx << 1) | lx) == 0)                /* zero */
 167                         return (_SVID_libm_err(x, x, 18));

 168                 if (hx < 0)                          /* negative */
 169                         return (_SVID_libm_err(x, x, 19));

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

 181         if (i >= 0x7f7 && i <= 0x806) {
 182                 /* 0.09375 (0x3fb80000) <= x < 24 (0x40380000) */
 183                 if (ix >= 0x3fec0000 && ix < 0x3ff20000) {
 184                         /* 0.875 <= x < 1.125 */
 185                         s = x - ONE;
 186                         z = s * s;

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

 189                         r = (A10 * s) * (A11 + s);
 190                         w = z * s;
 191                         return (LGH * s - (LGL * s - ((A1 * z) *
 192                                 ((A2 + (A3 * s) * (A4 + s)) + w * (A5 + s))) *
 193                                 (((A6 + s * (A7 + s)) + w * (A8 + s)) *
 194                                 ((A9 + r) + w * (A12 + s)))));
 195                 } else {
 196                         i = (ix - 0x3fb80000) >> 15;
 197                         tb = (double *)_TBL_log + (i + i + i);
 198                         s = (x - tb[0]) * tb[1];
 199                         return (LGH * tb[2] - (LGL * tb[2] - ((B1 * s) *
 200                                 (B2 + s * (B3 + s))) *
 201                                 (((B4 + s * B5) + (s * s) * (B6 + s)) *
 202                                 (B7 + s * (B8 + s)))));
 203                 }
 204         } else {
 205                 dn = (double)(n + ((ix >> 20) - 0x3ff));
 206                 dn1 = dn * LNAHI;
 207                 i = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
 208                 ((int *)&x)[HIWORD] = i;
 209                 i = (i - 0x3fb80000) >> 15;
 210                 tb = (double *)_TBL_log + (i + i + i);
 211                 s = (x - tb[0]) * tb[1];
 212                 dn = dn * LNALO + tb[2] * LG10V;
 213                 return (dn1 + (dn + ((B1 * s) *
 214                         (B2 + s * (B3 + s))) *
 215                         (((B4 + s * B5) + (s * s) * (B6 + s)) *
 216                         (B7 + s * (B8 + s)))));
 217         }
 218 }
   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 __log10 = log10
  32 
  33 
  34 /*
  35  * log10(x) = log(x)/log10
  36  *
  37  * Base on Table look-up algorithm with product polynomial
  38  * approximation for log(x).
  39  *
  40  * By K.C. Ng, Nov 29, 2004
  41  *
  42  * (a). For x in [1-0.125, 1+0.125], from log.c we have
  43  *      log(x) =  f + ((a1*f^2) *
  44  *                 ((a2 + (a3*f)*(a4+f)) + (f^3)*(a5+f))) *
  45  *                 (((a6 + f*(a7+f)) + (f^3)*(a8+f)) *
  46  *                 ((a9 + (a10*f)*(a11+f)) + (f^3)*(a12+f)))
  47  *      where f = x - 1.
  48  *      (i) modify a1 <- a1 / log10
  49  *      (ii) 1/log10 = 0.4342944819...
  50  *                   = 0.4375 - 0.003205518... (7 bit shift)
  51  *           Let lgv = 0.4375 - 1/log10, then
  52  *           lgv = 0.003205518096748172348871081083395...,
  53  *      (iii) f*0.4375 is exact because f has 3 trailing zero.


  70  *
  71  * (c). Otherwise, get "n", the exponent of x, and then normalize x to
  72  *      z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
  73  *      significant bits. Then
  74  *          log(x) = n*ln2 + log(Y[i]) + log(z/Y[i]).
  75  *          log10(x) = n*(ln2/ln10) + log10(z).
  76  *
  77  * Special cases:
  78  *      log10(x) is NaN with signal if x < 0 (including -INF) ;
  79  *      log10(+INF) is +INF; log10(0) is -INF with signal;
  80  *      log10(NaN) is that NaN with no signal.
  81  *
  82  * Maximum error observed: less than 0.89 ulp
  83  *
  84  * Constants:
  85  * The hexadecimal values are the intended ones for the following constants.
  86  * The decimal values may be used, provided that the compiler will convert
  87  * from decimal to binary accurately enough to produce the hexadecimal values
  88  * shown.
  89  */

  90 
  91 #include "libm.h"
  92 
  93 extern const double _TBL_log[];
  94 
  95 static const double P[] = {
  96 /* ONE   */
  97         1.0,
  98 /* TWO52 */ 4503599627370496.0,
  99 /* LNAHI */ 3.01029995607677847147e-01, /* 3FD34413 50900000 */
 100 /* LNALO */ 5.63033480667509769841e-11, /* 3DCEF3FD E623E256 */
 101 /* A1    */ -2.9142521960136582507385480707044582802184e-02,
 102 /* A2    */ 1.99628461483039965074226529395673424005508422852e+0000,
 103 /* A3    */ 2.26812367662950720159642514772713184356689453125e+0000,
 104 /* A4    */ -9.05030639084976384900471657601883634924888610840e-0001,
 105 /* A5    */ -1.48275767132434044270894446526654064655303955078e+0000,
 106 /* A6    */ 1.88158320939722756293122074566781520843505859375e+0000,
 107 /* A7    */ 1.83309386046986411145098827546462416648864746094e+0000,
 108 /* A8    */ 1.24847063988317086291601754055591300129890441895e+0000,
 109 /* A9    */ 1.98372421445537705508854742220137268304824829102e+0000,
 110 /* A10   */ -3.94711735767898475035764249696512706577777862549e-0001,
 111 /* A11   */ 3.07890395362954372160402272129431366920471191406e+0000,
 112 /* A12   */ -9.60099585275022149311041630426188930869102478027e-0001,
 113 /* B1    */ -5.4304894950350052960838096752491540286689e-02,
 114 /* B2    */ 1.87161713283355151891381127914642725337613123482e+0000,
 115 /* B3    */ -1.89082956295731507978530316904652863740921020508e+0000,
 116 /* B4    */ -2.50562891673640253387134180229622870683670043945e+0000,
 117 /* B5    */ 1.64822828085258366037635369139024987816810607910e+0000,


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

 211                 }
 212         } else {
 213                 dn = (double)(n + ((ix >> 20) - 0x3ff));
 214                 dn1 = dn * LNAHI;
 215                 i = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
 216                 ((int *)&x)[HIWORD] = i;
 217                 i = (i - 0x3fb80000) >> 15;
 218                 tb = (double *)_TBL_log + (i + i + i);
 219                 s = (x - tb[0]) * tb[1];
 220                 dn = dn * LNALO + tb[2] * LG10V;
 221                 return (dn1 + (dn + ((B1 * s) * (B2 + s * (B3 + s))) *

 222                     (((B4 + s * B5) + (s * s) * (B6 + s)) *
 223                     (B7 + s * (B8 + s)))));
 224         }
 225 }