Quelle  unity.c

/* unity.c
 *
 * Relative error approximations for function arguments near
 * unity.
 *
 *    log1p(x) = log(1+x)
 *    expm1(x) = exp(x) - 1
 *    cosm1(x) = cos(x) - 1
 *
 */

#include "mconf.h"

#ifdef ANSIPROT
extern int isnan(double);
extern int isfinite(double);
extern double log(double);
extern double polevl(double, void *, int);
extern double p1evl(double, void *, int);
extern double exp(double);
extern double cos(double);
#else
double log(), polevl(), p1evl(), exp(), cos();
int isnan(), isfinite();
#endif
extern double INFINITY;

/* log1p(x) = log(1 + x)  */

/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
 * 1/sqrt(2) <= x < sqrt(2)
 * Theoretical peak relative error = 2.32e-20
 */
static double LP[] = {
    4.5270000862445199635215E-5, 4.9854102823193375972212E-1,
    6.5787325942061044846969E0,  2.9911919328553073277375E1,
    6.0949667980987787057556E1,  5.7112963590585538103336E1,
    2.0039553499201281259648E1,
};
static double LQ[] = {
    /* 1.0000000000000000000000E0,*/
    1.5062909083469192043167E1, 8.3047565967967209469434E1,
    2.2176239823732856465394E2, 3.0909872225312059774938E2,
    2.1642788614495947685003E2, 6.0118660497603843919306E1,
};

#define SQRTH 0.70710678118654752440
#define SQRT2 1.41421356237309504880

double log1p(x) double x;
{
  double z;

  z = 1.0 + x;
  if ((z < SQRTH) || (z > SQRT2))
    return (log(z));
  z = x * x;
  z = -0.5 * z + x * (z * polevl(x, LP, 6) / p1evl(x, LQ, 6));
  return (x + z);
}

/* expm1(x) = exp(x) - 1  */

/*  e^x =  1 + 2x P(x^2)/( Q(x^2) - P(x^2) )
 * -0.5 <= x <= 0.5
 */

static double EP[3] = {
    1.2617719307481059087798E-4,
    3.0299440770744196129956E-2,
    9.9999999999999999991025E-1,
};
static double EQ[4] = {
    3.0019850513866445504159E-6,
    2.5244834034968410419224E-3,
    2.2726554820815502876593E-1,
    2.0000000000000000000897E0,
};

double expm1(x) double x;
{
  double r, xx;

#ifdef NANS
  if (isnan(x))
    return (x);
#endif
#ifdef INFINITIES
  if (x == INFINITY)
    return (INFINITY);
  if (x == -INFINITY)
    return (-1.0);
#endif
  if ((x < -0.5) || (x > 0.5))
    return (exp(x) - 1.0);
  xx = x * x;
  r = x * polevl(xx, EP, 2);
  r = r / (polevl(xx, EQ, 3) - r);
  return (r + r);
}

/* cosm1(x) = cos(x) - 1  */

static double coscof[7] = {
    4.7377507964246204691685E-14, -1.1470284843425359765671E-11,
    2.0876754287081521758361E-9,  -2.7557319214999787979814E-7,
    2.4801587301570552304991E-5,  -1.3888888888888872993737E-3,
    4.1666666666666666609054E-2,
};

extern double PIO4;

double cosm1(x) double x;
{
  double xx;

  if ((x < -PIO4) || (x > PIO4))
    return (cos(x) - 1.0);
  xx = x * x;
  xx = -0.5 * xx + xx * xx * polevl(xx, coscof, 6);
  return xx;
}