lib/libm/unity.c

   1 /*                                                      unity.c
   2  *
   3  * Relative error approximations for function arguments near
   4  * unity.
   5  *
   6  *    log1p(x) = log(1+x)
   7  *    expm1(x) = exp(x) - 1
   8  *    cosm1(x) = cos(x) - 1
   9  *
  10  */
  11
  12 #include "mconf.h"
  13
  14 #ifdef ANSIPROT
  15 extern int isnan (double);
  16 extern int isfinite (double);
  17 extern double log ( double );
  18 extern double polevl ( double, void *, int );
  19 extern double p1evl ( double, void *, int );
  20 extern double exp ( double );
  21 extern double cos ( double );
  22 #else
  23 double log(), polevl(), p1evl(), exp(), cos();
  24 int isnan(), isfinite();
  25 #endif
  26 extern double INFINITY;
  27
  28 /* log1p(x) = log(1 + x)  */
  29
  30 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  31  * 1/sqrt(2) <= x < sqrt(2)
  32  * Theoretical peak relative error = 2.32e-20
  33  */
  34 static double LP[] = {
  35  4.5270000862445199635215E-5,
  36  4.9854102823193375972212E-1,
  37  6.5787325942061044846969E0,
  38  2.9911919328553073277375E1,
  39  6.0949667980987787057556E1,
  40  5.7112963590585538103336E1,
  41  2.0039553499201281259648E1,
  42 };
  43 static double LQ[] = {
  44 /* 1.0000000000000000000000E0,*/
  45  1.5062909083469192043167E1,
  46  8.3047565967967209469434E1,
  47  2.2176239823732856465394E2,
  48  3.0909872225312059774938E2,
  49  2.1642788614495947685003E2,
  50  6.0118660497603843919306E1,
  51 };
  52
  53 #define SQRTH 0.70710678118654752440
  54 #define SQRT2 1.41421356237309504880
  55
  56 double log1p(x)
  57 double x;
  58 {
  59 double z;
  60
  61 z = 1.0 + x;
  62 if( (z < SQRTH) || (z > SQRT2) )
  63         return( log(z) );
  64 z = x*x;
  65 z = -0.5 * z + x * ( z * polevl( x, LP, 6 ) / p1evl( x, LQ, 6 ) );
  66 return (x + z);
  67 }
  68
  69
  70
  71 /* expm1(x) = exp(x) - 1  */
  72
  73 /*  e^x =  1 + 2x P(x^2)/( Q(x^2) - P(x^2) )
  74  * -0.5 <= x <= 0.5
  75  */
  76
  77 static double EP[3] = {
  78  1.2617719307481059087798E-4,
  79  3.0299440770744196129956E-2,
  80  9.9999999999999999991025E-1,
  81 };
  82 static double EQ[4] = {
  83  3.0019850513866445504159E-6,
  84  2.5244834034968410419224E-3,
  85  2.2726554820815502876593E-1,
  86  2.0000000000000000000897E0,
  87 };
  88
  89 double expm1(x)
  90 double x;
  91 {
  92 double r, xx;
  93
  94 #ifdef NANS
  95 if( isnan(x) )
  96         return(x);
  97 #endif
  98 #ifdef INFINITIES
  99 if( x == INFINITY )
 100         return(INFINITY);
 101 if( x == -INFINITY )
 102         return(-1.0);
 103 #endif
 104 if( (x < -0.5) || (x > 0.5) )
 105         return( exp(x) - 1.0 );
 106 xx = x * x;
 107 r = x * polevl( xx, EP, 2 );
 108 r = r/( polevl( xx, EQ, 3 ) - r );
 109 return (r + r);
 110 }
 111
 112
 113
 114 /* cosm1(x) = cos(x) - 1  */
 115
 116 static double coscof[7] = {
 117  4.7377507964246204691685E-14,
 118 -1.1470284843425359765671E-11,
 119  2.0876754287081521758361E-9,
 120 -2.7557319214999787979814E-7,
 121  2.4801587301570552304991E-5,
 122 -1.3888888888888872993737E-3,
 123  4.1666666666666666609054E-2,
 124 };
 125
 126 extern double PIO4;
 127
 128 double cosm1(x)
 129 double x;
 130 {
 131 double xx;
 132
 133 if( (x < -PIO4) || (x > PIO4) )
 134         return( cos(x) - 1.0 );
 135 xx = x * x;
 136 xx = -0.5*xx + xx * xx * polevl( xx, coscof, 6 );
 137 return xx;
 138 }