blob: 29c32a30d84e54d9b41ef3bcaeb378b7ed7b6cbe [file] [log] [blame]
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001/* Math module -- standard C math library functions, pi and e */
2
Christian Heimes53876d92008-04-19 00:31:39 +00003/* Here are some comments from Tim Peters, extracted from the
4 discussion attached to http://bugs.python.org/issue1640. They
5 describe the general aims of the math module with respect to
6 special values, IEEE-754 floating-point exceptions, and Python
7 exceptions.
8
9These are the "spirit of 754" rules:
10
111. If the mathematical result is a real number, but of magnitude too
12large to approximate by a machine float, overflow is signaled and the
13result is an infinity (with the appropriate sign).
14
152. If the mathematical result is a real number, but of magnitude too
16small to approximate by a machine float, underflow is signaled and the
17result is a zero (with the appropriate sign).
18
193. At a singularity (a value x such that the limit of f(y) as y
20approaches x exists and is an infinity), "divide by zero" is signaled
21and the result is an infinity (with the appropriate sign). This is
22complicated a little by that the left-side and right-side limits may
23not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0
24from the positive or negative directions. In that specific case, the
25sign of the zero determines the result of 1/0.
26
274. At a point where a function has no defined result in the extended
28reals (i.e., the reals plus an infinity or two), invalid operation is
29signaled and a NaN is returned.
30
31And these are what Python has historically /tried/ to do (but not
32always successfully, as platform libm behavior varies a lot):
33
34For #1, raise OverflowError.
35
36For #2, return a zero (with the appropriate sign if that happens by
37accident ;-)).
38
39For #3 and #4, raise ValueError. It may have made sense to raise
40Python's ZeroDivisionError in #3, but historically that's only been
41raised for division by zero and mod by zero.
42
43*/
44
45/*
46 In general, on an IEEE-754 platform the aim is to follow the C99
47 standard, including Annex 'F', whenever possible. Where the
48 standard recommends raising the 'divide-by-zero' or 'invalid'
49 floating-point exceptions, Python should raise a ValueError. Where
50 the standard recommends raising 'overflow', Python should raise an
51 OverflowError. In all other circumstances a value should be
52 returned.
53 */
54
Barry Warsaw8b43b191996-12-09 22:32:36 +000055#include "Python.h"
Mark Dickinson664b5112009-12-16 20:23:42 +000056#include "_math.h"
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000057
Christian Heimes969fe572008-01-25 11:23:10 +000058#ifdef _OSF_SOURCE
59/* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
60extern double copysign(double, double);
61#endif
62
Mark Dickinson12c4bdb2009-09-28 19:21:11 +000063/*
64 sin(pi*x), giving accurate results for all finite x (especially x
65 integral or close to an integer). This is here for use in the
66 reflection formula for the gamma function. It conforms to IEEE
67 754-2008 for finite arguments, but not for infinities or nans.
68*/
Tim Petersa40c7932001-09-05 22:36:56 +000069
Mark Dickinson12c4bdb2009-09-28 19:21:11 +000070static const double pi = 3.141592653589793238462643383279502884197;
Mark Dickinson45f992a2009-12-19 11:20:49 +000071static const double sqrtpi = 1.772453850905516027298167483341145182798;
Mark Dickinson9c91eb82010-07-07 16:17:31 +000072static const double logpi = 1.144729885849400174143427351353058711647;
Mark Dickinson12c4bdb2009-09-28 19:21:11 +000073
74static double
75sinpi(double x)
76{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +000077 double y, r;
78 int n;
79 /* this function should only ever be called for finite arguments */
80 assert(Py_IS_FINITE(x));
81 y = fmod(fabs(x), 2.0);
82 n = (int)round(2.0*y);
83 assert(0 <= n && n <= 4);
84 switch (n) {
85 case 0:
86 r = sin(pi*y);
87 break;
88 case 1:
89 r = cos(pi*(y-0.5));
90 break;
91 case 2:
92 /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give
93 -0.0 instead of 0.0 when y == 1.0. */
94 r = sin(pi*(1.0-y));
95 break;
96 case 3:
97 r = -cos(pi*(y-1.5));
98 break;
99 case 4:
100 r = sin(pi*(y-2.0));
101 break;
102 default:
103 assert(0); /* should never get here */
104 r = -1.23e200; /* silence gcc warning */
105 }
106 return copysign(1.0, x)*r;
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000107}
108
109/* Implementation of the real gamma function. In extensive but non-exhaustive
110 random tests, this function proved accurate to within <= 10 ulps across the
111 entire float domain. Note that accuracy may depend on the quality of the
112 system math functions, the pow function in particular. Special cases
113 follow C99 annex F. The parameters and method are tailored to platforms
114 whose double format is the IEEE 754 binary64 format.
115
116 Method: for x > 0.0 we use the Lanczos approximation with parameters N=13
117 and g=6.024680040776729583740234375; these parameters are amongst those
118 used by the Boost library. Following Boost (again), we re-express the
119 Lanczos sum as a rational function, and compute it that way. The
120 coefficients below were computed independently using MPFR, and have been
121 double-checked against the coefficients in the Boost source code.
122
123 For x < 0.0 we use the reflection formula.
124
125 There's one minor tweak that deserves explanation: Lanczos' formula for
126 Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x
127 values, x+g-0.5 can be represented exactly. However, in cases where it
128 can't be represented exactly the small error in x+g-0.5 can be magnified
129 significantly by the pow and exp calls, especially for large x. A cheap
130 correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error
131 involved in the computation of x+g-0.5 (that is, e = computed value of
132 x+g-0.5 - exact value of x+g-0.5). Here's the proof:
133
134 Correction factor
135 -----------------
136 Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754
137 double, and e is tiny. Then:
138
139 pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e)
140 = pow(y, x-0.5)/exp(y) * C,
141
142 where the correction_factor C is given by
143
144 C = pow(1-e/y, x-0.5) * exp(e)
145
146 Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so:
147
148 C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y
149
150 But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and
151
152 pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y),
153
154 Note that for accuracy, when computing r*C it's better to do
155
156 r + e*g/y*r;
157
158 than
159
160 r * (1 + e*g/y);
161
162 since the addition in the latter throws away most of the bits of
163 information in e*g/y.
164*/
165
166#define LANCZOS_N 13
167static const double lanczos_g = 6.024680040776729583740234375;
168static const double lanczos_g_minus_half = 5.524680040776729583740234375;
169static const double lanczos_num_coeffs[LANCZOS_N] = {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000170 23531376880.410759688572007674451636754734846804940,
171 42919803642.649098768957899047001988850926355848959,
172 35711959237.355668049440185451547166705960488635843,
173 17921034426.037209699919755754458931112671403265390,
174 6039542586.3520280050642916443072979210699388420708,
175 1439720407.3117216736632230727949123939715485786772,
176 248874557.86205415651146038641322942321632125127801,
177 31426415.585400194380614231628318205362874684987640,
178 2876370.6289353724412254090516208496135991145378768,
179 186056.26539522349504029498971604569928220784236328,
180 8071.6720023658162106380029022722506138218516325024,
181 210.82427775157934587250973392071336271166969580291,
182 2.5066282746310002701649081771338373386264310793408
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000183};
184
185/* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */
186static const double lanczos_den_coeffs[LANCZOS_N] = {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000187 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
188 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000189
190/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */
191#define NGAMMA_INTEGRAL 23
192static const double gamma_integral[NGAMMA_INTEGRAL] = {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000193 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
194 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
195 1307674368000.0, 20922789888000.0, 355687428096000.0,
196 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
197 51090942171709440000.0, 1124000727777607680000.0,
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000198};
199
200/* Lanczos' sum L_g(x), for positive x */
201
202static double
203lanczos_sum(double x)
204{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000205 double num = 0.0, den = 0.0;
206 int i;
207 assert(x > 0.0);
208 /* evaluate the rational function lanczos_sum(x). For large
209 x, the obvious algorithm risks overflow, so we instead
210 rescale the denominator and numerator of the rational
211 function by x**(1-LANCZOS_N) and treat this as a
212 rational function in 1/x. This also reduces the error for
213 larger x values. The choice of cutoff point (5.0 below) is
214 somewhat arbitrary; in tests, smaller cutoff values than
215 this resulted in lower accuracy. */
216 if (x < 5.0) {
217 for (i = LANCZOS_N; --i >= 0; ) {
218 num = num * x + lanczos_num_coeffs[i];
219 den = den * x + lanczos_den_coeffs[i];
220 }
221 }
222 else {
223 for (i = 0; i < LANCZOS_N; i++) {
224 num = num / x + lanczos_num_coeffs[i];
225 den = den / x + lanczos_den_coeffs[i];
226 }
227 }
228 return num/den;
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000229}
230
231static double
232m_tgamma(double x)
233{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000234 double absx, r, y, z, sqrtpow;
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000235
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000236 /* special cases */
237 if (!Py_IS_FINITE(x)) {
238 if (Py_IS_NAN(x) || x > 0.0)
239 return x; /* tgamma(nan) = nan, tgamma(inf) = inf */
240 else {
241 errno = EDOM;
242 return Py_NAN; /* tgamma(-inf) = nan, invalid */
243 }
244 }
245 if (x == 0.0) {
246 errno = EDOM;
247 return 1.0/x; /* tgamma(+-0.0) = +-inf, divide-by-zero */
248 }
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000249
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000250 /* integer arguments */
251 if (x == floor(x)) {
252 if (x < 0.0) {
253 errno = EDOM; /* tgamma(n) = nan, invalid for */
254 return Py_NAN; /* negative integers n */
255 }
256 if (x <= NGAMMA_INTEGRAL)
257 return gamma_integral[(int)x - 1];
258 }
259 absx = fabs(x);
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000260
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000261 /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */
262 if (absx < 1e-20) {
263 r = 1.0/x;
264 if (Py_IS_INFINITY(r))
265 errno = ERANGE;
266 return r;
267 }
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000268
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000269 /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for
270 x > 200, and underflows to +-0.0 for x < -200, not a negative
271 integer. */
272 if (absx > 200.0) {
273 if (x < 0.0) {
274 return 0.0/sinpi(x);
275 }
276 else {
277 errno = ERANGE;
278 return Py_HUGE_VAL;
279 }
280 }
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000281
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000282 y = absx + lanczos_g_minus_half;
283 /* compute error in sum */
284 if (absx > lanczos_g_minus_half) {
285 /* note: the correction can be foiled by an optimizing
286 compiler that (incorrectly) thinks that an expression like
287 a + b - a - b can be optimized to 0.0. This shouldn't
288 happen in a standards-conforming compiler. */
289 double q = y - absx;
290 z = q - lanczos_g_minus_half;
291 }
292 else {
293 double q = y - lanczos_g_minus_half;
294 z = q - absx;
295 }
296 z = z * lanczos_g / y;
297 if (x < 0.0) {
298 r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx);
299 r -= z * r;
300 if (absx < 140.0) {
301 r /= pow(y, absx - 0.5);
302 }
303 else {
304 sqrtpow = pow(y, absx / 2.0 - 0.25);
305 r /= sqrtpow;
306 r /= sqrtpow;
307 }
308 }
309 else {
310 r = lanczos_sum(absx) / exp(y);
311 r += z * r;
312 if (absx < 140.0) {
313 r *= pow(y, absx - 0.5);
314 }
315 else {
316 sqrtpow = pow(y, absx / 2.0 - 0.25);
317 r *= sqrtpow;
318 r *= sqrtpow;
319 }
320 }
321 if (Py_IS_INFINITY(r))
322 errno = ERANGE;
323 return r;
Guido van Rossum8832b621991-12-16 15:44:24 +0000324}
325
Christian Heimes53876d92008-04-19 00:31:39 +0000326/*
Mark Dickinson05d2e082009-12-11 20:17:17 +0000327 lgamma: natural log of the absolute value of the Gamma function.
328 For large arguments, Lanczos' formula works extremely well here.
329*/
330
331static double
332m_lgamma(double x)
333{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000334 double r, absx;
Mark Dickinson05d2e082009-12-11 20:17:17 +0000335
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000336 /* special cases */
337 if (!Py_IS_FINITE(x)) {
338 if (Py_IS_NAN(x))
339 return x; /* lgamma(nan) = nan */
340 else
341 return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */
342 }
Mark Dickinson05d2e082009-12-11 20:17:17 +0000343
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000344 /* integer arguments */
345 if (x == floor(x) && x <= 2.0) {
346 if (x <= 0.0) {
347 errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */
348 return Py_HUGE_VAL; /* integers n <= 0 */
349 }
350 else {
351 return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */
352 }
353 }
Mark Dickinson05d2e082009-12-11 20:17:17 +0000354
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000355 absx = fabs(x);
356 /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */
357 if (absx < 1e-20)
358 return -log(absx);
Mark Dickinson05d2e082009-12-11 20:17:17 +0000359
Mark Dickinson9c91eb82010-07-07 16:17:31 +0000360 /* Lanczos' formula. We could save a fraction of a ulp in accuracy by
361 having a second set of numerator coefficients for lanczos_sum that
362 absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g
363 subtraction below; it's probably not worth it. */
364 r = log(lanczos_sum(absx)) - lanczos_g;
365 r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1);
366 if (x < 0.0)
367 /* Use reflection formula to get value for negative x. */
368 r = logpi - log(fabs(sinpi(absx))) - log(absx) - r;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000369 if (Py_IS_INFINITY(r))
370 errno = ERANGE;
371 return r;
Mark Dickinson05d2e082009-12-11 20:17:17 +0000372}
373
Mark Dickinson45f992a2009-12-19 11:20:49 +0000374/*
375 Implementations of the error function erf(x) and the complementary error
376 function erfc(x).
377
378 Method: following 'Numerical Recipes' by Flannery, Press et. al. (2nd ed.,
379 Cambridge University Press), we use a series approximation for erf for
380 small x, and a continued fraction approximation for erfc(x) for larger x;
381 combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x),
382 this gives us erf(x) and erfc(x) for all x.
383
384 The series expansion used is:
385
386 erf(x) = x*exp(-x*x)/sqrt(pi) * [
387 2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...]
388
389 The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k).
390 This series converges well for smallish x, but slowly for larger x.
391
392 The continued fraction expansion used is:
393
394 erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - )
395 3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...]
396
397 after the first term, the general term has the form:
398
399 k*(k-0.5)/(2*k+0.5 + x**2 - ...).
400
401 This expansion converges fast for larger x, but convergence becomes
402 infinitely slow as x approaches 0.0. The (somewhat naive) continued
403 fraction evaluation algorithm used below also risks overflow for large x;
404 but for large x, erfc(x) == 0.0 to within machine precision. (For
405 example, erfc(30.0) is approximately 2.56e-393).
406
407 Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and
408 continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) <
409 ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the
410 numbers of terms to use for the relevant expansions. */
411
412#define ERF_SERIES_CUTOFF 1.5
413#define ERF_SERIES_TERMS 25
414#define ERFC_CONTFRAC_CUTOFF 30.0
415#define ERFC_CONTFRAC_TERMS 50
416
417/*
418 Error function, via power series.
419
420 Given a finite float x, return an approximation to erf(x).
421 Converges reasonably fast for small x.
422*/
423
424static double
425m_erf_series(double x)
426{
Mark Dickinsonbcdf9da2010-06-13 10:52:38 +0000427 double x2, acc, fk, result;
428 int i, saved_errno;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000429
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000430 x2 = x * x;
431 acc = 0.0;
432 fk = (double)ERF_SERIES_TERMS + 0.5;
433 for (i = 0; i < ERF_SERIES_TERMS; i++) {
434 acc = 2.0 + x2 * acc / fk;
435 fk -= 1.0;
436 }
Mark Dickinsonbcdf9da2010-06-13 10:52:38 +0000437 /* Make sure the exp call doesn't affect errno;
438 see m_erfc_contfrac for more. */
439 saved_errno = errno;
440 result = acc * x * exp(-x2) / sqrtpi;
441 errno = saved_errno;
442 return result;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000443}
444
445/*
446 Complementary error function, via continued fraction expansion.
447
448 Given a positive float x, return an approximation to erfc(x). Converges
449 reasonably fast for x large (say, x > 2.0), and should be safe from
450 overflow if x and nterms are not too large. On an IEEE 754 machine, with x
451 <= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller
452 than the smallest representable nonzero float. */
453
454static double
455m_erfc_contfrac(double x)
456{
Mark Dickinsonbcdf9da2010-06-13 10:52:38 +0000457 double x2, a, da, p, p_last, q, q_last, b, result;
458 int i, saved_errno;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000459
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000460 if (x >= ERFC_CONTFRAC_CUTOFF)
461 return 0.0;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000462
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000463 x2 = x*x;
464 a = 0.0;
465 da = 0.5;
466 p = 1.0; p_last = 0.0;
467 q = da + x2; q_last = 1.0;
468 for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) {
469 double temp;
470 a += da;
471 da += 2.0;
472 b = da + x2;
473 temp = p; p = b*p - a*p_last; p_last = temp;
474 temp = q; q = b*q - a*q_last; q_last = temp;
475 }
Mark Dickinsonbcdf9da2010-06-13 10:52:38 +0000476 /* Issue #8986: On some platforms, exp sets errno on underflow to zero;
477 save the current errno value so that we can restore it later. */
478 saved_errno = errno;
479 result = p / q * x * exp(-x2) / sqrtpi;
480 errno = saved_errno;
481 return result;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000482}
483
484/* Error function erf(x), for general x */
485
486static double
487m_erf(double x)
488{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000489 double absx, cf;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000490
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000491 if (Py_IS_NAN(x))
492 return x;
493 absx = fabs(x);
494 if (absx < ERF_SERIES_CUTOFF)
495 return m_erf_series(x);
496 else {
497 cf = m_erfc_contfrac(absx);
498 return x > 0.0 ? 1.0 - cf : cf - 1.0;
499 }
Mark Dickinson45f992a2009-12-19 11:20:49 +0000500}
501
502/* Complementary error function erfc(x), for general x. */
503
504static double
505m_erfc(double x)
506{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000507 double absx, cf;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000508
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000509 if (Py_IS_NAN(x))
510 return x;
511 absx = fabs(x);
512 if (absx < ERF_SERIES_CUTOFF)
513 return 1.0 - m_erf_series(x);
514 else {
515 cf = m_erfc_contfrac(absx);
516 return x > 0.0 ? cf : 2.0 - cf;
517 }
Mark Dickinson45f992a2009-12-19 11:20:49 +0000518}
Mark Dickinson05d2e082009-12-11 20:17:17 +0000519
520/*
Christian Heimese57950f2008-04-21 13:08:03 +0000521 wrapper for atan2 that deals directly with special cases before
522 delegating to the platform libm for the remaining cases. This
523 is necessary to get consistent behaviour across platforms.
524 Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
525 always follow C99.
526*/
527
528static double
529m_atan2(double y, double x)
530{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000531 if (Py_IS_NAN(x) || Py_IS_NAN(y))
532 return Py_NAN;
533 if (Py_IS_INFINITY(y)) {
534 if (Py_IS_INFINITY(x)) {
535 if (copysign(1., x) == 1.)
536 /* atan2(+-inf, +inf) == +-pi/4 */
537 return copysign(0.25*Py_MATH_PI, y);
538 else
539 /* atan2(+-inf, -inf) == +-pi*3/4 */
540 return copysign(0.75*Py_MATH_PI, y);
541 }
542 /* atan2(+-inf, x) == +-pi/2 for finite x */
543 return copysign(0.5*Py_MATH_PI, y);
544 }
545 if (Py_IS_INFINITY(x) || y == 0.) {
546 if (copysign(1., x) == 1.)
547 /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
548 return copysign(0., y);
549 else
550 /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
551 return copysign(Py_MATH_PI, y);
552 }
553 return atan2(y, x);
Christian Heimese57950f2008-04-21 13:08:03 +0000554}
555
556/*
Mark Dickinsone675f082008-12-11 21:56:00 +0000557 Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
558 log(-ve), log(NaN). Here are wrappers for log and log10 that deal with
559 special values directly, passing positive non-special values through to
560 the system log/log10.
561 */
562
563static double
564m_log(double x)
565{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000566 if (Py_IS_FINITE(x)) {
567 if (x > 0.0)
568 return log(x);
569 errno = EDOM;
570 if (x == 0.0)
571 return -Py_HUGE_VAL; /* log(0) = -inf */
572 else
573 return Py_NAN; /* log(-ve) = nan */
574 }
575 else if (Py_IS_NAN(x))
576 return x; /* log(nan) = nan */
577 else if (x > 0.0)
578 return x; /* log(inf) = inf */
579 else {
580 errno = EDOM;
581 return Py_NAN; /* log(-inf) = nan */
582 }
Mark Dickinsone675f082008-12-11 21:56:00 +0000583}
584
585static double
586m_log10(double x)
587{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000588 if (Py_IS_FINITE(x)) {
589 if (x > 0.0)
590 return log10(x);
591 errno = EDOM;
592 if (x == 0.0)
593 return -Py_HUGE_VAL; /* log10(0) = -inf */
594 else
595 return Py_NAN; /* log10(-ve) = nan */
596 }
597 else if (Py_IS_NAN(x))
598 return x; /* log10(nan) = nan */
599 else if (x > 0.0)
600 return x; /* log10(inf) = inf */
601 else {
602 errno = EDOM;
603 return Py_NAN; /* log10(-inf) = nan */
604 }
Mark Dickinsone675f082008-12-11 21:56:00 +0000605}
606
607
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000608/* Call is_error when errno != 0, and where x is the result libm
609 * returned. is_error will usually set up an exception and return
610 * true (1), but may return false (0) without setting up an exception.
611 */
612static int
613is_error(double x)
614{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000615 int result = 1; /* presumption of guilt */
616 assert(errno); /* non-zero errno is a precondition for calling */
617 if (errno == EDOM)
618 PyErr_SetString(PyExc_ValueError, "math domain error");
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000619
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000620 else if (errno == ERANGE) {
621 /* ANSI C generally requires libm functions to set ERANGE
622 * on overflow, but also generally *allows* them to set
623 * ERANGE on underflow too. There's no consistency about
624 * the latter across platforms.
625 * Alas, C99 never requires that errno be set.
626 * Here we suppress the underflow errors (libm functions
627 * should return a zero on underflow, and +- HUGE_VAL on
628 * overflow, so testing the result for zero suffices to
629 * distinguish the cases).
630 *
631 * On some platforms (Ubuntu/ia64) it seems that errno can be
632 * set to ERANGE for subnormal results that do *not* underflow
633 * to zero. So to be safe, we'll ignore ERANGE whenever the
634 * function result is less than one in absolute value.
635 */
636 if (fabs(x) < 1.0)
637 result = 0;
638 else
639 PyErr_SetString(PyExc_OverflowError,
640 "math range error");
641 }
642 else
643 /* Unexpected math error */
644 PyErr_SetFromErrno(PyExc_ValueError);
645 return result;
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000646}
647
Mark Dickinsone675f082008-12-11 21:56:00 +0000648/*
Christian Heimes53876d92008-04-19 00:31:39 +0000649 math_1 is used to wrap a libm function f that takes a double
650 arguments and returns a double.
651
652 The error reporting follows these rules, which are designed to do
653 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
654 platforms.
655
656 - a NaN result from non-NaN inputs causes ValueError to be raised
657 - an infinite result from finite inputs causes OverflowError to be
658 raised if can_overflow is 1, or raises ValueError if can_overflow
659 is 0.
660 - if the result is finite and errno == EDOM then ValueError is
661 raised
662 - if the result is finite and nonzero and errno == ERANGE then
663 OverflowError is raised
664
665 The last rule is used to catch overflow on platforms which follow
666 C89 but for which HUGE_VAL is not an infinity.
667
668 For the majority of one-argument functions these rules are enough
669 to ensure that Python's functions behave as specified in 'Annex F'
670 of the C99 standard, with the 'invalid' and 'divide-by-zero'
671 floating-point exceptions mapping to Python's ValueError and the
672 'overflow' floating-point exception mapping to OverflowError.
673 math_1 only works for functions that don't have singularities *and*
674 the possibility of overflow; fortunately, that covers everything we
675 care about right now.
676*/
677
Barry Warsaw8b43b191996-12-09 22:32:36 +0000678static PyObject *
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000679math_1_to_whatever(PyObject *arg, double (*func) (double),
Christian Heimes53876d92008-04-19 00:31:39 +0000680 PyObject *(*from_double_func) (double),
681 int can_overflow)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000682{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000683 double x, r;
684 x = PyFloat_AsDouble(arg);
685 if (x == -1.0 && PyErr_Occurred())
686 return NULL;
687 errno = 0;
688 PyFPE_START_PROTECT("in math_1", return 0);
689 r = (*func)(x);
690 PyFPE_END_PROTECT(r);
691 if (Py_IS_NAN(r) && !Py_IS_NAN(x)) {
692 PyErr_SetString(PyExc_ValueError,
693 "math domain error"); /* invalid arg */
694 return NULL;
695 }
696 if (Py_IS_INFINITY(r) && Py_IS_FINITE(x)) {
697 if (can_overflow)
698 PyErr_SetString(PyExc_OverflowError,
699 "math range error"); /* overflow */
700 else
701 PyErr_SetString(PyExc_ValueError,
702 "math domain error"); /* singularity */
703 return NULL;
704 }
705 if (Py_IS_FINITE(r) && errno && is_error(r))
706 /* this branch unnecessary on most platforms */
707 return NULL;
Mark Dickinsonde429622008-05-01 00:19:23 +0000708
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000709 return (*from_double_func)(r);
Christian Heimes53876d92008-04-19 00:31:39 +0000710}
711
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000712/* variant of math_1, to be used when the function being wrapped is known to
713 set errno properly (that is, errno = EDOM for invalid or divide-by-zero,
714 errno = ERANGE for overflow). */
715
716static PyObject *
717math_1a(PyObject *arg, double (*func) (double))
718{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000719 double x, r;
720 x = PyFloat_AsDouble(arg);
721 if (x == -1.0 && PyErr_Occurred())
722 return NULL;
723 errno = 0;
724 PyFPE_START_PROTECT("in math_1a", return 0);
725 r = (*func)(x);
726 PyFPE_END_PROTECT(r);
727 if (errno && is_error(r))
728 return NULL;
729 return PyFloat_FromDouble(r);
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000730}
731
Christian Heimes53876d92008-04-19 00:31:39 +0000732/*
733 math_2 is used to wrap a libm function f that takes two double
734 arguments and returns a double.
735
736 The error reporting follows these rules, which are designed to do
737 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
738 platforms.
739
740 - a NaN result from non-NaN inputs causes ValueError to be raised
741 - an infinite result from finite inputs causes OverflowError to be
742 raised.
743 - if the result is finite and errno == EDOM then ValueError is
744 raised
745 - if the result is finite and nonzero and errno == ERANGE then
746 OverflowError is raised
747
748 The last rule is used to catch overflow on platforms which follow
749 C89 but for which HUGE_VAL is not an infinity.
750
751 For most two-argument functions (copysign, fmod, hypot, atan2)
752 these rules are enough to ensure that Python's functions behave as
753 specified in 'Annex F' of the C99 standard, with the 'invalid' and
754 'divide-by-zero' floating-point exceptions mapping to Python's
755 ValueError and the 'overflow' floating-point exception mapping to
756 OverflowError.
757*/
758
759static PyObject *
760math_1(PyObject *arg, double (*func) (double), int can_overflow)
761{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000762 return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow);
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000763}
764
765static PyObject *
Christian Heimes53876d92008-04-19 00:31:39 +0000766math_1_to_int(PyObject *arg, double (*func) (double), int can_overflow)
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000767{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000768 return math_1_to_whatever(arg, func, PyLong_FromDouble, can_overflow);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000769}
770
Barry Warsaw8b43b191996-12-09 22:32:36 +0000771static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000772math_2(PyObject *args, double (*func) (double, double), char *funcname)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000773{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000774 PyObject *ox, *oy;
775 double x, y, r;
776 if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
777 return NULL;
778 x = PyFloat_AsDouble(ox);
779 y = PyFloat_AsDouble(oy);
780 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
781 return NULL;
782 errno = 0;
783 PyFPE_START_PROTECT("in math_2", return 0);
784 r = (*func)(x, y);
785 PyFPE_END_PROTECT(r);
786 if (Py_IS_NAN(r)) {
787 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
788 errno = EDOM;
789 else
790 errno = 0;
791 }
792 else if (Py_IS_INFINITY(r)) {
793 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
794 errno = ERANGE;
795 else
796 errno = 0;
797 }
798 if (errno && is_error(r))
799 return NULL;
800 else
801 return PyFloat_FromDouble(r);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000802}
803
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000804#define FUNC1(funcname, func, can_overflow, docstring) \
805 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
806 return math_1(args, func, can_overflow); \
807 }\
808 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000809
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000810#define FUNC1A(funcname, func, docstring) \
811 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
812 return math_1a(args, func); \
813 }\
814 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000815
Fred Drake40c48682000-07-03 18:11:56 +0000816#define FUNC2(funcname, func, docstring) \
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000817 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
818 return math_2(args, func, #funcname); \
819 }\
820 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000821
Christian Heimes53876d92008-04-19 00:31:39 +0000822FUNC1(acos, acos, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000823 "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
Mark Dickinsonf3718592009-12-21 15:27:41 +0000824FUNC1(acosh, m_acosh, 0,
Christian Heimes53876d92008-04-19 00:31:39 +0000825 "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
826FUNC1(asin, asin, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000827 "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
Mark Dickinsonf3718592009-12-21 15:27:41 +0000828FUNC1(asinh, m_asinh, 0,
Christian Heimes53876d92008-04-19 00:31:39 +0000829 "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
830FUNC1(atan, atan, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000831 "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
Christian Heimese57950f2008-04-21 13:08:03 +0000832FUNC2(atan2, m_atan2,
Tim Petersfe71f812001-08-07 22:10:00 +0000833 "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
834 "Unlike atan(y/x), the signs of both x and y are considered.")
Mark Dickinsonf3718592009-12-21 15:27:41 +0000835FUNC1(atanh, m_atanh, 0,
Christian Heimes53876d92008-04-19 00:31:39 +0000836 "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
Guido van Rossum13e05de2007-08-23 22:56:55 +0000837
838static PyObject * math_ceil(PyObject *self, PyObject *number) {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000839 static PyObject *ceil_str = NULL;
Mark Dickinson6d02d9c2010-07-02 16:05:15 +0000840 PyObject *method, *result;
Guido van Rossum13e05de2007-08-23 22:56:55 +0000841
Benjamin Petersonf751bc92010-07-02 13:46:42 +0000842 method = _PyObject_LookupSpecial(number, "__ceil__", &ceil_str);
843 if (method == NULL) {
844 if (PyErr_Occurred())
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000845 return NULL;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000846 return math_1_to_int(number, ceil, 0);
Benjamin Petersonf751bc92010-07-02 13:46:42 +0000847 }
Mark Dickinson6d02d9c2010-07-02 16:05:15 +0000848 result = PyObject_CallFunctionObjArgs(method, NULL);
849 Py_DECREF(method);
850 return result;
Guido van Rossum13e05de2007-08-23 22:56:55 +0000851}
852
853PyDoc_STRVAR(math_ceil_doc,
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000854 "ceil(x)\n\nReturn the ceiling of x as an int.\n"
855 "This is the smallest integral value >= x.");
Guido van Rossum13e05de2007-08-23 22:56:55 +0000856
Christian Heimes072c0f12008-01-03 23:01:04 +0000857FUNC2(copysign, copysign,
Benjamin Petersona0dfa822009-11-13 02:25:08 +0000858 "copysign(x, y)\n\nReturn x with the sign of y.")
Christian Heimes53876d92008-04-19 00:31:39 +0000859FUNC1(cos, cos, 0,
860 "cos(x)\n\nReturn the cosine of x (measured in radians).")
861FUNC1(cosh, cosh, 1,
862 "cosh(x)\n\nReturn the hyperbolic cosine of x.")
Mark Dickinson45f992a2009-12-19 11:20:49 +0000863FUNC1A(erf, m_erf,
864 "erf(x)\n\nError function at x.")
865FUNC1A(erfc, m_erfc,
866 "erfc(x)\n\nComplementary error function at x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000867FUNC1(exp, exp, 1,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000868 "exp(x)\n\nReturn e raised to the power of x.")
Mark Dickinson664b5112009-12-16 20:23:42 +0000869FUNC1(expm1, m_expm1, 1,
870 "expm1(x)\n\nReturn exp(x)-1.\n"
871 "This function avoids the loss of precision involved in the direct "
872 "evaluation of exp(x)-1 for small x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000873FUNC1(fabs, fabs, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000874 "fabs(x)\n\nReturn the absolute value of the float x.")
Guido van Rossum13e05de2007-08-23 22:56:55 +0000875
876static PyObject * math_floor(PyObject *self, PyObject *number) {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000877 static PyObject *floor_str = NULL;
Benjamin Petersonb0125892010-07-02 13:35:17 +0000878 PyObject *method, *result;
Guido van Rossum13e05de2007-08-23 22:56:55 +0000879
Benjamin Peterson8bb9cde2010-07-01 15:16:55 +0000880 method = _PyObject_LookupSpecial(number, "__floor__", &floor_str);
881 if (method == NULL) {
882 if (PyErr_Occurred())
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000883 return NULL;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000884 return math_1_to_int(number, floor, 0);
Benjamin Peterson8bb9cde2010-07-01 15:16:55 +0000885 }
Benjamin Petersonb0125892010-07-02 13:35:17 +0000886 result = PyObject_CallFunctionObjArgs(method, NULL);
887 Py_DECREF(method);
888 return result;
Guido van Rossum13e05de2007-08-23 22:56:55 +0000889}
890
891PyDoc_STRVAR(math_floor_doc,
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000892 "floor(x)\n\nReturn the floor of x as an int.\n"
893 "This is the largest integral value <= x.");
Guido van Rossum13e05de2007-08-23 22:56:55 +0000894
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000895FUNC1A(gamma, m_tgamma,
896 "gamma(x)\n\nGamma function at x.")
Mark Dickinson05d2e082009-12-11 20:17:17 +0000897FUNC1A(lgamma, m_lgamma,
898 "lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.")
Mark Dickinsonbe64d952010-07-07 16:21:29 +0000899FUNC1(log1p, m_log1p, 0,
Benjamin Petersona0dfa822009-11-13 02:25:08 +0000900 "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n"
901 "The result is computed in a way which is accurate for x near zero.")
Christian Heimes53876d92008-04-19 00:31:39 +0000902FUNC1(sin, sin, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000903 "sin(x)\n\nReturn the sine of x (measured in radians).")
Christian Heimes53876d92008-04-19 00:31:39 +0000904FUNC1(sinh, sinh, 1,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000905 "sinh(x)\n\nReturn the hyperbolic sine of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000906FUNC1(sqrt, sqrt, 0,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000907 "sqrt(x)\n\nReturn the square root of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000908FUNC1(tan, tan, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000909 "tan(x)\n\nReturn the tangent of x (measured in radians).")
Christian Heimes53876d92008-04-19 00:31:39 +0000910FUNC1(tanh, tanh, 0,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000911 "tanh(x)\n\nReturn the hyperbolic tangent of x.")
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000912
Benjamin Peterson2b7411d2008-05-26 17:36:47 +0000913/* Precision summation function as msum() by Raymond Hettinger in
914 <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
915 enhanced with the exact partials sum and roundoff from Mark
916 Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
917 See those links for more details, proofs and other references.
918
919 Note 1: IEEE 754R floating point semantics are assumed,
920 but the current implementation does not re-establish special
921 value semantics across iterations (i.e. handling -Inf + Inf).
922
923 Note 2: No provision is made for intermediate overflow handling;
Georg Brandlf78e02b2008-06-10 17:40:04 +0000924 therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while
Benjamin Peterson2b7411d2008-05-26 17:36:47 +0000925 sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
926 overflow of the first partial sum.
927
Benjamin Petersonfea6a942008-07-02 16:11:42 +0000928 Note 3: The intermediate values lo, yr, and hi are declared volatile so
929 aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
Georg Brandlf78e02b2008-06-10 17:40:04 +0000930 Also, the volatile declaration forces the values to be stored in memory as
931 regular doubles instead of extended long precision (80-bit) values. This
Benjamin Petersonfea6a942008-07-02 16:11:42 +0000932 prevents double rounding because any addition or subtraction of two doubles
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000933 can be resolved exactly into double-sized hi and lo values. As long as the
Georg Brandlf78e02b2008-06-10 17:40:04 +0000934 hi value gets forced into a double before yr and lo are computed, the extra
935 bits in downstream extended precision operations (x87 for example) will be
936 exactly zero and therefore can be losslessly stored back into a double,
937 thereby preventing double rounding.
Benjamin Peterson2b7411d2008-05-26 17:36:47 +0000938
939 Note 4: A similar implementation is in Modules/cmathmodule.c.
940 Be sure to update both when making changes.
941
Mark Dickinsonaa7633a2008-08-01 08:16:13 +0000942 Note 5: The signature of math.fsum() differs from __builtin__.sum()
Benjamin Peterson2b7411d2008-05-26 17:36:47 +0000943 because the start argument doesn't make sense in the context of
944 accurate summation. Since the partials table is collapsed before
945 returning a result, sum(seq2, start=sum(seq1)) may not equal the
946 accurate result returned by sum(itertools.chain(seq1, seq2)).
947*/
948
949#define NUM_PARTIALS 32 /* initial partials array size, on stack */
950
951/* Extend the partials array p[] by doubling its size. */
952static int /* non-zero on error */
Mark Dickinsonaa7633a2008-08-01 08:16:13 +0000953_fsum_realloc(double **p_ptr, Py_ssize_t n,
Benjamin Peterson2b7411d2008-05-26 17:36:47 +0000954 double *ps, Py_ssize_t *m_ptr)
955{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000956 void *v = NULL;
957 Py_ssize_t m = *m_ptr;
Benjamin Peterson2b7411d2008-05-26 17:36:47 +0000958
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000959 m += m; /* double */
960 if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
961 double *p = *p_ptr;
962 if (p == ps) {
963 v = PyMem_Malloc(sizeof(double) * m);
964 if (v != NULL)
965 memcpy(v, ps, sizeof(double) * n);
966 }
967 else
968 v = PyMem_Realloc(p, sizeof(double) * m);
969 }
970 if (v == NULL) { /* size overflow or no memory */
971 PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
972 return 1;
973 }
974 *p_ptr = (double*) v;
975 *m_ptr = m;
976 return 0;
Benjamin Peterson2b7411d2008-05-26 17:36:47 +0000977}
978
979/* Full precision summation of a sequence of floats.
980
981 def msum(iterable):
982 partials = [] # sorted, non-overlapping partial sums
983 for x in iterable:
Mark Dickinsonfdb0acc2010-06-25 20:22:24 +0000984 i = 0
985 for y in partials:
986 if abs(x) < abs(y):
987 x, y = y, x
988 hi = x + y
989 lo = y - (hi - x)
990 if lo:
991 partials[i] = lo
992 i += 1
993 x = hi
994 partials[i:] = [x]
Benjamin Peterson2b7411d2008-05-26 17:36:47 +0000995 return sum_exact(partials)
996
997 Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo
998 are exactly equal to x+y. The inner loop applies hi/lo summation to each
999 partial so that the list of partial sums remains exact.
1000
1001 Sum_exact() adds the partial sums exactly and correctly rounds the final
1002 result (using the round-half-to-even rule). The items in partials remain
1003 non-zero, non-special, non-overlapping and strictly increasing in
1004 magnitude, but possibly not all having the same sign.
1005
1006 Depends on IEEE 754 arithmetic guarantees and half-even rounding.
1007*/
1008
1009static PyObject*
Mark Dickinsonaa7633a2008-08-01 08:16:13 +00001010math_fsum(PyObject *self, PyObject *seq)
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001011{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001012 PyObject *item, *iter, *sum = NULL;
1013 Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
1014 double x, y, t, ps[NUM_PARTIALS], *p = ps;
1015 double xsave, special_sum = 0.0, inf_sum = 0.0;
1016 volatile double hi, yr, lo;
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001017
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001018 iter = PyObject_GetIter(seq);
1019 if (iter == NULL)
1020 return NULL;
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001021
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001022 PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001023
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001024 for(;;) { /* for x in iterable */
1025 assert(0 <= n && n <= m);
1026 assert((m == NUM_PARTIALS && p == ps) ||
1027 (m > NUM_PARTIALS && p != NULL));
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001028
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001029 item = PyIter_Next(iter);
1030 if (item == NULL) {
1031 if (PyErr_Occurred())
1032 goto _fsum_error;
1033 break;
1034 }
1035 x = PyFloat_AsDouble(item);
1036 Py_DECREF(item);
1037 if (PyErr_Occurred())
1038 goto _fsum_error;
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001039
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001040 xsave = x;
1041 for (i = j = 0; j < n; j++) { /* for y in partials */
1042 y = p[j];
1043 if (fabs(x) < fabs(y)) {
1044 t = x; x = y; y = t;
1045 }
1046 hi = x + y;
1047 yr = hi - x;
1048 lo = y - yr;
1049 if (lo != 0.0)
1050 p[i++] = lo;
1051 x = hi;
1052 }
Mark Dickinsonaa7633a2008-08-01 08:16:13 +00001053
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001054 n = i; /* ps[i:] = [x] */
1055 if (x != 0.0) {
1056 if (! Py_IS_FINITE(x)) {
1057 /* a nonfinite x could arise either as
1058 a result of intermediate overflow, or
1059 as a result of a nan or inf in the
1060 summands */
1061 if (Py_IS_FINITE(xsave)) {
1062 PyErr_SetString(PyExc_OverflowError,
1063 "intermediate overflow in fsum");
1064 goto _fsum_error;
1065 }
1066 if (Py_IS_INFINITY(xsave))
1067 inf_sum += xsave;
1068 special_sum += xsave;
1069 /* reset partials */
1070 n = 0;
1071 }
1072 else if (n >= m && _fsum_realloc(&p, n, ps, &m))
1073 goto _fsum_error;
1074 else
1075 p[n++] = x;
1076 }
1077 }
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001078
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001079 if (special_sum != 0.0) {
1080 if (Py_IS_NAN(inf_sum))
1081 PyErr_SetString(PyExc_ValueError,
1082 "-inf + inf in fsum");
1083 else
1084 sum = PyFloat_FromDouble(special_sum);
1085 goto _fsum_error;
1086 }
Mark Dickinsonaa7633a2008-08-01 08:16:13 +00001087
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001088 hi = 0.0;
1089 if (n > 0) {
1090 hi = p[--n];
1091 /* sum_exact(ps, hi) from the top, stop when the sum becomes
1092 inexact. */
1093 while (n > 0) {
1094 x = hi;
1095 y = p[--n];
1096 assert(fabs(y) < fabs(x));
1097 hi = x + y;
1098 yr = hi - x;
1099 lo = y - yr;
1100 if (lo != 0.0)
1101 break;
1102 }
1103 /* Make half-even rounding work across multiple partials.
1104 Needed so that sum([1e-16, 1, 1e16]) will round-up the last
1105 digit to two instead of down to zero (the 1e-16 makes the 1
1106 slightly closer to two). With a potential 1 ULP rounding
1107 error fixed-up, math.fsum() can guarantee commutativity. */
1108 if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
1109 (lo > 0.0 && p[n-1] > 0.0))) {
1110 y = lo * 2.0;
1111 x = hi + y;
1112 yr = x - hi;
1113 if (y == yr)
1114 hi = x;
1115 }
1116 }
1117 sum = PyFloat_FromDouble(hi);
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001118
Mark Dickinsonaa7633a2008-08-01 08:16:13 +00001119_fsum_error:
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001120 PyFPE_END_PROTECT(hi)
1121 Py_DECREF(iter);
1122 if (p != ps)
1123 PyMem_Free(p);
1124 return sum;
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001125}
1126
1127#undef NUM_PARTIALS
1128
Mark Dickinsonaa7633a2008-08-01 08:16:13 +00001129PyDoc_STRVAR(math_fsum_doc,
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001130"fsum(iterable)\n\n\
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001131Return an accurate floating point sum of values in the iterable.\n\
1132Assumes IEEE-754 floating point arithmetic.");
1133
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001134/* Return the smallest integer k such that n < 2**k, or 0 if n == 0.
1135 * Equivalent to floor(lg(x))+1. Also equivalent to: bitwidth_of_type -
1136 * count_leading_zero_bits(x)
1137 */
1138
1139/* XXX: This routine does more or less the same thing as
1140 * bits_in_digit() in Objects/longobject.c. Someday it would be nice to
1141 * consolidate them. On BSD, there's a library function called fls()
1142 * that we could use, and GCC provides __builtin_clz().
1143 */
1144
1145static unsigned long
1146bit_length(unsigned long n)
1147{
1148 unsigned long len = 0;
1149 while (n != 0) {
1150 ++len;
1151 n >>= 1;
1152 }
1153 return len;
1154}
1155
1156static unsigned long
1157count_set_bits(unsigned long n)
1158{
1159 unsigned long count = 0;
1160 while (n != 0) {
1161 ++count;
1162 n &= n - 1; /* clear least significant bit */
1163 }
1164 return count;
1165}
1166
1167/* Divide-and-conquer factorial algorithm
1168 *
1169 * Based on the formula and psuedo-code provided at:
1170 * http://www.luschny.de/math/factorial/binarysplitfact.html
1171 *
1172 * Faster algorithms exist, but they're more complicated and depend on
Ezio Melotti9527afd2010-07-08 15:03:02 +00001173 * a fast prime factorization algorithm.
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001174 *
1175 * Notes on the algorithm
1176 * ----------------------
1177 *
1178 * factorial(n) is written in the form 2**k * m, with m odd. k and m are
1179 * computed separately, and then combined using a left shift.
1180 *
1181 * The function factorial_odd_part computes the odd part m (i.e., the greatest
1182 * odd divisor) of factorial(n), using the formula:
1183 *
1184 * factorial_odd_part(n) =
1185 *
1186 * product_{i >= 0} product_{0 < j <= n / 2**i, j odd} j
1187 *
1188 * Example: factorial_odd_part(20) =
1189 *
1190 * (1) *
1191 * (1) *
1192 * (1 * 3 * 5) *
1193 * (1 * 3 * 5 * 7 * 9)
1194 * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1195 *
1196 * Here i goes from large to small: the first term corresponds to i=4 (any
1197 * larger i gives an empty product), and the last term corresponds to i=0.
1198 * Each term can be computed from the last by multiplying by the extra odd
1199 * numbers required: e.g., to get from the penultimate term to the last one,
1200 * we multiply by (11 * 13 * 15 * 17 * 19).
1201 *
1202 * To see a hint of why this formula works, here are the same numbers as above
1203 * but with the even parts (i.e., the appropriate powers of 2) included. For
1204 * each subterm in the product for i, we multiply that subterm by 2**i:
1205 *
1206 * factorial(20) =
1207 *
1208 * (16) *
1209 * (8) *
1210 * (4 * 12 * 20) *
1211 * (2 * 6 * 10 * 14 * 18) *
1212 * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1213 *
1214 * The factorial_partial_product function computes the product of all odd j in
1215 * range(start, stop) for given start and stop. It's used to compute the
1216 * partial products like (11 * 13 * 15 * 17 * 19) in the example above. It
1217 * operates recursively, repeatedly splitting the range into two roughly equal
1218 * pieces until the subranges are small enough to be computed using only C
1219 * integer arithmetic.
1220 *
1221 * The two-valuation k (i.e., the exponent of the largest power of 2 dividing
1222 * the factorial) is computed independently in the main math_factorial
1223 * function. By standard results, its value is:
1224 *
1225 * two_valuation = n//2 + n//4 + n//8 + ....
1226 *
1227 * It can be shown (e.g., by complete induction on n) that two_valuation is
1228 * equal to n - count_set_bits(n), where count_set_bits(n) gives the number of
1229 * '1'-bits in the binary expansion of n.
1230 */
1231
1232/* factorial_partial_product: Compute product(range(start, stop, 2)) using
1233 * divide and conquer. Assumes start and stop are odd and stop > start.
1234 * max_bits must be >= bit_length(stop - 2). */
1235
1236static PyObject *
1237factorial_partial_product(unsigned long start, unsigned long stop,
1238 unsigned long max_bits)
1239{
1240 unsigned long midpoint, num_operands;
1241 PyObject *left = NULL, *right = NULL, *result = NULL;
1242
1243 /* If the return value will fit an unsigned long, then we can
1244 * multiply in a tight, fast loop where each multiply is O(1).
1245 * Compute an upper bound on the number of bits required to store
1246 * the answer.
1247 *
1248 * Storing some integer z requires floor(lg(z))+1 bits, which is
1249 * conveniently the value returned by bit_length(z). The
1250 * product x*y will require at most
1251 * bit_length(x) + bit_length(y) bits to store, based
1252 * on the idea that lg product = lg x + lg y.
1253 *
1254 * We know that stop - 2 is the largest number to be multiplied. From
1255 * there, we have: bit_length(answer) <= num_operands *
1256 * bit_length(stop - 2)
1257 */
1258
1259 num_operands = (stop - start) / 2;
1260 /* The "num_operands <= 8 * SIZEOF_LONG" check guards against the
1261 * unlikely case of an overflow in num_operands * max_bits. */
1262 if (num_operands <= 8 * SIZEOF_LONG &&
1263 num_operands * max_bits <= 8 * SIZEOF_LONG) {
1264 unsigned long j, total;
1265 for (total = start, j = start + 2; j < stop; j += 2)
1266 total *= j;
1267 return PyLong_FromUnsignedLong(total);
1268 }
1269
1270 /* find midpoint of range(start, stop), rounded up to next odd number. */
1271 midpoint = (start + num_operands) | 1;
1272 left = factorial_partial_product(start, midpoint,
1273 bit_length(midpoint - 2));
1274 if (left == NULL)
1275 goto error;
1276 right = factorial_partial_product(midpoint, stop, max_bits);
1277 if (right == NULL)
1278 goto error;
1279 result = PyNumber_Multiply(left, right);
1280
1281 error:
1282 Py_XDECREF(left);
1283 Py_XDECREF(right);
1284 return result;
1285}
1286
1287/* factorial_odd_part: compute the odd part of factorial(n). */
1288
1289static PyObject *
1290factorial_odd_part(unsigned long n)
1291{
1292 long i;
1293 unsigned long v, lower, upper;
1294 PyObject *partial, *tmp, *inner, *outer;
1295
1296 inner = PyLong_FromLong(1);
1297 if (inner == NULL)
1298 return NULL;
1299 outer = inner;
1300 Py_INCREF(outer);
1301
1302 upper = 3;
1303 for (i = bit_length(n) - 2; i >= 0; i--) {
1304 v = n >> i;
1305 if (v <= 2)
1306 continue;
1307 lower = upper;
1308 /* (v + 1) | 1 = least odd integer strictly larger than n / 2**i */
1309 upper = (v + 1) | 1;
1310 /* Here inner is the product of all odd integers j in the range (0,
1311 n/2**(i+1)]. The factorial_partial_product call below gives the
1312 product of all odd integers j in the range (n/2**(i+1), n/2**i]. */
1313 partial = factorial_partial_product(lower, upper, bit_length(upper-2));
1314 /* inner *= partial */
1315 if (partial == NULL)
1316 goto error;
1317 tmp = PyNumber_Multiply(inner, partial);
1318 Py_DECREF(partial);
1319 if (tmp == NULL)
1320 goto error;
1321 Py_DECREF(inner);
1322 inner = tmp;
1323 /* Now inner is the product of all odd integers j in the range (0,
1324 n/2**i], giving the inner product in the formula above. */
1325
1326 /* outer *= inner; */
1327 tmp = PyNumber_Multiply(outer, inner);
1328 if (tmp == NULL)
1329 goto error;
1330 Py_DECREF(outer);
1331 outer = tmp;
1332 }
1333
1334 goto done;
1335
1336 error:
1337 Py_DECREF(outer);
1338 done:
1339 Py_DECREF(inner);
1340 return outer;
1341}
1342
1343/* Lookup table for small factorial values */
1344
1345static const unsigned long SmallFactorials[] = {
1346 1, 1, 2, 6, 24, 120, 720, 5040, 40320,
1347 362880, 3628800, 39916800, 479001600,
1348#if SIZEOF_LONG >= 8
1349 6227020800, 87178291200, 1307674368000,
1350 20922789888000, 355687428096000, 6402373705728000,
1351 121645100408832000, 2432902008176640000
1352#endif
1353};
1354
Barry Warsaw8b43b191996-12-09 22:32:36 +00001355static PyObject *
Georg Brandlc28e1fa2008-06-10 19:20:26 +00001356math_factorial(PyObject *self, PyObject *arg)
1357{
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001358 long x;
1359 PyObject *result, *odd_part, *two_valuation;
Georg Brandlc28e1fa2008-06-10 19:20:26 +00001360
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001361 if (PyFloat_Check(arg)) {
1362 PyObject *lx;
1363 double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);
1364 if (!(Py_IS_FINITE(dx) && dx == floor(dx))) {
1365 PyErr_SetString(PyExc_ValueError,
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001366 "factorial() only accepts integral values");
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001367 return NULL;
1368 }
1369 lx = PyLong_FromDouble(dx);
1370 if (lx == NULL)
1371 return NULL;
1372 x = PyLong_AsLong(lx);
1373 Py_DECREF(lx);
1374 }
1375 else
1376 x = PyLong_AsLong(arg);
Georg Brandlc28e1fa2008-06-10 19:20:26 +00001377
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001378 if (x == -1 && PyErr_Occurred())
1379 return NULL;
1380 if (x < 0) {
1381 PyErr_SetString(PyExc_ValueError,
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001382 "factorial() not defined for negative values");
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001383 return NULL;
1384 }
Georg Brandlc28e1fa2008-06-10 19:20:26 +00001385
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001386 /* use lookup table if x is small */
1387 if (x < (long)(sizeof(SmallFactorials)/sizeof(SmallFactorials[0])))
1388 return PyLong_FromUnsignedLong(SmallFactorials[x]);
1389
1390 /* else express in the form odd_part * 2**two_valuation, and compute as
1391 odd_part << two_valuation. */
1392 odd_part = factorial_odd_part(x);
1393 if (odd_part == NULL)
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001394 return NULL;
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001395 two_valuation = PyLong_FromLong(x - count_set_bits(x));
1396 if (two_valuation == NULL) {
1397 Py_DECREF(odd_part);
1398 return NULL;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001399 }
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001400 result = PyNumber_Lshift(odd_part, two_valuation);
1401 Py_DECREF(two_valuation);
1402 Py_DECREF(odd_part);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001403 return result;
Georg Brandlc28e1fa2008-06-10 19:20:26 +00001404}
1405
Benjamin Peterson6ebe78f2008-12-21 00:06:59 +00001406PyDoc_STRVAR(math_factorial_doc,
1407"factorial(x) -> Integral\n"
1408"\n"
1409"Find x!. Raise a ValueError if x is negative or non-integral.");
Georg Brandlc28e1fa2008-06-10 19:20:26 +00001410
1411static PyObject *
Christian Heimes400adb02008-02-01 08:12:03 +00001412math_trunc(PyObject *self, PyObject *number)
1413{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001414 static PyObject *trunc_str = NULL;
Benjamin Petersonb0125892010-07-02 13:35:17 +00001415 PyObject *trunc, *result;
Christian Heimes400adb02008-02-01 08:12:03 +00001416
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001417 if (Py_TYPE(number)->tp_dict == NULL) {
1418 if (PyType_Ready(Py_TYPE(number)) < 0)
1419 return NULL;
1420 }
Christian Heimes400adb02008-02-01 08:12:03 +00001421
Benjamin Peterson8bb9cde2010-07-01 15:16:55 +00001422 trunc = _PyObject_LookupSpecial(number, "__trunc__", &trunc_str);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001423 if (trunc == NULL) {
Benjamin Peterson8bb9cde2010-07-01 15:16:55 +00001424 if (!PyErr_Occurred())
1425 PyErr_Format(PyExc_TypeError,
1426 "type %.100s doesn't define __trunc__ method",
1427 Py_TYPE(number)->tp_name);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001428 return NULL;
1429 }
Benjamin Petersonb0125892010-07-02 13:35:17 +00001430 result = PyObject_CallFunctionObjArgs(trunc, NULL);
1431 Py_DECREF(trunc);
1432 return result;
Christian Heimes400adb02008-02-01 08:12:03 +00001433}
1434
1435PyDoc_STRVAR(math_trunc_doc,
1436"trunc(x:Real) -> Integral\n"
1437"\n"
Christian Heimes292d3512008-02-03 16:51:08 +00001438"Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
Christian Heimes400adb02008-02-01 08:12:03 +00001439
1440static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +00001441math_frexp(PyObject *self, PyObject *arg)
Guido van Rossumd18ad581991-10-24 14:57:21 +00001442{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001443 int i;
1444 double x = PyFloat_AsDouble(arg);
1445 if (x == -1.0 && PyErr_Occurred())
1446 return NULL;
1447 /* deal with special cases directly, to sidestep platform
1448 differences */
1449 if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
1450 i = 0;
1451 }
1452 else {
1453 PyFPE_START_PROTECT("in math_frexp", return 0);
1454 x = frexp(x, &i);
1455 PyFPE_END_PROTECT(x);
1456 }
1457 return Py_BuildValue("(di)", x, i);
Guido van Rossumd18ad581991-10-24 14:57:21 +00001458}
1459
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +00001460PyDoc_STRVAR(math_frexp_doc,
Tim Peters63c94532001-09-04 23:17:42 +00001461"frexp(x)\n"
1462"\n"
1463"Return the mantissa and exponent of x, as pair (m, e).\n"
1464"m is a float and e is an int, such that x = m * 2.**e.\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +00001465"If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0.");
Guido van Rossumc6e22901998-12-04 19:26:43 +00001466
Barry Warsaw8b43b191996-12-09 22:32:36 +00001467static PyObject *
Fred Drake40c48682000-07-03 18:11:56 +00001468math_ldexp(PyObject *self, PyObject *args)
Guido van Rossumd18ad581991-10-24 14:57:21 +00001469{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001470 double x, r;
1471 PyObject *oexp;
1472 long exp;
1473 int overflow;
1474 if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp))
1475 return NULL;
Alexandre Vassalotti6461e102008-05-15 22:09:29 +00001476
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001477 if (PyLong_Check(oexp)) {
1478 /* on overflow, replace exponent with either LONG_MAX
1479 or LONG_MIN, depending on the sign. */
1480 exp = PyLong_AsLongAndOverflow(oexp, &overflow);
1481 if (exp == -1 && PyErr_Occurred())
1482 return NULL;
1483 if (overflow)
1484 exp = overflow < 0 ? LONG_MIN : LONG_MAX;
1485 }
1486 else {
1487 PyErr_SetString(PyExc_TypeError,
1488 "Expected an int or long as second argument "
1489 "to ldexp.");
1490 return NULL;
1491 }
Alexandre Vassalotti6461e102008-05-15 22:09:29 +00001492
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001493 if (x == 0. || !Py_IS_FINITE(x)) {
1494 /* NaNs, zeros and infinities are returned unchanged */
1495 r = x;
1496 errno = 0;
1497 } else if (exp > INT_MAX) {
1498 /* overflow */
1499 r = copysign(Py_HUGE_VAL, x);
1500 errno = ERANGE;
1501 } else if (exp < INT_MIN) {
1502 /* underflow to +-0 */
1503 r = copysign(0., x);
1504 errno = 0;
1505 } else {
1506 errno = 0;
1507 PyFPE_START_PROTECT("in math_ldexp", return 0);
1508 r = ldexp(x, (int)exp);
1509 PyFPE_END_PROTECT(r);
1510 if (Py_IS_INFINITY(r))
1511 errno = ERANGE;
1512 }
Alexandre Vassalotti6461e102008-05-15 22:09:29 +00001513
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001514 if (errno && is_error(r))
1515 return NULL;
1516 return PyFloat_FromDouble(r);
Guido van Rossumd18ad581991-10-24 14:57:21 +00001517}
1518
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +00001519PyDoc_STRVAR(math_ldexp_doc,
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001520"ldexp(x, i)\n\n\
1521Return x * (2**i).");
Guido van Rossumc6e22901998-12-04 19:26:43 +00001522
Barry Warsaw8b43b191996-12-09 22:32:36 +00001523static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +00001524math_modf(PyObject *self, PyObject *arg)
Guido van Rossumd18ad581991-10-24 14:57:21 +00001525{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001526 double y, x = PyFloat_AsDouble(arg);
1527 if (x == -1.0 && PyErr_Occurred())
1528 return NULL;
1529 /* some platforms don't do the right thing for NaNs and
1530 infinities, so we take care of special cases directly. */
1531 if (!Py_IS_FINITE(x)) {
1532 if (Py_IS_INFINITY(x))
1533 return Py_BuildValue("(dd)", copysign(0., x), x);
1534 else if (Py_IS_NAN(x))
1535 return Py_BuildValue("(dd)", x, x);
1536 }
Christian Heimesa342c012008-04-20 21:01:16 +00001537
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001538 errno = 0;
1539 PyFPE_START_PROTECT("in math_modf", return 0);
1540 x = modf(x, &y);
1541 PyFPE_END_PROTECT(x);
1542 return Py_BuildValue("(dd)", x, y);
Guido van Rossumd18ad581991-10-24 14:57:21 +00001543}
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001544
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +00001545PyDoc_STRVAR(math_modf_doc,
Tim Peters63c94532001-09-04 23:17:42 +00001546"modf(x)\n"
1547"\n"
1548"Return the fractional and integer parts of x. Both results carry the sign\n"
Benjamin Peterson6ebe78f2008-12-21 00:06:59 +00001549"of x and are floats.");
Guido van Rossumc6e22901998-12-04 19:26:43 +00001550
Tim Peters78526162001-09-05 00:53:45 +00001551/* A decent logarithm is easy to compute even for huge longs, but libm can't
1552 do that by itself -- loghelper can. func is log or log10, and name is
Mark Dickinson6ecd9e52010-01-02 15:33:56 +00001553 "log" or "log10". Note that overflow of the result isn't possible: a long
1554 can contain no more than INT_MAX * SHIFT bits, so has value certainly less
1555 than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
Tim Peters78526162001-09-05 00:53:45 +00001556 small enough to fit in an IEEE single. log and log10 are even smaller.
Mark Dickinson6ecd9e52010-01-02 15:33:56 +00001557 However, intermediate overflow is possible for a long if the number of bits
1558 in that long is larger than PY_SSIZE_T_MAX. */
Tim Peters78526162001-09-05 00:53:45 +00001559
1560static PyObject*
Thomas Wouters89f507f2006-12-13 04:49:30 +00001561loghelper(PyObject* arg, double (*func)(double), char *funcname)
Tim Peters78526162001-09-05 00:53:45 +00001562{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001563 /* If it is long, do it ourselves. */
1564 if (PyLong_Check(arg)) {
Mark Dickinsonc6037172010-09-29 19:06:36 +00001565 double x, result;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001566 Py_ssize_t e;
Mark Dickinsonc6037172010-09-29 19:06:36 +00001567
1568 /* Negative or zero inputs give a ValueError. */
1569 if (Py_SIZE(arg) <= 0) {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001570 PyErr_SetString(PyExc_ValueError,
1571 "math domain error");
1572 return NULL;
1573 }
Mark Dickinsonfa41e602010-09-28 07:22:27 +00001574
Mark Dickinsonc6037172010-09-29 19:06:36 +00001575 x = PyLong_AsDouble(arg);
1576 if (x == -1.0 && PyErr_Occurred()) {
1577 if (!PyErr_ExceptionMatches(PyExc_OverflowError))
1578 return NULL;
1579 /* Here the conversion to double overflowed, but it's possible
1580 to compute the log anyway. Clear the exception and continue. */
1581 PyErr_Clear();
1582 x = _PyLong_Frexp((PyLongObject *)arg, &e);
1583 if (x == -1.0 && PyErr_Occurred())
1584 return NULL;
1585 /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */
1586 result = func(x) + func(2.0) * e;
1587 }
1588 else
1589 /* Successfully converted x to a double. */
1590 result = func(x);
1591 return PyFloat_FromDouble(result);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001592 }
Tim Peters78526162001-09-05 00:53:45 +00001593
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001594 /* Else let libm handle it by itself. */
1595 return math_1(arg, func, 0);
Tim Peters78526162001-09-05 00:53:45 +00001596}
1597
1598static PyObject *
1599math_log(PyObject *self, PyObject *args)
1600{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001601 PyObject *arg;
1602 PyObject *base = NULL;
1603 PyObject *num, *den;
1604 PyObject *ans;
Raymond Hettinger866964c2002-12-14 19:51:34 +00001605
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001606 if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
1607 return NULL;
Raymond Hettinger866964c2002-12-14 19:51:34 +00001608
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001609 num = loghelper(arg, m_log, "log");
1610 if (num == NULL || base == NULL)
1611 return num;
Raymond Hettinger866964c2002-12-14 19:51:34 +00001612
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001613 den = loghelper(base, m_log, "log");
1614 if (den == NULL) {
1615 Py_DECREF(num);
1616 return NULL;
1617 }
Raymond Hettinger866964c2002-12-14 19:51:34 +00001618
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001619 ans = PyNumber_TrueDivide(num, den);
1620 Py_DECREF(num);
1621 Py_DECREF(den);
1622 return ans;
Tim Peters78526162001-09-05 00:53:45 +00001623}
1624
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +00001625PyDoc_STRVAR(math_log_doc,
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001626"log(x[, base])\n\n\
1627Return the logarithm of x to the given base.\n\
Raymond Hettinger866964c2002-12-14 19:51:34 +00001628If the base not specified, returns the natural logarithm (base e) of x.");
Tim Peters78526162001-09-05 00:53:45 +00001629
1630static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +00001631math_log10(PyObject *self, PyObject *arg)
Tim Peters78526162001-09-05 00:53:45 +00001632{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001633 return loghelper(arg, m_log10, "log10");
Tim Peters78526162001-09-05 00:53:45 +00001634}
1635
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +00001636PyDoc_STRVAR(math_log10_doc,
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001637"log10(x)\n\nReturn the base 10 logarithm of x.");
Tim Peters78526162001-09-05 00:53:45 +00001638
Christian Heimes53876d92008-04-19 00:31:39 +00001639static PyObject *
1640math_fmod(PyObject *self, PyObject *args)
1641{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001642 PyObject *ox, *oy;
1643 double r, x, y;
1644 if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
1645 return NULL;
1646 x = PyFloat_AsDouble(ox);
1647 y = PyFloat_AsDouble(oy);
1648 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1649 return NULL;
1650 /* fmod(x, +/-Inf) returns x for finite x. */
1651 if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
1652 return PyFloat_FromDouble(x);
1653 errno = 0;
1654 PyFPE_START_PROTECT("in math_fmod", return 0);
1655 r = fmod(x, y);
1656 PyFPE_END_PROTECT(r);
1657 if (Py_IS_NAN(r)) {
1658 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1659 errno = EDOM;
1660 else
1661 errno = 0;
1662 }
1663 if (errno && is_error(r))
1664 return NULL;
1665 else
1666 return PyFloat_FromDouble(r);
Christian Heimes53876d92008-04-19 00:31:39 +00001667}
1668
1669PyDoc_STRVAR(math_fmod_doc,
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001670"fmod(x, y)\n\nReturn fmod(x, y), according to platform C."
Christian Heimes53876d92008-04-19 00:31:39 +00001671" x % y may differ.");
1672
1673static PyObject *
1674math_hypot(PyObject *self, PyObject *args)
1675{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001676 PyObject *ox, *oy;
1677 double r, x, y;
1678 if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
1679 return NULL;
1680 x = PyFloat_AsDouble(ox);
1681 y = PyFloat_AsDouble(oy);
1682 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1683 return NULL;
1684 /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
1685 if (Py_IS_INFINITY(x))
1686 return PyFloat_FromDouble(fabs(x));
1687 if (Py_IS_INFINITY(y))
1688 return PyFloat_FromDouble(fabs(y));
1689 errno = 0;
1690 PyFPE_START_PROTECT("in math_hypot", return 0);
1691 r = hypot(x, y);
1692 PyFPE_END_PROTECT(r);
1693 if (Py_IS_NAN(r)) {
1694 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1695 errno = EDOM;
1696 else
1697 errno = 0;
1698 }
1699 else if (Py_IS_INFINITY(r)) {
1700 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
1701 errno = ERANGE;
1702 else
1703 errno = 0;
1704 }
1705 if (errno && is_error(r))
1706 return NULL;
1707 else
1708 return PyFloat_FromDouble(r);
Christian Heimes53876d92008-04-19 00:31:39 +00001709}
1710
1711PyDoc_STRVAR(math_hypot_doc,
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001712"hypot(x, y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
Christian Heimes53876d92008-04-19 00:31:39 +00001713
1714/* pow can't use math_2, but needs its own wrapper: the problem is
1715 that an infinite result can arise either as a result of overflow
1716 (in which case OverflowError should be raised) or as a result of
1717 e.g. 0.**-5. (for which ValueError needs to be raised.)
1718*/
1719
1720static PyObject *
1721math_pow(PyObject *self, PyObject *args)
1722{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001723 PyObject *ox, *oy;
1724 double r, x, y;
1725 int odd_y;
Christian Heimes53876d92008-04-19 00:31:39 +00001726
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001727 if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
1728 return NULL;
1729 x = PyFloat_AsDouble(ox);
1730 y = PyFloat_AsDouble(oy);
1731 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1732 return NULL;
Christian Heimesa342c012008-04-20 21:01:16 +00001733
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001734 /* deal directly with IEEE specials, to cope with problems on various
1735 platforms whose semantics don't exactly match C99 */
1736 r = 0.; /* silence compiler warning */
1737 if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
1738 errno = 0;
1739 if (Py_IS_NAN(x))
1740 r = y == 0. ? 1. : x; /* NaN**0 = 1 */
1741 else if (Py_IS_NAN(y))
1742 r = x == 1. ? 1. : y; /* 1**NaN = 1 */
1743 else if (Py_IS_INFINITY(x)) {
1744 odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
1745 if (y > 0.)
1746 r = odd_y ? x : fabs(x);
1747 else if (y == 0.)
1748 r = 1.;
1749 else /* y < 0. */
1750 r = odd_y ? copysign(0., x) : 0.;
1751 }
1752 else if (Py_IS_INFINITY(y)) {
1753 if (fabs(x) == 1.0)
1754 r = 1.;
1755 else if (y > 0. && fabs(x) > 1.0)
1756 r = y;
1757 else if (y < 0. && fabs(x) < 1.0) {
1758 r = -y; /* result is +inf */
1759 if (x == 0.) /* 0**-inf: divide-by-zero */
1760 errno = EDOM;
1761 }
1762 else
1763 r = 0.;
1764 }
1765 }
1766 else {
1767 /* let libm handle finite**finite */
1768 errno = 0;
1769 PyFPE_START_PROTECT("in math_pow", return 0);
1770 r = pow(x, y);
1771 PyFPE_END_PROTECT(r);
1772 /* a NaN result should arise only from (-ve)**(finite
1773 non-integer); in this case we want to raise ValueError. */
1774 if (!Py_IS_FINITE(r)) {
1775 if (Py_IS_NAN(r)) {
1776 errno = EDOM;
1777 }
1778 /*
1779 an infinite result here arises either from:
1780 (A) (+/-0.)**negative (-> divide-by-zero)
1781 (B) overflow of x**y with x and y finite
1782 */
1783 else if (Py_IS_INFINITY(r)) {
1784 if (x == 0.)
1785 errno = EDOM;
1786 else
1787 errno = ERANGE;
1788 }
1789 }
1790 }
Christian Heimes53876d92008-04-19 00:31:39 +00001791
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001792 if (errno && is_error(r))
1793 return NULL;
1794 else
1795 return PyFloat_FromDouble(r);
Christian Heimes53876d92008-04-19 00:31:39 +00001796}
1797
1798PyDoc_STRVAR(math_pow_doc,
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001799"pow(x, y)\n\nReturn x**y (x to the power of y).");
Christian Heimes53876d92008-04-19 00:31:39 +00001800
Christian Heimes072c0f12008-01-03 23:01:04 +00001801static const double degToRad = Py_MATH_PI / 180.0;
1802static const double radToDeg = 180.0 / Py_MATH_PI;
Raymond Hettingerd6f22672002-05-13 03:56:10 +00001803
1804static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +00001805math_degrees(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +00001806{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001807 double x = PyFloat_AsDouble(arg);
1808 if (x == -1.0 && PyErr_Occurred())
1809 return NULL;
1810 return PyFloat_FromDouble(x * radToDeg);
Raymond Hettingerd6f22672002-05-13 03:56:10 +00001811}
1812
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +00001813PyDoc_STRVAR(math_degrees_doc,
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001814"degrees(x)\n\n\
1815Convert angle x from radians to degrees.");
Raymond Hettingerd6f22672002-05-13 03:56:10 +00001816
1817static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +00001818math_radians(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +00001819{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001820 double x = PyFloat_AsDouble(arg);
1821 if (x == -1.0 && PyErr_Occurred())
1822 return NULL;
1823 return PyFloat_FromDouble(x * degToRad);
Raymond Hettingerd6f22672002-05-13 03:56:10 +00001824}
1825
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +00001826PyDoc_STRVAR(math_radians_doc,
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001827"radians(x)\n\n\
1828Convert angle x from degrees to radians.");
Tim Peters78526162001-09-05 00:53:45 +00001829
Christian Heimes072c0f12008-01-03 23:01:04 +00001830static PyObject *
Mark Dickinson8e0c9962010-07-11 17:38:24 +00001831math_isfinite(PyObject *self, PyObject *arg)
1832{
1833 double x = PyFloat_AsDouble(arg);
1834 if (x == -1.0 && PyErr_Occurred())
1835 return NULL;
1836 return PyBool_FromLong((long)Py_IS_FINITE(x));
1837}
1838
1839PyDoc_STRVAR(math_isfinite_doc,
1840"isfinite(x) -> bool\n\n\
Mark Dickinson226f5442010-07-11 18:13:41 +00001841Return True if x is neither an infinity nor a NaN, and False otherwise.");
Mark Dickinson8e0c9962010-07-11 17:38:24 +00001842
1843static PyObject *
Christian Heimes072c0f12008-01-03 23:01:04 +00001844math_isnan(PyObject *self, PyObject *arg)
1845{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001846 double x = PyFloat_AsDouble(arg);
1847 if (x == -1.0 && PyErr_Occurred())
1848 return NULL;
1849 return PyBool_FromLong((long)Py_IS_NAN(x));
Christian Heimes072c0f12008-01-03 23:01:04 +00001850}
1851
1852PyDoc_STRVAR(math_isnan_doc,
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001853"isnan(x) -> bool\n\n\
Mark Dickinson226f5442010-07-11 18:13:41 +00001854Return True if x is a NaN (not a number), and False otherwise.");
Christian Heimes072c0f12008-01-03 23:01:04 +00001855
1856static PyObject *
1857math_isinf(PyObject *self, PyObject *arg)
1858{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001859 double x = PyFloat_AsDouble(arg);
1860 if (x == -1.0 && PyErr_Occurred())
1861 return NULL;
1862 return PyBool_FromLong((long)Py_IS_INFINITY(x));
Christian Heimes072c0f12008-01-03 23:01:04 +00001863}
1864
1865PyDoc_STRVAR(math_isinf_doc,
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001866"isinf(x) -> bool\n\n\
Mark Dickinson226f5442010-07-11 18:13:41 +00001867Return True if x is a positive or negative infinity, and False otherwise.");
Christian Heimes072c0f12008-01-03 23:01:04 +00001868
Barry Warsaw8b43b191996-12-09 22:32:36 +00001869static PyMethodDef math_methods[] = {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001870 {"acos", math_acos, METH_O, math_acos_doc},
1871 {"acosh", math_acosh, METH_O, math_acosh_doc},
1872 {"asin", math_asin, METH_O, math_asin_doc},
1873 {"asinh", math_asinh, METH_O, math_asinh_doc},
1874 {"atan", math_atan, METH_O, math_atan_doc},
1875 {"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
1876 {"atanh", math_atanh, METH_O, math_atanh_doc},
1877 {"ceil", math_ceil, METH_O, math_ceil_doc},
1878 {"copysign", math_copysign, METH_VARARGS, math_copysign_doc},
1879 {"cos", math_cos, METH_O, math_cos_doc},
1880 {"cosh", math_cosh, METH_O, math_cosh_doc},
1881 {"degrees", math_degrees, METH_O, math_degrees_doc},
1882 {"erf", math_erf, METH_O, math_erf_doc},
1883 {"erfc", math_erfc, METH_O, math_erfc_doc},
1884 {"exp", math_exp, METH_O, math_exp_doc},
1885 {"expm1", math_expm1, METH_O, math_expm1_doc},
1886 {"fabs", math_fabs, METH_O, math_fabs_doc},
1887 {"factorial", math_factorial, METH_O, math_factorial_doc},
1888 {"floor", math_floor, METH_O, math_floor_doc},
1889 {"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
1890 {"frexp", math_frexp, METH_O, math_frexp_doc},
1891 {"fsum", math_fsum, METH_O, math_fsum_doc},
1892 {"gamma", math_gamma, METH_O, math_gamma_doc},
1893 {"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
Mark Dickinson8e0c9962010-07-11 17:38:24 +00001894 {"isfinite", math_isfinite, METH_O, math_isfinite_doc},
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001895 {"isinf", math_isinf, METH_O, math_isinf_doc},
1896 {"isnan", math_isnan, METH_O, math_isnan_doc},
1897 {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
1898 {"lgamma", math_lgamma, METH_O, math_lgamma_doc},
1899 {"log", math_log, METH_VARARGS, math_log_doc},
1900 {"log1p", math_log1p, METH_O, math_log1p_doc},
1901 {"log10", math_log10, METH_O, math_log10_doc},
1902 {"modf", math_modf, METH_O, math_modf_doc},
1903 {"pow", math_pow, METH_VARARGS, math_pow_doc},
1904 {"radians", math_radians, METH_O, math_radians_doc},
1905 {"sin", math_sin, METH_O, math_sin_doc},
1906 {"sinh", math_sinh, METH_O, math_sinh_doc},
1907 {"sqrt", math_sqrt, METH_O, math_sqrt_doc},
1908 {"tan", math_tan, METH_O, math_tan_doc},
1909 {"tanh", math_tanh, METH_O, math_tanh_doc},
1910 {"trunc", math_trunc, METH_O, math_trunc_doc},
1911 {NULL, NULL} /* sentinel */
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001912};
1913
Guido van Rossumc6e22901998-12-04 19:26:43 +00001914
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +00001915PyDoc_STRVAR(module_doc,
Tim Peters63c94532001-09-04 23:17:42 +00001916"This module is always available. It provides access to the\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +00001917"mathematical functions defined by the C standard.");
Guido van Rossumc6e22901998-12-04 19:26:43 +00001918
Martin v. Löwis1a214512008-06-11 05:26:20 +00001919
1920static struct PyModuleDef mathmodule = {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001921 PyModuleDef_HEAD_INIT,
1922 "math",
1923 module_doc,
1924 -1,
1925 math_methods,
1926 NULL,
1927 NULL,
1928 NULL,
1929 NULL
Martin v. Löwis1a214512008-06-11 05:26:20 +00001930};
1931
Mark Hammondfe51c6d2002-08-02 02:27:13 +00001932PyMODINIT_FUNC
Martin v. Löwis1a214512008-06-11 05:26:20 +00001933PyInit_math(void)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001934{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001935 PyObject *m;
Tim Petersfe71f812001-08-07 22:10:00 +00001936
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001937 m = PyModule_Create(&mathmodule);
1938 if (m == NULL)
1939 goto finally;
Barry Warsawfc93f751996-12-17 00:47:03 +00001940
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001941 PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
1942 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
Barry Warsawfc93f751996-12-17 00:47:03 +00001943
Christian Heimes53876d92008-04-19 00:31:39 +00001944 finally:
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001945 return m;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001946}