]>
Commit | Line | Data |
---|---|---|
4710c53d | 1 | /* Math module -- standard C math library functions, pi and e */\r |
2 | \r | |
3 | /* Here are some comments from Tim Peters, extracted from the\r | |
4 | discussion attached to http://bugs.python.org/issue1640. They\r | |
5 | describe the general aims of the math module with respect to\r | |
6 | special values, IEEE-754 floating-point exceptions, and Python\r | |
7 | exceptions.\r | |
8 | \r | |
9 | These are the "spirit of 754" rules:\r | |
10 | \r | |
11 | 1. If the mathematical result is a real number, but of magnitude too\r | |
12 | large to approximate by a machine float, overflow is signaled and the\r | |
13 | result is an infinity (with the appropriate sign).\r | |
14 | \r | |
15 | 2. If the mathematical result is a real number, but of magnitude too\r | |
16 | small to approximate by a machine float, underflow is signaled and the\r | |
17 | result is a zero (with the appropriate sign).\r | |
18 | \r | |
19 | 3. At a singularity (a value x such that the limit of f(y) as y\r | |
20 | approaches x exists and is an infinity), "divide by zero" is signaled\r | |
21 | and the result is an infinity (with the appropriate sign). This is\r | |
22 | complicated a little by that the left-side and right-side limits may\r | |
23 | not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0\r | |
24 | from the positive or negative directions. In that specific case, the\r | |
25 | sign of the zero determines the result of 1/0.\r | |
26 | \r | |
27 | 4. At a point where a function has no defined result in the extended\r | |
28 | reals (i.e., the reals plus an infinity or two), invalid operation is\r | |
29 | signaled and a NaN is returned.\r | |
30 | \r | |
31 | And these are what Python has historically /tried/ to do (but not\r | |
32 | always successfully, as platform libm behavior varies a lot):\r | |
33 | \r | |
34 | For #1, raise OverflowError.\r | |
35 | \r | |
36 | For #2, return a zero (with the appropriate sign if that happens by\r | |
37 | accident ;-)).\r | |
38 | \r | |
39 | For #3 and #4, raise ValueError. It may have made sense to raise\r | |
40 | Python's ZeroDivisionError in #3, but historically that's only been\r | |
41 | raised for division by zero and mod by zero.\r | |
42 | \r | |
43 | */\r | |
44 | \r | |
45 | /*\r | |
46 | In general, on an IEEE-754 platform the aim is to follow the C99\r | |
47 | standard, including Annex 'F', whenever possible. Where the\r | |
48 | standard recommends raising the 'divide-by-zero' or 'invalid'\r | |
49 | floating-point exceptions, Python should raise a ValueError. Where\r | |
50 | the standard recommends raising 'overflow', Python should raise an\r | |
51 | OverflowError. In all other circumstances a value should be\r | |
52 | returned.\r | |
53 | */\r | |
54 | \r | |
55 | #include "Python.h"\r | |
56 | #include "_math.h"\r | |
57 | \r | |
58 | #ifdef _OSF_SOURCE\r | |
59 | /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */\r | |
60 | extern double copysign(double, double);\r | |
61 | #endif\r | |
62 | \r | |
63 | /*\r | |
64 | sin(pi*x), giving accurate results for all finite x (especially x\r | |
65 | integral or close to an integer). This is here for use in the\r | |
66 | reflection formula for the gamma function. It conforms to IEEE\r | |
67 | 754-2008 for finite arguments, but not for infinities or nans.\r | |
68 | */\r | |
69 | \r | |
70 | static const double pi = 3.141592653589793238462643383279502884197;\r | |
71 | static const double sqrtpi = 1.772453850905516027298167483341145182798;\r | |
72 | \r | |
73 | static double\r | |
74 | sinpi(double x)\r | |
75 | {\r | |
76 | double y, r;\r | |
77 | int n;\r | |
78 | /* this function should only ever be called for finite arguments */\r | |
79 | assert(Py_IS_FINITE(x));\r | |
80 | y = fmod(fabs(x), 2.0);\r | |
81 | n = (int)round(2.0*y);\r | |
82 | assert(0 <= n && n <= 4);\r | |
83 | switch (n) {\r | |
84 | case 0:\r | |
85 | r = sin(pi*y);\r | |
86 | break;\r | |
87 | case 1:\r | |
88 | r = cos(pi*(y-0.5));\r | |
89 | break;\r | |
90 | case 2:\r | |
91 | /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give\r | |
92 | -0.0 instead of 0.0 when y == 1.0. */\r | |
93 | r = sin(pi*(1.0-y));\r | |
94 | break;\r | |
95 | case 3:\r | |
96 | r = -cos(pi*(y-1.5));\r | |
97 | break;\r | |
98 | case 4:\r | |
99 | r = sin(pi*(y-2.0));\r | |
100 | break;\r | |
101 | default:\r | |
102 | assert(0); /* should never get here */\r | |
103 | r = -1.23e200; /* silence gcc warning */\r | |
104 | }\r | |
105 | return copysign(1.0, x)*r;\r | |
106 | }\r | |
107 | \r | |
108 | /* Implementation of the real gamma function. In extensive but non-exhaustive\r | |
109 | random tests, this function proved accurate to within <= 10 ulps across the\r | |
110 | entire float domain. Note that accuracy may depend on the quality of the\r | |
111 | system math functions, the pow function in particular. Special cases\r | |
112 | follow C99 annex F. The parameters and method are tailored to platforms\r | |
113 | whose double format is the IEEE 754 binary64 format.\r | |
114 | \r | |
115 | Method: for x > 0.0 we use the Lanczos approximation with parameters N=13\r | |
116 | and g=6.024680040776729583740234375; these parameters are amongst those\r | |
117 | used by the Boost library. Following Boost (again), we re-express the\r | |
118 | Lanczos sum as a rational function, and compute it that way. The\r | |
119 | coefficients below were computed independently using MPFR, and have been\r | |
120 | double-checked against the coefficients in the Boost source code.\r | |
121 | \r | |
122 | For x < 0.0 we use the reflection formula.\r | |
123 | \r | |
124 | There's one minor tweak that deserves explanation: Lanczos' formula for\r | |
125 | Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x\r | |
126 | values, x+g-0.5 can be represented exactly. However, in cases where it\r | |
127 | can't be represented exactly the small error in x+g-0.5 can be magnified\r | |
128 | significantly by the pow and exp calls, especially for large x. A cheap\r | |
129 | correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error\r | |
130 | involved in the computation of x+g-0.5 (that is, e = computed value of\r | |
131 | x+g-0.5 - exact value of x+g-0.5). Here's the proof:\r | |
132 | \r | |
133 | Correction factor\r | |
134 | -----------------\r | |
135 | Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754\r | |
136 | double, and e is tiny. Then:\r | |
137 | \r | |
138 | pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e)\r | |
139 | = pow(y, x-0.5)/exp(y) * C,\r | |
140 | \r | |
141 | where the correction_factor C is given by\r | |
142 | \r | |
143 | C = pow(1-e/y, x-0.5) * exp(e)\r | |
144 | \r | |
145 | Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so:\r | |
146 | \r | |
147 | C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y\r | |
148 | \r | |
149 | But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and\r | |
150 | \r | |
151 | pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y),\r | |
152 | \r | |
153 | Note that for accuracy, when computing r*C it's better to do\r | |
154 | \r | |
155 | r + e*g/y*r;\r | |
156 | \r | |
157 | than\r | |
158 | \r | |
159 | r * (1 + e*g/y);\r | |
160 | \r | |
161 | since the addition in the latter throws away most of the bits of\r | |
162 | information in e*g/y.\r | |
163 | */\r | |
164 | \r | |
165 | #define LANCZOS_N 13\r | |
166 | static const double lanczos_g = 6.024680040776729583740234375;\r | |
167 | static const double lanczos_g_minus_half = 5.524680040776729583740234375;\r | |
168 | static const double lanczos_num_coeffs[LANCZOS_N] = {\r | |
169 | 23531376880.410759688572007674451636754734846804940,\r | |
170 | 42919803642.649098768957899047001988850926355848959,\r | |
171 | 35711959237.355668049440185451547166705960488635843,\r | |
172 | 17921034426.037209699919755754458931112671403265390,\r | |
173 | 6039542586.3520280050642916443072979210699388420708,\r | |
174 | 1439720407.3117216736632230727949123939715485786772,\r | |
175 | 248874557.86205415651146038641322942321632125127801,\r | |
176 | 31426415.585400194380614231628318205362874684987640,\r | |
177 | 2876370.6289353724412254090516208496135991145378768,\r | |
178 | 186056.26539522349504029498971604569928220784236328,\r | |
179 | 8071.6720023658162106380029022722506138218516325024,\r | |
180 | 210.82427775157934587250973392071336271166969580291,\r | |
181 | 2.5066282746310002701649081771338373386264310793408\r | |
182 | };\r | |
183 | \r | |
184 | /* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */\r | |
185 | static const double lanczos_den_coeffs[LANCZOS_N] = {\r | |
186 | 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,\r | |
187 | 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};\r | |
188 | \r | |
189 | /* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */\r | |
190 | #define NGAMMA_INTEGRAL 23\r | |
191 | static const double gamma_integral[NGAMMA_INTEGRAL] = {\r | |
192 | 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,\r | |
193 | 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,\r | |
194 | 1307674368000.0, 20922789888000.0, 355687428096000.0,\r | |
195 | 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,\r | |
196 | 51090942171709440000.0, 1124000727777607680000.0,\r | |
197 | };\r | |
198 | \r | |
199 | /* Lanczos' sum L_g(x), for positive x */\r | |
200 | \r | |
201 | static double\r | |
202 | lanczos_sum(double x)\r | |
203 | {\r | |
204 | double num = 0.0, den = 0.0;\r | |
205 | int i;\r | |
206 | assert(x > 0.0);\r | |
207 | /* evaluate the rational function lanczos_sum(x). For large\r | |
208 | x, the obvious algorithm risks overflow, so we instead\r | |
209 | rescale the denominator and numerator of the rational\r | |
210 | function by x**(1-LANCZOS_N) and treat this as a\r | |
211 | rational function in 1/x. This also reduces the error for\r | |
212 | larger x values. The choice of cutoff point (5.0 below) is\r | |
213 | somewhat arbitrary; in tests, smaller cutoff values than\r | |
214 | this resulted in lower accuracy. */\r | |
215 | if (x < 5.0) {\r | |
216 | for (i = LANCZOS_N; --i >= 0; ) {\r | |
217 | num = num * x + lanczos_num_coeffs[i];\r | |
218 | den = den * x + lanczos_den_coeffs[i];\r | |
219 | }\r | |
220 | }\r | |
221 | else {\r | |
222 | for (i = 0; i < LANCZOS_N; i++) {\r | |
223 | num = num / x + lanczos_num_coeffs[i];\r | |
224 | den = den / x + lanczos_den_coeffs[i];\r | |
225 | }\r | |
226 | }\r | |
227 | return num/den;\r | |
228 | }\r | |
229 | \r | |
230 | static double\r | |
231 | m_tgamma(double x)\r | |
232 | {\r | |
233 | double absx, r, y, z, sqrtpow;\r | |
234 | \r | |
235 | /* special cases */\r | |
236 | if (!Py_IS_FINITE(x)) {\r | |
237 | if (Py_IS_NAN(x) || x > 0.0)\r | |
238 | return x; /* tgamma(nan) = nan, tgamma(inf) = inf */\r | |
239 | else {\r | |
240 | errno = EDOM;\r | |
241 | return Py_NAN; /* tgamma(-inf) = nan, invalid */\r | |
242 | }\r | |
243 | }\r | |
244 | if (x == 0.0) {\r | |
245 | errno = EDOM;\r | |
246 | return 1.0/x; /* tgamma(+-0.0) = +-inf, divide-by-zero */\r | |
247 | }\r | |
248 | \r | |
249 | /* integer arguments */\r | |
250 | if (x == floor(x)) {\r | |
251 | if (x < 0.0) {\r | |
252 | errno = EDOM; /* tgamma(n) = nan, invalid for */\r | |
253 | return Py_NAN; /* negative integers n */\r | |
254 | }\r | |
255 | if (x <= NGAMMA_INTEGRAL)\r | |
256 | return gamma_integral[(int)x - 1];\r | |
257 | }\r | |
258 | absx = fabs(x);\r | |
259 | \r | |
260 | /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */\r | |
261 | if (absx < 1e-20) {\r | |
262 | r = 1.0/x;\r | |
263 | if (Py_IS_INFINITY(r))\r | |
264 | errno = ERANGE;\r | |
265 | return r;\r | |
266 | }\r | |
267 | \r | |
268 | /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for\r | |
269 | x > 200, and underflows to +-0.0 for x < -200, not a negative\r | |
270 | integer. */\r | |
271 | if (absx > 200.0) {\r | |
272 | if (x < 0.0) {\r | |
273 | return 0.0/sinpi(x);\r | |
274 | }\r | |
275 | else {\r | |
276 | errno = ERANGE;\r | |
277 | return Py_HUGE_VAL;\r | |
278 | }\r | |
279 | }\r | |
280 | \r | |
281 | y = absx + lanczos_g_minus_half;\r | |
282 | /* compute error in sum */\r | |
283 | if (absx > lanczos_g_minus_half) {\r | |
284 | /* note: the correction can be foiled by an optimizing\r | |
285 | compiler that (incorrectly) thinks that an expression like\r | |
286 | a + b - a - b can be optimized to 0.0. This shouldn't\r | |
287 | happen in a standards-conforming compiler. */\r | |
288 | double q = y - absx;\r | |
289 | z = q - lanczos_g_minus_half;\r | |
290 | }\r | |
291 | else {\r | |
292 | double q = y - lanczos_g_minus_half;\r | |
293 | z = q - absx;\r | |
294 | }\r | |
295 | z = z * lanczos_g / y;\r | |
296 | if (x < 0.0) {\r | |
297 | r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx);\r | |
298 | r -= z * r;\r | |
299 | if (absx < 140.0) {\r | |
300 | r /= pow(y, absx - 0.5);\r | |
301 | }\r | |
302 | else {\r | |
303 | sqrtpow = pow(y, absx / 2.0 - 0.25);\r | |
304 | r /= sqrtpow;\r | |
305 | r /= sqrtpow;\r | |
306 | }\r | |
307 | }\r | |
308 | else {\r | |
309 | r = lanczos_sum(absx) / exp(y);\r | |
310 | r += z * r;\r | |
311 | if (absx < 140.0) {\r | |
312 | r *= pow(y, absx - 0.5);\r | |
313 | }\r | |
314 | else {\r | |
315 | sqrtpow = pow(y, absx / 2.0 - 0.25);\r | |
316 | r *= sqrtpow;\r | |
317 | r *= sqrtpow;\r | |
318 | }\r | |
319 | }\r | |
320 | if (Py_IS_INFINITY(r))\r | |
321 | errno = ERANGE;\r | |
322 | return r;\r | |
323 | }\r | |
324 | \r | |
325 | /*\r | |
326 | lgamma: natural log of the absolute value of the Gamma function.\r | |
327 | For large arguments, Lanczos' formula works extremely well here.\r | |
328 | */\r | |
329 | \r | |
330 | static double\r | |
331 | m_lgamma(double x)\r | |
332 | {\r | |
333 | double r, absx;\r | |
334 | \r | |
335 | /* special cases */\r | |
336 | if (!Py_IS_FINITE(x)) {\r | |
337 | if (Py_IS_NAN(x))\r | |
338 | return x; /* lgamma(nan) = nan */\r | |
339 | else\r | |
340 | return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */\r | |
341 | }\r | |
342 | \r | |
343 | /* integer arguments */\r | |
344 | if (x == floor(x) && x <= 2.0) {\r | |
345 | if (x <= 0.0) {\r | |
346 | errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */\r | |
347 | return Py_HUGE_VAL; /* integers n <= 0 */\r | |
348 | }\r | |
349 | else {\r | |
350 | return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */\r | |
351 | }\r | |
352 | }\r | |
353 | \r | |
354 | absx = fabs(x);\r | |
355 | /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */\r | |
356 | if (absx < 1e-20)\r | |
357 | return -log(absx);\r | |
358 | \r | |
359 | /* Lanczos' formula */\r | |
360 | if (x > 0.0) {\r | |
361 | /* we could save a fraction of a ulp in accuracy by having a\r | |
362 | second set of numerator coefficients for lanczos_sum that\r | |
363 | absorbed the exp(-lanczos_g) term, and throwing out the\r | |
364 | lanczos_g subtraction below; it's probably not worth it. */\r | |
365 | r = log(lanczos_sum(x)) - lanczos_g +\r | |
366 | (x-0.5)*(log(x+lanczos_g-0.5)-1);\r | |
367 | }\r | |
368 | else {\r | |
369 | r = log(pi) - log(fabs(sinpi(absx))) - log(absx) -\r | |
370 | (log(lanczos_sum(absx)) - lanczos_g +\r | |
371 | (absx-0.5)*(log(absx+lanczos_g-0.5)-1));\r | |
372 | }\r | |
373 | if (Py_IS_INFINITY(r))\r | |
374 | errno = ERANGE;\r | |
375 | return r;\r | |
376 | }\r | |
377 | \r | |
378 | /*\r | |
379 | Implementations of the error function erf(x) and the complementary error\r | |
380 | function erfc(x).\r | |
381 | \r | |
382 | Method: following 'Numerical Recipes' by Flannery, Press et. al. (2nd ed.,\r | |
383 | Cambridge University Press), we use a series approximation for erf for\r | |
384 | small x, and a continued fraction approximation for erfc(x) for larger x;\r | |
385 | combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x),\r | |
386 | this gives us erf(x) and erfc(x) for all x.\r | |
387 | \r | |
388 | The series expansion used is:\r | |
389 | \r | |
390 | erf(x) = x*exp(-x*x)/sqrt(pi) * [\r | |
391 | 2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...]\r | |
392 | \r | |
393 | The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k).\r | |
394 | This series converges well for smallish x, but slowly for larger x.\r | |
395 | \r | |
396 | The continued fraction expansion used is:\r | |
397 | \r | |
398 | erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - )\r | |
399 | 3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...]\r | |
400 | \r | |
401 | after the first term, the general term has the form:\r | |
402 | \r | |
403 | k*(k-0.5)/(2*k+0.5 + x**2 - ...).\r | |
404 | \r | |
405 | This expansion converges fast for larger x, but convergence becomes\r | |
406 | infinitely slow as x approaches 0.0. The (somewhat naive) continued\r | |
407 | fraction evaluation algorithm used below also risks overflow for large x;\r | |
408 | but for large x, erfc(x) == 0.0 to within machine precision. (For\r | |
409 | example, erfc(30.0) is approximately 2.56e-393).\r | |
410 | \r | |
411 | Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and\r | |
412 | continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) <\r | |
413 | ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the\r | |
414 | numbers of terms to use for the relevant expansions. */\r | |
415 | \r | |
416 | #define ERF_SERIES_CUTOFF 1.5\r | |
417 | #define ERF_SERIES_TERMS 25\r | |
418 | #define ERFC_CONTFRAC_CUTOFF 30.0\r | |
419 | #define ERFC_CONTFRAC_TERMS 50\r | |
420 | \r | |
421 | /*\r | |
422 | Error function, via power series.\r | |
423 | \r | |
424 | Given a finite float x, return an approximation to erf(x).\r | |
425 | Converges reasonably fast for small x.\r | |
426 | */\r | |
427 | \r | |
428 | static double\r | |
429 | m_erf_series(double x)\r | |
430 | {\r | |
431 | double x2, acc, fk, result;\r | |
432 | int i, saved_errno;\r | |
433 | \r | |
434 | x2 = x * x;\r | |
435 | acc = 0.0;\r | |
436 | fk = (double)ERF_SERIES_TERMS + 0.5;\r | |
437 | for (i = 0; i < ERF_SERIES_TERMS; i++) {\r | |
438 | acc = 2.0 + x2 * acc / fk;\r | |
439 | fk -= 1.0;\r | |
440 | }\r | |
441 | /* Make sure the exp call doesn't affect errno;\r | |
442 | see m_erfc_contfrac for more. */\r | |
443 | saved_errno = errno;\r | |
444 | result = acc * x * exp(-x2) / sqrtpi;\r | |
445 | errno = saved_errno;\r | |
446 | return result;\r | |
447 | }\r | |
448 | \r | |
449 | /*\r | |
450 | Complementary error function, via continued fraction expansion.\r | |
451 | \r | |
452 | Given a positive float x, return an approximation to erfc(x). Converges\r | |
453 | reasonably fast for x large (say, x > 2.0), and should be safe from\r | |
454 | overflow if x and nterms are not too large. On an IEEE 754 machine, with x\r | |
455 | <= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller\r | |
456 | than the smallest representable nonzero float. */\r | |
457 | \r | |
458 | static double\r | |
459 | m_erfc_contfrac(double x)\r | |
460 | {\r | |
461 | double x2, a, da, p, p_last, q, q_last, b, result;\r | |
462 | int i, saved_errno;\r | |
463 | \r | |
464 | if (x >= ERFC_CONTFRAC_CUTOFF)\r | |
465 | return 0.0;\r | |
466 | \r | |
467 | x2 = x*x;\r | |
468 | a = 0.0;\r | |
469 | da = 0.5;\r | |
470 | p = 1.0; p_last = 0.0;\r | |
471 | q = da + x2; q_last = 1.0;\r | |
472 | for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) {\r | |
473 | double temp;\r | |
474 | a += da;\r | |
475 | da += 2.0;\r | |
476 | b = da + x2;\r | |
477 | temp = p; p = b*p - a*p_last; p_last = temp;\r | |
478 | temp = q; q = b*q - a*q_last; q_last = temp;\r | |
479 | }\r | |
480 | /* Issue #8986: On some platforms, exp sets errno on underflow to zero;\r | |
481 | save the current errno value so that we can restore it later. */\r | |
482 | saved_errno = errno;\r | |
483 | result = p / q * x * exp(-x2) / sqrtpi;\r | |
484 | errno = saved_errno;\r | |
485 | return result;\r | |
486 | }\r | |
487 | \r | |
488 | /* Error function erf(x), for general x */\r | |
489 | \r | |
490 | static double\r | |
491 | m_erf(double x)\r | |
492 | {\r | |
493 | double absx, cf;\r | |
494 | \r | |
495 | if (Py_IS_NAN(x))\r | |
496 | return x;\r | |
497 | absx = fabs(x);\r | |
498 | if (absx < ERF_SERIES_CUTOFF)\r | |
499 | return m_erf_series(x);\r | |
500 | else {\r | |
501 | cf = m_erfc_contfrac(absx);\r | |
502 | return x > 0.0 ? 1.0 - cf : cf - 1.0;\r | |
503 | }\r | |
504 | }\r | |
505 | \r | |
506 | /* Complementary error function erfc(x), for general x. */\r | |
507 | \r | |
508 | static double\r | |
509 | m_erfc(double x)\r | |
510 | {\r | |
511 | double absx, cf;\r | |
512 | \r | |
513 | if (Py_IS_NAN(x))\r | |
514 | return x;\r | |
515 | absx = fabs(x);\r | |
516 | if (absx < ERF_SERIES_CUTOFF)\r | |
517 | return 1.0 - m_erf_series(x);\r | |
518 | else {\r | |
519 | cf = m_erfc_contfrac(absx);\r | |
520 | return x > 0.0 ? cf : 2.0 - cf;\r | |
521 | }\r | |
522 | }\r | |
523 | \r | |
524 | /*\r | |
525 | wrapper for atan2 that deals directly with special cases before\r | |
526 | delegating to the platform libm for the remaining cases. This\r | |
527 | is necessary to get consistent behaviour across platforms.\r | |
528 | Windows, FreeBSD and alpha Tru64 are amongst platforms that don't\r | |
529 | always follow C99.\r | |
530 | */\r | |
531 | \r | |
532 | static double\r | |
533 | m_atan2(double y, double x)\r | |
534 | {\r | |
535 | if (Py_IS_NAN(x) || Py_IS_NAN(y))\r | |
536 | return Py_NAN;\r | |
537 | if (Py_IS_INFINITY(y)) {\r | |
538 | if (Py_IS_INFINITY(x)) {\r | |
539 | if (copysign(1., x) == 1.)\r | |
540 | /* atan2(+-inf, +inf) == +-pi/4 */\r | |
541 | return copysign(0.25*Py_MATH_PI, y);\r | |
542 | else\r | |
543 | /* atan2(+-inf, -inf) == +-pi*3/4 */\r | |
544 | return copysign(0.75*Py_MATH_PI, y);\r | |
545 | }\r | |
546 | /* atan2(+-inf, x) == +-pi/2 for finite x */\r | |
547 | return copysign(0.5*Py_MATH_PI, y);\r | |
548 | }\r | |
549 | if (Py_IS_INFINITY(x) || y == 0.) {\r | |
550 | if (copysign(1., x) == 1.)\r | |
551 | /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */\r | |
552 | return copysign(0., y);\r | |
553 | else\r | |
554 | /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */\r | |
555 | return copysign(Py_MATH_PI, y);\r | |
556 | }\r | |
557 | return atan2(y, x);\r | |
558 | }\r | |
559 | \r | |
560 | /*\r | |
561 | Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),\r | |
562 | log(-ve), log(NaN). Here are wrappers for log and log10 that deal with\r | |
563 | special values directly, passing positive non-special values through to\r | |
564 | the system log/log10.\r | |
565 | */\r | |
566 | \r | |
567 | static double\r | |
568 | m_log(double x)\r | |
569 | {\r | |
570 | if (Py_IS_FINITE(x)) {\r | |
571 | if (x > 0.0)\r | |
572 | return log(x);\r | |
573 | errno = EDOM;\r | |
574 | if (x == 0.0)\r | |
575 | return -Py_HUGE_VAL; /* log(0) = -inf */\r | |
576 | else\r | |
577 | return Py_NAN; /* log(-ve) = nan */\r | |
578 | }\r | |
579 | else if (Py_IS_NAN(x))\r | |
580 | return x; /* log(nan) = nan */\r | |
581 | else if (x > 0.0)\r | |
582 | return x; /* log(inf) = inf */\r | |
583 | else {\r | |
584 | errno = EDOM;\r | |
585 | return Py_NAN; /* log(-inf) = nan */\r | |
586 | }\r | |
587 | }\r | |
588 | \r | |
589 | static double\r | |
590 | m_log10(double x)\r | |
591 | {\r | |
592 | if (Py_IS_FINITE(x)) {\r | |
593 | if (x > 0.0)\r | |
594 | return log10(x);\r | |
595 | errno = EDOM;\r | |
596 | if (x == 0.0)\r | |
597 | return -Py_HUGE_VAL; /* log10(0) = -inf */\r | |
598 | else\r | |
599 | return Py_NAN; /* log10(-ve) = nan */\r | |
600 | }\r | |
601 | else if (Py_IS_NAN(x))\r | |
602 | return x; /* log10(nan) = nan */\r | |
603 | else if (x > 0.0)\r | |
604 | return x; /* log10(inf) = inf */\r | |
605 | else {\r | |
606 | errno = EDOM;\r | |
607 | return Py_NAN; /* log10(-inf) = nan */\r | |
608 | }\r | |
609 | }\r | |
610 | \r | |
611 | \r | |
612 | /* Call is_error when errno != 0, and where x is the result libm\r | |
613 | * returned. is_error will usually set up an exception and return\r | |
614 | * true (1), but may return false (0) without setting up an exception.\r | |
615 | */\r | |
616 | static int\r | |
617 | is_error(double x)\r | |
618 | {\r | |
619 | int result = 1; /* presumption of guilt */\r | |
620 | assert(errno); /* non-zero errno is a precondition for calling */\r | |
621 | if (errno == EDOM)\r | |
622 | PyErr_SetString(PyExc_ValueError, "math domain error");\r | |
623 | \r | |
624 | else if (errno == ERANGE) {\r | |
625 | /* ANSI C generally requires libm functions to set ERANGE\r | |
626 | * on overflow, but also generally *allows* them to set\r | |
627 | * ERANGE on underflow too. There's no consistency about\r | |
628 | * the latter across platforms.\r | |
629 | * Alas, C99 never requires that errno be set.\r | |
630 | * Here we suppress the underflow errors (libm functions\r | |
631 | * should return a zero on underflow, and +- HUGE_VAL on\r | |
632 | * overflow, so testing the result for zero suffices to\r | |
633 | * distinguish the cases).\r | |
634 | *\r | |
635 | * On some platforms (Ubuntu/ia64) it seems that errno can be\r | |
636 | * set to ERANGE for subnormal results that do *not* underflow\r | |
637 | * to zero. So to be safe, we'll ignore ERANGE whenever the\r | |
638 | * function result is less than one in absolute value.\r | |
639 | */\r | |
640 | if (fabs(x) < 1.0)\r | |
641 | result = 0;\r | |
642 | else\r | |
643 | PyErr_SetString(PyExc_OverflowError,\r | |
644 | "math range error");\r | |
645 | }\r | |
646 | else\r | |
647 | /* Unexpected math error */\r | |
648 | PyErr_SetFromErrno(PyExc_ValueError);\r | |
649 | return result;\r | |
650 | }\r | |
651 | \r | |
652 | /*\r | |
653 | math_1 is used to wrap a libm function f that takes a double\r | |
654 | arguments and returns a double.\r | |
655 | \r | |
656 | The error reporting follows these rules, which are designed to do\r | |
657 | the right thing on C89/C99 platforms and IEEE 754/non IEEE 754\r | |
658 | platforms.\r | |
659 | \r | |
660 | - a NaN result from non-NaN inputs causes ValueError to be raised\r | |
661 | - an infinite result from finite inputs causes OverflowError to be\r | |
662 | raised if can_overflow is 1, or raises ValueError if can_overflow\r | |
663 | is 0.\r | |
664 | - if the result is finite and errno == EDOM then ValueError is\r | |
665 | raised\r | |
666 | - if the result is finite and nonzero and errno == ERANGE then\r | |
667 | OverflowError is raised\r | |
668 | \r | |
669 | The last rule is used to catch overflow on platforms which follow\r | |
670 | C89 but for which HUGE_VAL is not an infinity.\r | |
671 | \r | |
672 | For the majority of one-argument functions these rules are enough\r | |
673 | to ensure that Python's functions behave as specified in 'Annex F'\r | |
674 | of the C99 standard, with the 'invalid' and 'divide-by-zero'\r | |
675 | floating-point exceptions mapping to Python's ValueError and the\r | |
676 | 'overflow' floating-point exception mapping to OverflowError.\r | |
677 | math_1 only works for functions that don't have singularities *and*\r | |
678 | the possibility of overflow; fortunately, that covers everything we\r | |
679 | care about right now.\r | |
680 | */\r | |
681 | \r | |
682 | static PyObject *\r | |
683 | math_1(PyObject *arg, double (*func) (double), int can_overflow)\r | |
684 | {\r | |
685 | double x, r;\r | |
686 | x = PyFloat_AsDouble(arg);\r | |
687 | if (x == -1.0 && PyErr_Occurred())\r | |
688 | return NULL;\r | |
689 | errno = 0;\r | |
690 | PyFPE_START_PROTECT("in math_1", return 0);\r | |
691 | r = (*func)(x);\r | |
692 | PyFPE_END_PROTECT(r);\r | |
693 | if (Py_IS_NAN(r)) {\r | |
694 | if (!Py_IS_NAN(x))\r | |
695 | errno = EDOM;\r | |
696 | else\r | |
697 | errno = 0;\r | |
698 | }\r | |
699 | else if (Py_IS_INFINITY(r)) {\r | |
700 | if (Py_IS_FINITE(x))\r | |
701 | errno = can_overflow ? ERANGE : EDOM;\r | |
702 | else\r | |
703 | errno = 0;\r | |
704 | }\r | |
705 | if (errno && is_error(r))\r | |
706 | return NULL;\r | |
707 | else\r | |
708 | return PyFloat_FromDouble(r);\r | |
709 | }\r | |
710 | \r | |
711 | /* variant of math_1, to be used when the function being wrapped is known to\r | |
712 | set errno properly (that is, errno = EDOM for invalid or divide-by-zero,\r | |
713 | errno = ERANGE for overflow). */\r | |
714 | \r | |
715 | static PyObject *\r | |
716 | math_1a(PyObject *arg, double (*func) (double))\r | |
717 | {\r | |
718 | double x, r;\r | |
719 | x = PyFloat_AsDouble(arg);\r | |
720 | if (x == -1.0 && PyErr_Occurred())\r | |
721 | return NULL;\r | |
722 | errno = 0;\r | |
723 | PyFPE_START_PROTECT("in math_1a", return 0);\r | |
724 | r = (*func)(x);\r | |
725 | PyFPE_END_PROTECT(r);\r | |
726 | if (errno && is_error(r))\r | |
727 | return NULL;\r | |
728 | return PyFloat_FromDouble(r);\r | |
729 | }\r | |
730 | \r | |
731 | /*\r | |
732 | math_2 is used to wrap a libm function f that takes two double\r | |
733 | arguments and returns a double.\r | |
734 | \r | |
735 | The error reporting follows these rules, which are designed to do\r | |
736 | the right thing on C89/C99 platforms and IEEE 754/non IEEE 754\r | |
737 | platforms.\r | |
738 | \r | |
739 | - a NaN result from non-NaN inputs causes ValueError to be raised\r | |
740 | - an infinite result from finite inputs causes OverflowError to be\r | |
741 | raised.\r | |
742 | - if the result is finite and errno == EDOM then ValueError is\r | |
743 | raised\r | |
744 | - if the result is finite and nonzero and errno == ERANGE then\r | |
745 | OverflowError is raised\r | |
746 | \r | |
747 | The last rule is used to catch overflow on platforms which follow\r | |
748 | C89 but for which HUGE_VAL is not an infinity.\r | |
749 | \r | |
750 | For most two-argument functions (copysign, fmod, hypot, atan2)\r | |
751 | these rules are enough to ensure that Python's functions behave as\r | |
752 | specified in 'Annex F' of the C99 standard, with the 'invalid' and\r | |
753 | 'divide-by-zero' floating-point exceptions mapping to Python's\r | |
754 | ValueError and the 'overflow' floating-point exception mapping to\r | |
755 | OverflowError.\r | |
756 | */\r | |
757 | \r | |
758 | static PyObject *\r | |
759 | math_2(PyObject *args, double (*func) (double, double), char *funcname)\r | |
760 | {\r | |
761 | PyObject *ox, *oy;\r | |
762 | double x, y, r;\r | |
763 | if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))\r | |
764 | return NULL;\r | |
765 | x = PyFloat_AsDouble(ox);\r | |
766 | y = PyFloat_AsDouble(oy);\r | |
767 | if ((x == -1.0 || y == -1.0) && PyErr_Occurred())\r | |
768 | return NULL;\r | |
769 | errno = 0;\r | |
770 | PyFPE_START_PROTECT("in math_2", return 0);\r | |
771 | r = (*func)(x, y);\r | |
772 | PyFPE_END_PROTECT(r);\r | |
773 | if (Py_IS_NAN(r)) {\r | |
774 | if (!Py_IS_NAN(x) && !Py_IS_NAN(y))\r | |
775 | errno = EDOM;\r | |
776 | else\r | |
777 | errno = 0;\r | |
778 | }\r | |
779 | else if (Py_IS_INFINITY(r)) {\r | |
780 | if (Py_IS_FINITE(x) && Py_IS_FINITE(y))\r | |
781 | errno = ERANGE;\r | |
782 | else\r | |
783 | errno = 0;\r | |
784 | }\r | |
785 | if (errno && is_error(r))\r | |
786 | return NULL;\r | |
787 | else\r | |
788 | return PyFloat_FromDouble(r);\r | |
789 | }\r | |
790 | \r | |
791 | #define FUNC1(funcname, func, can_overflow, docstring) \\r | |
792 | static PyObject * math_##funcname(PyObject *self, PyObject *args) { \\r | |
793 | return math_1(args, func, can_overflow); \\r | |
794 | }\\r | |
795 | PyDoc_STRVAR(math_##funcname##_doc, docstring);\r | |
796 | \r | |
797 | #define FUNC1A(funcname, func, docstring) \\r | |
798 | static PyObject * math_##funcname(PyObject *self, PyObject *args) { \\r | |
799 | return math_1a(args, func); \\r | |
800 | }\\r | |
801 | PyDoc_STRVAR(math_##funcname##_doc, docstring);\r | |
802 | \r | |
803 | #define FUNC2(funcname, func, docstring) \\r | |
804 | static PyObject * math_##funcname(PyObject *self, PyObject *args) { \\r | |
805 | return math_2(args, func, #funcname); \\r | |
806 | }\\r | |
807 | PyDoc_STRVAR(math_##funcname##_doc, docstring);\r | |
808 | \r | |
809 | FUNC1(acos, acos, 0,\r | |
810 | "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")\r | |
811 | FUNC1(acosh, m_acosh, 0,\r | |
812 | "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")\r | |
813 | FUNC1(asin, asin, 0,\r | |
814 | "asin(x)\n\nReturn the arc sine (measured in radians) of x.")\r | |
815 | FUNC1(asinh, m_asinh, 0,\r | |
816 | "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")\r | |
817 | FUNC1(atan, atan, 0,\r | |
818 | "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")\r | |
819 | FUNC2(atan2, m_atan2,\r | |
820 | "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"\r | |
821 | "Unlike atan(y/x), the signs of both x and y are considered.")\r | |
822 | FUNC1(atanh, m_atanh, 0,\r | |
823 | "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")\r | |
824 | FUNC1(ceil, ceil, 0,\r | |
825 | "ceil(x)\n\nReturn the ceiling of x as a float.\n"\r | |
826 | "This is the smallest integral value >= x.")\r | |
827 | FUNC2(copysign, copysign,\r | |
828 | "copysign(x, y)\n\nReturn x with the sign of y.")\r | |
829 | FUNC1(cos, cos, 0,\r | |
830 | "cos(x)\n\nReturn the cosine of x (measured in radians).")\r | |
831 | FUNC1(cosh, cosh, 1,\r | |
832 | "cosh(x)\n\nReturn the hyperbolic cosine of x.")\r | |
833 | FUNC1A(erf, m_erf,\r | |
834 | "erf(x)\n\nError function at x.")\r | |
835 | FUNC1A(erfc, m_erfc,\r | |
836 | "erfc(x)\n\nComplementary error function at x.")\r | |
837 | FUNC1(exp, exp, 1,\r | |
838 | "exp(x)\n\nReturn e raised to the power of x.")\r | |
839 | FUNC1(expm1, m_expm1, 1,\r | |
840 | "expm1(x)\n\nReturn exp(x)-1.\n"\r | |
841 | "This function avoids the loss of precision involved in the direct "\r | |
842 | "evaluation of exp(x)-1 for small x.")\r | |
843 | FUNC1(fabs, fabs, 0,\r | |
844 | "fabs(x)\n\nReturn the absolute value of the float x.")\r | |
845 | FUNC1(floor, floor, 0,\r | |
846 | "floor(x)\n\nReturn the floor of x as a float.\n"\r | |
847 | "This is the largest integral value <= x.")\r | |
848 | FUNC1A(gamma, m_tgamma,\r | |
849 | "gamma(x)\n\nGamma function at x.")\r | |
850 | FUNC1A(lgamma, m_lgamma,\r | |
851 | "lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.")\r | |
852 | FUNC1(log1p, m_log1p, 1,\r | |
853 | "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n"\r | |
854 | "The result is computed in a way which is accurate for x near zero.")\r | |
855 | FUNC1(sin, sin, 0,\r | |
856 | "sin(x)\n\nReturn the sine of x (measured in radians).")\r | |
857 | FUNC1(sinh, sinh, 1,\r | |
858 | "sinh(x)\n\nReturn the hyperbolic sine of x.")\r | |
859 | FUNC1(sqrt, sqrt, 0,\r | |
860 | "sqrt(x)\n\nReturn the square root of x.")\r | |
861 | FUNC1(tan, tan, 0,\r | |
862 | "tan(x)\n\nReturn the tangent of x (measured in radians).")\r | |
863 | FUNC1(tanh, tanh, 0,\r | |
864 | "tanh(x)\n\nReturn the hyperbolic tangent of x.")\r | |
865 | \r | |
866 | /* Precision summation function as msum() by Raymond Hettinger in\r | |
867 | <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,\r | |
868 | enhanced with the exact partials sum and roundoff from Mark\r | |
869 | Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.\r | |
870 | See those links for more details, proofs and other references.\r | |
871 | \r | |
872 | Note 1: IEEE 754R floating point semantics are assumed,\r | |
873 | but the current implementation does not re-establish special\r | |
874 | value semantics across iterations (i.e. handling -Inf + Inf).\r | |
875 | \r | |
876 | Note 2: No provision is made for intermediate overflow handling;\r | |
877 | therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while\r | |
878 | sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the\r | |
879 | overflow of the first partial sum.\r | |
880 | \r | |
881 | Note 3: The intermediate values lo, yr, and hi are declared volatile so\r | |
882 | aggressive compilers won't algebraically reduce lo to always be exactly 0.0.\r | |
883 | Also, the volatile declaration forces the values to be stored in memory as\r | |
884 | regular doubles instead of extended long precision (80-bit) values. This\r | |
885 | prevents double rounding because any addition or subtraction of two doubles\r | |
886 | can be resolved exactly into double-sized hi and lo values. As long as the\r | |
887 | hi value gets forced into a double before yr and lo are computed, the extra\r | |
888 | bits in downstream extended precision operations (x87 for example) will be\r | |
889 | exactly zero and therefore can be losslessly stored back into a double,\r | |
890 | thereby preventing double rounding.\r | |
891 | \r | |
892 | Note 4: A similar implementation is in Modules/cmathmodule.c.\r | |
893 | Be sure to update both when making changes.\r | |
894 | \r | |
895 | Note 5: The signature of math.fsum() differs from __builtin__.sum()\r | |
896 | because the start argument doesn't make sense in the context of\r | |
897 | accurate summation. Since the partials table is collapsed before\r | |
898 | returning a result, sum(seq2, start=sum(seq1)) may not equal the\r | |
899 | accurate result returned by sum(itertools.chain(seq1, seq2)).\r | |
900 | */\r | |
901 | \r | |
902 | #define NUM_PARTIALS 32 /* initial partials array size, on stack */\r | |
903 | \r | |
904 | /* Extend the partials array p[] by doubling its size. */\r | |
905 | static int /* non-zero on error */\r | |
906 | _fsum_realloc(double **p_ptr, Py_ssize_t n,\r | |
907 | double *ps, Py_ssize_t *m_ptr)\r | |
908 | {\r | |
909 | void *v = NULL;\r | |
910 | Py_ssize_t m = *m_ptr;\r | |
911 | \r | |
912 | m += m; /* double */\r | |
913 | if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {\r | |
914 | double *p = *p_ptr;\r | |
915 | if (p == ps) {\r | |
916 | v = PyMem_Malloc(sizeof(double) * m);\r | |
917 | if (v != NULL)\r | |
918 | memcpy(v, ps, sizeof(double) * n);\r | |
919 | }\r | |
920 | else\r | |
921 | v = PyMem_Realloc(p, sizeof(double) * m);\r | |
922 | }\r | |
923 | if (v == NULL) { /* size overflow or no memory */\r | |
924 | PyErr_SetString(PyExc_MemoryError, "math.fsum partials");\r | |
925 | return 1;\r | |
926 | }\r | |
927 | *p_ptr = (double*) v;\r | |
928 | *m_ptr = m;\r | |
929 | return 0;\r | |
930 | }\r | |
931 | \r | |
932 | /* Full precision summation of a sequence of floats.\r | |
933 | \r | |
934 | def msum(iterable):\r | |
935 | partials = [] # sorted, non-overlapping partial sums\r | |
936 | for x in iterable:\r | |
937 | i = 0\r | |
938 | for y in partials:\r | |
939 | if abs(x) < abs(y):\r | |
940 | x, y = y, x\r | |
941 | hi = x + y\r | |
942 | lo = y - (hi - x)\r | |
943 | if lo:\r | |
944 | partials[i] = lo\r | |
945 | i += 1\r | |
946 | x = hi\r | |
947 | partials[i:] = [x]\r | |
948 | return sum_exact(partials)\r | |
949 | \r | |
950 | Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo\r | |
951 | are exactly equal to x+y. The inner loop applies hi/lo summation to each\r | |
952 | partial so that the list of partial sums remains exact.\r | |
953 | \r | |
954 | Sum_exact() adds the partial sums exactly and correctly rounds the final\r | |
955 | result (using the round-half-to-even rule). The items in partials remain\r | |
956 | non-zero, non-special, non-overlapping and strictly increasing in\r | |
957 | magnitude, but possibly not all having the same sign.\r | |
958 | \r | |
959 | Depends on IEEE 754 arithmetic guarantees and half-even rounding.\r | |
960 | */\r | |
961 | \r | |
962 | static PyObject*\r | |
963 | math_fsum(PyObject *self, PyObject *seq)\r | |
964 | {\r | |
965 | PyObject *item, *iter, *sum = NULL;\r | |
966 | Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;\r | |
967 | double x, y, t, ps[NUM_PARTIALS], *p = ps;\r | |
968 | double xsave, special_sum = 0.0, inf_sum = 0.0;\r | |
969 | volatile double hi, yr, lo;\r | |
970 | \r | |
971 | iter = PyObject_GetIter(seq);\r | |
972 | if (iter == NULL)\r | |
973 | return NULL;\r | |
974 | \r | |
975 | PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)\r | |
976 | \r | |
977 | for(;;) { /* for x in iterable */\r | |
978 | assert(0 <= n && n <= m);\r | |
979 | assert((m == NUM_PARTIALS && p == ps) ||\r | |
980 | (m > NUM_PARTIALS && p != NULL));\r | |
981 | \r | |
982 | item = PyIter_Next(iter);\r | |
983 | if (item == NULL) {\r | |
984 | if (PyErr_Occurred())\r | |
985 | goto _fsum_error;\r | |
986 | break;\r | |
987 | }\r | |
988 | x = PyFloat_AsDouble(item);\r | |
989 | Py_DECREF(item);\r | |
990 | if (PyErr_Occurred())\r | |
991 | goto _fsum_error;\r | |
992 | \r | |
993 | xsave = x;\r | |
994 | for (i = j = 0; j < n; j++) { /* for y in partials */\r | |
995 | y = p[j];\r | |
996 | if (fabs(x) < fabs(y)) {\r | |
997 | t = x; x = y; y = t;\r | |
998 | }\r | |
999 | hi = x + y;\r | |
1000 | yr = hi - x;\r | |
1001 | lo = y - yr;\r | |
1002 | if (lo != 0.0)\r | |
1003 | p[i++] = lo;\r | |
1004 | x = hi;\r | |
1005 | }\r | |
1006 | \r | |
1007 | n = i; /* ps[i:] = [x] */\r | |
1008 | if (x != 0.0) {\r | |
1009 | if (! Py_IS_FINITE(x)) {\r | |
1010 | /* a nonfinite x could arise either as\r | |
1011 | a result of intermediate overflow, or\r | |
1012 | as a result of a nan or inf in the\r | |
1013 | summands */\r | |
1014 | if (Py_IS_FINITE(xsave)) {\r | |
1015 | PyErr_SetString(PyExc_OverflowError,\r | |
1016 | "intermediate overflow in fsum");\r | |
1017 | goto _fsum_error;\r | |
1018 | }\r | |
1019 | if (Py_IS_INFINITY(xsave))\r | |
1020 | inf_sum += xsave;\r | |
1021 | special_sum += xsave;\r | |
1022 | /* reset partials */\r | |
1023 | n = 0;\r | |
1024 | }\r | |
1025 | else if (n >= m && _fsum_realloc(&p, n, ps, &m))\r | |
1026 | goto _fsum_error;\r | |
1027 | else\r | |
1028 | p[n++] = x;\r | |
1029 | }\r | |
1030 | }\r | |
1031 | \r | |
1032 | if (special_sum != 0.0) {\r | |
1033 | if (Py_IS_NAN(inf_sum))\r | |
1034 | PyErr_SetString(PyExc_ValueError,\r | |
1035 | "-inf + inf in fsum");\r | |
1036 | else\r | |
1037 | sum = PyFloat_FromDouble(special_sum);\r | |
1038 | goto _fsum_error;\r | |
1039 | }\r | |
1040 | \r | |
1041 | hi = 0.0;\r | |
1042 | if (n > 0) {\r | |
1043 | hi = p[--n];\r | |
1044 | /* sum_exact(ps, hi) from the top, stop when the sum becomes\r | |
1045 | inexact. */\r | |
1046 | while (n > 0) {\r | |
1047 | x = hi;\r | |
1048 | y = p[--n];\r | |
1049 | assert(fabs(y) < fabs(x));\r | |
1050 | hi = x + y;\r | |
1051 | yr = hi - x;\r | |
1052 | lo = y - yr;\r | |
1053 | if (lo != 0.0)\r | |
1054 | break;\r | |
1055 | }\r | |
1056 | /* Make half-even rounding work across multiple partials.\r | |
1057 | Needed so that sum([1e-16, 1, 1e16]) will round-up the last\r | |
1058 | digit to two instead of down to zero (the 1e-16 makes the 1\r | |
1059 | slightly closer to two). With a potential 1 ULP rounding\r | |
1060 | error fixed-up, math.fsum() can guarantee commutativity. */\r | |
1061 | if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||\r | |
1062 | (lo > 0.0 && p[n-1] > 0.0))) {\r | |
1063 | y = lo * 2.0;\r | |
1064 | x = hi + y;\r | |
1065 | yr = x - hi;\r | |
1066 | if (y == yr)\r | |
1067 | hi = x;\r | |
1068 | }\r | |
1069 | }\r | |
1070 | sum = PyFloat_FromDouble(hi);\r | |
1071 | \r | |
1072 | _fsum_error:\r | |
1073 | PyFPE_END_PROTECT(hi)\r | |
1074 | Py_DECREF(iter);\r | |
1075 | if (p != ps)\r | |
1076 | PyMem_Free(p);\r | |
1077 | return sum;\r | |
1078 | }\r | |
1079 | \r | |
1080 | #undef NUM_PARTIALS\r | |
1081 | \r | |
1082 | PyDoc_STRVAR(math_fsum_doc,\r | |
1083 | "fsum(iterable)\n\n\\r | |
1084 | Return an accurate floating point sum of values in the iterable.\n\\r | |
1085 | Assumes IEEE-754 floating point arithmetic.");\r | |
1086 | \r | |
1087 | static PyObject *\r | |
1088 | math_factorial(PyObject *self, PyObject *arg)\r | |
1089 | {\r | |
1090 | long i, x;\r | |
1091 | PyObject *result, *iobj, *newresult;\r | |
1092 | \r | |
1093 | if (PyFloat_Check(arg)) {\r | |
1094 | PyObject *lx;\r | |
1095 | double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);\r | |
1096 | if (!(Py_IS_FINITE(dx) && dx == floor(dx))) {\r | |
1097 | PyErr_SetString(PyExc_ValueError,\r | |
1098 | "factorial() only accepts integral values");\r | |
1099 | return NULL;\r | |
1100 | }\r | |
1101 | lx = PyLong_FromDouble(dx);\r | |
1102 | if (lx == NULL)\r | |
1103 | return NULL;\r | |
1104 | x = PyLong_AsLong(lx);\r | |
1105 | Py_DECREF(lx);\r | |
1106 | }\r | |
1107 | else\r | |
1108 | x = PyInt_AsLong(arg);\r | |
1109 | \r | |
1110 | if (x == -1 && PyErr_Occurred())\r | |
1111 | return NULL;\r | |
1112 | if (x < 0) {\r | |
1113 | PyErr_SetString(PyExc_ValueError,\r | |
1114 | "factorial() not defined for negative values");\r | |
1115 | return NULL;\r | |
1116 | }\r | |
1117 | \r | |
1118 | result = (PyObject *)PyInt_FromLong(1);\r | |
1119 | if (result == NULL)\r | |
1120 | return NULL;\r | |
1121 | for (i=1 ; i<=x ; i++) {\r | |
1122 | iobj = (PyObject *)PyInt_FromLong(i);\r | |
1123 | if (iobj == NULL)\r | |
1124 | goto error;\r | |
1125 | newresult = PyNumber_Multiply(result, iobj);\r | |
1126 | Py_DECREF(iobj);\r | |
1127 | if (newresult == NULL)\r | |
1128 | goto error;\r | |
1129 | Py_DECREF(result);\r | |
1130 | result = newresult;\r | |
1131 | }\r | |
1132 | return result;\r | |
1133 | \r | |
1134 | error:\r | |
1135 | Py_DECREF(result);\r | |
1136 | return NULL;\r | |
1137 | }\r | |
1138 | \r | |
1139 | PyDoc_STRVAR(math_factorial_doc,\r | |
1140 | "factorial(x) -> Integral\n"\r | |
1141 | "\n"\r | |
1142 | "Find x!. Raise a ValueError if x is negative or non-integral.");\r | |
1143 | \r | |
1144 | static PyObject *\r | |
1145 | math_trunc(PyObject *self, PyObject *number)\r | |
1146 | {\r | |
1147 | return PyObject_CallMethod(number, "__trunc__", NULL);\r | |
1148 | }\r | |
1149 | \r | |
1150 | PyDoc_STRVAR(math_trunc_doc,\r | |
1151 | "trunc(x:Real) -> Integral\n"\r | |
1152 | "\n"\r | |
1153 | "Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");\r | |
1154 | \r | |
1155 | static PyObject *\r | |
1156 | math_frexp(PyObject *self, PyObject *arg)\r | |
1157 | {\r | |
1158 | int i;\r | |
1159 | double x = PyFloat_AsDouble(arg);\r | |
1160 | if (x == -1.0 && PyErr_Occurred())\r | |
1161 | return NULL;\r | |
1162 | /* deal with special cases directly, to sidestep platform\r | |
1163 | differences */\r | |
1164 | if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {\r | |
1165 | i = 0;\r | |
1166 | }\r | |
1167 | else {\r | |
1168 | PyFPE_START_PROTECT("in math_frexp", return 0);\r | |
1169 | x = frexp(x, &i);\r | |
1170 | PyFPE_END_PROTECT(x);\r | |
1171 | }\r | |
1172 | return Py_BuildValue("(di)", x, i);\r | |
1173 | }\r | |
1174 | \r | |
1175 | PyDoc_STRVAR(math_frexp_doc,\r | |
1176 | "frexp(x)\n"\r | |
1177 | "\n"\r | |
1178 | "Return the mantissa and exponent of x, as pair (m, e).\n"\r | |
1179 | "m is a float and e is an int, such that x = m * 2.**e.\n"\r | |
1180 | "If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0.");\r | |
1181 | \r | |
1182 | static PyObject *\r | |
1183 | math_ldexp(PyObject *self, PyObject *args)\r | |
1184 | {\r | |
1185 | double x, r;\r | |
1186 | PyObject *oexp;\r | |
1187 | long exp;\r | |
1188 | int overflow;\r | |
1189 | if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp))\r | |
1190 | return NULL;\r | |
1191 | \r | |
1192 | if (PyLong_Check(oexp) || PyInt_Check(oexp)) {\r | |
1193 | /* on overflow, replace exponent with either LONG_MAX\r | |
1194 | or LONG_MIN, depending on the sign. */\r | |
1195 | exp = PyLong_AsLongAndOverflow(oexp, &overflow);\r | |
1196 | if (exp == -1 && PyErr_Occurred())\r | |
1197 | return NULL;\r | |
1198 | if (overflow)\r | |
1199 | exp = overflow < 0 ? LONG_MIN : LONG_MAX;\r | |
1200 | }\r | |
1201 | else {\r | |
1202 | PyErr_SetString(PyExc_TypeError,\r | |
1203 | "Expected an int or long as second argument "\r | |
1204 | "to ldexp.");\r | |
1205 | return NULL;\r | |
1206 | }\r | |
1207 | \r | |
1208 | if (x == 0. || !Py_IS_FINITE(x)) {\r | |
1209 | /* NaNs, zeros and infinities are returned unchanged */\r | |
1210 | r = x;\r | |
1211 | errno = 0;\r | |
1212 | } else if (exp > INT_MAX) {\r | |
1213 | /* overflow */\r | |
1214 | r = copysign(Py_HUGE_VAL, x);\r | |
1215 | errno = ERANGE;\r | |
1216 | } else if (exp < INT_MIN) {\r | |
1217 | /* underflow to +-0 */\r | |
1218 | r = copysign(0., x);\r | |
1219 | errno = 0;\r | |
1220 | } else {\r | |
1221 | errno = 0;\r | |
1222 | PyFPE_START_PROTECT("in math_ldexp", return 0);\r | |
1223 | r = ldexp(x, (int)exp);\r | |
1224 | PyFPE_END_PROTECT(r);\r | |
1225 | if (Py_IS_INFINITY(r))\r | |
1226 | errno = ERANGE;\r | |
1227 | }\r | |
1228 | \r | |
1229 | if (errno && is_error(r))\r | |
1230 | return NULL;\r | |
1231 | return PyFloat_FromDouble(r);\r | |
1232 | }\r | |
1233 | \r | |
1234 | PyDoc_STRVAR(math_ldexp_doc,\r | |
1235 | "ldexp(x, i)\n\n\\r | |
1236 | Return x * (2**i).");\r | |
1237 | \r | |
1238 | static PyObject *\r | |
1239 | math_modf(PyObject *self, PyObject *arg)\r | |
1240 | {\r | |
1241 | double y, x = PyFloat_AsDouble(arg);\r | |
1242 | if (x == -1.0 && PyErr_Occurred())\r | |
1243 | return NULL;\r | |
1244 | /* some platforms don't do the right thing for NaNs and\r | |
1245 | infinities, so we take care of special cases directly. */\r | |
1246 | if (!Py_IS_FINITE(x)) {\r | |
1247 | if (Py_IS_INFINITY(x))\r | |
1248 | return Py_BuildValue("(dd)", copysign(0., x), x);\r | |
1249 | else if (Py_IS_NAN(x))\r | |
1250 | return Py_BuildValue("(dd)", x, x);\r | |
1251 | }\r | |
1252 | \r | |
1253 | errno = 0;\r | |
1254 | PyFPE_START_PROTECT("in math_modf", return 0);\r | |
1255 | x = modf(x, &y);\r | |
1256 | PyFPE_END_PROTECT(x);\r | |
1257 | return Py_BuildValue("(dd)", x, y);\r | |
1258 | }\r | |
1259 | \r | |
1260 | PyDoc_STRVAR(math_modf_doc,\r | |
1261 | "modf(x)\n"\r | |
1262 | "\n"\r | |
1263 | "Return the fractional and integer parts of x. Both results carry the sign\n"\r | |
1264 | "of x and are floats.");\r | |
1265 | \r | |
1266 | /* A decent logarithm is easy to compute even for huge longs, but libm can't\r | |
1267 | do that by itself -- loghelper can. func is log or log10, and name is\r | |
1268 | "log" or "log10". Note that overflow of the result isn't possible: a long\r | |
1269 | can contain no more than INT_MAX * SHIFT bits, so has value certainly less\r | |
1270 | than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is\r | |
1271 | small enough to fit in an IEEE single. log and log10 are even smaller.\r | |
1272 | However, intermediate overflow is possible for a long if the number of bits\r | |
1273 | in that long is larger than PY_SSIZE_T_MAX. */\r | |
1274 | \r | |
1275 | static PyObject*\r | |
1276 | loghelper(PyObject* arg, double (*func)(double), char *funcname)\r | |
1277 | {\r | |
1278 | /* If it is long, do it ourselves. */\r | |
1279 | if (PyLong_Check(arg)) {\r | |
1280 | double x;\r | |
1281 | Py_ssize_t e;\r | |
1282 | x = _PyLong_Frexp((PyLongObject *)arg, &e);\r | |
1283 | if (x == -1.0 && PyErr_Occurred())\r | |
1284 | return NULL;\r | |
1285 | if (x <= 0.0) {\r | |
1286 | PyErr_SetString(PyExc_ValueError,\r | |
1287 | "math domain error");\r | |
1288 | return NULL;\r | |
1289 | }\r | |
1290 | /* Special case for log(1), to make sure we get an\r | |
1291 | exact result there. */\r | |
1292 | if (e == 1 && x == 0.5)\r | |
1293 | return PyFloat_FromDouble(0.0);\r | |
1294 | /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */\r | |
1295 | x = func(x) + func(2.0) * e;\r | |
1296 | return PyFloat_FromDouble(x);\r | |
1297 | }\r | |
1298 | \r | |
1299 | /* Else let libm handle it by itself. */\r | |
1300 | return math_1(arg, func, 0);\r | |
1301 | }\r | |
1302 | \r | |
1303 | static PyObject *\r | |
1304 | math_log(PyObject *self, PyObject *args)\r | |
1305 | {\r | |
1306 | PyObject *arg;\r | |
1307 | PyObject *base = NULL;\r | |
1308 | PyObject *num, *den;\r | |
1309 | PyObject *ans;\r | |
1310 | \r | |
1311 | if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))\r | |
1312 | return NULL;\r | |
1313 | \r | |
1314 | num = loghelper(arg, m_log, "log");\r | |
1315 | if (num == NULL || base == NULL)\r | |
1316 | return num;\r | |
1317 | \r | |
1318 | den = loghelper(base, m_log, "log");\r | |
1319 | if (den == NULL) {\r | |
1320 | Py_DECREF(num);\r | |
1321 | return NULL;\r | |
1322 | }\r | |
1323 | \r | |
1324 | ans = PyNumber_Divide(num, den);\r | |
1325 | Py_DECREF(num);\r | |
1326 | Py_DECREF(den);\r | |
1327 | return ans;\r | |
1328 | }\r | |
1329 | \r | |
1330 | PyDoc_STRVAR(math_log_doc,\r | |
1331 | "log(x[, base])\n\n\\r | |
1332 | Return the logarithm of x to the given base.\n\\r | |
1333 | If the base not specified, returns the natural logarithm (base e) of x.");\r | |
1334 | \r | |
1335 | static PyObject *\r | |
1336 | math_log10(PyObject *self, PyObject *arg)\r | |
1337 | {\r | |
1338 | return loghelper(arg, m_log10, "log10");\r | |
1339 | }\r | |
1340 | \r | |
1341 | PyDoc_STRVAR(math_log10_doc,\r | |
1342 | "log10(x)\n\nReturn the base 10 logarithm of x.");\r | |
1343 | \r | |
1344 | static PyObject *\r | |
1345 | math_fmod(PyObject *self, PyObject *args)\r | |
1346 | {\r | |
1347 | PyObject *ox, *oy;\r | |
1348 | double r, x, y;\r | |
1349 | if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))\r | |
1350 | return NULL;\r | |
1351 | x = PyFloat_AsDouble(ox);\r | |
1352 | y = PyFloat_AsDouble(oy);\r | |
1353 | if ((x == -1.0 || y == -1.0) && PyErr_Occurred())\r | |
1354 | return NULL;\r | |
1355 | /* fmod(x, +/-Inf) returns x for finite x. */\r | |
1356 | if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))\r | |
1357 | return PyFloat_FromDouble(x);\r | |
1358 | errno = 0;\r | |
1359 | PyFPE_START_PROTECT("in math_fmod", return 0);\r | |
1360 | r = fmod(x, y);\r | |
1361 | PyFPE_END_PROTECT(r);\r | |
1362 | if (Py_IS_NAN(r)) {\r | |
1363 | if (!Py_IS_NAN(x) && !Py_IS_NAN(y))\r | |
1364 | errno = EDOM;\r | |
1365 | else\r | |
1366 | errno = 0;\r | |
1367 | }\r | |
1368 | if (errno && is_error(r))\r | |
1369 | return NULL;\r | |
1370 | else\r | |
1371 | return PyFloat_FromDouble(r);\r | |
1372 | }\r | |
1373 | \r | |
1374 | PyDoc_STRVAR(math_fmod_doc,\r | |
1375 | "fmod(x, y)\n\nReturn fmod(x, y), according to platform C."\r | |
1376 | " x % y may differ.");\r | |
1377 | \r | |
1378 | static PyObject *\r | |
1379 | math_hypot(PyObject *self, PyObject *args)\r | |
1380 | {\r | |
1381 | PyObject *ox, *oy;\r | |
1382 | double r, x, y;\r | |
1383 | if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))\r | |
1384 | return NULL;\r | |
1385 | x = PyFloat_AsDouble(ox);\r | |
1386 | y = PyFloat_AsDouble(oy);\r | |
1387 | if ((x == -1.0 || y == -1.0) && PyErr_Occurred())\r | |
1388 | return NULL;\r | |
1389 | /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */\r | |
1390 | if (Py_IS_INFINITY(x))\r | |
1391 | return PyFloat_FromDouble(fabs(x));\r | |
1392 | if (Py_IS_INFINITY(y))\r | |
1393 | return PyFloat_FromDouble(fabs(y));\r | |
1394 | errno = 0;\r | |
1395 | PyFPE_START_PROTECT("in math_hypot", return 0);\r | |
1396 | r = hypot(x, y);\r | |
1397 | PyFPE_END_PROTECT(r);\r | |
1398 | if (Py_IS_NAN(r)) {\r | |
1399 | if (!Py_IS_NAN(x) && !Py_IS_NAN(y))\r | |
1400 | errno = EDOM;\r | |
1401 | else\r | |
1402 | errno = 0;\r | |
1403 | }\r | |
1404 | else if (Py_IS_INFINITY(r)) {\r | |
1405 | if (Py_IS_FINITE(x) && Py_IS_FINITE(y))\r | |
1406 | errno = ERANGE;\r | |
1407 | else\r | |
1408 | errno = 0;\r | |
1409 | }\r | |
1410 | if (errno && is_error(r))\r | |
1411 | return NULL;\r | |
1412 | else\r | |
1413 | return PyFloat_FromDouble(r);\r | |
1414 | }\r | |
1415 | \r | |
1416 | PyDoc_STRVAR(math_hypot_doc,\r | |
1417 | "hypot(x, y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");\r | |
1418 | \r | |
1419 | /* pow can't use math_2, but needs its own wrapper: the problem is\r | |
1420 | that an infinite result can arise either as a result of overflow\r | |
1421 | (in which case OverflowError should be raised) or as a result of\r | |
1422 | e.g. 0.**-5. (for which ValueError needs to be raised.)\r | |
1423 | */\r | |
1424 | \r | |
1425 | static PyObject *\r | |
1426 | math_pow(PyObject *self, PyObject *args)\r | |
1427 | {\r | |
1428 | PyObject *ox, *oy;\r | |
1429 | double r, x, y;\r | |
1430 | int odd_y;\r | |
1431 | \r | |
1432 | if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))\r | |
1433 | return NULL;\r | |
1434 | x = PyFloat_AsDouble(ox);\r | |
1435 | y = PyFloat_AsDouble(oy);\r | |
1436 | if ((x == -1.0 || y == -1.0) && PyErr_Occurred())\r | |
1437 | return NULL;\r | |
1438 | \r | |
1439 | /* deal directly with IEEE specials, to cope with problems on various\r | |
1440 | platforms whose semantics don't exactly match C99 */\r | |
1441 | r = 0.; /* silence compiler warning */\r | |
1442 | if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {\r | |
1443 | errno = 0;\r | |
1444 | if (Py_IS_NAN(x))\r | |
1445 | r = y == 0. ? 1. : x; /* NaN**0 = 1 */\r | |
1446 | else if (Py_IS_NAN(y))\r | |
1447 | r = x == 1. ? 1. : y; /* 1**NaN = 1 */\r | |
1448 | else if (Py_IS_INFINITY(x)) {\r | |
1449 | odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;\r | |
1450 | if (y > 0.)\r | |
1451 | r = odd_y ? x : fabs(x);\r | |
1452 | else if (y == 0.)\r | |
1453 | r = 1.;\r | |
1454 | else /* y < 0. */\r | |
1455 | r = odd_y ? copysign(0., x) : 0.;\r | |
1456 | }\r | |
1457 | else if (Py_IS_INFINITY(y)) {\r | |
1458 | if (fabs(x) == 1.0)\r | |
1459 | r = 1.;\r | |
1460 | else if (y > 0. && fabs(x) > 1.0)\r | |
1461 | r = y;\r | |
1462 | else if (y < 0. && fabs(x) < 1.0) {\r | |
1463 | r = -y; /* result is +inf */\r | |
1464 | if (x == 0.) /* 0**-inf: divide-by-zero */\r | |
1465 | errno = EDOM;\r | |
1466 | }\r | |
1467 | else\r | |
1468 | r = 0.;\r | |
1469 | }\r | |
1470 | }\r | |
1471 | else {\r | |
1472 | /* let libm handle finite**finite */\r | |
1473 | errno = 0;\r | |
1474 | PyFPE_START_PROTECT("in math_pow", return 0);\r | |
1475 | r = pow(x, y);\r | |
1476 | PyFPE_END_PROTECT(r);\r | |
1477 | /* a NaN result should arise only from (-ve)**(finite\r | |
1478 | non-integer); in this case we want to raise ValueError. */\r | |
1479 | if (!Py_IS_FINITE(r)) {\r | |
1480 | if (Py_IS_NAN(r)) {\r | |
1481 | errno = EDOM;\r | |
1482 | }\r | |
1483 | /*\r | |
1484 | an infinite result here arises either from:\r | |
1485 | (A) (+/-0.)**negative (-> divide-by-zero)\r | |
1486 | (B) overflow of x**y with x and y finite\r | |
1487 | */\r | |
1488 | else if (Py_IS_INFINITY(r)) {\r | |
1489 | if (x == 0.)\r | |
1490 | errno = EDOM;\r | |
1491 | else\r | |
1492 | errno = ERANGE;\r | |
1493 | }\r | |
1494 | }\r | |
1495 | }\r | |
1496 | \r | |
1497 | if (errno && is_error(r))\r | |
1498 | return NULL;\r | |
1499 | else\r | |
1500 | return PyFloat_FromDouble(r);\r | |
1501 | }\r | |
1502 | \r | |
1503 | PyDoc_STRVAR(math_pow_doc,\r | |
1504 | "pow(x, y)\n\nReturn x**y (x to the power of y).");\r | |
1505 | \r | |
1506 | static const double degToRad = Py_MATH_PI / 180.0;\r | |
1507 | static const double radToDeg = 180.0 / Py_MATH_PI;\r | |
1508 | \r | |
1509 | static PyObject *\r | |
1510 | math_degrees(PyObject *self, PyObject *arg)\r | |
1511 | {\r | |
1512 | double x = PyFloat_AsDouble(arg);\r | |
1513 | if (x == -1.0 && PyErr_Occurred())\r | |
1514 | return NULL;\r | |
1515 | return PyFloat_FromDouble(x * radToDeg);\r | |
1516 | }\r | |
1517 | \r | |
1518 | PyDoc_STRVAR(math_degrees_doc,\r | |
1519 | "degrees(x)\n\n\\r | |
1520 | Convert angle x from radians to degrees.");\r | |
1521 | \r | |
1522 | static PyObject *\r | |
1523 | math_radians(PyObject *self, PyObject *arg)\r | |
1524 | {\r | |
1525 | double x = PyFloat_AsDouble(arg);\r | |
1526 | if (x == -1.0 && PyErr_Occurred())\r | |
1527 | return NULL;\r | |
1528 | return PyFloat_FromDouble(x * degToRad);\r | |
1529 | }\r | |
1530 | \r | |
1531 | PyDoc_STRVAR(math_radians_doc,\r | |
1532 | "radians(x)\n\n\\r | |
1533 | Convert angle x from degrees to radians.");\r | |
1534 | \r | |
1535 | static PyObject *\r | |
1536 | math_isnan(PyObject *self, PyObject *arg)\r | |
1537 | {\r | |
1538 | double x = PyFloat_AsDouble(arg);\r | |
1539 | if (x == -1.0 && PyErr_Occurred())\r | |
1540 | return NULL;\r | |
1541 | return PyBool_FromLong((long)Py_IS_NAN(x));\r | |
1542 | }\r | |
1543 | \r | |
1544 | PyDoc_STRVAR(math_isnan_doc,\r | |
1545 | "isnan(x) -> bool\n\n\\r | |
1546 | Check if float x is not a number (NaN).");\r | |
1547 | \r | |
1548 | static PyObject *\r | |
1549 | math_isinf(PyObject *self, PyObject *arg)\r | |
1550 | {\r | |
1551 | double x = PyFloat_AsDouble(arg);\r | |
1552 | if (x == -1.0 && PyErr_Occurred())\r | |
1553 | return NULL;\r | |
1554 | return PyBool_FromLong((long)Py_IS_INFINITY(x));\r | |
1555 | }\r | |
1556 | \r | |
1557 | PyDoc_STRVAR(math_isinf_doc,\r | |
1558 | "isinf(x) -> bool\n\n\\r | |
1559 | Check if float x is infinite (positive or negative).");\r | |
1560 | \r | |
1561 | static PyMethodDef math_methods[] = {\r | |
1562 | {"acos", math_acos, METH_O, math_acos_doc},\r | |
1563 | {"acosh", math_acosh, METH_O, math_acosh_doc},\r | |
1564 | {"asin", math_asin, METH_O, math_asin_doc},\r | |
1565 | {"asinh", math_asinh, METH_O, math_asinh_doc},\r | |
1566 | {"atan", math_atan, METH_O, math_atan_doc},\r | |
1567 | {"atan2", math_atan2, METH_VARARGS, math_atan2_doc},\r | |
1568 | {"atanh", math_atanh, METH_O, math_atanh_doc},\r | |
1569 | {"ceil", math_ceil, METH_O, math_ceil_doc},\r | |
1570 | {"copysign", math_copysign, METH_VARARGS, math_copysign_doc},\r | |
1571 | {"cos", math_cos, METH_O, math_cos_doc},\r | |
1572 | {"cosh", math_cosh, METH_O, math_cosh_doc},\r | |
1573 | {"degrees", math_degrees, METH_O, math_degrees_doc},\r | |
1574 | {"erf", math_erf, METH_O, math_erf_doc},\r | |
1575 | {"erfc", math_erfc, METH_O, math_erfc_doc},\r | |
1576 | {"exp", math_exp, METH_O, math_exp_doc},\r | |
1577 | {"expm1", math_expm1, METH_O, math_expm1_doc},\r | |
1578 | {"fabs", math_fabs, METH_O, math_fabs_doc},\r | |
1579 | {"factorial", math_factorial, METH_O, math_factorial_doc},\r | |
1580 | {"floor", math_floor, METH_O, math_floor_doc},\r | |
1581 | {"fmod", math_fmod, METH_VARARGS, math_fmod_doc},\r | |
1582 | {"frexp", math_frexp, METH_O, math_frexp_doc},\r | |
1583 | {"fsum", math_fsum, METH_O, math_fsum_doc},\r | |
1584 | {"gamma", math_gamma, METH_O, math_gamma_doc},\r | |
1585 | {"hypot", math_hypot, METH_VARARGS, math_hypot_doc},\r | |
1586 | {"isinf", math_isinf, METH_O, math_isinf_doc},\r | |
1587 | {"isnan", math_isnan, METH_O, math_isnan_doc},\r | |
1588 | {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},\r | |
1589 | {"lgamma", math_lgamma, METH_O, math_lgamma_doc},\r | |
1590 | {"log", math_log, METH_VARARGS, math_log_doc},\r | |
1591 | {"log1p", math_log1p, METH_O, math_log1p_doc},\r | |
1592 | {"log10", math_log10, METH_O, math_log10_doc},\r | |
1593 | {"modf", math_modf, METH_O, math_modf_doc},\r | |
1594 | {"pow", math_pow, METH_VARARGS, math_pow_doc},\r | |
1595 | {"radians", math_radians, METH_O, math_radians_doc},\r | |
1596 | {"sin", math_sin, METH_O, math_sin_doc},\r | |
1597 | {"sinh", math_sinh, METH_O, math_sinh_doc},\r | |
1598 | {"sqrt", math_sqrt, METH_O, math_sqrt_doc},\r | |
1599 | {"tan", math_tan, METH_O, math_tan_doc},\r | |
1600 | {"tanh", math_tanh, METH_O, math_tanh_doc},\r | |
1601 | {"trunc", math_trunc, METH_O, math_trunc_doc},\r | |
1602 | {NULL, NULL} /* sentinel */\r | |
1603 | };\r | |
1604 | \r | |
1605 | \r | |
1606 | PyDoc_STRVAR(module_doc,\r | |
1607 | "This module is always available. It provides access to the\n"\r | |
1608 | "mathematical functions defined by the C standard.");\r | |
1609 | \r | |
1610 | PyMODINIT_FUNC\r | |
1611 | initmath(void)\r | |
1612 | {\r | |
1613 | PyObject *m;\r | |
1614 | \r | |
1615 | m = Py_InitModule3("math", math_methods, module_doc);\r | |
1616 | if (m == NULL)\r | |
1617 | goto finally;\r | |
1618 | \r | |
1619 | PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));\r | |
1620 | PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));\r | |
1621 | \r | |
1622 | finally:\r | |
1623 | return;\r | |
1624 | }\r |