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 * lgam1p(x) = lgam(1+x)
10 *
11 */
12
13 /* Scipy changes:
14 * - 06-10-2016: added lgam1p
15 */
16
17 #include "mconf.h"
18
19 extern double MACHEP;
20
21
22
23 /* log1p(x) = log(1 + x) */
24
25 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
26 * 1/sqrt(2) <= x < sqrt(2)
27 * Theoretical peak relative error = 2.32e-20
28 */
29 static const double LP[] = {
30 4.5270000862445199635215E-5,
31 4.9854102823193375972212E-1,
32 6.5787325942061044846969E0,
33 2.9911919328553073277375E1,
34 6.0949667980987787057556E1,
35 5.7112963590585538103336E1,
36 2.0039553499201281259648E1,
37 };
38
39 static const double LQ[] = {
40 /* 1.0000000000000000000000E0, */
41 1.5062909083469192043167E1,
42 8.3047565967967209469434E1,
43 2.2176239823732856465394E2,
44 3.0909872225312059774938E2,
45 2.1642788614495947685003E2,
46 6.0118660497603843919306E1,
47 };
48
log1p(double x)49 double log1p(double x)
50 {
51 double z;
52
53 z = 1.0 + x;
54 if ((z < NPY_SQRT1_2) || (z > NPY_SQRT2))
55 return (log(z));
56 z = x * x;
57 z = -0.5 * z + x * (z * polevl(x, LP, 6) / p1evl(x, LQ, 6));
58 return (x + z);
59 }
60
61
62 /* log(1 + x) - x */
log1pmx(double x)63 double log1pmx(double x)
64 {
65 if (fabs(x) < 0.5) {
66 int n;
67 double xfac = x;
68 double term;
69 double res = 0;
70
71 for(n = 2; n < MAXITER; n++) {
72 xfac *= -x;
73 term = xfac / n;
74 res += term;
75 if (fabs(term) < MACHEP * fabs(res)) {
76 break;
77 }
78 }
79 return res;
80 }
81 else {
82 return log1p(x) - x;
83 }
84 }
85
86
87 /* expm1(x) = exp(x) - 1 */
88
89 /* e^x = 1 + 2x P(x^2)/( Q(x^2) - P(x^2) )
90 * -0.5 <= x <= 0.5
91 */
92
93 static double EP[3] = {
94 1.2617719307481059087798E-4,
95 3.0299440770744196129956E-2,
96 9.9999999999999999991025E-1,
97 };
98
99 static double EQ[4] = {
100 3.0019850513866445504159E-6,
101 2.5244834034968410419224E-3,
102 2.2726554820815502876593E-1,
103 2.0000000000000000000897E0,
104 };
105
expm1(double x)106 double expm1(double x)
107 {
108 double r, xx;
109
110 if (!cephes_isfinite(x)) {
111 if (cephes_isnan(x)) {
112 return x;
113 }
114 else if (x > 0) {
115 return x;
116 }
117 else {
118 return -1.0;
119 }
120
121 }
122 if ((x < -0.5) || (x > 0.5))
123 return (exp(x) - 1.0);
124 xx = x * x;
125 r = x * polevl(xx, EP, 2);
126 r = r / (polevl(xx, EQ, 3) - r);
127 return (r + r);
128 }
129
130
131
132 /* cosm1(x) = cos(x) - 1 */
133
134 static double coscof[7] = {
135 4.7377507964246204691685E-14,
136 -1.1470284843425359765671E-11,
137 2.0876754287081521758361E-9,
138 -2.7557319214999787979814E-7,
139 2.4801587301570552304991E-5,
140 -1.3888888888888872993737E-3,
141 4.1666666666666666609054E-2,
142 };
143
cosm1(double x)144 double cosm1(double x)
145 {
146 double xx;
147
148 if ((x < -NPY_PI_4) || (x > NPY_PI_4))
149 return (cos(x) - 1.0);
150 xx = x * x;
151 xx = -0.5 * xx + xx * xx * polevl(xx, coscof, 6);
152 return xx;
153 }
154
155
156 /* Compute lgam(x + 1) around x = 0 using its Taylor series. */
lgam1p_taylor(double x)157 static double lgam1p_taylor(double x)
158 {
159 int n;
160 double xfac, coeff, res;
161
162 if (x == 0) {
163 return 0;
164 }
165 res = -NPY_EULER * x;
166 xfac = -x;
167 for (n = 2; n < 42; n++) {
168 xfac *= -x;
169 coeff = zeta(n, 1) * xfac / n;
170 res += coeff;
171 if (fabs(coeff) < MACHEP * fabs(res)) {
172 break;
173 }
174 }
175
176 return res;
177 }
178
179
180 /* Compute lgam(x + 1). */
lgam1p(double x)181 double lgam1p(double x)
182 {
183 if (fabs(x) <= 0.5) {
184 return lgam1p_taylor(x);
185 } else if (fabs(x - 1) < 0.5) {
186 return log(x) + lgam1p_taylor(x - 1);
187 } else {
188 return lgam(x + 1);
189 }
190 }
191