lib/libm/log2f.c

   1 /*                                                      log2f.c
   2  *
   3  *      Base 2 logarithm
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, log2f();
  10  *
  11  * y = log2f( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns the base 2 logarithm of x.
  18  *
  19  * The argument is separated into its exponent and fractional
  20  * parts.  If the exponent is between -1 and +1, the base e
  21  * logarithm of the fraction is approximated by
  22  *
  23  *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
  24  *
  25  * Otherwise, setting  z = 2(x-1)/x+1),
  26  *
  27  *     log(x) = z + z**3 P(z)/Q(z).
  28  *
  29  *
  30  *
  31  * ACCURACY:
  32  *
  33  *                      Relative error:
  34  * arithmetic   domain     # trials      peak         rms
  35  *    IEEE      exp(+-88)   100000      1.1e-7      2.4e-8
  36  *    IEEE      0.5, 2.0    100000      1.1e-7      3.0e-8
  37  *
  38  * In the tests over the interval [exp(+-88)], the logarithms
  39  * of the random arguments were uniformly distributed.
  40  *
  41  * ERROR MESSAGES:
  42  *
  43  * log singularity:  x = 0; returns MINLOGF/log(2)
  44  * log domain:       x < 0; returns MINLOGF/log(2)
  45  */
  46 \f
  47 /*
  48 Cephes Math Library Release 2.2:  June, 1992
  49 Copyright 1984, 1992 by Stephen L. Moshier
  50 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  51 */
  52
  53 #include "mconf.h"
  54 static char fname[] = {"log2"};
  55
  56 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)
  57  * 1/sqrt(2) <= x < sqrt(2)
  58  */
  59
  60 static float P[] = {
  61  7.0376836292E-2,
  62 -1.1514610310E-1,
  63  1.1676998740E-1,
  64 -1.2420140846E-1,
  65  1.4249322787E-1,
  66 -1.6668057665E-1,
  67  2.0000714765E-1,
  68 -2.4999993993E-1,
  69  3.3333331174E-1
  70 };
  71
  72 #define LOG2EA 0.44269504088896340735992
  73 #define SQRTH 0.70710678118654752440
  74 extern float MINLOGF, LOGE2F;
  75
  76 #ifdef ANSIC
  77 float frexpf(float, int *), polevlf(float, float *, int);
  78
  79 float log2f(float xx)
  80 #else
  81 float frexpf(), polevlf();
  82
  83 float log2f(xx)
  84 double xx;
  85 #endif
  86 {
  87 float x, y, z;
  88 int e;
  89
  90 x = xx;
  91 /* Test for domain */
  92 if( x <= 0.0 )
  93         {
  94         if( x == 0.0 )
  95                 mtherr( fname, SING );
  96         else
  97                 mtherr( fname, DOMAIN );
  98         return( MINLOGF/LOGE2F );
  99         }
 100
 101 /* separate mantissa from exponent */
 102 x = frexpf( x, &e );
 103
 104
 105 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 106
 107 if( x < SQRTH )
 108         {
 109         e -= 1;
 110         x = 2.0*x - 1.0;
 111         }
 112 else
 113         {
 114         x = x - 1.0;
 115         }
 116
 117 z = x*x;
 118 y = x * ( z * polevlf( x, P, 8 ) );
 119 y = y - 0.5 * z;   /*  y - 0.5 * x**2  */
 120
 121
 122 /* Multiply log of fraction by log2(e)
 123  * and base 2 exponent by 1
 124  *
 125  * ***CAUTION***
 126  *
 127  * This sequence of operations is critical and it may
 128  * be horribly defeated by some compiler optimizers.
 129  */
 130 z = y * LOG2EA;
 131 z += x * LOG2EA;
 132 z += y;
 133 z += x;
 134 z += (float )e;
 135 return( z );
 136 }