newlib/libm/mathfp/sf_logarithm.c

   1
   2 /* @(#)z_logarithmf.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  * Logarithm
  11  *
  12  * Input:
  13  *   x - floating point value
  14  *   ten - indicates base ten numbers
  15  *
  16  * Output:
  17  *   logarithm of x
  18  *
  19  * Description:
  20  *   This routine calculates logarithms.
  21  *
  22  *****************************************************************/
  23
  24 #include "fdlibm.h"
  25 #include "zmath.h"
  26
  27 static const float a[] = { -0.5527074855 };
  28 static const float b[] = { -0.6632718214e+1 };
  29 static const float C1 = 0.693145752;
  30 static const float C2 = 1.428606820e-06;
  31 static const float C3 = 0.4342944819;
  32
  33 float
  34 logarithmf (float x,
  35         int ten)
  36 {
  37   int N;
  38   float f, w, z;
  39
  40   /* Check for domain/range errors here. */
  41   if (x == 0.0)
  42     {
  43       errno = ERANGE;
  44       return (-z_infinity_f.f);
  45     }
  46   else if (x < 0.0)
  47     {
  48       errno = EDOM;
  49       return (z_notanum_f.f);
  50     }
  51   else if (!isfinite(x))
  52     {
  53       if (isnanf(x))
  54         return (z_notanum_f.f);
  55       else
  56         return (z_infinity_f.f);
  57     }
  58
  59   /* Get the exponent and mantissa where x = f * 2^N. */
  60   f = frexpf (x, &N);
  61
  62   z = f - 0.5;
  63
  64   if (f > __SQRT_HALF)
  65     z = (z - 0.5) / (f * 0.5 + 0.5);
  66   else
  67     {
  68       N--;
  69       z /= (z * 0.5 + 0.5);
  70     }
  71   w = z * z;
  72
  73   /* Use Newton's method with 4 terms. */
  74   z += z * w * (a[0]) / ((w + 1.0) * w + b[0]);
  75
  76   if (N != 0)
  77     z = (N * C2 + z) + N * C1;
  78
  79   if (ten)
  80     z *= C3;
  81
  82   return (z);
  83 }