newlib/libm/mathfp/s_logarithm.c

   1
   2 /* @(#)z_logarithm.c 1.0 98/08/13 */
   3 /******************************************************************
   4  * The following routines are coded directly from the algorithms
   5  * and coefficients given in "Software Manual for the Elementary
   6  * Functions" by William J. Cody, Jr. and William Waite, Prentice
   7  * Hall, 1980.
   8  ******************************************************************/
   9
  10 /*
  11 FUNCTION
  12        <<log>>, <<logf>>, <<log10>>, <<log10f>>, <<logarithm>>, <<logarithmf>>---natural or base 10 logarithms
  13
  14 INDEX
  15     log
  16 INDEX
  17     logf
  18 INDEX
  19     log10
  20 INDEX
  21     log10f
  22
  23 SYNOPSIS
  24        #include <math.h>
  25        double log(double <[x]>);
  26        float logf(float <[x]>);
  27        double log10(double <[x]>);
  28        float log10f(float <[x]>);
  29
  30 DESCRIPTION
  31 Return the natural or base 10 logarithm of <[x]>, that is, its logarithm base e
  32 (where e is the base of the natural system of logarithms, 2.71828@dots{}) or
  33 base 10.
  34 <<log>> and <<logf>> are identical save for the return and argument types.
  35 <<log10>> and <<log10f>> are identical save for the return and argument types.
  36
  37 RETURNS
  38 Normally, returns the calculated value.  When <[x]> is zero, the
  39 returned value is <<-HUGE_VAL>> and <<errno>> is set to <<ERANGE>>.
  40 When <[x]> is negative, the returned value is <<-HUGE_VAL>> and
  41 <<errno>> is set to <<EDOM>>.
  42
  43 PORTABILITY
  44 <<log>> is ANSI. <<logf>> is an extension.
  45
  46 <<log10>> is ANSI. <<log10f>> is an extension.
  47 */
  48
  49
  50 /******************************************************************
  51  * Logarithm
  52  *
  53  * Input:
  54  *   x - floating point value
  55  *   ten - indicates base ten numbers
  56  *
  57  * Output:
  58  *   logarithm of x
  59  *
  60  * Description:
  61  *   This routine calculates logarithms.
  62  *
  63  *****************************************************************/
  64
  65 #include "fdlibm.h"
  66 #include "zmath.h"
  67
  68 #ifndef _DOUBLE_IS_32BITS
  69
  70 static const double a[] = { -0.64124943423745581147e+02,
  71                              0.16383943563021534222e+02,
  72                             -0.78956112887481257267 };
  73 static const double b[] = { -0.76949932108494879777e+03,
  74                              0.31203222091924532844e+03,
  75                             -0.35667977739034646171e+02 };
  76 static const double C1 =  22713.0 / 32768.0;
  77 static const double C2 =  1.428606820309417232e-06;
  78 static const double C3 =  0.43429448190325182765;
  79
  80 double
  81 logarithm (double x,
  82         int ten)
  83 {
  84   int N;
  85   double f, w, z;
  86
  87   /* Check for range and domain errors here. */
  88   if (x == 0.0)
  89     {
  90       errno = ERANGE;
  91       return (-z_infinity.d);
  92     }
  93   else if (x < 0.0)
  94     {
  95       errno = EDOM;
  96       return (z_notanum.d);
  97     }
  98   else if (!isfinite(x))
  99     {
 100       if (isnan(x))
 101         return (z_notanum.d);
 102       else
 103         return (z_infinity.d);
 104     }
 105
 106   /* Get the exponent and mantissa where x = f * 2^N. */
 107   f = frexp (x, &N);
 108
 109   z = f - 0.5;
 110
 111   if (f > __SQRT_HALF)
 112     z = (z - 0.5) / (f * 0.5 + 0.5);
 113   else
 114     {
 115       N--;
 116       z /= (z * 0.5 + 0.5);
 117     }
 118   w = z * z;
 119
 120   /* Use Newton's method with 4 terms. */
 121   z += z * w * ((a[2] * w + a[1]) * w + a[0]) / (((w + b[2]) * w + b[1]) * w + b[0]);
 122
 123   if (N != 0)
 124     z = (N * C2 + z) + N * C1;
 125
 126   if (ten)
 127     z *= C3;
 128
 129   return (z);
 130 }
 131
 132 #endif /* _DOUBLE_IS_32BITS */