lib/libm/log10.c

   1 /*                                                      log10.c
   2  *
   3  *      Common logarithm
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, log10();
  10  *
  11  * y = log10( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns logarithm to the base 10 of x.
  18  *
  19  * The argument is separated into its exponent and fractional
  20  * parts.  The logarithm of the fraction is approximated by
  21  *
  22  *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
  23  *
  24  *
  25  *
  26  * ACCURACY:
  27  *
  28  *                      Relative error:
  29  * arithmetic   domain     # trials      peak         rms
  30  *    IEEE      0.5, 2.0     30000      1.5e-16     5.0e-17
  31  *    IEEE      0, MAXNUM    30000      1.4e-16     4.8e-17
  32  *    DEC       1, MAXNUM    50000      2.5e-17     6.0e-18
  33  *
  34  * In the tests over the interval [1, MAXNUM], the logarithms
  35  * of the random arguments were uniformly distributed over
  36  * [0, MAXLOG].
  37  *
  38  * ERROR MESSAGES:
  39  *
  40  * log10 singularity:  x = 0; returns -INFINITY
  41  * log10 domain:       x < 0; returns NAN
  42  */
  43 \f
  44 /*
  45 Cephes Math Library Release 2.8:  June, 2000
  46 Copyright 1984, 1995, 2000 by Stephen L. Moshier
  47 */
  48
  49 #include "mconf.h"
  50 static char fname[] = {"log10"};
  51
  52 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  53  * 1/sqrt(2) <= x < sqrt(2)
  54  */
  55 #ifdef UNK
  56 static double P[] = {
  57   4.58482948458143443514E-5,
  58   4.98531067254050724270E-1,
  59   6.56312093769992875930E0,
  60   2.97877425097986925891E1,
  61   6.06127134467767258030E1,
  62   5.67349287391754285487E1,
  63   1.98892446572874072159E1
  64 };
  65 static double Q[] = {
  66 /* 1.00000000000000000000E0, */
  67   1.50314182634250003249E1,
  68   8.27410449222435217021E1,
  69   2.20664384982121929218E2,
  70   3.07254189979530058263E2,
  71   2.14955586696422947765E2,
  72   5.96677339718622216300E1
  73 };
  74 #endif
  75
  76 #ifdef DEC
  77 static unsigned short P[] = {
  78 0034500,0046473,0051374,0135174,
  79 0037777,0037566,0145712,0150321,
  80 0040722,0002426,0031543,0123107,
  81 0041356,0046513,0170752,0004346,
  82 0041562,0071553,0023536,0163343,
  83 0041542,0170221,0024316,0114216,
  84 0041237,0016454,0046611,0104602
  85 };
  86 static unsigned short Q[] = {
  87 /*0040200,0000000,0000000,0000000,*/
  88 0041160,0100260,0067736,0102424,
  89 0041645,0075552,0036563,0147072,
  90 0042134,0125025,0021132,0025320,
  91 0042231,0120211,0046030,0103271,
  92 0042126,0172241,0052151,0120426,
  93 0041556,0125702,0072116,0047103
  94 };
  95 #endif
  96
  97 #ifdef IBMPC
  98 static unsigned short P[] = {
  99 0x974f,0x6a5f,0x09a7,0x3f08,
 100 0x5a1a,0xd979,0xe7ee,0x3fdf,
 101 0x74c9,0xc66c,0x40a2,0x401a,
 102 0x411d,0x7e3d,0xc9a9,0x403d,
 103 0xdcdc,0x64eb,0x4e6d,0x404e,
 104 0xd312,0x2519,0x5e12,0x404c,
 105 0x3130,0x89b1,0xe3a5,0x4033
 106 };
 107 static unsigned short Q[] = {
 108 /*0x0000,0x0000,0x0000,0x3ff0,*/
 109 0xd0a2,0x0dfb,0x1016,0x402e,
 110 0x79c7,0x47ae,0xaf6d,0x4054,
 111 0x455a,0xa44b,0x9542,0x406b,
 112 0x10d7,0x2983,0x3411,0x4073,
 113 0x3423,0x2a8d,0xde94,0x406a,
 114 0xc9c8,0x4e89,0xd578,0x404d
 115 };
 116 #endif
 117
 118 #ifdef MIEEE
 119 static unsigned short P[] = {
 120 0x3f08,0x09a7,0x6a5f,0x974f,
 121 0x3fdf,0xe7ee,0xd979,0x5a1a,
 122 0x401a,0x40a2,0xc66c,0x74c9,
 123 0x403d,0xc9a9,0x7e3d,0x411d,
 124 0x404e,0x4e6d,0x64eb,0xdcdc,
 125 0x404c,0x5e12,0x2519,0xd312,
 126 0x4033,0xe3a5,0x89b1,0x3130
 127 };
 128 static unsigned short Q[] = {
 129 0x402e,0x1016,0x0dfb,0xd0a2,
 130 0x4054,0xaf6d,0x47ae,0x79c7,
 131 0x406b,0x9542,0xa44b,0x455a,
 132 0x4073,0x3411,0x2983,0x10d7,
 133 0x406a,0xde94,0x2a8d,0x3423,
 134 0x404d,0xd578,0x4e89,0xc9c8
 135 };
 136 #endif
 137
 138 #define SQRTH 0.70710678118654752440
 139 #define L102A 3.0078125E-1
 140 #define L102B 2.48745663981195213739E-4
 141 #define L10EA 4.3359375E-1
 142 #define L10EB 7.00731903251827651129E-4
 143
 144 #ifdef ANSIPROT
 145 extern double frexp ( double, int * );
 146 extern double ldexp ( double, int );
 147 extern double polevl ( double, void *, int );
 148 extern double p1evl ( double, void *, int );
 149 extern int isnan ( double );
 150 extern int isfinite ( double );
 151 #else
 152 double frexp(), ldexp(), polevl(), p1evl();
 153 int isnan(), isfinite();
 154 #endif
 155 extern double LOGE2, SQRT2, INFINITY, NAN;
 156
 157 double log10(x)
 158 double x;
 159 {
 160 VOLATILE double z;
 161 double y;
 162 #ifdef DEC
 163 short *q;
 164 #endif
 165 int e;
 166
 167 #ifdef NANS
 168 if( isnan(x) )
 169         return(x);
 170 #endif
 171 #ifdef INFINITIES
 172 if( x == INFINITY )
 173         return(x);
 174 #endif
 175 /* Test for domain */
 176 if( x <= 0.0 )
 177         {
 178         if( x == 0.0 )
 179                 {
 180                 mtherr( fname, SING );
 181                 return( -INFINITY );
 182                 }
 183         else
 184                 {
 185                 mtherr( fname, DOMAIN );
 186                 return( NAN );
 187                 }
 188         }
 189
 190 /* separate mantissa from exponent */
 191
 192 #ifdef DEC
 193 q = (short *)&x;
 194 e = *q;                 /* short containing exponent */
 195 e = ((e >> 7) & 0377) - 0200;   /* the exponent */
 196 *q &= 0177;     /* strip exponent from x */
 197 *q |= 040000;   /* x now between 0.5 and 1 */
 198 #endif
 199
 200 #ifdef IBMPC
 201 x = frexp( x, &e );
 202 /*
 203 q = (short *)&x;
 204 q += 3;
 205 e = *q;
 206 e = ((e >> 4) & 0x0fff) - 0x3fe;
 207 *q &= 0x0f;
 208 *q |= 0x3fe0;
 209 */
 210 #endif
 211
 212 /* Equivalent C language standard library function: */
 213 #ifdef UNK
 214 x = frexp( x, &e );
 215 #endif
 216
 217 #ifdef MIEEE
 218 x = frexp( x, &e );
 219 #endif
 220
 221 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 222
 223 if( x < SQRTH )
 224         {
 225         e -= 1;
 226         x = ldexp( x, 1 ) - 1.0; /*  2x - 1  */
 227         }
 228 else
 229         {
 230         x = x - 1.0;
 231         }
 232
 233
 234 /* rational form */
 235 z = x*x;
 236 y = x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) );
 237 y = y - ldexp( z, -1 );   /*  y - 0.5 * x**2  */
 238
 239 /* multiply log of fraction by log10(e)
 240  * and base 2 exponent by log10(2)
 241  */
 242 z = (x + y) * L10EB;  /* accumulate terms in order of size */
 243 z += y * L10EA;
 244 z += x * L10EA;
 245 z += e * L102B;
 246 z += e * L102A;
 247
 248
 249 return( z );
 250 }