Quelle  log.c

/* log.c
 *
 * Natural logarithm
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, log();
 *
 * y = log( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the base e (2.718...) logarithm of x.
 *
 * The argument is separated into its exponent and fractional
 * parts.  If the exponent is between -1 and +1, the logarithm
 * of the fraction is approximated by
 *
 *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
 *
 * Otherwise, setting  z = 2(x-1)/x+1),
 *
 *     log(x) = z + z**3 P(z)/Q(z).
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0.5, 2.0    150000      1.44e-16    5.06e-17
 *    IEEE      +-MAXNUM    30000       1.20e-16    4.78e-17
 *    DEC       0, 10       170000      1.8e-17     6.3e-18
 *
 * In the tests over the interval [+-MAXNUM], the logarithms
 * of the random arguments were uniformly distributed over
 * [0, MAXLOG].
 *
 * ERROR MESSAGES:
 *
 * log singularity:  x = 0; returns -INFINITY
 * log domain:       x < 0; returns NAN
 */

/*
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/

#include "mconf.h"
static char fname[] = {"log"};

/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
 * 1/sqrt(2) <= x < sqrt(2)
 */
#ifdef UNK
static double P[] = {
    1.01875663804580931796E-4, 4.97494994976747001425E-1,
    4.70579119878881725854E0,  1.44989225341610930846E1,
    1.79368678507819816313E1,  7.70838733755885391666E0,
};
static double Q[] = {
    /* 1.00000000000000000000E0, */
    1.12873587189167450590E1, 4.52279145837532221105E1,
    8.29875266912776603211E1, 7.11544750618563894466E1,
    2.31251620126765340583E1,
};
#endif

#ifdef DEC
static unsigned short P[] = {
    0037777, 0127270, 0162547, 0057274, 0041001, 0054665, 0164317, 0005341,
    0041451, 0034104, 0031640, 0105773, 0041677, 0011276, 0123617, 0160135,
    0041701, 0126603, 0053215, 0117250, 0041420, 0115777, 0135206, 0030232,
};
static unsigned short Q[] = {
    /*0040200,0000000,0000000,0000000,*/
    0041220, 0144332, 0045272, 0174241, 0041742, 0164566, 0035720, 0130431,
    0042246, 0126327, 0166065, 0116357, 0042372, 0033420, 0157525, 0124560,
    0042271, 0167002, 0066537, 0172303, 0041730, 0164777, 0113711, 0044407,
};
#endif

#ifdef IBMPC
static unsigned short P[] = {
    0x1bb0, 0x93c3, 0xb4c2, 0x3f1a, 0x52f2, 0x3f56, 0xd6f5, 0x3fdf,
    0x6911, 0xed92, 0xd2ba, 0x4012, 0xeb2e, 0xc63e, 0xff72, 0x402c,
    0xc84d, 0x924b, 0xefd6, 0x4031, 0xdcf8, 0x7d7e, 0xd563, 0x401e,
};
static unsigned short Q[] = {
    /*0x0000,0x0000,0x0000,0x3ff0,*/
    0xef8e, 0xae97, 0x9320, 0x4026, 0xc033, 0x4e19, 0x9d2c,
    0x4046, 0xbdbd, 0xa326, 0xbf33, 0x4054, 0xae21, 0xeb5e,
    0xc9e2, 0x4051, 0x25b2, 0x9e1f, 0x200a, 0x4037,
};
#endif

#ifdef MIEEE
static unsigned short P[] = {
    0x3f1a, 0xb4c2, 0x93c3, 0x1bb0, 0x3fdf, 0xd6f5, 0x3f56, 0x52f2,
    0x4012, 0xd2ba, 0xed92, 0x6911, 0x402c, 0xff72, 0xc63e, 0xeb2e,
    0x4031, 0xefd6, 0x924b, 0xc84d, 0x401e, 0xd563, 0x7d7e, 0xdcf8,
};
static unsigned short Q[] = {
    /*0x3ff0,0x0000,0x0000,0x0000,*/
    0x4026, 0x9320, 0xae97, 0xef8e, 0x4046, 0x9d2c, 0x4e19,
    0xc033, 0x4054, 0xbf33, 0xa326, 0xbdbd, 0x4051, 0xc9e2,
    0xeb5e, 0xae21, 0x4037, 0x200a, 0x9e1f, 0x25b2,
};
#endif

/* Coefficients for log(x) = z + z**3 P(z)/Q(z),
 * where z = 2(x-1)/(x+1)
 * 1/sqrt(2) <= x < sqrt(2)
 */

#ifdef UNK
static double R[3] = {
    -7.89580278884799154124E-1,
    1.63866645699558079767E1,
    -6.41409952958715622951E1,
};
static double S[3] = {
    /* 1.00000000000000000000E0,*/
    -3.56722798256324312549E1,
    3.12093766372244180303E2,
    -7.69691943550460008604E2,
};
#endif
#ifdef DEC
static unsigned short R[12] = {
    0140112, 0020756, 0161540, 0072035, 0041203, 0013743,
    0114023, 0155527, 0141600, 0044060, 0104421, 0050400,
};
static unsigned short S[12] = {
    /*0040200,0000000,0000000,0000000,*/
    0141416, 0130152, 0017543, 0064122, 0042234, 0006000,
    0104527, 0020155, 0142500, 0066110, 0146631, 0174731,
};
#endif
#ifdef IBMPC
static unsigned short R[12] = {
    0x0e84, 0xdc6c, 0x443d, 0xbfe9, 0x7b6b, 0x7302,
    0x62fc, 0x4030, 0x2a20, 0x1122, 0x0906, 0xc050,
};
static unsigned short S[12] = {
    /*0x0000,0x0000,0x0000,0x3ff0,*/
    0x6d0a, 0x43ec, 0xd60d, 0xc041, 0xe40e, 0x112a,
    0x8180, 0x4073, 0x3f3b, 0x19b3, 0x0d89, 0xc088,
};
#endif
#ifdef MIEEE
static unsigned short R[12] = {
    0xbfe9, 0x443d, 0xdc6c, 0x0e84, 0x4030, 0x62fc,
    0x7302, 0x7b6b, 0xc050, 0x0906, 0x1122, 0x2a20,
};
static unsigned short S[12] = {
    /*0x3ff0,0x0000,0x0000,0x0000,*/
    0xc041, 0xd60d, 0x43ec, 0x6d0a, 0x4073, 0x8180,
    0x112a, 0xe40e, 0xc088, 0x0d89, 0x19b3, 0x3f3b,
};
#endif

#ifdef ANSIPROT
extern double frexp(double, int *);
extern double ldexp(double, int);
extern double polevl(double, void *, int);
extern double p1evl(double, void *, int);
extern int isnan(double);
extern int isfinite(double);
#else
double frexp(), ldexp(), polevl(), p1evl();
int isnan(), isfinite();
#endif
#define SQRTH 0.70710678118654752440
extern double INFINITY, NAN;

double log(x) double x;
{
  int e;
#ifdef DEC
  short *q;
#endif
  double y, z;

#ifdef NANS
  if (isnan(x))
    return (x);
#endif
#ifdef INFINITIES
  if (x == INFINITY)
    return (x);
#endif
  /* Test for domain */
  if (x <= 0.0) {
    if (x == 0.0) {
      mtherr(fname, SING);
      return (-INFINITY);
    } else {
      mtherr(fname, DOMAIN);
      return (NAN);
    }
  }

  /* separate mantissa from exponent */

#ifdef DEC
  q = (short *)&x;
  e = *q;                       /* short containing exponent */
  e = ((e >> 7) & 0377) - 0200; /* the exponent */
  *q &= 0177;                   /* strip exponent from x */
  *q |= 040000;                 /* x now between 0.5 and 1 */
#endif

/* Note, frexp is used so that denormal numbers
 * will be handled properly.
 */
#ifdef IBMPC
  x = frexp(x, &e);
/*
q = (short *)&x;
q += 3;
e = *q;
e = ((e >> 4) & 0x0fff) - 0x3fe;
*q &= 0x0f;
*q |= 0x3fe0;
*/
#endif

/* Equivalent C language standard library function: */
#ifdef UNK
  x = frexp(x, &e);
#endif

#ifdef MIEEE
  x = frexp(x, &e);
#endif

  /* logarithm using log(x) = z + z**3 P(z)/Q(z),
   * where z = 2(x-1)/x+1)
   */

  if ((e > 2) || (e < -2)) {
    if (x < SQRTH) { /* 2( 2x-1 )/( 2x+1 ) */
      e -= 1;
      z = x - 0.5;
      y = 0.5 * z + 0.5;
    } else { /*  2 (x-1)/(x+1)   */
      z = x - 0.5;
      z -= 0.5;
      y = 0.5 * x + 0.5;
    }

    x = z / y;

    /* rational form */
    z = x * x;
    z = x * (z * polevl(z, R, 2) / p1evl(z, S, 3));
    y = e;
    z = z - y * 2.121944400546905827679e-4;
    z = z + x;
    z = z + e * 0.693359375;
    goto ldone;
  }

  /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */

  if (x < SQRTH) {
    e -= 1;
    x = ldexp(x, 1) - 1.0; /*  2x - 1  */
  } else {
    x = x - 1.0;
  }

  /* rational form */
  z = x * x;
#if DEC
  y = x * (z * polevl(x, P, 5) / p1evl(x, Q, 6));
#else
  y = x * (z * polevl(x, P, 5) / p1evl(x, Q, 5));
#endif
  if (e)
    y = y - e * 2.121944400546905827679e-4;
  y = y - ldexp(z, -1); /*  y - 0.5 * z  */
  z = x + y;
  if (e)
    z = z + e * 0.693359375;

ldone:

  return (z);
}