1 /* Definitions of some C99 math library functions, for those platforms
2 that don't implement these functions already. */
3
4 #include "Python.h"
5 #include <float.h>
6 #include "_math.h"
7
8 /* The following copyright notice applies to the original
9 implementations of acosh, asinh and atanh. */
10
11 /*
12 * ====================================================
13 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
14 *
15 * Developed at SunPro, a Sun Microsystems, Inc. business.
16 * Permission to use, copy, modify, and distribute this
17 * software is freely granted, provided that this notice
18 * is preserved.
19 * ====================================================
20 */
21
22 static const double ln2 = 6.93147180559945286227E-01;
23 static const double two_pow_m28 = 3.7252902984619141E-09; /* 2**-28 */
24 static const double two_pow_p28 = 268435456.0; /* 2**28 */
25 #ifndef Py_NAN
26 static const double zero = 0.0;
27 #endif
28
29 /* acosh(x)
30 * Method :
31 * Based on
32 * acosh(x) = log [ x + sqrt(x*x-1) ]
33 * we have
34 * acosh(x) := log(x)+ln2, if x is large; else
35 * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
36 * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
37 *
38 * Special cases:
39 * acosh(x) is NaN with signal if x<1.
40 * acosh(NaN) is NaN without signal.
41 */
42
43 double
_Py_acosh(double x)44 _Py_acosh(double x)
45 {
46 if (Py_IS_NAN(x)) {
47 return x+x;
48 }
49 if (x < 1.) { /* x < 1; return a signaling NaN */
50 errno = EDOM;
51 #ifdef Py_NAN
52 return Py_NAN;
53 #else
54 return (x-x)/(x-x);
55 #endif
56 }
57 else if (x >= two_pow_p28) { /* x > 2**28 */
58 if (Py_IS_INFINITY(x)) {
59 return x+x;
60 }
61 else {
62 return log(x)+ln2; /* acosh(huge)=log(2x) */
63 }
64 }
65 else if (x == 1.) {
66 return 0.0; /* acosh(1) = 0 */
67 }
68 else if (x > 2.) { /* 2 < x < 2**28 */
69 double t = x*x;
70 return log(2.0*x - 1.0 / (x + sqrt(t - 1.0)));
71 }
72 else { /* 1 < x <= 2 */
73 double t = x - 1.0;
74 return m_log1p(t + sqrt(2.0*t + t*t));
75 }
76 }
77
78
79 /* asinh(x)
80 * Method :
81 * Based on
82 * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
83 * we have
84 * asinh(x) := x if 1+x*x=1,
85 * := sign(x)*(log(x)+ln2)) for large |x|, else
86 * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
87 * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
88 */
89
90 double
_Py_asinh(double x)91 _Py_asinh(double x)
92 {
93 double w;
94 double absx = fabs(x);
95
96 if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) {
97 return x+x;
98 }
99 if (absx < two_pow_m28) { /* |x| < 2**-28 */
100 return x; /* return x inexact except 0 */
101 }
102 if (absx > two_pow_p28) { /* |x| > 2**28 */
103 w = log(absx)+ln2;
104 }
105 else if (absx > 2.0) { /* 2 < |x| < 2**28 */
106 w = log(2.0*absx + 1.0 / (sqrt(x*x + 1.0) + absx));
107 }
108 else { /* 2**-28 <= |x| < 2= */
109 double t = x*x;
110 w = m_log1p(absx + t / (1.0 + sqrt(1.0 + t)));
111 }
112 return copysign(w, x);
113
114 }
115
116 /* atanh(x)
117 * Method :
118 * 1.Reduced x to positive by atanh(-x) = -atanh(x)
119 * 2.For x>=0.5
120 * 1 2x x
121 * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * -------)
122 * 2 1 - x 1 - x
123 *
124 * For x<0.5
125 * atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
126 *
127 * Special cases:
128 * atanh(x) is NaN if |x| >= 1 with signal;
129 * atanh(NaN) is that NaN with no signal;
130 *
131 */
132
133 double
_Py_atanh(double x)134 _Py_atanh(double x)
135 {
136 double absx;
137 double t;
138
139 if (Py_IS_NAN(x)) {
140 return x+x;
141 }
142 absx = fabs(x);
143 if (absx >= 1.) { /* |x| >= 1 */
144 errno = EDOM;
145 #ifdef Py_NAN
146 return Py_NAN;
147 #else
148 return x/zero;
149 #endif
150 }
151 if (absx < two_pow_m28) { /* |x| < 2**-28 */
152 return x;
153 }
154 if (absx < 0.5) { /* |x| < 0.5 */
155 t = absx+absx;
156 t = 0.5 * m_log1p(t + t*absx / (1.0 - absx));
157 }
158 else { /* 0.5 <= |x| <= 1.0 */
159 t = 0.5 * m_log1p((absx + absx) / (1.0 - absx));
160 }
161 return copysign(t, x);
162 }
163
164 /* Mathematically, expm1(x) = exp(x) - 1. The expm1 function is designed
165 to avoid the significant loss of precision that arises from direct
166 evaluation of the expression exp(x) - 1, for x near 0. */
167
168 double
_Py_expm1(double x)169 _Py_expm1(double x)
170 {
171 /* For abs(x) >= log(2), it's safe to evaluate exp(x) - 1 directly; this
172 also works fine for infinities and nans.
173
174 For smaller x, we can use a method due to Kahan that achieves close to
175 full accuracy.
176 */
177
178 if (fabs(x) < 0.7) {
179 double u;
180 u = exp(x);
181 if (u == 1.0)
182 return x;
183 else
184 return (u - 1.0) * x / log(u);
185 }
186 else
187 return exp(x) - 1.0;
188 }
189
190 /* log1p(x) = log(1+x). The log1p function is designed to avoid the
191 significant loss of precision that arises from direct evaluation when x is
192 small. */
193
194 #ifdef HAVE_LOG1P
195
196 double
_Py_log1p(double x)197 _Py_log1p(double x)
198 {
199 /* Some platforms supply a log1p function but don't respect the sign of
200 zero: log1p(-0.0) gives 0.0 instead of the correct result of -0.0.
201
202 To save fiddling with configure tests and platform checks, we handle the
203 special case of zero input directly on all platforms.
204 */
205 if (x == 0.0) {
206 return x;
207 }
208 else {
209 return log1p(x);
210 }
211 }
212
213 #else
214
215 double
_Py_log1p(double x)216 _Py_log1p(double x)
217 {
218 /* For x small, we use the following approach. Let y be the nearest float
219 to 1+x, then
220
221 1+x = y * (1 - (y-1-x)/y)
222
223 so log(1+x) = log(y) + log(1-(y-1-x)/y). Since (y-1-x)/y is tiny, the
224 second term is well approximated by (y-1-x)/y. If abs(x) >=
225 DBL_EPSILON/2 or the rounding-mode is some form of round-to-nearest
226 then y-1-x will be exactly representable, and is computed exactly by
227 (y-1)-x.
228
229 If abs(x) < DBL_EPSILON/2 and the rounding mode is not known to be
230 round-to-nearest then this method is slightly dangerous: 1+x could be
231 rounded up to 1+DBL_EPSILON instead of down to 1, and in that case
232 y-1-x will not be exactly representable any more and the result can be
233 off by many ulps. But this is easily fixed: for a floating-point
234 number |x| < DBL_EPSILON/2., the closest floating-point number to
235 log(1+x) is exactly x.
236 */
237
238 double y;
239 if (fabs(x) < DBL_EPSILON/2.) {
240 return x;
241 }
242 else if (-0.5 <= x && x <= 1.) {
243 /* WARNING: it's possible that an overeager compiler
244 will incorrectly optimize the following two lines
245 to the equivalent of "return log(1.+x)". If this
246 happens, then results from log1p will be inaccurate
247 for small x. */
248 y = 1.+x;
249 return log(y)-((y-1.)-x)/y;
250 }
251 else {
252 /* NaNs and infinities should end up here */
253 return log(1.+x);
254 }
255 }
256
257 #endif /* ifdef HAVE_LOG1P */
258