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