lib/libm/log10f.c

   1 /*                                                      log10f.c
   2  *
   3  *      Common logarithm
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, log10f();
  10  *
  11  * y = log10f( 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).
  23  *
  24  *
  25  *
  26  * ACCURACY:
  27  *
  28  *                      Relative error:
  29  * arithmetic   domain     # trials      peak         rms
  30  *    IEEE      0.5, 2.0    100000      1.3e-7      3.4e-8
  31  *    IEEE      0, MAXNUMF  100000      1.3e-7      2.6e-8
  32  *
  33  * In the tests over the interval [0, MAXNUM], the logarithms
  34  * of the random arguments were uniformly distributed over
  35  * [-MAXL10, MAXL10].
  36  *
  37  * ERROR MESSAGES:
  38  *
  39  * log10f singularity:  x = 0; returns -MAXL10
  40  * log10f domain:       x < 0; returns -MAXL10
  41  * MAXL10 = 38.230809449325611792
  42  */
  43 \f
  44 /*
  45 Cephes Math Library Release 2.1:  December, 1988
  46 Copyright 1984, 1987, 1988 by Stephen L. Moshier
  47 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  48 */
  49
  50 #include "mconf.h"
  51 static char fname[] = {"log10"};
  52
  53 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  54  * 1/sqrt(2) <= x < sqrt(2)
  55  */
  56 static float P[] = {
  57  7.0376836292E-2,
  58 -1.1514610310E-1,
  59  1.1676998740E-1,
  60 -1.2420140846E-1,
  61  1.4249322787E-1,
  62 -1.6668057665E-1,
  63  2.0000714765E-1,
  64 -2.4999993993E-1,
  65  3.3333331174E-1
  66 };
  67
  68
  69 #define SQRTH 0.70710678118654752440
  70 #define L102A 3.0078125E-1
  71 #define L102B 2.48745663981195213739E-4
  72 #define L10EA 4.3359375E-1
  73 #define L10EB 7.00731903251827651129E-4
  74
  75 static float MAXL10 = 38.230809449325611792;
  76
  77 #ifdef ANSIC
  78 float frexpf(float, int *), polevlf(float, float *, int);
  79
  80 float log10f(float xx)
  81 #else
  82 float frexpf(), polevlf();
  83
  84 float log10f(xx)
  85 double xx;
  86 #endif
  87 {
  88 float x, y, z;
  89 int e;
  90
  91 x = xx;
  92 /* Test for domain */
  93 if( x <= 0.0 )
  94         {
  95         if( x == 0.0 )
  96                 mtherr( fname, SING );
  97         else
  98                 mtherr( fname, DOMAIN );
  99         return( -MAXL10 );
 100         }
 101
 102 /* separate mantissa from exponent */
 103
 104 x = frexpf( x, &e );
 105
 106 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x) */
 107
 108 if( x < SQRTH )
 109         {
 110         e -= 1;
 111         x = 2.0*x - 1.0;
 112         }
 113 else
 114         {
 115         x = x - 1.0;
 116         }
 117
 118
 119 /* rational form */
 120 z = x*x;
 121 y = x * ( z * polevlf( x, P, 8 ) );
 122 y = y - 0.5 * z;   /*  y - 0.5 * x**2  */
 123
 124 /* multiply log of fraction by log10(e)
 125  * and base 2 exponent by log10(2)
 126  */
 127 z = (x + y) * L10EB;  /* accumulate terms in order of size */
 128 z += y * L10EA;
 129 z += x * L10EA;
 130 x = e;
 131 z += x * L102B;
 132 z += x * L102A;
 133
 134
 135 return( z );
 136 }