1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
/* unity.c
*
* Relative error approximations for function arguments near
* unity.
*
* md_log1p(x) = md_log(1+x)
* expm1(x) = md_exp(x) - 1
* cosm1(x) = md_cos(x) - 1
*
*/
#include "mconf.h"
#ifdef ANSIPROT
extern int isnan (double);
extern int isfinite (double);
extern double md_log ( double );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double md_exp ( double );
extern double md_cos ( double );
#else
double md_log(), polevl(), p1evl(), md_exp(), md_cos();
int isnan(), isfinite();
#endif
extern double INFINITY;
/* md_log1p(x) = md_log(1 + x) */
/* Coefficients for md_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 md_log1p(x)
double x;
{
double z;
z = 1.0 + x;
if( (z < SQRTH) || (z > SQRT2) )
return( md_log(z) );
z = x*x;
z = -0.5 * z + x * ( z * polevl( x, LP, 6 ) / p1evl( x, LQ, 6 ) );
return (x + z);
}
/* expm1(x) = md_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( md_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) = md_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( md_cos(x) - 1.0 );
xx = x * x;
xx = -0.5*xx + xx * xx * polevl( xx, coscof, 6 );
return xx;
}
|
