1
/* Math module -- standard C math library functions, pi and e */
3
/* 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
9
These are the "spirit of 754" rules:
11
1. If the mathematical result is a real number, but of magnitude too
12
large to approximate by a machine float, overflow is signaled and the
13
result is an infinity (with the appropriate sign).
15
2. If the mathematical result is a real number, but of magnitude too
16
small to approximate by a machine float, underflow is signaled and the
17
result is a zero (with the appropriate sign).
19
3. At a singularity (a value x such that the limit of f(y) as y
20
approaches x exists and is an infinity), "divide by zero" is signaled
21
and the result is an infinity (with the appropriate sign). This is
22
complicated a little by that the left-side and right-side limits may
23
not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0
24
from the positive or negative directions. In that specific case, the
25
sign of the zero determines the result of 1/0.
27
4. At a point where a function has no defined result in the extended
28
reals (i.e., the reals plus an infinity or two), invalid operation is
29
signaled and a NaN is returned.
31
And these are what Python has historically /tried/ to do (but not
32
always successfully, as platform libm behavior varies a lot):
34
For #1, raise OverflowError.
36
For #2, return a zero (with the appropriate sign if that happens by
39
For #3 and #4, raise ValueError. It may have made sense to raise
40
Python's ZeroDivisionError in #3, but historically that's only been
41
raised for division by zero and mod by zero.
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
59
sin(pi*x), giving accurate results for all finite x (especially x
60
integral or close to an integer). This is here for use in the
61
reflection formula for the gamma function. It conforms to IEEE
62
754-2008 for finite arguments, but not for infinities or nans.
65
static const double pi = 3.141592653589793238462643383279502884197;
66
static const double sqrtpi = 1.772453850905516027298167483341145182798;
67
static const double logpi = 1.144729885849400174143427351353058711647;
74
/* this function should only ever be called for finite arguments */
75
assert(Py_IS_FINITE(x));
76
y = fmod(fabs(x), 2.0);
77
n = (int)round(2.0*y);
78
assert(0 <= n && n <= 4);
87
/* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give
88
-0.0 instead of 0.0 when y == 1.0. */
98
assert(0); /* should never get here */
99
r = -1.23e200; /* silence gcc warning */
101
return copysign(1.0, x)*r;
104
/* Implementation of the real gamma function. In extensive but non-exhaustive
105
random tests, this function proved accurate to within <= 10 ulps across the
106
entire float domain. Note that accuracy may depend on the quality of the
107
system math functions, the pow function in particular. Special cases
108
follow C99 annex F. The parameters and method are tailored to platforms
109
whose double format is the IEEE 754 binary64 format.
111
Method: for x > 0.0 we use the Lanczos approximation with parameters N=13
112
and g=6.024680040776729583740234375; these parameters are amongst those
113
used by the Boost library. Following Boost (again), we re-express the
114
Lanczos sum as a rational function, and compute it that way. The
115
coefficients below were computed independently using MPFR, and have been
116
double-checked against the coefficients in the Boost source code.
118
For x < 0.0 we use the reflection formula.
120
There's one minor tweak that deserves explanation: Lanczos' formula for
121
Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x
122
values, x+g-0.5 can be represented exactly. However, in cases where it
123
can't be represented exactly the small error in x+g-0.5 can be magnified
124
significantly by the pow and exp calls, especially for large x. A cheap
125
correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error
126
involved in the computation of x+g-0.5 (that is, e = computed value of
127
x+g-0.5 - exact value of x+g-0.5). Here's the proof:
131
Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754
132
double, and e is tiny. Then:
134
pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e)
135
= pow(y, x-0.5)/exp(y) * C,
137
where the correction_factor C is given by
139
C = pow(1-e/y, x-0.5) * exp(e)
141
Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so:
143
C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y
145
But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and
147
pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y),
149
Note that for accuracy, when computing r*C it's better to do
157
since the addition in the latter throws away most of the bits of
158
information in e*g/y.
162
static const double lanczos_g = 6.024680040776729583740234375;
163
static const double lanczos_g_minus_half = 5.524680040776729583740234375;
164
static const double lanczos_num_coeffs[LANCZOS_N] = {
165
23531376880.410759688572007674451636754734846804940,
166
42919803642.649098768957899047001988850926355848959,
167
35711959237.355668049440185451547166705960488635843,
168
17921034426.037209699919755754458931112671403265390,
169
6039542586.3520280050642916443072979210699388420708,
170
1439720407.3117216736632230727949123939715485786772,
171
248874557.86205415651146038641322942321632125127801,
172
31426415.585400194380614231628318205362874684987640,
173
2876370.6289353724412254090516208496135991145378768,
174
186056.26539522349504029498971604569928220784236328,
175
8071.6720023658162106380029022722506138218516325024,
176
210.82427775157934587250973392071336271166969580291,
177
2.5066282746310002701649081771338373386264310793408
180
/* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */
181
static const double lanczos_den_coeffs[LANCZOS_N] = {
182
0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
183
13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};
185
/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */
186
#define NGAMMA_INTEGRAL 23
187
static const double gamma_integral[NGAMMA_INTEGRAL] = {
188
1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
189
3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
190
1307674368000.0, 20922789888000.0, 355687428096000.0,
191
6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
192
51090942171709440000.0, 1124000727777607680000.0,
195
/* Lanczos' sum L_g(x), for positive x */
198
lanczos_sum(double x)
200
double num = 0.0, den = 0.0;
203
/* evaluate the rational function lanczos_sum(x). For large
204
x, the obvious algorithm risks overflow, so we instead
205
rescale the denominator and numerator of the rational
206
function by x**(1-LANCZOS_N) and treat this as a
207
rational function in 1/x. This also reduces the error for
208
larger x values. The choice of cutoff point (5.0 below) is
209
somewhat arbitrary; in tests, smaller cutoff values than
210
this resulted in lower accuracy. */
212
for (i = LANCZOS_N; --i >= 0; ) {
213
num = num * x + lanczos_num_coeffs[i];
214
den = den * x + lanczos_den_coeffs[i];
218
for (i = 0; i < LANCZOS_N; i++) {
219
num = num / x + lanczos_num_coeffs[i];
220
den = den / x + lanczos_den_coeffs[i];
229
double absx, r, y, z, sqrtpow;
232
if (!Py_IS_FINITE(x)) {
233
if (Py_IS_NAN(x) || x > 0.0)
234
return x; /* tgamma(nan) = nan, tgamma(inf) = inf */
237
return Py_NAN; /* tgamma(-inf) = nan, invalid */
242
/* tgamma(+-0.0) = +-inf, divide-by-zero */
243
return copysign(Py_HUGE_VAL, x);
246
/* integer arguments */
249
errno = EDOM; /* tgamma(n) = nan, invalid for */
250
return Py_NAN; /* negative integers n */
252
if (x <= NGAMMA_INTEGRAL)
253
return gamma_integral[(int)x - 1];
257
/* tiny arguments: tgamma(x) ~ 1/x for x near 0 */
260
if (Py_IS_INFINITY(r))
265
/* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for
266
x > 200, and underflows to +-0.0 for x < -200, not a negative
278
y = absx + lanczos_g_minus_half;
279
/* compute error in sum */
280
if (absx > lanczos_g_minus_half) {
281
/* note: the correction can be foiled by an optimizing
282
compiler that (incorrectly) thinks that an expression like
283
a + b - a - b can be optimized to 0.0. This shouldn't
284
happen in a standards-conforming compiler. */
286
z = q - lanczos_g_minus_half;
289
double q = y - lanczos_g_minus_half;
292
z = z * lanczos_g / y;
294
r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx);
297
r /= pow(y, absx - 0.5);
300
sqrtpow = pow(y, absx / 2.0 - 0.25);
306
r = lanczos_sum(absx) / exp(y);
309
r *= pow(y, absx - 0.5);
312
sqrtpow = pow(y, absx / 2.0 - 0.25);
317
if (Py_IS_INFINITY(r))
323
lgamma: natural log of the absolute value of the Gamma function.
324
For large arguments, Lanczos' formula works extremely well here.
333
if (!Py_IS_FINITE(x)) {
335
return x; /* lgamma(nan) = nan */
337
return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */
340
/* integer arguments */
341
if (x == floor(x) && x <= 2.0) {
343
errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */
344
return Py_HUGE_VAL; /* integers n <= 0 */
347
return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */
352
/* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */
356
/* Lanczos' formula. We could save a fraction of a ulp in accuracy by
357
having a second set of numerator coefficients for lanczos_sum that
358
absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g
359
subtraction below; it's probably not worth it. */
360
r = log(lanczos_sum(absx)) - lanczos_g;
361
r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1);
363
/* Use reflection formula to get value for negative x. */
364
r = logpi - log(fabs(sinpi(absx))) - log(absx) - r;
365
if (Py_IS_INFINITY(r))
371
Implementations of the error function erf(x) and the complementary error
374
Method: following 'Numerical Recipes' by Flannery, Press et. al. (2nd ed.,
375
Cambridge University Press), we use a series approximation for erf for
376
small x, and a continued fraction approximation for erfc(x) for larger x;
377
combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x),
378
this gives us erf(x) and erfc(x) for all x.
380
The series expansion used is:
382
erf(x) = x*exp(-x*x)/sqrt(pi) * [
383
2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...]
385
The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k).
386
This series converges well for smallish x, but slowly for larger x.
388
The continued fraction expansion used is:
390
erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - )
391
3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...]
393
after the first term, the general term has the form:
395
k*(k-0.5)/(2*k+0.5 + x**2 - ...).
397
This expansion converges fast for larger x, but convergence becomes
398
infinitely slow as x approaches 0.0. The (somewhat naive) continued
399
fraction evaluation algorithm used below also risks overflow for large x;
400
but for large x, erfc(x) == 0.0 to within machine precision. (For
401
example, erfc(30.0) is approximately 2.56e-393).
403
Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and
404
continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) <
405
ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the
406
numbers of terms to use for the relevant expansions. */
408
#define ERF_SERIES_CUTOFF 1.5
409
#define ERF_SERIES_TERMS 25
410
#define ERFC_CONTFRAC_CUTOFF 30.0
411
#define ERFC_CONTFRAC_TERMS 50
414
Error function, via power series.
416
Given a finite float x, return an approximation to erf(x).
417
Converges reasonably fast for small x.
421
m_erf_series(double x)
423
double x2, acc, fk, result;
428
fk = (double)ERF_SERIES_TERMS + 0.5;
429
for (i = 0; i < ERF_SERIES_TERMS; i++) {
430
acc = 2.0 + x2 * acc / fk;
433
/* Make sure the exp call doesn't affect errno;
434
see m_erfc_contfrac for more. */
436
result = acc * x * exp(-x2) / sqrtpi;
442
Complementary error function, via continued fraction expansion.
444
Given a positive float x, return an approximation to erfc(x). Converges
445
reasonably fast for x large (say, x > 2.0), and should be safe from
446
overflow if x and nterms are not too large. On an IEEE 754 machine, with x
447
<= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller
448
than the smallest representable nonzero float. */
451
m_erfc_contfrac(double x)
453
double x2, a, da, p, p_last, q, q_last, b, result;
456
if (x >= ERFC_CONTFRAC_CUTOFF)
462
p = 1.0; p_last = 0.0;
463
q = da + x2; q_last = 1.0;
464
for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) {
469
temp = p; p = b*p - a*p_last; p_last = temp;
470
temp = q; q = b*q - a*q_last; q_last = temp;
472
/* Issue #8986: On some platforms, exp sets errno on underflow to zero;
473
save the current errno value so that we can restore it later. */
475
result = p / q * x * exp(-x2) / sqrtpi;
480
/* Error function erf(x), for general x */
490
if (absx < ERF_SERIES_CUTOFF)
491
return m_erf_series(x);
493
cf = m_erfc_contfrac(absx);
494
return x > 0.0 ? 1.0 - cf : cf - 1.0;
498
/* Complementary error function erfc(x), for general x. */
508
if (absx < ERF_SERIES_CUTOFF)
509
return 1.0 - m_erf_series(x);
511
cf = m_erfc_contfrac(absx);
512
return x > 0.0 ? cf : 2.0 - cf;
517
wrapper for atan2 that deals directly with special cases before
518
delegating to the platform libm for the remaining cases. This
519
is necessary to get consistent behaviour across platforms.
520
Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
525
m_atan2(double y, double x)
527
if (Py_IS_NAN(x) || Py_IS_NAN(y))
529
if (Py_IS_INFINITY(y)) {
530
if (Py_IS_INFINITY(x)) {
531
if (copysign(1., x) == 1.)
532
/* atan2(+-inf, +inf) == +-pi/4 */
533
return copysign(0.25*Py_MATH_PI, y);
535
/* atan2(+-inf, -inf) == +-pi*3/4 */
536
return copysign(0.75*Py_MATH_PI, y);
538
/* atan2(+-inf, x) == +-pi/2 for finite x */
539
return copysign(0.5*Py_MATH_PI, y);
541
if (Py_IS_INFINITY(x) || y == 0.) {
542
if (copysign(1., x) == 1.)
543
/* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
544
return copysign(0., y);
546
/* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
547
return copysign(Py_MATH_PI, y);
553
Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
554
log(-ve), log(NaN). Here are wrappers for log and log10 that deal with
555
special values directly, passing positive non-special values through to
556
the system log/log10.
562
if (Py_IS_FINITE(x)) {
567
return -Py_HUGE_VAL; /* log(0) = -inf */
569
return Py_NAN; /* log(-ve) = nan */
571
else if (Py_IS_NAN(x))
572
return x; /* log(nan) = nan */
574
return x; /* log(inf) = inf */
577
return Py_NAN; /* log(-inf) = nan */
584
Uses an algorithm that should:
586
(a) produce exact results for powers of 2, and
587
(b) give a monotonic log2 (for positive finite floats),
588
assuming that the system log is monotonic.
594
if (!Py_IS_FINITE(x)) {
596
return x; /* log2(nan) = nan */
598
return x; /* log2(+inf) = +inf */
601
return Py_NAN; /* log2(-inf) = nan, invalid-operation */
612
/* We want log2(m * 2**e) == log(m) / log(2) + e. Care is needed when
613
* x is just greater than 1.0: in that case e is 1, log(m) is negative,
614
* and we get significant cancellation error from the addition of
615
* log(m) / log(2) to e. The slight rewrite of the expression below
616
* avoids this problem.
619
return log(2.0 * m) / log(2.0) + (e - 1);
622
return log(m) / log(2.0) + e;
628
return -Py_HUGE_VAL; /* log2(0) = -inf, divide-by-zero */
632
return Py_NAN; /* log2(-inf) = nan, invalid-operation */
639
if (Py_IS_FINITE(x)) {
644
return -Py_HUGE_VAL; /* log10(0) = -inf */
646
return Py_NAN; /* log10(-ve) = nan */
648
else if (Py_IS_NAN(x))
649
return x; /* log10(nan) = nan */
651
return x; /* log10(inf) = inf */
654
return Py_NAN; /* log10(-inf) = nan */
659
/* Call is_error when errno != 0, and where x is the result libm
660
* returned. is_error will usually set up an exception and return
661
* true (1), but may return false (0) without setting up an exception.
666
int result = 1; /* presumption of guilt */
667
assert(errno); /* non-zero errno is a precondition for calling */
669
PyErr_SetString(PyExc_ValueError, "math domain error");
671
else if (errno == ERANGE) {
672
/* ANSI C generally requires libm functions to set ERANGE
673
* on overflow, but also generally *allows* them to set
674
* ERANGE on underflow too. There's no consistency about
675
* the latter across platforms.
676
* Alas, C99 never requires that errno be set.
677
* Here we suppress the underflow errors (libm functions
678
* should return a zero on underflow, and +- HUGE_VAL on
679
* overflow, so testing the result for zero suffices to
680
* distinguish the cases).
682
* On some platforms (Ubuntu/ia64) it seems that errno can be
683
* set to ERANGE for subnormal results that do *not* underflow
684
* to zero. So to be safe, we'll ignore ERANGE whenever the
685
* function result is less than one in absolute value.
690
PyErr_SetString(PyExc_OverflowError,
694
/* Unexpected math error */
695
PyErr_SetFromErrno(PyExc_ValueError);
700
math_1 is used to wrap a libm function f that takes a double
701
arguments and returns a double.
703
The error reporting follows these rules, which are designed to do
704
the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
707
- a NaN result from non-NaN inputs causes ValueError to be raised
708
- an infinite result from finite inputs causes OverflowError to be
709
raised if can_overflow is 1, or raises ValueError if can_overflow
711
- if the result is finite and errno == EDOM then ValueError is
713
- if the result is finite and nonzero and errno == ERANGE then
714
OverflowError is raised
716
The last rule is used to catch overflow on platforms which follow
717
C89 but for which HUGE_VAL is not an infinity.
719
For the majority of one-argument functions these rules are enough
720
to ensure that Python's functions behave as specified in 'Annex F'
721
of the C99 standard, with the 'invalid' and 'divide-by-zero'
722
floating-point exceptions mapping to Python's ValueError and the
723
'overflow' floating-point exception mapping to OverflowError.
724
math_1 only works for functions that don't have singularities *and*
725
the possibility of overflow; fortunately, that covers everything we
726
care about right now.
730
math_1_to_whatever(PyObject *arg, double (*func) (double),
731
PyObject *(*from_double_func) (double),
735
x = PyFloat_AsDouble(arg);
736
if (x == -1.0 && PyErr_Occurred())
739
PyFPE_START_PROTECT("in math_1", return 0);
741
PyFPE_END_PROTECT(r);
742
if (Py_IS_NAN(r) && !Py_IS_NAN(x)) {
743
PyErr_SetString(PyExc_ValueError,
744
"math domain error"); /* invalid arg */
747
if (Py_IS_INFINITY(r) && Py_IS_FINITE(x)) {
749
PyErr_SetString(PyExc_OverflowError,
750
"math range error"); /* overflow */
752
PyErr_SetString(PyExc_ValueError,
753
"math domain error"); /* singularity */
756
if (Py_IS_FINITE(r) && errno && is_error(r))
757
/* this branch unnecessary on most platforms */
760
return (*from_double_func)(r);
763
/* variant of math_1, to be used when the function being wrapped is known to
764
set errno properly (that is, errno = EDOM for invalid or divide-by-zero,
765
errno = ERANGE for overflow). */
768
math_1a(PyObject *arg, double (*func) (double))
771
x = PyFloat_AsDouble(arg);
772
if (x == -1.0 && PyErr_Occurred())
775
PyFPE_START_PROTECT("in math_1a", return 0);
777
PyFPE_END_PROTECT(r);
778
if (errno && is_error(r))
780
return PyFloat_FromDouble(r);
784
math_2 is used to wrap a libm function f that takes two double
785
arguments and returns a double.
787
The error reporting follows these rules, which are designed to do
788
the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
791
- a NaN result from non-NaN inputs causes ValueError to be raised
792
- an infinite result from finite inputs causes OverflowError to be
794
- if the result is finite and errno == EDOM then ValueError is
796
- if the result is finite and nonzero and errno == ERANGE then
797
OverflowError is raised
799
The last rule is used to catch overflow on platforms which follow
800
C89 but for which HUGE_VAL is not an infinity.
802
For most two-argument functions (copysign, fmod, hypot, atan2)
803
these rules are enough to ensure that Python's functions behave as
804
specified in 'Annex F' of the C99 standard, with the 'invalid' and
805
'divide-by-zero' floating-point exceptions mapping to Python's
806
ValueError and the 'overflow' floating-point exception mapping to
811
math_1(PyObject *arg, double (*func) (double), int can_overflow)
813
return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow);
817
math_1_to_int(PyObject *arg, double (*func) (double), int can_overflow)
819
return math_1_to_whatever(arg, func, PyLong_FromDouble, can_overflow);
823
math_2(PyObject *args, double (*func) (double, double), char *funcname)
827
if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
829
x = PyFloat_AsDouble(ox);
830
y = PyFloat_AsDouble(oy);
831
if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
834
PyFPE_START_PROTECT("in math_2", return 0);
836
PyFPE_END_PROTECT(r);
838
if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
843
else if (Py_IS_INFINITY(r)) {
844
if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
849
if (errno && is_error(r))
852
return PyFloat_FromDouble(r);
855
#define FUNC1(funcname, func, can_overflow, docstring) \
856
static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
857
return math_1(args, func, can_overflow); \
859
PyDoc_STRVAR(math_##funcname##_doc, docstring);
861
#define FUNC1A(funcname, func, docstring) \
862
static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
863
return math_1a(args, func); \
865
PyDoc_STRVAR(math_##funcname##_doc, docstring);
867
#define FUNC2(funcname, func, docstring) \
868
static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
869
return math_2(args, func, #funcname); \
871
PyDoc_STRVAR(math_##funcname##_doc, docstring);
874
"acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
875
FUNC1(acosh, m_acosh, 0,
876
"acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
878
"asin(x)\n\nReturn the arc sine (measured in radians) of x.")
879
FUNC1(asinh, m_asinh, 0,
880
"asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
882
"atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
883
FUNC2(atan2, m_atan2,
884
"atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
885
"Unlike atan(y/x), the signs of both x and y are considered.")
886
FUNC1(atanh, m_atanh, 0,
887
"atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
889
static PyObject * math_ceil(PyObject *self, PyObject *number) {
890
_Py_IDENTIFIER(__ceil__);
891
PyObject *method, *result;
893
method = _PyObject_LookupSpecial(number, &PyId___ceil__);
894
if (method == NULL) {
895
if (PyErr_Occurred())
897
return math_1_to_int(number, ceil, 0);
899
result = PyObject_CallFunctionObjArgs(method, NULL);
904
PyDoc_STRVAR(math_ceil_doc,
905
"ceil(x)\n\nReturn the ceiling of x as an int.\n"
906
"This is the smallest integral value >= x.");
908
FUNC2(copysign, copysign,
909
"copysign(x, y)\n\nReturn x with the sign of y.")
911
"cos(x)\n\nReturn the cosine of x (measured in radians).")
913
"cosh(x)\n\nReturn the hyperbolic cosine of x.")
915
"erf(x)\n\nError function at x.")
917
"erfc(x)\n\nComplementary error function at x.")
919
"exp(x)\n\nReturn e raised to the power of x.")
920
FUNC1(expm1, m_expm1, 1,
921
"expm1(x)\n\nReturn exp(x)-1.\n"
922
"This function avoids the loss of precision involved in the direct "
923
"evaluation of exp(x)-1 for small x.")
925
"fabs(x)\n\nReturn the absolute value of the float x.")
927
static PyObject * math_floor(PyObject *self, PyObject *number) {
928
_Py_IDENTIFIER(__floor__);
929
PyObject *method, *result;
931
method = _PyObject_LookupSpecial(number, &PyId___floor__);
932
if (method == NULL) {
933
if (PyErr_Occurred())
935
return math_1_to_int(number, floor, 0);
937
result = PyObject_CallFunctionObjArgs(method, NULL);
942
PyDoc_STRVAR(math_floor_doc,
943
"floor(x)\n\nReturn the floor of x as an int.\n"
944
"This is the largest integral value <= x.");
946
FUNC1A(gamma, m_tgamma,
947
"gamma(x)\n\nGamma function at x.")
948
FUNC1A(lgamma, m_lgamma,
949
"lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.")
950
FUNC1(log1p, m_log1p, 0,
951
"log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n"
952
"The result is computed in a way which is accurate for x near zero.")
954
"sin(x)\n\nReturn the sine of x (measured in radians).")
956
"sinh(x)\n\nReturn the hyperbolic sine of x.")
958
"sqrt(x)\n\nReturn the square root of x.")
960
"tan(x)\n\nReturn the tangent of x (measured in radians).")
962
"tanh(x)\n\nReturn the hyperbolic tangent of x.")
964
/* Precision summation function as msum() by Raymond Hettinger in
965
<http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
966
enhanced with the exact partials sum and roundoff from Mark
967
Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
968
See those links for more details, proofs and other references.
970
Note 1: IEEE 754R floating point semantics are assumed,
971
but the current implementation does not re-establish special
972
value semantics across iterations (i.e. handling -Inf + Inf).
974
Note 2: No provision is made for intermediate overflow handling;
975
therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while
976
sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
977
overflow of the first partial sum.
979
Note 3: The intermediate values lo, yr, and hi are declared volatile so
980
aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
981
Also, the volatile declaration forces the values to be stored in memory as
982
regular doubles instead of extended long precision (80-bit) values. This
983
prevents double rounding because any addition or subtraction of two doubles
984
can be resolved exactly into double-sized hi and lo values. As long as the
985
hi value gets forced into a double before yr and lo are computed, the extra
986
bits in downstream extended precision operations (x87 for example) will be
987
exactly zero and therefore can be losslessly stored back into a double,
988
thereby preventing double rounding.
990
Note 4: A similar implementation is in Modules/cmathmodule.c.
991
Be sure to update both when making changes.
993
Note 5: The signature of math.fsum() differs from __builtin__.sum()
994
because the start argument doesn't make sense in the context of
995
accurate summation. Since the partials table is collapsed before
996
returning a result, sum(seq2, start=sum(seq1)) may not equal the
997
accurate result returned by sum(itertools.chain(seq1, seq2)).
1000
#define NUM_PARTIALS 32 /* initial partials array size, on stack */
1002
/* Extend the partials array p[] by doubling its size. */
1003
static int /* non-zero on error */
1004
_fsum_realloc(double **p_ptr, Py_ssize_t n,
1005
double *ps, Py_ssize_t *m_ptr)
1008
Py_ssize_t m = *m_ptr;
1010
m += m; /* double */
1011
if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
1014
v = PyMem_Malloc(sizeof(double) * m);
1016
memcpy(v, ps, sizeof(double) * n);
1019
v = PyMem_Realloc(p, sizeof(double) * m);
1021
if (v == NULL) { /* size overflow or no memory */
1022
PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
1025
*p_ptr = (double*) v;
1030
/* Full precision summation of a sequence of floats.
1033
partials = [] # sorted, non-overlapping partial sums
1046
return sum_exact(partials)
1048
Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo
1049
are exactly equal to x+y. The inner loop applies hi/lo summation to each
1050
partial so that the list of partial sums remains exact.
1052
Sum_exact() adds the partial sums exactly and correctly rounds the final
1053
result (using the round-half-to-even rule). The items in partials remain
1054
non-zero, non-special, non-overlapping and strictly increasing in
1055
magnitude, but possibly not all having the same sign.
1057
Depends on IEEE 754 arithmetic guarantees and half-even rounding.
1061
math_fsum(PyObject *self, PyObject *seq)
1063
PyObject *item, *iter, *sum = NULL;
1064
Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
1065
double x, y, t, ps[NUM_PARTIALS], *p = ps;
1066
double xsave, special_sum = 0.0, inf_sum = 0.0;
1067
volatile double hi, yr, lo;
1069
iter = PyObject_GetIter(seq);
1073
PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
1075
for(;;) { /* for x in iterable */
1076
assert(0 <= n && n <= m);
1077
assert((m == NUM_PARTIALS && p == ps) ||
1078
(m > NUM_PARTIALS && p != NULL));
1080
item = PyIter_Next(iter);
1082
if (PyErr_Occurred())
1086
x = PyFloat_AsDouble(item);
1088
if (PyErr_Occurred())
1092
for (i = j = 0; j < n; j++) { /* for y in partials */
1094
if (fabs(x) < fabs(y)) {
1095
t = x; x = y; y = t;
1105
n = i; /* ps[i:] = [x] */
1107
if (! Py_IS_FINITE(x)) {
1108
/* a nonfinite x could arise either as
1109
a result of intermediate overflow, or
1110
as a result of a nan or inf in the
1112
if (Py_IS_FINITE(xsave)) {
1113
PyErr_SetString(PyExc_OverflowError,
1114
"intermediate overflow in fsum");
1117
if (Py_IS_INFINITY(xsave))
1119
special_sum += xsave;
1120
/* reset partials */
1123
else if (n >= m && _fsum_realloc(&p, n, ps, &m))
1130
if (special_sum != 0.0) {
1131
if (Py_IS_NAN(inf_sum))
1132
PyErr_SetString(PyExc_ValueError,
1133
"-inf + inf in fsum");
1135
sum = PyFloat_FromDouble(special_sum);
1142
/* sum_exact(ps, hi) from the top, stop when the sum becomes
1147
assert(fabs(y) < fabs(x));
1154
/* Make half-even rounding work across multiple partials.
1155
Needed so that sum([1e-16, 1, 1e16]) will round-up the last
1156
digit to two instead of down to zero (the 1e-16 makes the 1
1157
slightly closer to two). With a potential 1 ULP rounding
1158
error fixed-up, math.fsum() can guarantee commutativity. */
1159
if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
1160
(lo > 0.0 && p[n-1] > 0.0))) {
1168
sum = PyFloat_FromDouble(hi);
1171
PyFPE_END_PROTECT(hi)
1180
PyDoc_STRVAR(math_fsum_doc,
1181
"fsum(iterable)\n\n\
1182
Return an accurate floating point sum of values in the iterable.\n\
1183
Assumes IEEE-754 floating point arithmetic.");
1185
/* Return the smallest integer k such that n < 2**k, or 0 if n == 0.
1186
* Equivalent to floor(lg(x))+1. Also equivalent to: bitwidth_of_type -
1187
* count_leading_zero_bits(x)
1190
/* XXX: This routine does more or less the same thing as
1191
* bits_in_digit() in Objects/longobject.c. Someday it would be nice to
1192
* consolidate them. On BSD, there's a library function called fls()
1193
* that we could use, and GCC provides __builtin_clz().
1196
static unsigned long
1197
bit_length(unsigned long n)
1199
unsigned long len = 0;
1207
static unsigned long
1208
count_set_bits(unsigned long n)
1210
unsigned long count = 0;
1213
n &= n - 1; /* clear least significant bit */
1218
/* Divide-and-conquer factorial algorithm
1220
* Based on the formula and psuedo-code provided at:
1221
* http://www.luschny.de/math/factorial/binarysplitfact.html
1223
* Faster algorithms exist, but they're more complicated and depend on
1224
* a fast prime factorization algorithm.
1226
* Notes on the algorithm
1227
* ----------------------
1229
* factorial(n) is written in the form 2**k * m, with m odd. k and m are
1230
* computed separately, and then combined using a left shift.
1232
* The function factorial_odd_part computes the odd part m (i.e., the greatest
1233
* odd divisor) of factorial(n), using the formula:
1235
* factorial_odd_part(n) =
1237
* product_{i >= 0} product_{0 < j <= n / 2**i, j odd} j
1239
* Example: factorial_odd_part(20) =
1244
* (1 * 3 * 5 * 7 * 9)
1245
* (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1247
* Here i goes from large to small: the first term corresponds to i=4 (any
1248
* larger i gives an empty product), and the last term corresponds to i=0.
1249
* Each term can be computed from the last by multiplying by the extra odd
1250
* numbers required: e.g., to get from the penultimate term to the last one,
1251
* we multiply by (11 * 13 * 15 * 17 * 19).
1253
* To see a hint of why this formula works, here are the same numbers as above
1254
* but with the even parts (i.e., the appropriate powers of 2) included. For
1255
* each subterm in the product for i, we multiply that subterm by 2**i:
1262
* (2 * 6 * 10 * 14 * 18) *
1263
* (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1265
* The factorial_partial_product function computes the product of all odd j in
1266
* range(start, stop) for given start and stop. It's used to compute the
1267
* partial products like (11 * 13 * 15 * 17 * 19) in the example above. It
1268
* operates recursively, repeatedly splitting the range into two roughly equal
1269
* pieces until the subranges are small enough to be computed using only C
1270
* integer arithmetic.
1272
* The two-valuation k (i.e., the exponent of the largest power of 2 dividing
1273
* the factorial) is computed independently in the main math_factorial
1274
* function. By standard results, its value is:
1276
* two_valuation = n//2 + n//4 + n//8 + ....
1278
* It can be shown (e.g., by complete induction on n) that two_valuation is
1279
* equal to n - count_set_bits(n), where count_set_bits(n) gives the number of
1280
* '1'-bits in the binary expansion of n.
1283
/* factorial_partial_product: Compute product(range(start, stop, 2)) using
1284
* divide and conquer. Assumes start and stop are odd and stop > start.
1285
* max_bits must be >= bit_length(stop - 2). */
1288
factorial_partial_product(unsigned long start, unsigned long stop,
1289
unsigned long max_bits)
1291
unsigned long midpoint, num_operands;
1292
PyObject *left = NULL, *right = NULL, *result = NULL;
1294
/* If the return value will fit an unsigned long, then we can
1295
* multiply in a tight, fast loop where each multiply is O(1).
1296
* Compute an upper bound on the number of bits required to store
1299
* Storing some integer z requires floor(lg(z))+1 bits, which is
1300
* conveniently the value returned by bit_length(z). The
1301
* product x*y will require at most
1302
* bit_length(x) + bit_length(y) bits to store, based
1303
* on the idea that lg product = lg x + lg y.
1305
* We know that stop - 2 is the largest number to be multiplied. From
1306
* there, we have: bit_length(answer) <= num_operands *
1307
* bit_length(stop - 2)
1310
num_operands = (stop - start) / 2;
1311
/* The "num_operands <= 8 * SIZEOF_LONG" check guards against the
1312
* unlikely case of an overflow in num_operands * max_bits. */
1313
if (num_operands <= 8 * SIZEOF_LONG &&
1314
num_operands * max_bits <= 8 * SIZEOF_LONG) {
1315
unsigned long j, total;
1316
for (total = start, j = start + 2; j < stop; j += 2)
1318
return PyLong_FromUnsignedLong(total);
1321
/* find midpoint of range(start, stop), rounded up to next odd number. */
1322
midpoint = (start + num_operands) | 1;
1323
left = factorial_partial_product(start, midpoint,
1324
bit_length(midpoint - 2));
1327
right = factorial_partial_product(midpoint, stop, max_bits);
1330
result = PyNumber_Multiply(left, right);
1338
/* factorial_odd_part: compute the odd part of factorial(n). */
1341
factorial_odd_part(unsigned long n)
1344
unsigned long v, lower, upper;
1345
PyObject *partial, *tmp, *inner, *outer;
1347
inner = PyLong_FromLong(1);
1354
for (i = bit_length(n) - 2; i >= 0; i--) {
1359
/* (v + 1) | 1 = least odd integer strictly larger than n / 2**i */
1360
upper = (v + 1) | 1;
1361
/* Here inner is the product of all odd integers j in the range (0,
1362
n/2**(i+1)]. The factorial_partial_product call below gives the
1363
product of all odd integers j in the range (n/2**(i+1), n/2**i]. */
1364
partial = factorial_partial_product(lower, upper, bit_length(upper-2));
1365
/* inner *= partial */
1366
if (partial == NULL)
1368
tmp = PyNumber_Multiply(inner, partial);
1374
/* Now inner is the product of all odd integers j in the range (0,
1375
n/2**i], giving the inner product in the formula above. */
1377
/* outer *= inner; */
1378
tmp = PyNumber_Multiply(outer, inner);
1393
/* Lookup table for small factorial values */
1395
static const unsigned long SmallFactorials[] = {
1396
1, 1, 2, 6, 24, 120, 720, 5040, 40320,
1397
362880, 3628800, 39916800, 479001600,
1398
#if SIZEOF_LONG >= 8
1399
6227020800, 87178291200, 1307674368000,
1400
20922789888000, 355687428096000, 6402373705728000,
1401
121645100408832000, 2432902008176640000
1406
math_factorial(PyObject *self, PyObject *arg)
1409
PyObject *result, *odd_part, *two_valuation;
1411
if (PyFloat_Check(arg)) {
1413
double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);
1414
if (!(Py_IS_FINITE(dx) && dx == floor(dx))) {
1415
PyErr_SetString(PyExc_ValueError,
1416
"factorial() only accepts integral values");
1419
lx = PyLong_FromDouble(dx);
1422
x = PyLong_AsLong(lx);
1426
x = PyLong_AsLong(arg);
1428
if (x == -1 && PyErr_Occurred())
1431
PyErr_SetString(PyExc_ValueError,
1432
"factorial() not defined for negative values");
1436
/* use lookup table if x is small */
1437
if (x < (long)Py_ARRAY_LENGTH(SmallFactorials))
1438
return PyLong_FromUnsignedLong(SmallFactorials[x]);
1440
/* else express in the form odd_part * 2**two_valuation, and compute as
1441
odd_part << two_valuation. */
1442
odd_part = factorial_odd_part(x);
1443
if (odd_part == NULL)
1445
two_valuation = PyLong_FromLong(x - count_set_bits(x));
1446
if (two_valuation == NULL) {
1447
Py_DECREF(odd_part);
1450
result = PyNumber_Lshift(odd_part, two_valuation);
1451
Py_DECREF(two_valuation);
1452
Py_DECREF(odd_part);
1456
PyDoc_STRVAR(math_factorial_doc,
1457
"factorial(x) -> Integral\n"
1459
"Find x!. Raise a ValueError if x is negative or non-integral.");
1462
math_trunc(PyObject *self, PyObject *number)
1464
_Py_IDENTIFIER(__trunc__);
1465
PyObject *trunc, *result;
1467
if (Py_TYPE(number)->tp_dict == NULL) {
1468
if (PyType_Ready(Py_TYPE(number)) < 0)
1472
trunc = _PyObject_LookupSpecial(number, &PyId___trunc__);
1473
if (trunc == NULL) {
1474
if (!PyErr_Occurred())
1475
PyErr_Format(PyExc_TypeError,
1476
"type %.100s doesn't define __trunc__ method",
1477
Py_TYPE(number)->tp_name);
1480
result = PyObject_CallFunctionObjArgs(trunc, NULL);
1485
PyDoc_STRVAR(math_trunc_doc,
1486
"trunc(x:Real) -> Integral\n"
1488
"Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
1491
math_frexp(PyObject *self, PyObject *arg)
1494
double x = PyFloat_AsDouble(arg);
1495
if (x == -1.0 && PyErr_Occurred())
1497
/* deal with special cases directly, to sidestep platform
1499
if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
1503
PyFPE_START_PROTECT("in math_frexp", return 0);
1505
PyFPE_END_PROTECT(x);
1507
return Py_BuildValue("(di)", x, i);
1510
PyDoc_STRVAR(math_frexp_doc,
1513
"Return the mantissa and exponent of x, as pair (m, e).\n"
1514
"m is a float and e is an int, such that x = m * 2.**e.\n"
1515
"If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0.");
1518
math_ldexp(PyObject *self, PyObject *args)
1524
if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp))
1527
if (PyLong_Check(oexp)) {
1528
/* on overflow, replace exponent with either LONG_MAX
1529
or LONG_MIN, depending on the sign. */
1530
exp = PyLong_AsLongAndOverflow(oexp, &overflow);
1531
if (exp == -1 && PyErr_Occurred())
1534
exp = overflow < 0 ? LONG_MIN : LONG_MAX;
1537
PyErr_SetString(PyExc_TypeError,
1538
"Expected an int as second argument to ldexp.");
1542
if (x == 0. || !Py_IS_FINITE(x)) {
1543
/* NaNs, zeros and infinities are returned unchanged */
1546
} else if (exp > INT_MAX) {
1548
r = copysign(Py_HUGE_VAL, x);
1550
} else if (exp < INT_MIN) {
1551
/* underflow to +-0 */
1552
r = copysign(0., x);
1556
PyFPE_START_PROTECT("in math_ldexp", return 0);
1557
r = ldexp(x, (int)exp);
1558
PyFPE_END_PROTECT(r);
1559
if (Py_IS_INFINITY(r))
1563
if (errno && is_error(r))
1565
return PyFloat_FromDouble(r);
1568
PyDoc_STRVAR(math_ldexp_doc,
1570
Return x * (2**i).");
1573
math_modf(PyObject *self, PyObject *arg)
1575
double y, x = PyFloat_AsDouble(arg);
1576
if (x == -1.0 && PyErr_Occurred())
1578
/* some platforms don't do the right thing for NaNs and
1579
infinities, so we take care of special cases directly. */
1580
if (!Py_IS_FINITE(x)) {
1581
if (Py_IS_INFINITY(x))
1582
return Py_BuildValue("(dd)", copysign(0., x), x);
1583
else if (Py_IS_NAN(x))
1584
return Py_BuildValue("(dd)", x, x);
1588
PyFPE_START_PROTECT("in math_modf", return 0);
1590
PyFPE_END_PROTECT(x);
1591
return Py_BuildValue("(dd)", x, y);
1594
PyDoc_STRVAR(math_modf_doc,
1597
"Return the fractional and integer parts of x. Both results carry the sign\n"
1598
"of x and are floats.");
1600
/* A decent logarithm is easy to compute even for huge ints, but libm can't
1601
do that by itself -- loghelper can. func is log or log10, and name is
1602
"log" or "log10". Note that overflow of the result isn't possible: an int
1603
can contain no more than INT_MAX * SHIFT bits, so has value certainly less
1604
than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
1605
small enough to fit in an IEEE single. log and log10 are even smaller.
1606
However, intermediate overflow is possible for an int if the number of bits
1607
in that int is larger than PY_SSIZE_T_MAX. */
1610
loghelper(PyObject* arg, double (*func)(double), char *funcname)
1612
/* If it is int, do it ourselves. */
1613
if (PyLong_Check(arg)) {
1617
/* Negative or zero inputs give a ValueError. */
1618
if (Py_SIZE(arg) <= 0) {
1619
PyErr_SetString(PyExc_ValueError,
1620
"math domain error");
1624
x = PyLong_AsDouble(arg);
1625
if (x == -1.0 && PyErr_Occurred()) {
1626
if (!PyErr_ExceptionMatches(PyExc_OverflowError))
1628
/* Here the conversion to double overflowed, but it's possible
1629
to compute the log anyway. Clear the exception and continue. */
1631
x = _PyLong_Frexp((PyLongObject *)arg, &e);
1632
if (x == -1.0 && PyErr_Occurred())
1634
/* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */
1635
result = func(x) + func(2.0) * e;
1638
/* Successfully converted x to a double. */
1640
return PyFloat_FromDouble(result);
1643
/* Else let libm handle it by itself. */
1644
return math_1(arg, func, 0);
1648
math_log(PyObject *self, PyObject *args)
1651
PyObject *base = NULL;
1652
PyObject *num, *den;
1655
if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
1658
num = loghelper(arg, m_log, "log");
1659
if (num == NULL || base == NULL)
1662
den = loghelper(base, m_log, "log");
1668
ans = PyNumber_TrueDivide(num, den);
1674
PyDoc_STRVAR(math_log_doc,
1675
"log(x[, base])\n\n\
1676
Return the logarithm of x to the given base.\n\
1677
If the base not specified, returns the natural logarithm (base e) of x.");
1680
math_log2(PyObject *self, PyObject *arg)
1682
return loghelper(arg, m_log2, "log2");
1685
PyDoc_STRVAR(math_log2_doc,
1686
"log2(x)\n\nReturn the base 2 logarithm of x.");
1689
math_log10(PyObject *self, PyObject *arg)
1691
return loghelper(arg, m_log10, "log10");
1694
PyDoc_STRVAR(math_log10_doc,
1695
"log10(x)\n\nReturn the base 10 logarithm of x.");
1698
math_fmod(PyObject *self, PyObject *args)
1702
if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
1704
x = PyFloat_AsDouble(ox);
1705
y = PyFloat_AsDouble(oy);
1706
if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1708
/* fmod(x, +/-Inf) returns x for finite x. */
1709
if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
1710
return PyFloat_FromDouble(x);
1712
PyFPE_START_PROTECT("in math_fmod", return 0);
1714
PyFPE_END_PROTECT(r);
1716
if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1721
if (errno && is_error(r))
1724
return PyFloat_FromDouble(r);
1727
PyDoc_STRVAR(math_fmod_doc,
1728
"fmod(x, y)\n\nReturn fmod(x, y), according to platform C."
1729
" x % y may differ.");
1732
math_hypot(PyObject *self, PyObject *args)
1736
if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
1738
x = PyFloat_AsDouble(ox);
1739
y = PyFloat_AsDouble(oy);
1740
if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1742
/* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
1743
if (Py_IS_INFINITY(x))
1744
return PyFloat_FromDouble(fabs(x));
1745
if (Py_IS_INFINITY(y))
1746
return PyFloat_FromDouble(fabs(y));
1748
PyFPE_START_PROTECT("in math_hypot", return 0);
1750
PyFPE_END_PROTECT(r);
1752
if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1757
else if (Py_IS_INFINITY(r)) {
1758
if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
1763
if (errno && is_error(r))
1766
return PyFloat_FromDouble(r);
1769
PyDoc_STRVAR(math_hypot_doc,
1770
"hypot(x, y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
1772
/* pow can't use math_2, but needs its own wrapper: the problem is
1773
that an infinite result can arise either as a result of overflow
1774
(in which case OverflowError should be raised) or as a result of
1775
e.g. 0.**-5. (for which ValueError needs to be raised.)
1779
math_pow(PyObject *self, PyObject *args)
1785
if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
1787
x = PyFloat_AsDouble(ox);
1788
y = PyFloat_AsDouble(oy);
1789
if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1792
/* deal directly with IEEE specials, to cope with problems on various
1793
platforms whose semantics don't exactly match C99 */
1794
r = 0.; /* silence compiler warning */
1795
if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
1798
r = y == 0. ? 1. : x; /* NaN**0 = 1 */
1799
else if (Py_IS_NAN(y))
1800
r = x == 1. ? 1. : y; /* 1**NaN = 1 */
1801
else if (Py_IS_INFINITY(x)) {
1802
odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
1804
r = odd_y ? x : fabs(x);
1808
r = odd_y ? copysign(0., x) : 0.;
1810
else if (Py_IS_INFINITY(y)) {
1813
else if (y > 0. && fabs(x) > 1.0)
1815
else if (y < 0. && fabs(x) < 1.0) {
1816
r = -y; /* result is +inf */
1817
if (x == 0.) /* 0**-inf: divide-by-zero */
1825
/* let libm handle finite**finite */
1827
PyFPE_START_PROTECT("in math_pow", return 0);
1829
PyFPE_END_PROTECT(r);
1830
/* a NaN result should arise only from (-ve)**(finite
1831
non-integer); in this case we want to raise ValueError. */
1832
if (!Py_IS_FINITE(r)) {
1837
an infinite result here arises either from:
1838
(A) (+/-0.)**negative (-> divide-by-zero)
1839
(B) overflow of x**y with x and y finite
1841
else if (Py_IS_INFINITY(r)) {
1850
if (errno && is_error(r))
1853
return PyFloat_FromDouble(r);
1856
PyDoc_STRVAR(math_pow_doc,
1857
"pow(x, y)\n\nReturn x**y (x to the power of y).");
1859
static const double degToRad = Py_MATH_PI / 180.0;
1860
static const double radToDeg = 180.0 / Py_MATH_PI;
1863
math_degrees(PyObject *self, PyObject *arg)
1865
double x = PyFloat_AsDouble(arg);
1866
if (x == -1.0 && PyErr_Occurred())
1868
return PyFloat_FromDouble(x * radToDeg);
1871
PyDoc_STRVAR(math_degrees_doc,
1873
Convert angle x from radians to degrees.");
1876
math_radians(PyObject *self, PyObject *arg)
1878
double x = PyFloat_AsDouble(arg);
1879
if (x == -1.0 && PyErr_Occurred())
1881
return PyFloat_FromDouble(x * degToRad);
1884
PyDoc_STRVAR(math_radians_doc,
1886
Convert angle x from degrees to radians.");
1889
math_isfinite(PyObject *self, PyObject *arg)
1891
double x = PyFloat_AsDouble(arg);
1892
if (x == -1.0 && PyErr_Occurred())
1894
return PyBool_FromLong((long)Py_IS_FINITE(x));
1897
PyDoc_STRVAR(math_isfinite_doc,
1898
"isfinite(x) -> bool\n\n\
1899
Return True if x is neither an infinity nor a NaN, and False otherwise.");
1902
math_isnan(PyObject *self, PyObject *arg)
1904
double x = PyFloat_AsDouble(arg);
1905
if (x == -1.0 && PyErr_Occurred())
1907
return PyBool_FromLong((long)Py_IS_NAN(x));
1910
PyDoc_STRVAR(math_isnan_doc,
1911
"isnan(x) -> bool\n\n\
1912
Return True if x is a NaN (not a number), and False otherwise.");
1915
math_isinf(PyObject *self, PyObject *arg)
1917
double x = PyFloat_AsDouble(arg);
1918
if (x == -1.0 && PyErr_Occurred())
1920
return PyBool_FromLong((long)Py_IS_INFINITY(x));
1923
PyDoc_STRVAR(math_isinf_doc,
1924
"isinf(x) -> bool\n\n\
1925
Return True if x is a positive or negative infinity, and False otherwise.");
1927
static PyMethodDef math_methods[] = {
1928
{"acos", math_acos, METH_O, math_acos_doc},
1929
{"acosh", math_acosh, METH_O, math_acosh_doc},
1930
{"asin", math_asin, METH_O, math_asin_doc},
1931
{"asinh", math_asinh, METH_O, math_asinh_doc},
1932
{"atan", math_atan, METH_O, math_atan_doc},
1933
{"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
1934
{"atanh", math_atanh, METH_O, math_atanh_doc},
1935
{"ceil", math_ceil, METH_O, math_ceil_doc},
1936
{"copysign", math_copysign, METH_VARARGS, math_copysign_doc},
1937
{"cos", math_cos, METH_O, math_cos_doc},
1938
{"cosh", math_cosh, METH_O, math_cosh_doc},
1939
{"degrees", math_degrees, METH_O, math_degrees_doc},
1940
{"erf", math_erf, METH_O, math_erf_doc},
1941
{"erfc", math_erfc, METH_O, math_erfc_doc},
1942
{"exp", math_exp, METH_O, math_exp_doc},
1943
{"expm1", math_expm1, METH_O, math_expm1_doc},
1944
{"fabs", math_fabs, METH_O, math_fabs_doc},
1945
{"factorial", math_factorial, METH_O, math_factorial_doc},
1946
{"floor", math_floor, METH_O, math_floor_doc},
1947
{"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
1948
{"frexp", math_frexp, METH_O, math_frexp_doc},
1949
{"fsum", math_fsum, METH_O, math_fsum_doc},
1950
{"gamma", math_gamma, METH_O, math_gamma_doc},
1951
{"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
1952
{"isfinite", math_isfinite, METH_O, math_isfinite_doc},
1953
{"isinf", math_isinf, METH_O, math_isinf_doc},
1954
{"isnan", math_isnan, METH_O, math_isnan_doc},
1955
{"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
1956
{"lgamma", math_lgamma, METH_O, math_lgamma_doc},
1957
{"log", math_log, METH_VARARGS, math_log_doc},
1958
{"log1p", math_log1p, METH_O, math_log1p_doc},
1959
{"log10", math_log10, METH_O, math_log10_doc},
1960
{"log2", math_log2, METH_O, math_log2_doc},
1961
{"modf", math_modf, METH_O, math_modf_doc},
1962
{"pow", math_pow, METH_VARARGS, math_pow_doc},
1963
{"radians", math_radians, METH_O, math_radians_doc},
1964
{"sin", math_sin, METH_O, math_sin_doc},
1965
{"sinh", math_sinh, METH_O, math_sinh_doc},
1966
{"sqrt", math_sqrt, METH_O, math_sqrt_doc},
1967
{"tan", math_tan, METH_O, math_tan_doc},
1968
{"tanh", math_tanh, METH_O, math_tanh_doc},
1969
{"trunc", math_trunc, METH_O, math_trunc_doc},
1970
{NULL, NULL} /* sentinel */
1974
PyDoc_STRVAR(module_doc,
1975
"This module is always available. It provides access to the\n"
1976
"mathematical functions defined by the C standard.");
1979
static struct PyModuleDef mathmodule = {
1980
PyModuleDef_HEAD_INIT,
1996
m = PyModule_Create(&mathmodule);
2000
PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
2001
PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
added
removed
Modules/mathmodule.c
