libm/float/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 <math.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 float frexpf(float, int *), polevlf(float, float *, int);
  78
  79 float log10f(float xx)
  80 {
  81 float x, y, z;
  82 int e;
  83
  84 x = xx;
  85 /* Test for domain */
  86 if( x <= 0.0 )
  87         {
  88         if( x == 0.0 )
  89                 mtherr( fname, SING );
  90         else
  91                 mtherr( fname, DOMAIN );
  92         return( -MAXL10 );
  93         }
  94
  95 /* separate mantissa from exponent */
  96
  97 x = frexpf( x, &e );
  98
  99 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x) */
 100
 101 if( x < SQRTH )
 102         {
 103         e -= 1;
 104         x = 2.0*x - 1.0;
 105         }
 106 else
 107         {
 108         x = x - 1.0;
 109         }
 110
 111
 112 /* rational form */
 113 z = x*x;
 114 y = x * ( z * polevlf( x, P, 8 ) );
 115 y = y - 0.5 * z;   /*  y - 0.5 * x**2  */
 116
 117 /* multiply log of fraction by log10(e)
 118  * and base 2 exponent by log10(2)
 119  */
 120 z = (x + y) * L10EB;  /* accumulate terms in order of size */
 121 z += y * L10EA;
 122 z += x * L10EA;
 123 x = e;
 124 z += x * L102B;
 125 z += x * L102A;
 126
 127
 128 return( z );
 129 }