1/* Math module -- standard C math library functions, pi and e */
2
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
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
55#include "Python.h"
56#include "pycore_bitutils.h" // _Py_bit_length()
57#include "pycore_dtoa.h"
58#include "pycore_long.h" // _PyLong_GetZero()
59#include "_math.h"
60
61#include "clinic/mathmodule.c.h"
62
63/*[clinic input]
64module math
65[clinic start generated code]*/
66/*[clinic end generated code: output=da39a3ee5e6b4b0d input=76bc7002685dd942]*/
67
68
69/*
70 sin(pi*x), giving accurate results for all finite x (especially x
71 integral or close to an integer). This is here for use in the
72 reflection formula for the gamma function. It conforms to IEEE
73 754-2008 for finite arguments, but not for infinities or nans.
74*/
75
76static const double pi = 3.141592653589793238462643383279502884197;
77static const double logpi = 1.144729885849400174143427351353058711647;
78#if !defined(HAVE_ERF) || !defined(HAVE_ERFC)
79static const double sqrtpi = 1.772453850905516027298167483341145182798;
80#endif /* !defined(HAVE_ERF) || !defined(HAVE_ERFC) */
81
82
83/* Version of PyFloat_AsDouble() with in-line fast paths
84 for exact floats and integers. Gives a substantial
85 speed improvement for extracting float arguments.
86*/
87
88#define ASSIGN_DOUBLE(target_var, obj, error_label) \
89 if (PyFloat_CheckExact(obj)) { \
90 target_var = PyFloat_AS_DOUBLE(obj); \
91 } \
92 else if (PyLong_CheckExact(obj)) { \
93 target_var = PyLong_AsDouble(obj); \
94 if (target_var == -1.0 && PyErr_Occurred()) { \
95 goto error_label; \
96 } \
97 } \
98 else { \
99 target_var = PyFloat_AsDouble(obj); \
100 if (target_var == -1.0 && PyErr_Occurred()) { \
101 goto error_label; \
102 } \
103 }
104
105static double
106m_sinpi(double x)
107{
108 double y, r;
109 int n;
110 /* this function should only ever be called for finite arguments */
111 assert(Py_IS_FINITE(x));
112 y = fmod(fabs(x), 2.0);
113 n = (int)round(2.0*y);
114 assert(0 <= n && n <= 4);
115 switch (n) {
116 case 0:
117 r = sin(pi*y);
118 break;
119 case 1:
120 r = cos(pi*(y-0.5));
121 break;
122 case 2:
123 /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give
124 -0.0 instead of 0.0 when y == 1.0. */
125 r = sin(pi*(1.0-y));
126 break;
127 case 3:
128 r = -cos(pi*(y-1.5));
129 break;
130 case 4:
131 r = sin(pi*(y-2.0));
132 break;
133 default:
134 Py_UNREACHABLE();
135 }
136 return copysign(1.0, x)*r;
137}
138
139/* Implementation of the real gamma function. In extensive but non-exhaustive
140 random tests, this function proved accurate to within <= 10 ulps across the
141 entire float domain. Note that accuracy may depend on the quality of the
142 system math functions, the pow function in particular. Special cases
143 follow C99 annex F. The parameters and method are tailored to platforms
144 whose double format is the IEEE 754 binary64 format.
145
146 Method: for x > 0.0 we use the Lanczos approximation with parameters N=13
147 and g=6.024680040776729583740234375; these parameters are amongst those
148 used by the Boost library. Following Boost (again), we re-express the
149 Lanczos sum as a rational function, and compute it that way. The
150 coefficients below were computed independently using MPFR, and have been
151 double-checked against the coefficients in the Boost source code.
152
153 For x < 0.0 we use the reflection formula.
154
155 There's one minor tweak that deserves explanation: Lanczos' formula for
156 Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x
157 values, x+g-0.5 can be represented exactly. However, in cases where it
158 can't be represented exactly the small error in x+g-0.5 can be magnified
159 significantly by the pow and exp calls, especially for large x. A cheap
160 correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error
161 involved in the computation of x+g-0.5 (that is, e = computed value of
162 x+g-0.5 - exact value of x+g-0.5). Here's the proof:
163
164 Correction factor
165 -----------------
166 Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754
167 double, and e is tiny. Then:
168
169 pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e)
170 = pow(y, x-0.5)/exp(y) * C,
171
172 where the correction_factor C is given by
173
174 C = pow(1-e/y, x-0.5) * exp(e)
175
176 Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so:
177
178 C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y
179
180 But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and
181
182 pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y),
183
184 Note that for accuracy, when computing r*C it's better to do
185
186 r + e*g/y*r;
187
188 than
189
190 r * (1 + e*g/y);
191
192 since the addition in the latter throws away most of the bits of
193 information in e*g/y.
194*/
195
196#define LANCZOS_N 13
197static const double lanczos_g = 6.024680040776729583740234375;
198static const double lanczos_g_minus_half = 5.524680040776729583740234375;
199static const double lanczos_num_coeffs[LANCZOS_N] = {
200 23531376880.410759688572007674451636754734846804940,
201 42919803642.649098768957899047001988850926355848959,
202 35711959237.355668049440185451547166705960488635843,
203 17921034426.037209699919755754458931112671403265390,
204 6039542586.3520280050642916443072979210699388420708,
205 1439720407.3117216736632230727949123939715485786772,
206 248874557.86205415651146038641322942321632125127801,
207 31426415.585400194380614231628318205362874684987640,
208 2876370.6289353724412254090516208496135991145378768,
209 186056.26539522349504029498971604569928220784236328,
210 8071.6720023658162106380029022722506138218516325024,
211 210.82427775157934587250973392071336271166969580291,
212 2.5066282746310002701649081771338373386264310793408
213};
214
215/* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */
216static const double lanczos_den_coeffs[LANCZOS_N] = {
217 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
218 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};
219
220/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */
221#define NGAMMA_INTEGRAL 23
222static const double gamma_integral[NGAMMA_INTEGRAL] = {
223 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
224 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
225 1307674368000.0, 20922789888000.0, 355687428096000.0,
226 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
227 51090942171709440000.0, 1124000727777607680000.0,
228};
229
230/* Lanczos' sum L_g(x), for positive x */
231
232static double
233lanczos_sum(double x)
234{
235 double num = 0.0, den = 0.0;
236 int i;
237 assert(x > 0.0);
238 /* evaluate the rational function lanczos_sum(x). For large
239 x, the obvious algorithm risks overflow, so we instead
240 rescale the denominator and numerator of the rational
241 function by x**(1-LANCZOS_N) and treat this as a
242 rational function in 1/x. This also reduces the error for
243 larger x values. The choice of cutoff point (5.0 below) is
244 somewhat arbitrary; in tests, smaller cutoff values than
245 this resulted in lower accuracy. */
246 if (x < 5.0) {
247 for (i = LANCZOS_N; --i >= 0; ) {
248 num = num * x + lanczos_num_coeffs[i];
249 den = den * x + lanczos_den_coeffs[i];
250 }
251 }
252 else {
253 for (i = 0; i < LANCZOS_N; i++) {
254 num = num / x + lanczos_num_coeffs[i];
255 den = den / x + lanczos_den_coeffs[i];
256 }
257 }
258 return num/den;
259}
260
261/* Constant for +infinity, generated in the same way as float('inf'). */
262
263static double
264m_inf(void)
265{
266#ifndef PY_NO_SHORT_FLOAT_REPR
267 return _Py_dg_infinity(0);
268#else
269 return Py_HUGE_VAL;
270#endif
271}
272
273/* Constant nan value, generated in the same way as float('nan'). */
274/* We don't currently assume that Py_NAN is defined everywhere. */
275
276#if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN)
277
278static double
279m_nan(void)
280{
281#ifndef PY_NO_SHORT_FLOAT_REPR
282 return _Py_dg_stdnan(0);
283#else
284 return Py_NAN;
285#endif
286}
287
288#endif
289
290static double
291m_tgamma(double x)
292{
293 double absx, r, y, z, sqrtpow;
294
295 /* special cases */
296 if (!Py_IS_FINITE(x)) {
297 if (Py_IS_NAN(x) || x > 0.0)
298 return x; /* tgamma(nan) = nan, tgamma(inf) = inf */
299 else {
300 errno = EDOM;
301 return Py_NAN; /* tgamma(-inf) = nan, invalid */
302 }
303 }
304 if (x == 0.0) {
305 errno = EDOM;
306 /* tgamma(+-0.0) = +-inf, divide-by-zero */
307 return copysign(Py_HUGE_VAL, x);
308 }
309
310 /* integer arguments */
311 if (x == floor(x)) {
312 if (x < 0.0) {
313 errno = EDOM; /* tgamma(n) = nan, invalid for */
314 return Py_NAN; /* negative integers n */
315 }
316 if (x <= NGAMMA_INTEGRAL)
317 return gamma_integral[(int)x - 1];
318 }
319 absx = fabs(x);
320
321 /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */
322 if (absx < 1e-20) {
323 r = 1.0/x;
324 if (Py_IS_INFINITY(r))
325 errno = ERANGE;
326 return r;
327 }
328
329 /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for
330 x > 200, and underflows to +-0.0 for x < -200, not a negative
331 integer. */
332 if (absx > 200.0) {
333 if (x < 0.0) {
334 return 0.0/m_sinpi(x);
335 }
336 else {
337 errno = ERANGE;
338 return Py_HUGE_VAL;
339 }
340 }
341
342 y = absx + lanczos_g_minus_half;
343 /* compute error in sum */
344 if (absx > lanczos_g_minus_half) {
345 /* note: the correction can be foiled by an optimizing
346 compiler that (incorrectly) thinks that an expression like
347 a + b - a - b can be optimized to 0.0. This shouldn't
348 happen in a standards-conforming compiler. */
349 double q = y - absx;
350 z = q - lanczos_g_minus_half;
351 }
352 else {
353 double q = y - lanczos_g_minus_half;
354 z = q - absx;
355 }
356 z = z * lanczos_g / y;
357 if (x < 0.0) {
358 r = -pi / m_sinpi(absx) / absx * exp(y) / lanczos_sum(absx);
359 r -= z * r;
360 if (absx < 140.0) {
361 r /= pow(y, absx - 0.5);
362 }
363 else {
364 sqrtpow = pow(y, absx / 2.0 - 0.25);
365 r /= sqrtpow;
366 r /= sqrtpow;
367 }
368 }
369 else {
370 r = lanczos_sum(absx) / exp(y);
371 r += z * r;
372 if (absx < 140.0) {
373 r *= pow(y, absx - 0.5);
374 }
375 else {
376 sqrtpow = pow(y, absx / 2.0 - 0.25);
377 r *= sqrtpow;
378 r *= sqrtpow;
379 }
380 }
381 if (Py_IS_INFINITY(r))
382 errno = ERANGE;
383 return r;
384}
385
386/*
387 lgamma: natural log of the absolute value of the Gamma function.
388 For large arguments, Lanczos' formula works extremely well here.
389*/
390
391static double
392m_lgamma(double x)
393{
394 double r;
395 double absx;
396
397 /* special cases */
398 if (!Py_IS_FINITE(x)) {
399 if (Py_IS_NAN(x))
400 return x; /* lgamma(nan) = nan */
401 else
402 return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */
403 }
404
405 /* integer arguments */
406 if (x == floor(x) && x <= 2.0) {
407 if (x <= 0.0) {
408 errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */
409 return Py_HUGE_VAL; /* integers n <= 0 */
410 }
411 else {
412 return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */
413 }
414 }
415
416 absx = fabs(x);
417 /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */
418 if (absx < 1e-20)
419 return -log(absx);
420
421 /* Lanczos' formula. We could save a fraction of a ulp in accuracy by
422 having a second set of numerator coefficients for lanczos_sum that
423 absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g
424 subtraction below; it's probably not worth it. */
425 r = log(lanczos_sum(absx)) - lanczos_g;
426 r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1);
427 if (x < 0.0)
428 /* Use reflection formula to get value for negative x. */
429 r = logpi - log(fabs(m_sinpi(absx))) - log(absx) - r;
430 if (Py_IS_INFINITY(r))
431 errno = ERANGE;
432 return r;
433}
434
435#if !defined(HAVE_ERF) || !defined(HAVE_ERFC)
436
437/*
438 Implementations of the error function erf(x) and the complementary error
439 function erfc(x).
440
441 Method: we use a series approximation for erf for small x, and a continued
442 fraction approximation for erfc(x) for larger x;
443 combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x),
444 this gives us erf(x) and erfc(x) for all x.
445
446 The series expansion used is:
447
448 erf(x) = x*exp(-x*x)/sqrt(pi) * [
449 2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...]
450
451 The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k).
452 This series converges well for smallish x, but slowly for larger x.
453
454 The continued fraction expansion used is:
455
456 erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - )
457 3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...]
458
459 after the first term, the general term has the form:
460
461 k*(k-0.5)/(2*k+0.5 + x**2 - ...).
462
463 This expansion converges fast for larger x, but convergence becomes
464 infinitely slow as x approaches 0.0. The (somewhat naive) continued
465 fraction evaluation algorithm used below also risks overflow for large x;
466 but for large x, erfc(x) == 0.0 to within machine precision. (For
467 example, erfc(30.0) is approximately 2.56e-393).
468
469 Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and
470 continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) <
471 ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the
472 numbers of terms to use for the relevant expansions. */
473
474#define ERF_SERIES_CUTOFF 1.5
475#define ERF_SERIES_TERMS 25
476#define ERFC_CONTFRAC_CUTOFF 30.0
477#define ERFC_CONTFRAC_TERMS 50
478
479/*
480 Error function, via power series.
481
482 Given a finite float x, return an approximation to erf(x).
483 Converges reasonably fast for small x.
484*/
485
486static double
487m_erf_series(double x)
488{
489 double x2, acc, fk, result;
490 int i, saved_errno;
491
492 x2 = x * x;
493 acc = 0.0;
494 fk = (double)ERF_SERIES_TERMS + 0.5;
495 for (i = 0; i < ERF_SERIES_TERMS; i++) {
496 acc = 2.0 + x2 * acc / fk;
497 fk -= 1.0;
498 }
499 /* Make sure the exp call doesn't affect errno;
500 see m_erfc_contfrac for more. */
501 saved_errno = errno;
502 result = acc * x * exp(-x2) / sqrtpi;
503 errno = saved_errno;
504 return result;
505}
506
507/*
508 Complementary error function, via continued fraction expansion.
509
510 Given a positive float x, return an approximation to erfc(x). Converges
511 reasonably fast for x large (say, x > 2.0), and should be safe from
512 overflow if x and nterms are not too large. On an IEEE 754 machine, with x
513 <= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller
514 than the smallest representable nonzero float. */
515
516static double
517m_erfc_contfrac(double x)
518{
519 double x2, a, da, p, p_last, q, q_last, b, result;
520 int i, saved_errno;
521
522 if (x >= ERFC_CONTFRAC_CUTOFF)
523 return 0.0;
524
525 x2 = x*x;
526 a = 0.0;
527 da = 0.5;
528 p = 1.0; p_last = 0.0;
529 q = da + x2; q_last = 1.0;
530 for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) {
531 double temp;
532 a += da;
533 da += 2.0;
534 b = da + x2;
535 temp = p; p = b*p - a*p_last; p_last = temp;
536 temp = q; q = b*q - a*q_last; q_last = temp;
537 }
538 /* Issue #8986: On some platforms, exp sets errno on underflow to zero;
539 save the current errno value so that we can restore it later. */
540 saved_errno = errno;
541 result = p / q * x * exp(-x2) / sqrtpi;
542 errno = saved_errno;
543 return result;
544}
545
546#endif /* !defined(HAVE_ERF) || !defined(HAVE_ERFC) */
547
548/* Error function erf(x), for general x */
549
550static double
551m_erf(double x)
552{
553#ifdef HAVE_ERF
554 return erf(x);
555#else
556 double absx, cf;
557
558 if (Py_IS_NAN(x))
559 return x;
560 absx = fabs(x);
561 if (absx < ERF_SERIES_CUTOFF)
562 return m_erf_series(x);
563 else {
564 cf = m_erfc_contfrac(absx);
565 return x > 0.0 ? 1.0 - cf : cf - 1.0;
566 }
567#endif
568}
569
570/* Complementary error function erfc(x), for general x. */
571
572static double
573m_erfc(double x)
574{
575#ifdef HAVE_ERFC
576 return erfc(x);
577#else
578 double absx, cf;
579
580 if (Py_IS_NAN(x))
581 return x;
582 absx = fabs(x);
583 if (absx < ERF_SERIES_CUTOFF)
584 return 1.0 - m_erf_series(x);
585 else {
586 cf = m_erfc_contfrac(absx);
587 return x > 0.0 ? cf : 2.0 - cf;
588 }
589#endif
590}
591
592/*
593 wrapper for atan2 that deals directly with special cases before
594 delegating to the platform libm for the remaining cases. This
595 is necessary to get consistent behaviour across platforms.
596 Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
597 always follow C99.
598*/
599
600static double
601m_atan2(double y, double x)
602{
603 if (Py_IS_NAN(x) || Py_IS_NAN(y))
604 return Py_NAN;
605 if (Py_IS_INFINITY(y)) {
606 if (Py_IS_INFINITY(x)) {
607 if (copysign(1., x) == 1.)
608 /* atan2(+-inf, +inf) == +-pi/4 */
609 return copysign(0.25*Py_MATH_PI, y);
610 else
611 /* atan2(+-inf, -inf) == +-pi*3/4 */
612 return copysign(0.75*Py_MATH_PI, y);
613 }
614 /* atan2(+-inf, x) == +-pi/2 for finite x */
615 return copysign(0.5*Py_MATH_PI, y);
616 }
617 if (Py_IS_INFINITY(x) || y == 0.) {
618 if (copysign(1., x) == 1.)
619 /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
620 return copysign(0., y);
621 else
622 /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
623 return copysign(Py_MATH_PI, y);
624 }
625 return atan2(y, x);
626}
627
628
629/* IEEE 754-style remainder operation: x - n*y where n*y is the nearest
630 multiple of y to x, taking n even in the case of a tie. Assuming an IEEE 754
631 binary floating-point format, the result is always exact. */
632
633static double
634m_remainder(double x, double y)
635{
636 /* Deal with most common case first. */
637 if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) {
638 double absx, absy, c, m, r;
639
640 if (y == 0.0) {
641 return Py_NAN;
642 }
643
644 absx = fabs(x);
645 absy = fabs(y);
646 m = fmod(absx, absy);
647
648 /*
649 Warning: some subtlety here. What we *want* to know at this point is
650 whether the remainder m is less than, equal to, or greater than half
651 of absy. However, we can't do that comparison directly because we
652 can't be sure that 0.5*absy is representable (the multiplication
653 might incur precision loss due to underflow). So instead we compare
654 m with the complement c = absy - m: m < 0.5*absy if and only if m <
655 c, and so on. The catch is that absy - m might also not be
656 representable, but it turns out that it doesn't matter:
657
658 - if m > 0.5*absy then absy - m is exactly representable, by
659 Sterbenz's lemma, so m > c
660 - if m == 0.5*absy then again absy - m is exactly representable
661 and m == c
662 - if m < 0.5*absy then either (i) 0.5*absy is exactly representable,
663 in which case 0.5*absy < absy - m, so 0.5*absy <= c and hence m <
664 c, or (ii) absy is tiny, either subnormal or in the lowest normal
665 binade. Then absy - m is exactly representable and again m < c.
666 */
667
668 c = absy - m;
669 if (m < c) {
670 r = m;
671 }
672 else if (m > c) {
673 r = -c;
674 }
675 else {
676 /*
677 Here absx is exactly halfway between two multiples of absy,
678 and we need to choose the even multiple. x now has the form
679
680 absx = n * absy + m
681
682 for some integer n (recalling that m = 0.5*absy at this point).
683 If n is even we want to return m; if n is odd, we need to
684 return -m.
685
686 So
687
688 0.5 * (absx - m) = (n/2) * absy
689
690 and now reducing modulo absy gives us:
691
692 | m, if n is odd
693 fmod(0.5 * (absx - m), absy) = |
694 | 0, if n is even
695
696 Now m - 2.0 * fmod(...) gives the desired result: m
697 if n is even, -m if m is odd.
698
699 Note that all steps in fmod(0.5 * (absx - m), absy)
700 will be computed exactly, with no rounding error
701 introduced.
702 */
703 assert(m == c);
704 r = m - 2.0 * fmod(0.5 * (absx - m), absy);
705 }
706 return copysign(1.0, x) * r;
707 }
708
709 /* Special values. */
710 if (Py_IS_NAN(x)) {
711 return x;
712 }
713 if (Py_IS_NAN(y)) {
714 return y;
715 }
716 if (Py_IS_INFINITY(x)) {
717 return Py_NAN;
718 }
719 assert(Py_IS_INFINITY(y));
720 return x;
721}
722
723
724/*
725 Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
726 log(-ve), log(NaN). Here are wrappers for log and log10 that deal with
727 special values directly, passing positive non-special values through to
728 the system log/log10.
729 */
730
731static double
732m_log(double x)
733{
734 if (Py_IS_FINITE(x)) {
735 if (x > 0.0)
736 return log(x);
737 errno = EDOM;
738 if (x == 0.0)
739 return -Py_HUGE_VAL; /* log(0) = -inf */
740 else
741 return Py_NAN; /* log(-ve) = nan */
742 }
743 else if (Py_IS_NAN(x))
744 return x; /* log(nan) = nan */
745 else if (x > 0.0)
746 return x; /* log(inf) = inf */
747 else {
748 errno = EDOM;
749 return Py_NAN; /* log(-inf) = nan */
750 }
751}
752
753/*
754 log2: log to base 2.
755
756 Uses an algorithm that should:
757
758 (a) produce exact results for powers of 2, and
759 (b) give a monotonic log2 (for positive finite floats),
760 assuming that the system log is monotonic.
761*/
762
763static double
764m_log2(double x)
765{
766 if (!Py_IS_FINITE(x)) {
767 if (Py_IS_NAN(x))
768 return x; /* log2(nan) = nan */
769 else if (x > 0.0)
770 return x; /* log2(+inf) = +inf */
771 else {
772 errno = EDOM;
773 return Py_NAN; /* log2(-inf) = nan, invalid-operation */
774 }
775 }
776
777 if (x > 0.0) {
778#ifdef HAVE_LOG2
779 return log2(x);
780#else
781 double m;
782 int e;
783 m = frexp(x, &e);
784 /* We want log2(m * 2**e) == log(m) / log(2) + e. Care is needed when
785 * x is just greater than 1.0: in that case e is 1, log(m) is negative,
786 * and we get significant cancellation error from the addition of
787 * log(m) / log(2) to e. The slight rewrite of the expression below
788 * avoids this problem.
789 */
790 if (x >= 1.0) {
791 return log(2.0 * m) / log(2.0) + (e - 1);
792 }
793 else {
794 return log(m) / log(2.0) + e;
795 }
796#endif
797 }
798 else if (x == 0.0) {
799 errno = EDOM;
800 return -Py_HUGE_VAL; /* log2(0) = -inf, divide-by-zero */
801 }
802 else {
803 errno = EDOM;
804 return Py_NAN; /* log2(-inf) = nan, invalid-operation */
805 }
806}
807
808static double
809m_log10(double x)
810{
811 if (Py_IS_FINITE(x)) {
812 if (x > 0.0)
813 return log10(x);
814 errno = EDOM;
815 if (x == 0.0)
816 return -Py_HUGE_VAL; /* log10(0) = -inf */
817 else
818 return Py_NAN; /* log10(-ve) = nan */
819 }
820 else if (Py_IS_NAN(x))
821 return x; /* log10(nan) = nan */
822 else if (x > 0.0)
823 return x; /* log10(inf) = inf */
824 else {
825 errno = EDOM;
826 return Py_NAN; /* log10(-inf) = nan */
827 }
828}
829
830
831static PyObject *
832math_gcd(PyObject *module, PyObject * const *args, Py_ssize_t nargs)
833{
834 PyObject *res, *x;
835 Py_ssize_t i;
836
837 if (nargs == 0) {
838 return PyLong_FromLong(0);
839 }
840 res = PyNumber_Index(args[0]);
841 if (res == NULL) {
842 return NULL;
843 }
844 if (nargs == 1) {
845 Py_SETREF(res, PyNumber_Absolute(res));
846 return res;
847 }
848
849 PyObject *one = _PyLong_GetOne(); // borrowed ref
850 for (i = 1; i < nargs; i++) {
851 x = _PyNumber_Index(args[i]);
852 if (x == NULL) {
853 Py_DECREF(res);
854 return NULL;
855 }
856 if (res == one) {
857 /* Fast path: just check arguments.
858 It is okay to use identity comparison here. */
859 Py_DECREF(x);
860 continue;
861 }
862 Py_SETREF(res, _PyLong_GCD(res, x));
863 Py_DECREF(x);
864 if (res == NULL) {
865 return NULL;
866 }
867 }
868 return res;
869}
870
871PyDoc_STRVAR(math_gcd_doc,
872"gcd($module, *integers)\n"
873"--\n"
874"\n"
875"Greatest Common Divisor.");
876
877
878static PyObject *
879long_lcm(PyObject *a, PyObject *b)
880{
881 PyObject *g, *m, *f, *ab;
882
883 if (Py_SIZE(a) == 0 || Py_SIZE(b) == 0) {
884 return PyLong_FromLong(0);
885 }
886 g = _PyLong_GCD(a, b);
887 if (g == NULL) {
888 return NULL;
889 }
890 f = PyNumber_FloorDivide(a, g);
891 Py_DECREF(g);
892 if (f == NULL) {
893 return NULL;
894 }
895 m = PyNumber_Multiply(f, b);
896 Py_DECREF(f);
897 if (m == NULL) {
898 return NULL;
899 }
900 ab = PyNumber_Absolute(m);
901 Py_DECREF(m);
902 return ab;
903}
904
905
906static PyObject *
907math_lcm(PyObject *module, PyObject * const *args, Py_ssize_t nargs)
908{
909 PyObject *res, *x;
910 Py_ssize_t i;
911
912 if (nargs == 0) {
913 return PyLong_FromLong(1);
914 }
915 res = PyNumber_Index(args[0]);
916 if (res == NULL) {
917 return NULL;
918 }
919 if (nargs == 1) {
920 Py_SETREF(res, PyNumber_Absolute(res));
921 return res;
922 }
923
924 PyObject *zero = _PyLong_GetZero(); // borrowed ref
925 for (i = 1; i < nargs; i++) {
926 x = PyNumber_Index(args[i]);
927 if (x == NULL) {
928 Py_DECREF(res);
929 return NULL;
930 }
931 if (res == zero) {
932 /* Fast path: just check arguments.
933 It is okay to use identity comparison here. */
934 Py_DECREF(x);
935 continue;
936 }
937 Py_SETREF(res, long_lcm(res, x));
938 Py_DECREF(x);
939 if (res == NULL) {
940 return NULL;
941 }
942 }
943 return res;
944}
945
946
947PyDoc_STRVAR(math_lcm_doc,
948"lcm($module, *integers)\n"
949"--\n"
950"\n"
951"Least Common Multiple.");
952
953
954/* Call is_error when errno != 0, and where x is the result libm
955 * returned. is_error will usually set up an exception and return
956 * true (1), but may return false (0) without setting up an exception.
957 */
958static int
959is_error(double x)
960{
961 int result = 1; /* presumption of guilt */
962 assert(errno); /* non-zero errno is a precondition for calling */
963 if (errno == EDOM)
964 PyErr_SetString(PyExc_ValueError, "math domain error");
965
966 else if (errno == ERANGE) {
967 /* ANSI C generally requires libm functions to set ERANGE
968 * on overflow, but also generally *allows* them to set
969 * ERANGE on underflow too. There's no consistency about
970 * the latter across platforms.
971 * Alas, C99 never requires that errno be set.
972 * Here we suppress the underflow errors (libm functions
973 * should return a zero on underflow, and +- HUGE_VAL on
974 * overflow, so testing the result for zero suffices to
975 * distinguish the cases).
976 *
977 * On some platforms (Ubuntu/ia64) it seems that errno can be
978 * set to ERANGE for subnormal results that do *not* underflow
979 * to zero. So to be safe, we'll ignore ERANGE whenever the
980 * function result is less than 1.5 in absolute value.
981 *
982 * bpo-46018: Changed to 1.5 to ensure underflows in expm1()
983 * are correctly detected, since the function may underflow
984 * toward -1.0 rather than 0.0.
985 */
986 if (fabs(x) < 1.5)
987 result = 0;
988 else
989 PyErr_SetString(PyExc_OverflowError,
990 "math range error");
991 }
992 else
993 /* Unexpected math error */
994 PyErr_SetFromErrno(PyExc_ValueError);
995 return result;
996}
997
998/*
999 math_1 is used to wrap a libm function f that takes a double
1000 argument and returns a double.
1001
1002 The error reporting follows these rules, which are designed to do
1003 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
1004 platforms.
1005
1006 - a NaN result from non-NaN inputs causes ValueError to be raised
1007 - an infinite result from finite inputs causes OverflowError to be
1008 raised if can_overflow is 1, or raises ValueError if can_overflow
1009 is 0.
1010 - if the result is finite and errno == EDOM then ValueError is
1011 raised
1012 - if the result is finite and nonzero and errno == ERANGE then
1013 OverflowError is raised
1014
1015 The last rule is used to catch overflow on platforms which follow
1016 C89 but for which HUGE_VAL is not an infinity.
1017
1018 For the majority of one-argument functions these rules are enough
1019 to ensure that Python's functions behave as specified in 'Annex F'
1020 of the C99 standard, with the 'invalid' and 'divide-by-zero'
1021 floating-point exceptions mapping to Python's ValueError and the
1022 'overflow' floating-point exception mapping to OverflowError.
1023 math_1 only works for functions that don't have singularities *and*
1024 the possibility of overflow; fortunately, that covers everything we
1025 care about right now.
1026*/
1027
1028static PyObject *
1029math_1_to_whatever(PyObject *arg, double (*func) (double),
1030 PyObject *(*from_double_func) (double),
1031 int can_overflow)
1032{
1033 double x, r;
1034 x = PyFloat_AsDouble(arg);
1035 if (x == -1.0 && PyErr_Occurred())
1036 return NULL;
1037 errno = 0;
1038 r = (*func)(x);
1039 if (Py_IS_NAN(r) && !Py_IS_NAN(x)) {
1040 PyErr_SetString(PyExc_ValueError,
1041 "math domain error"); /* invalid arg */
1042 return NULL;
1043 }
1044 if (Py_IS_INFINITY(r) && Py_IS_FINITE(x)) {
1045 if (can_overflow)
1046 PyErr_SetString(PyExc_OverflowError,
1047 "math range error"); /* overflow */
1048 else
1049 PyErr_SetString(PyExc_ValueError,
1050 "math domain error"); /* singularity */
1051 return NULL;
1052 }
1053 if (Py_IS_FINITE(r) && errno && is_error(r))
1054 /* this branch unnecessary on most platforms */
1055 return NULL;
1056
1057 return (*from_double_func)(r);
1058}
1059
1060/* variant of math_1, to be used when the function being wrapped is known to
1061 set errno properly (that is, errno = EDOM for invalid or divide-by-zero,
1062 errno = ERANGE for overflow). */
1063
1064static PyObject *
1065math_1a(PyObject *arg, double (*func) (double))
1066{
1067 double x, r;
1068 x = PyFloat_AsDouble(arg);
1069 if (x == -1.0 && PyErr_Occurred())
1070 return NULL;
1071 errno = 0;
1072 r = (*func)(x);
1073 if (errno && is_error(r))
1074 return NULL;
1075 return PyFloat_FromDouble(r);
1076}
1077
1078/*
1079 math_2 is used to wrap a libm function f that takes two double
1080 arguments and returns a double.
1081
1082 The error reporting follows these rules, which are designed to do
1083 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
1084 platforms.
1085
1086 - a NaN result from non-NaN inputs causes ValueError to be raised
1087 - an infinite result from finite inputs causes OverflowError to be
1088 raised.
1089 - if the result is finite and errno == EDOM then ValueError is
1090 raised
1091 - if the result is finite and nonzero and errno == ERANGE then
1092 OverflowError is raised
1093
1094 The last rule is used to catch overflow on platforms which follow
1095 C89 but for which HUGE_VAL is not an infinity.
1096
1097 For most two-argument functions (copysign, fmod, hypot, atan2)
1098 these rules are enough to ensure that Python's functions behave as
1099 specified in 'Annex F' of the C99 standard, with the 'invalid' and
1100 'divide-by-zero' floating-point exceptions mapping to Python's
1101 ValueError and the 'overflow' floating-point exception mapping to
1102 OverflowError.
1103*/
1104
1105static PyObject *
1106math_1(PyObject *arg, double (*func) (double), int can_overflow)
1107{
1108 return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow);
1109}
1110
1111static PyObject *
1112math_2(PyObject *const *args, Py_ssize_t nargs,
1113 double (*func) (double, double), const char *funcname)
1114{
1115 double x, y, r;
1116 if (!_PyArg_CheckPositional(funcname, nargs, 2, 2))
1117 return NULL;
1118 x = PyFloat_AsDouble(args[0]);
1119 if (x == -1.0 && PyErr_Occurred()) {
1120 return NULL;
1121 }
1122 y = PyFloat_AsDouble(args[1]);
1123 if (y == -1.0 && PyErr_Occurred()) {
1124 return NULL;
1125 }
1126 errno = 0;
1127 r = (*func)(x, y);
1128 if (Py_IS_NAN(r)) {
1129 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1130 errno = EDOM;
1131 else
1132 errno = 0;
1133 }
1134 else if (Py_IS_INFINITY(r)) {
1135 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
1136 errno = ERANGE;
1137 else
1138 errno = 0;
1139 }
1140 if (errno && is_error(r))
1141 return NULL;
1142 else
1143 return PyFloat_FromDouble(r);
1144}
1145
1146#define FUNC1(funcname, func, can_overflow, docstring) \
1147 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
1148 return math_1(args, func, can_overflow); \
1149 }\
1150 PyDoc_STRVAR(math_##funcname##_doc, docstring);
1151
1152#define FUNC1A(funcname, func, docstring) \
1153 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
1154 return math_1a(args, func); \
1155 }\
1156 PyDoc_STRVAR(math_##funcname##_doc, docstring);
1157
1158#define FUNC2(funcname, func, docstring) \
1159 static PyObject * math_##funcname(PyObject *self, PyObject *const *args, Py_ssize_t nargs) { \
1160 return math_2(args, nargs, func, #funcname); \
1161 }\
1162 PyDoc_STRVAR(math_##funcname##_doc, docstring);
1163
1164FUNC1(acos, acos, 0,
1165 "acos($module, x, /)\n--\n\n"
1166 "Return the arc cosine (measured in radians) of x.\n\n"
1167 "The result is between 0 and pi.")
1168FUNC1(acosh, m_acosh, 0,
1169 "acosh($module, x, /)\n--\n\n"
1170 "Return the inverse hyperbolic cosine of x.")
1171FUNC1(asin, asin, 0,
1172 "asin($module, x, /)\n--\n\n"
1173 "Return the arc sine (measured in radians) of x.\n\n"
1174 "The result is between -pi/2 and pi/2.")
1175FUNC1(asinh, m_asinh, 0,
1176 "asinh($module, x, /)\n--\n\n"
1177 "Return the inverse hyperbolic sine of x.")
1178FUNC1(atan, atan, 0,
1179 "atan($module, x, /)\n--\n\n"
1180 "Return the arc tangent (measured in radians) of x.\n\n"
1181 "The result is between -pi/2 and pi/2.")
1182FUNC2(atan2, m_atan2,
1183 "atan2($module, y, x, /)\n--\n\n"
1184 "Return the arc tangent (measured in radians) of y/x.\n\n"
1185 "Unlike atan(y/x), the signs of both x and y are considered.")
1186FUNC1(atanh, m_atanh, 0,
1187 "atanh($module, x, /)\n--\n\n"
1188 "Return the inverse hyperbolic tangent of x.")
1189
1190/*[clinic input]
1191math.ceil
1192
1193 x as number: object
1194 /
1195
1196Return the ceiling of x as an Integral.
1197
1198This is the smallest integer >= x.
1199[clinic start generated code]*/
1200
1201static PyObject *
1202math_ceil(PyObject *module, PyObject *number)
1203/*[clinic end generated code: output=6c3b8a78bc201c67 input=2725352806399cab]*/
1204{
1205 _Py_IDENTIFIER(__ceil__);
1206
1207 if (!PyFloat_CheckExact(number)) {
1208 PyObject *method = _PyObject_LookupSpecial(number, &PyId___ceil__);
1209 if (method != NULL) {
1210 PyObject *result = _PyObject_CallNoArg(method);
1211 Py_DECREF(method);
1212 return result;
1213 }
1214 if (PyErr_Occurred())
1215 return NULL;
1216 }
1217 double x = PyFloat_AsDouble(number);
1218 if (x == -1.0 && PyErr_Occurred())
1219 return NULL;
1220
1221 return PyLong_FromDouble(ceil(x));
1222}
1223
1224FUNC2(copysign, copysign,
1225 "copysign($module, x, y, /)\n--\n\n"
1226 "Return a float with the magnitude (absolute value) of x but the sign of y.\n\n"
1227 "On platforms that support signed zeros, copysign(1.0, -0.0)\n"
1228 "returns -1.0.\n")
1229FUNC1(cos, cos, 0,
1230 "cos($module, x, /)\n--\n\n"
1231 "Return the cosine of x (measured in radians).")
1232FUNC1(cosh, cosh, 1,
1233 "cosh($module, x, /)\n--\n\n"
1234 "Return the hyperbolic cosine of x.")
1235FUNC1A(erf, m_erf,
1236 "erf($module, x, /)\n--\n\n"
1237 "Error function at x.")
1238FUNC1A(erfc, m_erfc,
1239 "erfc($module, x, /)\n--\n\n"
1240 "Complementary error function at x.")
1241FUNC1(exp, exp, 1,
1242 "exp($module, x, /)\n--\n\n"
1243 "Return e raised to the power of x.")
1244FUNC1(expm1, m_expm1, 1,
1245 "expm1($module, x, /)\n--\n\n"
1246 "Return exp(x)-1.\n\n"
1247 "This function avoids the loss of precision involved in the direct "
1248 "evaluation of exp(x)-1 for small x.")
1249FUNC1(fabs, fabs, 0,
1250 "fabs($module, x, /)\n--\n\n"
1251 "Return the absolute value of the float x.")
1252
1253/*[clinic input]
1254math.floor
1255
1256 x as number: object
1257 /
1258
1259Return the floor of x as an Integral.
1260
1261This is the largest integer <= x.
1262[clinic start generated code]*/
1263
1264static PyObject *
1265math_floor(PyObject *module, PyObject *number)
1266/*[clinic end generated code: output=c6a65c4884884b8a input=63af6b5d7ebcc3d6]*/
1267{
1268 double x;
1269
1270 _Py_IDENTIFIER(__floor__);
1271
1272 if (PyFloat_CheckExact(number)) {
1273 x = PyFloat_AS_DOUBLE(number);
1274 }
1275 else
1276 {
1277 PyObject *method = _PyObject_LookupSpecial(number, &PyId___floor__);
1278 if (method != NULL) {
1279 PyObject *result = _PyObject_CallNoArg(method);
1280 Py_DECREF(method);
1281 return result;
1282 }
1283 if (PyErr_Occurred())
1284 return NULL;
1285 x = PyFloat_AsDouble(number);
1286 if (x == -1.0 && PyErr_Occurred())
1287 return NULL;
1288 }
1289 return PyLong_FromDouble(floor(x));
1290}
1291
1292FUNC1A(gamma, m_tgamma,
1293 "gamma($module, x, /)\n--\n\n"
1294 "Gamma function at x.")
1295FUNC1A(lgamma, m_lgamma,
1296 "lgamma($module, x, /)\n--\n\n"
1297 "Natural logarithm of absolute value of Gamma function at x.")
1298FUNC1(log1p, m_log1p, 0,
1299 "log1p($module, x, /)\n--\n\n"
1300 "Return the natural logarithm of 1+x (base e).\n\n"
1301 "The result is computed in a way which is accurate for x near zero.")
1302FUNC2(remainder, m_remainder,
1303 "remainder($module, x, y, /)\n--\n\n"
1304 "Difference between x and the closest integer multiple of y.\n\n"
1305 "Return x - n*y where n*y is the closest integer multiple of y.\n"
1306 "In the case where x is exactly halfway between two multiples of\n"
1307 "y, the nearest even value of n is used. The result is always exact.")
1308FUNC1(sin, sin, 0,
1309 "sin($module, x, /)\n--\n\n"
1310 "Return the sine of x (measured in radians).")
1311FUNC1(sinh, sinh, 1,
1312 "sinh($module, x, /)\n--\n\n"
1313 "Return the hyperbolic sine of x.")
1314FUNC1(sqrt, sqrt, 0,
1315 "sqrt($module, x, /)\n--\n\n"
1316 "Return the square root of x.")
1317FUNC1(tan, tan, 0,
1318 "tan($module, x, /)\n--\n\n"
1319 "Return the tangent of x (measured in radians).")
1320FUNC1(tanh, tanh, 0,
1321 "tanh($module, x, /)\n--\n\n"
1322 "Return the hyperbolic tangent of x.")
1323
1324/* Precision summation function as msum() by Raymond Hettinger in
1325 <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
1326 enhanced with the exact partials sum and roundoff from Mark
1327 Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
1328 See those links for more details, proofs and other references.
1329
1330 Note 1: IEEE 754R floating point semantics are assumed,
1331 but the current implementation does not re-establish special
1332 value semantics across iterations (i.e. handling -Inf + Inf).
1333
1334 Note 2: No provision is made for intermediate overflow handling;
1335 therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while
1336 sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
1337 overflow of the first partial sum.
1338
1339 Note 3: The intermediate values lo, yr, and hi are declared volatile so
1340 aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
1341 Also, the volatile declaration forces the values to be stored in memory as
1342 regular doubles instead of extended long precision (80-bit) values. This
1343 prevents double rounding because any addition or subtraction of two doubles
1344 can be resolved exactly into double-sized hi and lo values. As long as the
1345 hi value gets forced into a double before yr and lo are computed, the extra
1346 bits in downstream extended precision operations (x87 for example) will be
1347 exactly zero and therefore can be losslessly stored back into a double,
1348 thereby preventing double rounding.
1349
1350 Note 4: A similar implementation is in Modules/cmathmodule.c.
1351 Be sure to update both when making changes.
1352
1353 Note 5: The signature of math.fsum() differs from builtins.sum()
1354 because the start argument doesn't make sense in the context of
1355 accurate summation. Since the partials table is collapsed before
1356 returning a result, sum(seq2, start=sum(seq1)) may not equal the
1357 accurate result returned by sum(itertools.chain(seq1, seq2)).
1358*/
1359
1360#define NUM_PARTIALS 32 /* initial partials array size, on stack */
1361
1362/* Extend the partials array p[] by doubling its size. */
1363static int /* non-zero on error */
1364_fsum_realloc(double **p_ptr, Py_ssize_t n,
1365 double *ps, Py_ssize_t *m_ptr)
1366{
1367 void *v = NULL;
1368 Py_ssize_t m = *m_ptr;
1369
1370 m += m; /* double */
1371 if (n < m && (size_t)m < ((size_t)PY_SSIZE_T_MAX / sizeof(double))) {
1372 double *p = *p_ptr;
1373 if (p == ps) {
1374 v = PyMem_Malloc(sizeof(double) * m);
1375 if (v != NULL)
1376 memcpy(v, ps, sizeof(double) * n);
1377 }
1378 else
1379 v = PyMem_Realloc(p, sizeof(double) * m);
1380 }
1381 if (v == NULL) { /* size overflow or no memory */
1382 PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
1383 return 1;
1384 }
1385 *p_ptr = (double*) v;
1386 *m_ptr = m;
1387 return 0;
1388}
1389
1390/* Full precision summation of a sequence of floats.
1391
1392 def msum(iterable):
1393 partials = [] # sorted, non-overlapping partial sums
1394 for x in iterable:
1395 i = 0
1396 for y in partials:
1397 if abs(x) < abs(y):
1398 x, y = y, x
1399 hi = x + y
1400 lo = y - (hi - x)
1401 if lo:
1402 partials[i] = lo
1403 i += 1
1404 x = hi
1405 partials[i:] = [x]
1406 return sum_exact(partials)
1407
1408 Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo
1409 are exactly equal to x+y. The inner loop applies hi/lo summation to each
1410 partial so that the list of partial sums remains exact.
1411
1412 Sum_exact() adds the partial sums exactly and correctly rounds the final
1413 result (using the round-half-to-even rule). The items in partials remain
1414 non-zero, non-special, non-overlapping and strictly increasing in
1415 magnitude, but possibly not all having the same sign.
1416
1417 Depends on IEEE 754 arithmetic guarantees and half-even rounding.
1418*/
1419
1420/*[clinic input]
1421math.fsum
1422
1423 seq: object
1424 /
1425
1426Return an accurate floating point sum of values in the iterable seq.
1427
1428Assumes IEEE-754 floating point arithmetic.
1429[clinic start generated code]*/
1430
1431static PyObject *
1432math_fsum(PyObject *module, PyObject *seq)
1433/*[clinic end generated code: output=ba5c672b87fe34fc input=c51b7d8caf6f6e82]*/
1434{
1435 PyObject *item, *iter, *sum = NULL;
1436 Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
1437 double x, y, t, ps[NUM_PARTIALS], *p = ps;
1438 double xsave, special_sum = 0.0, inf_sum = 0.0;
1439 volatile double hi, yr, lo;
1440
1441 iter = PyObject_GetIter(seq);
1442 if (iter == NULL)
1443 return NULL;
1444
1445 for(;;) { /* for x in iterable */
1446 assert(0 <= n && n <= m);
1447 assert((m == NUM_PARTIALS && p == ps) ||
1448 (m > NUM_PARTIALS && p != NULL));
1449
1450 item = PyIter_Next(iter);
1451 if (item == NULL) {
1452 if (PyErr_Occurred())
1453 goto _fsum_error;
1454 break;
1455 }
1456 ASSIGN_DOUBLE(x, item, error_with_item);
1457 Py_DECREF(item);
1458
1459 xsave = x;
1460 for (i = j = 0; j < n; j++) { /* for y in partials */
1461 y = p[j];
1462 if (fabs(x) < fabs(y)) {
1463 t = x; x = y; y = t;
1464 }
1465 hi = x + y;
1466 yr = hi - x;
1467 lo = y - yr;
1468 if (lo != 0.0)
1469 p[i++] = lo;
1470 x = hi;
1471 }
1472
1473 n = i; /* ps[i:] = [x] */
1474 if (x != 0.0) {
1475 if (! Py_IS_FINITE(x)) {
1476 /* a nonfinite x could arise either as
1477 a result of intermediate overflow, or
1478 as a result of a nan or inf in the
1479 summands */
1480 if (Py_IS_FINITE(xsave)) {
1481 PyErr_SetString(PyExc_OverflowError,
1482 "intermediate overflow in fsum");
1483 goto _fsum_error;
1484 }
1485 if (Py_IS_INFINITY(xsave))
1486 inf_sum += xsave;
1487 special_sum += xsave;
1488 /* reset partials */
1489 n = 0;
1490 }
1491 else if (n >= m && _fsum_realloc(&p, n, ps, &m))
1492 goto _fsum_error;
1493 else
1494 p[n++] = x;
1495 }
1496 }
1497
1498 if (special_sum != 0.0) {
1499 if (Py_IS_NAN(inf_sum))
1500 PyErr_SetString(PyExc_ValueError,
1501 "-inf + inf in fsum");
1502 else
1503 sum = PyFloat_FromDouble(special_sum);
1504 goto _fsum_error;
1505 }
1506
1507 hi = 0.0;
1508 if (n > 0) {
1509 hi = p[--n];
1510 /* sum_exact(ps, hi) from the top, stop when the sum becomes
1511 inexact. */
1512 while (n > 0) {
1513 x = hi;
1514 y = p[--n];
1515 assert(fabs(y) < fabs(x));
1516 hi = x + y;
1517 yr = hi - x;
1518 lo = y - yr;
1519 if (lo != 0.0)
1520 break;
1521 }
1522 /* Make half-even rounding work across multiple partials.
1523 Needed so that sum([1e-16, 1, 1e16]) will round-up the last
1524 digit to two instead of down to zero (the 1e-16 makes the 1
1525 slightly closer to two). With a potential 1 ULP rounding
1526 error fixed-up, math.fsum() can guarantee commutativity. */
1527 if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
1528 (lo > 0.0 && p[n-1] > 0.0))) {
1529 y = lo * 2.0;
1530 x = hi + y;
1531 yr = x - hi;
1532 if (y == yr)
1533 hi = x;
1534 }
1535 }
1536 sum = PyFloat_FromDouble(hi);
1537
1538 _fsum_error:
1539 Py_DECREF(iter);
1540 if (p != ps)
1541 PyMem_Free(p);
1542 return sum;
1543
1544 error_with_item:
1545 Py_DECREF(item);
1546 goto _fsum_error;
1547}
1548
1549#undef NUM_PARTIALS
1550
1551
1552static unsigned long
1553count_set_bits(unsigned long n)
1554{
1555 unsigned long count = 0;
1556 while (n != 0) {
1557 ++count;
1558 n &= n - 1; /* clear least significant bit */
1559 }
1560 return count;
1561}
1562
1563/* Integer square root
1564
1565Given a nonnegative integer `n`, we want to compute the largest integer
1566`a` for which `a * a <= n`, or equivalently the integer part of the exact
1567square root of `n`.
1568
1569We use an adaptive-precision pure-integer version of Newton's iteration. Given
1570a positive integer `n`, the algorithm produces at each iteration an integer
1571approximation `a` to the square root of `n >> s` for some even integer `s`,
1572with `s` decreasing as the iterations progress. On the final iteration, `s` is
1573zero and we have an approximation to the square root of `n` itself.
1574
1575At every step, the approximation `a` is strictly within 1.0 of the true square
1576root, so we have
1577
1578 (a - 1)**2 < (n >> s) < (a + 1)**2
1579
1580After the final iteration, a check-and-correct step is needed to determine
1581whether `a` or `a - 1` gives the desired integer square root of `n`.
1582
1583The algorithm is remarkable in its simplicity. There's no need for a
1584per-iteration check-and-correct step, and termination is straightforward: the
1585number of iterations is known in advance (it's exactly `floor(log2(log2(n)))`
1586for `n > 1`). The only tricky part of the correctness proof is in establishing
1587that the bound `(a - 1)**2 < (n >> s) < (a + 1)**2` is maintained from one
1588iteration to the next. A sketch of the proof of this is given below.
1589
1590In addition to the proof sketch, a formal, computer-verified proof
1591of correctness (using Lean) of an equivalent recursive algorithm can be found
1592here:
1593
1594 https://github.com/mdickinson/snippets/blob/master/proofs/isqrt/src/isqrt.lean
1595
1596
1597Here's Python code equivalent to the C implementation below:
1598
1599 def isqrt(n):
1600 """
1601 Return the integer part of the square root of the input.
1602 """
1603 n = operator.index(n)
1604
1605 if n < 0:
1606 raise ValueError("isqrt() argument must be nonnegative")
1607 if n == 0:
1608 return 0
1609
1610 c = (n.bit_length() - 1) // 2
1611 a = 1
1612 d = 0
1613 for s in reversed(range(c.bit_length())):
1614 # Loop invariant: (a-1)**2 < (n >> 2*(c - d)) < (a+1)**2
1615 e = d
1616 d = c >> s
1617 a = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a
1618
1619 return a - (a*a > n)
1620
1621
1622Sketch of proof of correctness
1623------------------------------
1624
1625The delicate part of the correctness proof is showing that the loop invariant
1626is preserved from one iteration to the next. That is, just before the line
1627
1628 a = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a
1629
1630is executed in the above code, we know that
1631
1632 (1) (a - 1)**2 < (n >> 2*(c - e)) < (a + 1)**2.
1633
1634(since `e` is always the value of `d` from the previous iteration). We must
1635prove that after that line is executed, we have
1636
1637 (a - 1)**2 < (n >> 2*(c - d)) < (a + 1)**2
1638
1639To facilitate the proof, we make some changes of notation. Write `m` for
1640`n >> 2*(c-d)`, and write `b` for the new value of `a`, so
1641
1642 b = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a
1643
1644or equivalently:
1645
1646 (2) b = (a << d - e - 1) + (m >> d - e + 1) // a
1647
1648Then we can rewrite (1) as:
1649
1650 (3) (a - 1)**2 < (m >> 2*(d - e)) < (a + 1)**2
1651
1652and we must show that (b - 1)**2 < m < (b + 1)**2.
1653
1654From this point on, we switch to mathematical notation, so `/` means exact
1655division rather than integer division and `^` is used for exponentiation. We
1656use the `√` symbol for the exact square root. In (3), we can remove the
1657implicit floor operation to give:
1658
1659 (4) (a - 1)^2 < m / 4^(d - e) < (a + 1)^2
1660
1661Taking square roots throughout (4), scaling by `2^(d-e)`, and rearranging gives
1662
1663 (5) 0 <= | 2^(d-e)a - √m | < 2^(d-e)
1664
1665Squaring and dividing through by `2^(d-e+1) a` gives
1666
1667 (6) 0 <= 2^(d-e-1) a + m / (2^(d-e+1) a) - √m < 2^(d-e-1) / a
1668
1669We'll show below that `2^(d-e-1) <= a`. Given that, we can replace the
1670right-hand side of (6) with `1`, and now replacing the central
1671term `m / (2^(d-e+1) a)` with its floor in (6) gives
1672
1673 (7) -1 < 2^(d-e-1) a + m // 2^(d-e+1) a - √m < 1
1674
1675Or equivalently, from (2):
1676
1677 (7) -1 < b - √m < 1
1678
1679and rearranging gives that `(b-1)^2 < m < (b+1)^2`, which is what we needed
1680to prove.
1681
1682We're not quite done: we still have to prove the inequality `2^(d - e - 1) <=
1683a` that was used to get line (7) above. From the definition of `c`, we have
1684`4^c <= n`, which implies
1685
1686 (8) 4^d <= m
1687
1688also, since `e == d >> 1`, `d` is at most `2e + 1`, from which it follows
1689that `2d - 2e - 1 <= d` and hence that
1690
1691 (9) 4^(2d - 2e - 1) <= m
1692
1693Dividing both sides by `4^(d - e)` gives
1694
1695 (10) 4^(d - e - 1) <= m / 4^(d - e)
1696
1697But we know from (4) that `m / 4^(d-e) < (a + 1)^2`, hence
1698
1699 (11) 4^(d - e - 1) < (a + 1)^2
1700
1701Now taking square roots of both sides and observing that both `2^(d-e-1)` and
1702`a` are integers gives `2^(d - e - 1) <= a`, which is what we needed. This
1703completes the proof sketch.
1704
1705*/
1706
1707
1708/* Approximate square root of a large 64-bit integer.
1709
1710 Given `n` satisfying `2**62 <= n < 2**64`, return `a`
1711 satisfying `(a - 1)**2 < n < (a + 1)**2`. */
1712
1713static uint64_t
1714_approximate_isqrt(uint64_t n)
1715{
1716 uint32_t u = 1U + (n >> 62);
1717 u = (u << 1) + (n >> 59) / u;
1718 u = (u << 3) + (n >> 53) / u;
1719 u = (u << 7) + (n >> 41) / u;
1720 return (u << 15) + (n >> 17) / u;
1721}
1722
1723/*[clinic input]
1724math.isqrt
1725
1726 n: object
1727 /
1728
1729Return the integer part of the square root of the input.
1730[clinic start generated code]*/
1731
1732static PyObject *
1733math_isqrt(PyObject *module, PyObject *n)
1734/*[clinic end generated code: output=35a6f7f980beab26 input=5b6e7ae4fa6c43d6]*/
1735{
1736 int a_too_large, c_bit_length;
1737 size_t c, d;
1738 uint64_t m, u;
1739 PyObject *a = NULL, *b;
1740
1741 n = _PyNumber_Index(n);
1742 if (n == NULL) {
1743 return NULL;
1744 }
1745
1746 if (_PyLong_Sign(n) < 0) {
1747 PyErr_SetString(
1748 PyExc_ValueError,
1749 "isqrt() argument must be nonnegative");
1750 goto error;
1751 }
1752 if (_PyLong_Sign(n) == 0) {
1753 Py_DECREF(n);
1754 return PyLong_FromLong(0);
1755 }
1756
1757 /* c = (n.bit_length() - 1) // 2 */
1758 c = _PyLong_NumBits(n);
1759 if (c == (size_t)(-1)) {
1760 goto error;
1761 }
1762 c = (c - 1U) / 2U;
1763
1764 /* Fast path: if c <= 31 then n < 2**64 and we can compute directly with a
1765 fast, almost branch-free algorithm. In the final correction, we use `u*u
1766 - 1 >= m` instead of the simpler `u*u > m` in order to get the correct
1767 result in the corner case where `u=2**32`. */
1768 if (c <= 31U) {
1769 m = (uint64_t)PyLong_AsUnsignedLongLong(n);
1770 Py_DECREF(n);
1771 if (m == (uint64_t)(-1) && PyErr_Occurred()) {
1772 return NULL;
1773 }
1774 u = _approximate_isqrt(m << (62U - 2U*c)) >> (31U - c);
1775 u -= u * u - 1U >= m;
1776 return PyLong_FromUnsignedLongLong((unsigned long long)u);
1777 }
1778
1779 /* Slow path: n >= 2**64. We perform the first five iterations in C integer
1780 arithmetic, then switch to using Python long integers. */
1781
1782 /* From n >= 2**64 it follows that c.bit_length() >= 6. */
1783 c_bit_length = 6;
1784 while ((c >> c_bit_length) > 0U) {
1785 ++c_bit_length;
1786 }
1787
1788 /* Initialise d and a. */
1789 d = c >> (c_bit_length - 5);
1790 b = _PyLong_Rshift(n, 2U*c - 62U);
1791 if (b == NULL) {
1792 goto error;
1793 }
1794 m = (uint64_t)PyLong_AsUnsignedLongLong(b);
1795 Py_DECREF(b);
1796 if (m == (uint64_t)(-1) && PyErr_Occurred()) {
1797 goto error;
1798 }
1799 u = _approximate_isqrt(m) >> (31U - d);
1800 a = PyLong_FromUnsignedLongLong((unsigned long long)u);
1801 if (a == NULL) {
1802 goto error;
1803 }
1804
1805 for (int s = c_bit_length - 6; s >= 0; --s) {
1806 PyObject *q;
1807 size_t e = d;
1808
1809 d = c >> s;
1810
1811 /* q = (n >> 2*c - e - d + 1) // a */
1812 q = _PyLong_Rshift(n, 2U*c - d - e + 1U);
1813 if (q == NULL) {
1814 goto error;
1815 }
1816 Py_SETREF(q, PyNumber_FloorDivide(q, a));
1817 if (q == NULL) {
1818 goto error;
1819 }
1820
1821 /* a = (a << d - 1 - e) + q */
1822 Py_SETREF(a, _PyLong_Lshift(a, d - 1U - e));
1823 if (a == NULL) {
1824 Py_DECREF(q);
1825 goto error;
1826 }
1827 Py_SETREF(a, PyNumber_Add(a, q));
1828 Py_DECREF(q);
1829 if (a == NULL) {
1830 goto error;
1831 }
1832 }
1833
1834 /* The correct result is either a or a - 1. Figure out which, and
1835 decrement a if necessary. */
1836
1837 /* a_too_large = n < a * a */
1838 b = PyNumber_Multiply(a, a);
1839 if (b == NULL) {
1840 goto error;
1841 }
1842 a_too_large = PyObject_RichCompareBool(n, b, Py_LT);
1843 Py_DECREF(b);
1844 if (a_too_large == -1) {
1845 goto error;
1846 }
1847
1848 if (a_too_large) {
1849 Py_SETREF(a, PyNumber_Subtract(a, _PyLong_GetOne()));
1850 }
1851 Py_DECREF(n);
1852 return a;
1853
1854 error:
1855 Py_XDECREF(a);
1856 Py_DECREF(n);
1857 return NULL;
1858}
1859
1860/* Divide-and-conquer factorial algorithm
1861 *
1862 * Based on the formula and pseudo-code provided at:
1863 * http://www.luschny.de/math/factorial/binarysplitfact.html
1864 *
1865 * Faster algorithms exist, but they're more complicated and depend on
1866 * a fast prime factorization algorithm.
1867 *
1868 * Notes on the algorithm
1869 * ----------------------
1870 *
1871 * factorial(n) is written in the form 2**k * m, with m odd. k and m are
1872 * computed separately, and then combined using a left shift.
1873 *
1874 * The function factorial_odd_part computes the odd part m (i.e., the greatest
1875 * odd divisor) of factorial(n), using the formula:
1876 *
1877 * factorial_odd_part(n) =
1878 *
1879 * product_{i >= 0} product_{0 < j <= n / 2**i, j odd} j
1880 *
1881 * Example: factorial_odd_part(20) =
1882 *
1883 * (1) *
1884 * (1) *
1885 * (1 * 3 * 5) *
1886 * (1 * 3 * 5 * 7 * 9) *
1887 * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1888 *
1889 * Here i goes from large to small: the first term corresponds to i=4 (any
1890 * larger i gives an empty product), and the last term corresponds to i=0.
1891 * Each term can be computed from the last by multiplying by the extra odd
1892 * numbers required: e.g., to get from the penultimate term to the last one,
1893 * we multiply by (11 * 13 * 15 * 17 * 19).
1894 *
1895 * To see a hint of why this formula works, here are the same numbers as above
1896 * but with the even parts (i.e., the appropriate powers of 2) included. For
1897 * each subterm in the product for i, we multiply that subterm by 2**i:
1898 *
1899 * factorial(20) =
1900 *
1901 * (16) *
1902 * (8) *
1903 * (4 * 12 * 20) *
1904 * (2 * 6 * 10 * 14 * 18) *
1905 * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1906 *
1907 * The factorial_partial_product function computes the product of all odd j in
1908 * range(start, stop) for given start and stop. It's used to compute the
1909 * partial products like (11 * 13 * 15 * 17 * 19) in the example above. It
1910 * operates recursively, repeatedly splitting the range into two roughly equal
1911 * pieces until the subranges are small enough to be computed using only C
1912 * integer arithmetic.
1913 *
1914 * The two-valuation k (i.e., the exponent of the largest power of 2 dividing
1915 * the factorial) is computed independently in the main math_factorial
1916 * function. By standard results, its value is:
1917 *
1918 * two_valuation = n//2 + n//4 + n//8 + ....
1919 *
1920 * It can be shown (e.g., by complete induction on n) that two_valuation is
1921 * equal to n - count_set_bits(n), where count_set_bits(n) gives the number of
1922 * '1'-bits in the binary expansion of n.
1923 */
1924
1925/* factorial_partial_product: Compute product(range(start, stop, 2)) using
1926 * divide and conquer. Assumes start and stop are odd and stop > start.
1927 * max_bits must be >= bit_length(stop - 2). */
1928
1929static PyObject *
1930factorial_partial_product(unsigned long start, unsigned long stop,
1931 unsigned long max_bits)
1932{
1933 unsigned long midpoint, num_operands;
1934 PyObject *left = NULL, *right = NULL, *result = NULL;
1935
1936 /* If the return value will fit an unsigned long, then we can
1937 * multiply in a tight, fast loop where each multiply is O(1).
1938 * Compute an upper bound on the number of bits required to store
1939 * the answer.
1940 *
1941 * Storing some integer z requires floor(lg(z))+1 bits, which is
1942 * conveniently the value returned by bit_length(z). The
1943 * product x*y will require at most
1944 * bit_length(x) + bit_length(y) bits to store, based
1945 * on the idea that lg product = lg x + lg y.
1946 *
1947 * We know that stop - 2 is the largest number to be multiplied. From
1948 * there, we have: bit_length(answer) <= num_operands *
1949 * bit_length(stop - 2)
1950 */
1951
1952 num_operands = (stop - start) / 2;
1953 /* The "num_operands <= 8 * SIZEOF_LONG" check guards against the
1954 * unlikely case of an overflow in num_operands * max_bits. */
1955 if (num_operands <= 8 * SIZEOF_LONG &&
1956 num_operands * max_bits <= 8 * SIZEOF_LONG) {
1957 unsigned long j, total;
1958 for (total = start, j = start + 2; j < stop; j += 2)
1959 total *= j;
1960 return PyLong_FromUnsignedLong(total);
1961 }
1962
1963 /* find midpoint of range(start, stop), rounded up to next odd number. */
1964 midpoint = (start + num_operands) | 1;
1965 left = factorial_partial_product(start, midpoint,
1966 _Py_bit_length(midpoint - 2));
1967 if (left == NULL)
1968 goto error;
1969 right = factorial_partial_product(midpoint, stop, max_bits);
1970 if (right == NULL)
1971 goto error;
1972 result = PyNumber_Multiply(left, right);
1973
1974 error:
1975 Py_XDECREF(left);
1976 Py_XDECREF(right);
1977 return result;
1978}
1979
1980/* factorial_odd_part: compute the odd part of factorial(n). */
1981
1982static PyObject *
1983factorial_odd_part(unsigned long n)
1984{
1985 long i;
1986 unsigned long v, lower, upper;
1987 PyObject *partial, *tmp, *inner, *outer;
1988
1989 inner = PyLong_FromLong(1);
1990 if (inner == NULL)
1991 return NULL;
1992 outer = inner;
1993 Py_INCREF(outer);
1994
1995 upper = 3;
1996 for (i = _Py_bit_length(n) - 2; i >= 0; i--) {
1997 v = n >> i;
1998 if (v <= 2)
1999 continue;
2000 lower = upper;
2001 /* (v + 1) | 1 = least odd integer strictly larger than n / 2**i */
2002 upper = (v + 1) | 1;
2003 /* Here inner is the product of all odd integers j in the range (0,
2004 n/2**(i+1)]. The factorial_partial_product call below gives the
2005 product of all odd integers j in the range (n/2**(i+1), n/2**i]. */
2006 partial = factorial_partial_product(lower, upper, _Py_bit_length(upper-2));
2007 /* inner *= partial */
2008 if (partial == NULL)
2009 goto error;
2010 tmp = PyNumber_Multiply(inner, partial);
2011 Py_DECREF(partial);
2012 if (tmp == NULL)
2013 goto error;
2014 Py_DECREF(inner);
2015 inner = tmp;
2016 /* Now inner is the product of all odd integers j in the range (0,
2017 n/2**i], giving the inner product in the formula above. */
2018
2019 /* outer *= inner; */
2020 tmp = PyNumber_Multiply(outer, inner);
2021 if (tmp == NULL)
2022 goto error;
2023 Py_DECREF(outer);
2024 outer = tmp;
2025 }
2026 Py_DECREF(inner);
2027 return outer;
2028
2029 error:
2030 Py_DECREF(outer);
2031 Py_DECREF(inner);
2032 return NULL;
2033}
2034
2035
2036/* Lookup table for small factorial values */
2037
2038static const unsigned long SmallFactorials[] = {
2039 1, 1, 2, 6, 24, 120, 720, 5040, 40320,
2040 362880, 3628800, 39916800, 479001600,
2041#if SIZEOF_LONG >= 8
2042 6227020800, 87178291200, 1307674368000,
2043 20922789888000, 355687428096000, 6402373705728000,
2044 121645100408832000, 2432902008176640000
2045#endif
2046};
2047
2048/*[clinic input]
2049math.factorial
2050
2051 x as arg: object
2052 /
2053
2054Find x!.
2055
2056Raise a ValueError if x is negative or non-integral.
2057[clinic start generated code]*/
2058
2059static PyObject *
2060math_factorial(PyObject *module, PyObject *arg)
2061/*[clinic end generated code: output=6686f26fae00e9ca input=6d1c8105c0d91fb4]*/
2062{
2063 long x, two_valuation;
2064 int overflow;
2065 PyObject *result, *odd_part;
2066
2067 x = PyLong_AsLongAndOverflow(arg, &overflow);
2068 if (x == -1 && PyErr_Occurred()) {
2069 return NULL;
2070 }
2071 else if (overflow == 1) {
2072 PyErr_Format(PyExc_OverflowError,
2073 "factorial() argument should not exceed %ld",
2074 LONG_MAX);
2075 return NULL;
2076 }
2077 else if (overflow == -1 || x < 0) {
2078 PyErr_SetString(PyExc_ValueError,
2079 "factorial() not defined for negative values");
2080 return NULL;
2081 }
2082
2083 /* use lookup table if x is small */
2084 if (x < (long)Py_ARRAY_LENGTH(SmallFactorials))
2085 return PyLong_FromUnsignedLong(SmallFactorials[x]);
2086
2087 /* else express in the form odd_part * 2**two_valuation, and compute as
2088 odd_part << two_valuation. */
2089 odd_part = factorial_odd_part(x);
2090 if (odd_part == NULL)
2091 return NULL;
2092 two_valuation = x - count_set_bits(x);
2093 result = _PyLong_Lshift(odd_part, two_valuation);
2094 Py_DECREF(odd_part);
2095 return result;
2096}
2097
2098
2099/*[clinic input]
2100math.trunc
2101
2102 x: object
2103 /
2104
2105Truncates the Real x to the nearest Integral toward 0.
2106
2107Uses the __trunc__ magic method.
2108[clinic start generated code]*/
2109
2110static PyObject *
2111math_trunc(PyObject *module, PyObject *x)
2112/*[clinic end generated code: output=34b9697b707e1031 input=2168b34e0a09134d]*/
2113{
2114 _Py_IDENTIFIER(__trunc__);
2115 PyObject *trunc, *result;
2116
2117 if (PyFloat_CheckExact(x)) {
2118 return PyFloat_Type.tp_as_number->nb_int(x);
2119 }
2120
2121 if (Py_TYPE(x)->tp_dict == NULL) {
2122 if (PyType_Ready(Py_TYPE(x)) < 0)
2123 return NULL;
2124 }
2125
2126 trunc = _PyObject_LookupSpecial(x, &PyId___trunc__);
2127 if (trunc == NULL) {
2128 if (!PyErr_Occurred())
2129 PyErr_Format(PyExc_TypeError,
2130 "type %.100s doesn't define __trunc__ method",
2131 Py_TYPE(x)->tp_name);
2132 return NULL;
2133 }
2134 result = _PyObject_CallNoArg(trunc);
2135 Py_DECREF(trunc);
2136 return result;
2137}
2138
2139
2140/*[clinic input]
2141math.frexp
2142
2143 x: double
2144 /
2145
2146Return the mantissa and exponent of x, as pair (m, e).
2147
2148m is a float and e is an int, such that x = m * 2.**e.
2149If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0.
2150[clinic start generated code]*/
2151
2152static PyObject *
2153math_frexp_impl(PyObject *module, double x)
2154/*[clinic end generated code: output=03e30d252a15ad4a input=96251c9e208bc6e9]*/
2155{
2156 int i;
2157 /* deal with special cases directly, to sidestep platform
2158 differences */
2159 if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
2160 i = 0;
2161 }
2162 else {
2163 x = frexp(x, &i);
2164 }
2165 return Py_BuildValue("(di)", x, i);
2166}
2167
2168
2169/*[clinic input]
2170math.ldexp
2171
2172 x: double
2173 i: object
2174 /
2175
2176Return x * (2**i).
2177
2178This is essentially the inverse of frexp().
2179[clinic start generated code]*/
2180
2181static PyObject *
2182math_ldexp_impl(PyObject *module, double x, PyObject *i)
2183/*[clinic end generated code: output=b6892f3c2df9cc6a input=17d5970c1a40a8c1]*/
2184{
2185 double r;
2186 long exp;
2187 int overflow;
2188
2189 if (PyLong_Check(i)) {
2190 /* on overflow, replace exponent with either LONG_MAX
2191 or LONG_MIN, depending on the sign. */
2192 exp = PyLong_AsLongAndOverflow(i, &overflow);
2193 if (exp == -1 && PyErr_Occurred())
2194 return NULL;
2195 if (overflow)
2196 exp = overflow < 0 ? LONG_MIN : LONG_MAX;
2197 }
2198 else {
2199 PyErr_SetString(PyExc_TypeError,
2200 "Expected an int as second argument to ldexp.");
2201 return NULL;
2202 }
2203
2204 if (x == 0. || !Py_IS_FINITE(x)) {
2205 /* NaNs, zeros and infinities are returned unchanged */
2206 r = x;
2207 errno = 0;
2208 } else if (exp > INT_MAX) {
2209 /* overflow */
2210 r = copysign(Py_HUGE_VAL, x);
2211 errno = ERANGE;
2212 } else if (exp < INT_MIN) {
2213 /* underflow to +-0 */
2214 r = copysign(0., x);
2215 errno = 0;
2216 } else {
2217 errno = 0;
2218 r = ldexp(x, (int)exp);
2219 if (Py_IS_INFINITY(r))
2220 errno = ERANGE;
2221 }
2222
2223 if (errno && is_error(r))
2224 return NULL;
2225 return PyFloat_FromDouble(r);
2226}
2227
2228
2229/*[clinic input]
2230math.modf
2231
2232 x: double
2233 /
2234
2235Return the fractional and integer parts of x.
2236
2237Both results carry the sign of x and are floats.
2238[clinic start generated code]*/
2239
2240static PyObject *
2241math_modf_impl(PyObject *module, double x)
2242/*[clinic end generated code: output=90cee0260014c3c0 input=b4cfb6786afd9035]*/
2243{
2244 double y;
2245 /* some platforms don't do the right thing for NaNs and
2246 infinities, so we take care of special cases directly. */
2247 if (!Py_IS_FINITE(x)) {
2248 if (Py_IS_INFINITY(x))
2249 return Py_BuildValue("(dd)", copysign(0., x), x);
2250 else if (Py_IS_NAN(x))
2251 return Py_BuildValue("(dd)", x, x);
2252 }
2253
2254 errno = 0;
2255 x = modf(x, &y);
2256 return Py_BuildValue("(dd)", x, y);
2257}
2258
2259
2260/* A decent logarithm is easy to compute even for huge ints, but libm can't
2261 do that by itself -- loghelper can. func is log or log10, and name is
2262 "log" or "log10". Note that overflow of the result isn't possible: an int
2263 can contain no more than INT_MAX * SHIFT bits, so has value certainly less
2264 than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
2265 small enough to fit in an IEEE single. log and log10 are even smaller.
2266 However, intermediate overflow is possible for an int if the number of bits
2267 in that int is larger than PY_SSIZE_T_MAX. */
2268
2269static PyObject*
2270loghelper(PyObject* arg, double (*func)(double), const char *funcname)
2271{
2272 /* If it is int, do it ourselves. */
2273 if (PyLong_Check(arg)) {
2274 double x, result;
2275 Py_ssize_t e;
2276
2277 /* Negative or zero inputs give a ValueError. */
2278 if (Py_SIZE(arg) <= 0) {
2279 PyErr_SetString(PyExc_ValueError,
2280 "math domain error");
2281 return NULL;
2282 }
2283
2284 x = PyLong_AsDouble(arg);
2285 if (x == -1.0 && PyErr_Occurred()) {
2286 if (!PyErr_ExceptionMatches(PyExc_OverflowError))
2287 return NULL;
2288 /* Here the conversion to double overflowed, but it's possible
2289 to compute the log anyway. Clear the exception and continue. */
2290 PyErr_Clear();
2291 x = _PyLong_Frexp((PyLongObject *)arg, &e);
2292 if (x == -1.0 && PyErr_Occurred())
2293 return NULL;
2294 /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */
2295 result = func(x) + func(2.0) * e;
2296 }
2297 else
2298 /* Successfully converted x to a double. */
2299 result = func(x);
2300 return PyFloat_FromDouble(result);
2301 }
2302
2303 /* Else let libm handle it by itself. */
2304 return math_1(arg, func, 0);
2305}
2306
2307
2308/*[clinic input]
2309math.log
2310
2311 x: object
2312 [
2313 base: object(c_default="NULL") = math.e
2314 ]
2315 /
2316
2317Return the logarithm of x to the given base.
2318
2319If the base not specified, returns the natural logarithm (base e) of x.
2320[clinic start generated code]*/
2321
2322static PyObject *
2323math_log_impl(PyObject *module, PyObject *x, int group_right_1,
2324 PyObject *base)
2325/*[clinic end generated code: output=7b5a39e526b73fc9 input=0f62d5726cbfebbd]*/
2326{
2327 PyObject *num, *den;
2328 PyObject *ans;
2329
2330 num = loghelper(x, m_log, "log");
2331 if (num == NULL || base == NULL)
2332 return num;
2333
2334 den = loghelper(base, m_log, "log");
2335 if (den == NULL) {
2336 Py_DECREF(num);
2337 return NULL;
2338 }
2339
2340 ans = PyNumber_TrueDivide(num, den);
2341 Py_DECREF(num);
2342 Py_DECREF(den);
2343 return ans;
2344}
2345
2346
2347/*[clinic input]
2348math.log2
2349
2350 x: object
2351 /
2352
2353Return the base 2 logarithm of x.
2354[clinic start generated code]*/
2355
2356static PyObject *
2357math_log2(PyObject *module, PyObject *x)
2358/*[clinic end generated code: output=5425899a4d5d6acb input=08321262bae4f39b]*/
2359{
2360 return loghelper(x, m_log2, "log2");
2361}
2362
2363
2364/*[clinic input]
2365math.log10
2366
2367 x: object
2368 /
2369
2370Return the base 10 logarithm of x.
2371[clinic start generated code]*/
2372
2373static PyObject *
2374math_log10(PyObject *module, PyObject *x)
2375/*[clinic end generated code: output=be72a64617df9c6f input=b2469d02c6469e53]*/
2376{
2377 return loghelper(x, m_log10, "log10");
2378}
2379
2380
2381/*[clinic input]
2382math.fmod
2383
2384 x: double
2385 y: double
2386 /
2387
2388Return fmod(x, y), according to platform C.
2389
2390x % y may differ.
2391[clinic start generated code]*/
2392
2393static PyObject *
2394math_fmod_impl(PyObject *module, double x, double y)
2395/*[clinic end generated code: output=7559d794343a27b5 input=4f84caa8cfc26a03]*/
2396{
2397 double r;
2398 /* fmod(x, +/-Inf) returns x for finite x. */
2399 if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
2400 return PyFloat_FromDouble(x);
2401 errno = 0;
2402 r = fmod(x, y);
2403 if (Py_IS_NAN(r)) {
2404 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
2405 errno = EDOM;
2406 else
2407 errno = 0;
2408 }
2409 if (errno && is_error(r))
2410 return NULL;
2411 else
2412 return PyFloat_FromDouble(r);
2413}
2414
2415/*
2416Given a *vec* of values, compute the vector norm:
2417
2418 sqrt(sum(x ** 2 for x in vec))
2419
2420The *max* variable should be equal to the largest fabs(x).
2421The *n* variable is the length of *vec*.
2422If n==0, then *max* should be 0.0.
2423If an infinity is present in the vec, *max* should be INF.
2424The *found_nan* variable indicates whether some member of
2425the *vec* is a NaN.
2426
2427To avoid overflow/underflow and to achieve high accuracy giving results
2428that are almost always correctly rounded, four techniques are used:
2429
2430* lossless scaling using a power-of-two scaling factor
2431* accurate squaring using Veltkamp-Dekker splitting [1]
2432* compensated summation using a variant of the Neumaier algorithm [2]
2433* differential correction of the square root [3]
2434
2435The usual presentation of the Neumaier summation algorithm has an
2436expensive branch depending on which operand has the larger
2437magnitude. We avoid this cost by arranging the calculation so that
2438fabs(csum) is always as large as fabs(x).
2439
2440To establish the invariant, *csum* is initialized to 1.0 which is
2441always larger than x**2 after scaling or after division by *max*.
2442After the loop is finished, the initial 1.0 is subtracted out for a
2443net zero effect on the final sum. Since *csum* will be greater than
24441.0, the subtraction of 1.0 will not cause fractional digits to be
2445dropped from *csum*.
2446
2447To get the full benefit from compensated summation, the largest
2448addend should be in the range: 0.5 <= |x| <= 1.0. Accordingly,
2449scaling or division by *max* should not be skipped even if not
2450otherwise needed to prevent overflow or loss of precision.
2451
2452The assertion that hi*hi <= 1.0 is a bit subtle. Each vector element
2453gets scaled to a magnitude below 1.0. The Veltkamp-Dekker splitting
2454algorithm gives a *hi* value that is correctly rounded to half
2455precision. When a value at or below 1.0 is correctly rounded, it
2456never goes above 1.0. And when values at or below 1.0 are squared,
2457they remain at or below 1.0, thus preserving the summation invariant.
2458
2459Another interesting assertion is that csum+lo*lo == csum. In the loop,
2460each scaled vector element has a magnitude less than 1.0. After the
2461Veltkamp split, *lo* has a maximum value of 2**-27. So the maximum
2462value of *lo* squared is 2**-54. The value of ulp(1.0)/2.0 is 2**-53.
2463Given that csum >= 1.0, we have:
2464 lo**2 <= 2**-54 < 2**-53 == 1/2*ulp(1.0) <= ulp(csum)/2
2465Since lo**2 is less than 1/2 ulp(csum), we have csum+lo*lo == csum.
2466
2467To minimize loss of information during the accumulation of fractional
2468values, each term has a separate accumulator. This also breaks up
2469sequential dependencies in the inner loop so the CPU can maximize
2470floating point throughput. [4] On a 2.6 GHz Haswell, adding one
2471dimension has an incremental cost of only 5ns -- for example when
2472moving from hypot(x,y) to hypot(x,y,z).
2473
2474The square root differential correction is needed because a
2475correctly rounded square root of a correctly rounded sum of
2476squares can still be off by as much as one ulp.
2477
2478The differential correction starts with a value *x* that is
2479the difference between the square of *h*, the possibly inaccurately
2480rounded square root, and the accurately computed sum of squares.
2481The correction is the first order term of the Maclaurin series
2482expansion of sqrt(h**2 + x) == h + x/(2*h) + O(x**2). [5]
2483
2484Essentially, this differential correction is equivalent to one
2485refinement step in Newton's divide-and-average square root
2486algorithm, effectively doubling the number of accurate bits.
2487This technique is used in Dekker's SQRT2 algorithm and again in
2488Borges' ALGORITHM 4 and 5.
2489
2490Without proof for all cases, hypot() cannot claim to be always
2491correctly rounded. However for n <= 1000, prior to the final addition
2492that rounds the overall result, the internal accuracy of "h" together
2493with its correction of "x / (2.0 * h)" is at least 100 bits. [6]
2494Also, hypot() was tested against a Decimal implementation with
2495prec=300. After 100 million trials, no incorrectly rounded examples
2496were found. In addition, perfect commutativity (all permutations are
2497exactly equal) was verified for 1 billion random inputs with n=5. [7]
2498
2499References:
2500
25011. Veltkamp-Dekker splitting: http://csclub.uwaterloo.ca/~pbarfuss/dekker1971.pdf
25022. Compensated summation: http://www.ti3.tu-harburg.de/paper/rump/Ru08b.pdf
25033. Square root differential correction: https://arxiv.org/pdf/1904.09481.pdf
25044. Data dependency graph: https://bugs.python.org/file49439/hypot.png
25055. https://www.wolframalpha.com/input/?i=Maclaurin+series+sqrt%28h**2+%2B+x%29+at+x%3D0
25066. Analysis of internal accuracy: https://bugs.python.org/file49484/best_frac.py
25077. Commutativity test: https://bugs.python.org/file49448/test_hypot_commutativity.py
2508
2509*/
2510
2511static inline double
2512vector_norm(Py_ssize_t n, double *vec, double max, int found_nan)
2513{
2514 const double T27 = 134217729.0; /* ldexp(1.0, 27) + 1.0) */
2515 double x, scale, oldcsum, csum = 1.0, frac1 = 0.0, frac2 = 0.0, frac3 = 0.0;
2516 double t, hi, lo, h;
2517 int max_e;
2518 Py_ssize_t i;
2519
2520 if (Py_IS_INFINITY(max)) {
2521 return max;
2522 }
2523 if (found_nan) {
2524 return Py_NAN;
2525 }
2526 if (max == 0.0 || n <= 1) {
2527 return max;
2528 }
2529 frexp(max, &max_e);
2530 if (max_e >= -1023) {
2531 scale = ldexp(1.0, -max_e);
2532 assert(max * scale >= 0.5);
2533 assert(max * scale < 1.0);
2534 for (i=0 ; i < n ; i++) {
2535 x = vec[i];
2536 assert(Py_IS_FINITE(x) && fabs(x) <= max);
2537
2538 x *= scale;
2539 assert(fabs(x) < 1.0);
2540
2541 t = x * T27;
2542 hi = t - (t - x);
2543 lo = x - hi;
2544 assert(hi + lo == x);
2545
2546 x = hi * hi;
2547 assert(x <= 1.0);
2548 assert(fabs(csum) >= fabs(x));
2549 oldcsum = csum;
2550 csum += x;
2551 frac1 += (oldcsum - csum) + x;
2552
2553 x = 2.0 * hi * lo;
2554 assert(fabs(csum) >= fabs(x));
2555 oldcsum = csum;
2556 csum += x;
2557 frac2 += (oldcsum - csum) + x;
2558
2559 assert(csum + lo * lo == csum);
2560 frac3 += lo * lo;
2561 }
2562 h = sqrt(csum - 1.0 + (frac1 + frac2 + frac3));
2563
2564 x = h;
2565 t = x * T27;
2566 hi = t - (t - x);
2567 lo = x - hi;
2568 assert (hi + lo == x);
2569
2570 x = -hi * hi;
2571 assert(fabs(csum) >= fabs(x));
2572 oldcsum = csum;
2573 csum += x;
2574 frac1 += (oldcsum - csum) + x;
2575
2576 x = -2.0 * hi * lo;
2577 assert(fabs(csum) >= fabs(x));
2578 oldcsum = csum;
2579 csum += x;
2580 frac2 += (oldcsum - csum) + x;
2581
2582 x = -lo * lo;
2583 assert(fabs(csum) >= fabs(x));
2584 oldcsum = csum;
2585 csum += x;
2586 frac3 += (oldcsum - csum) + x;
2587
2588 x = csum - 1.0 + (frac1 + frac2 + frac3);
2589 return (h + x / (2.0 * h)) / scale;
2590 }
2591 /* When max_e < -1023, ldexp(1.0, -max_e) overflows.
2592 So instead of multiplying by a scale, we just divide by *max*.
2593 */
2594 for (i=0 ; i < n ; i++) {
2595 x = vec[i];
2596 assert(Py_IS_FINITE(x) && fabs(x) <= max);
2597 x /= max;
2598 x = x*x;
2599 assert(x <= 1.0);
2600 assert(fabs(csum) >= fabs(x));
2601 oldcsum = csum;
2602 csum += x;
2603 frac1 += (oldcsum - csum) + x;
2604 }
2605 return max * sqrt(csum - 1.0 + frac1);
2606}
2607
2608#define NUM_STACK_ELEMS 16
2609
2610/*[clinic input]
2611math.dist
2612
2613 p: object
2614 q: object
2615 /
2616
2617Return the Euclidean distance between two points p and q.
2618
2619The points should be specified as sequences (or iterables) of
2620coordinates. Both inputs must have the same dimension.
2621
2622Roughly equivalent to:
2623 sqrt(sum((px - qx) ** 2.0 for px, qx in zip(p, q)))
2624[clinic start generated code]*/
2625
2626static PyObject *
2627math_dist_impl(PyObject *module, PyObject *p, PyObject *q)
2628/*[clinic end generated code: output=56bd9538d06bbcfe input=74e85e1b6092e68e]*/
2629{
2630 PyObject *item;
2631 double max = 0.0;
2632 double x, px, qx, result;
2633 Py_ssize_t i, m, n;
2634 int found_nan = 0, p_allocated = 0, q_allocated = 0;
2635 double diffs_on_stack[NUM_STACK_ELEMS];
2636 double *diffs = diffs_on_stack;
2637
2638 if (!PyTuple_Check(p)) {
2639 p = PySequence_Tuple(p);
2640 if (p == NULL) {
2641 return NULL;
2642 }
2643 p_allocated = 1;
2644 }
2645 if (!PyTuple_Check(q)) {
2646 q = PySequence_Tuple(q);
2647 if (q == NULL) {
2648 if (p_allocated) {
2649 Py_DECREF(p);
2650 }
2651 return NULL;
2652 }
2653 q_allocated = 1;
2654 }
2655
2656 m = PyTuple_GET_SIZE(p);
2657 n = PyTuple_GET_SIZE(q);
2658 if (m != n) {
2659 PyErr_SetString(PyExc_ValueError,
2660 "both points must have the same number of dimensions");
2661 return NULL;
2662
2663 }
2664 if (n > NUM_STACK_ELEMS) {
2665 diffs = (double *) PyObject_Malloc(n * sizeof(double));
2666 if (diffs == NULL) {
2667 return PyErr_NoMemory();
2668 }
2669 }
2670 for (i=0 ; i<n ; i++) {
2671 item = PyTuple_GET_ITEM(p, i);
2672 ASSIGN_DOUBLE(px, item, error_exit);
2673 item = PyTuple_GET_ITEM(q, i);
2674 ASSIGN_DOUBLE(qx, item, error_exit);
2675 x = fabs(px - qx);
2676 diffs[i] = x;
2677 found_nan |= Py_IS_NAN(x);
2678 if (x > max) {
2679 max = x;
2680 }
2681 }
2682 result = vector_norm(n, diffs, max, found_nan);
2683 if (diffs != diffs_on_stack) {
2684 PyObject_Free(diffs);
2685 }
2686 if (p_allocated) {
2687 Py_DECREF(p);
2688 }
2689 if (q_allocated) {
2690 Py_DECREF(q);
2691 }
2692 return PyFloat_FromDouble(result);
2693
2694 error_exit:
2695 if (diffs != diffs_on_stack) {
2696 PyObject_Free(diffs);
2697 }
2698 if (p_allocated) {
2699 Py_DECREF(p);
2700 }
2701 if (q_allocated) {
2702 Py_DECREF(q);
2703 }
2704 return NULL;
2705}
2706
2707/* AC: cannot convert yet, waiting for *args support */
2708static PyObject *
2709math_hypot(PyObject *self, PyObject *const *args, Py_ssize_t nargs)
2710{
2711 Py_ssize_t i;
2712 PyObject *item;
2713 double max = 0.0;
2714 double x, result;
2715 int found_nan = 0;
2716 double coord_on_stack[NUM_STACK_ELEMS];
2717 double *coordinates = coord_on_stack;
2718
2719 if (nargs > NUM_STACK_ELEMS) {
2720 coordinates = (double *) PyObject_Malloc(nargs * sizeof(double));
2721 if (coordinates == NULL) {
2722 return PyErr_NoMemory();
2723 }
2724 }
2725 for (i = 0; i < nargs; i++) {
2726 item = args[i];
2727 ASSIGN_DOUBLE(x, item, error_exit);
2728 x = fabs(x);
2729 coordinates[i] = x;
2730 found_nan |= Py_IS_NAN(x);
2731 if (x > max) {
2732 max = x;
2733 }
2734 }
2735 result = vector_norm(nargs, coordinates, max, found_nan);
2736 if (coordinates != coord_on_stack) {
2737 PyObject_Free(coordinates);
2738 }
2739 return PyFloat_FromDouble(result);
2740
2741 error_exit:
2742 if (coordinates != coord_on_stack) {
2743 PyObject_Free(coordinates);
2744 }
2745 return NULL;
2746}
2747
2748#undef NUM_STACK_ELEMS
2749
2750PyDoc_STRVAR(math_hypot_doc,
2751 "hypot(*coordinates) -> value\n\n\
2752Multidimensional Euclidean distance from the origin to a point.\n\
2753\n\
2754Roughly equivalent to:\n\
2755 sqrt(sum(x**2 for x in coordinates))\n\
2756\n\
2757For a two dimensional point (x, y), gives the hypotenuse\n\
2758using the Pythagorean theorem: sqrt(x*x + y*y).\n\
2759\n\
2760For example, the hypotenuse of a 3/4/5 right triangle is:\n\
2761\n\
2762 >>> hypot(3.0, 4.0)\n\
2763 5.0\n\
2764");
2765
2766/* pow can't use math_2, but needs its own wrapper: the problem is
2767 that an infinite result can arise either as a result of overflow
2768 (in which case OverflowError should be raised) or as a result of
2769 e.g. 0.**-5. (for which ValueError needs to be raised.)
2770*/
2771
2772/*[clinic input]
2773math.pow
2774
2775 x: double
2776 y: double
2777 /
2778
2779Return x**y (x to the power of y).
2780[clinic start generated code]*/
2781
2782static PyObject *
2783math_pow_impl(PyObject *module, double x, double y)
2784/*[clinic end generated code: output=fff93e65abccd6b0 input=c26f1f6075088bfd]*/
2785{
2786 double r;
2787 int odd_y;
2788
2789 /* deal directly with IEEE specials, to cope with problems on various
2790 platforms whose semantics don't exactly match C99 */
2791 r = 0.; /* silence compiler warning */
2792 if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
2793 errno = 0;
2794 if (Py_IS_NAN(x))
2795 r = y == 0. ? 1. : x; /* NaN**0 = 1 */
2796 else if (Py_IS_NAN(y))
2797 r = x == 1. ? 1. : y; /* 1**NaN = 1 */
2798 else if (Py_IS_INFINITY(x)) {
2799 odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
2800 if (y > 0.)
2801 r = odd_y ? x : fabs(x);
2802 else if (y == 0.)
2803 r = 1.;
2804 else /* y < 0. */
2805 r = odd_y ? copysign(0., x) : 0.;
2806 }
2807 else if (Py_IS_INFINITY(y)) {
2808 if (fabs(x) == 1.0)
2809 r = 1.;
2810 else if (y > 0. && fabs(x) > 1.0)
2811 r = y;
2812 else if (y < 0. && fabs(x) < 1.0) {
2813 r = -y; /* result is +inf */
2814 if (x == 0.) /* 0**-inf: divide-by-zero */
2815 errno = EDOM;
2816 }
2817 else
2818 r = 0.;
2819 }
2820 }
2821 else {
2822 /* let libm handle finite**finite */
2823 errno = 0;
2824 r = pow(x, y);
2825 /* a NaN result should arise only from (-ve)**(finite
2826 non-integer); in this case we want to raise ValueError. */
2827 if (!Py_IS_FINITE(r)) {
2828 if (Py_IS_NAN(r)) {
2829 errno = EDOM;
2830 }
2831 /*
2832 an infinite result here arises either from:
2833 (A) (+/-0.)**negative (-> divide-by-zero)
2834 (B) overflow of x**y with x and y finite
2835 */
2836 else if (Py_IS_INFINITY(r)) {
2837 if (x == 0.)
2838 errno = EDOM;
2839 else
2840 errno = ERANGE;
2841 }
2842 }
2843 }
2844
2845 if (errno && is_error(r))
2846 return NULL;
2847 else
2848 return PyFloat_FromDouble(r);
2849}
2850
2851
2852static const double degToRad = Py_MATH_PI / 180.0;
2853static const double radToDeg = 180.0 / Py_MATH_PI;
2854
2855/*[clinic input]
2856math.degrees
2857
2858 x: double
2859 /
2860
2861Convert angle x from radians to degrees.
2862[clinic start generated code]*/
2863
2864static PyObject *
2865math_degrees_impl(PyObject *module, double x)
2866/*[clinic end generated code: output=7fea78b294acd12f input=81e016555d6e3660]*/
2867{
2868 return PyFloat_FromDouble(x * radToDeg);
2869}
2870
2871
2872/*[clinic input]
2873math.radians
2874
2875 x: double
2876 /
2877
2878Convert angle x from degrees to radians.
2879[clinic start generated code]*/
2880
2881static PyObject *
2882math_radians_impl(PyObject *module, double x)
2883/*[clinic end generated code: output=34daa47caf9b1590 input=91626fc489fe3d63]*/
2884{
2885 return PyFloat_FromDouble(x * degToRad);
2886}
2887
2888
2889/*[clinic input]
2890math.isfinite
2891
2892 x: double
2893 /
2894
2895Return True if x is neither an infinity nor a NaN, and False otherwise.
2896[clinic start generated code]*/
2897
2898static PyObject *
2899math_isfinite_impl(PyObject *module, double x)
2900/*[clinic end generated code: output=8ba1f396440c9901 input=46967d254812e54a]*/
2901{
2902 return PyBool_FromLong((long)Py_IS_FINITE(x));
2903}
2904
2905
2906/*[clinic input]
2907math.isnan
2908
2909 x: double
2910 /
2911
2912Return True if x is a NaN (not a number), and False otherwise.
2913[clinic start generated code]*/
2914
2915static PyObject *
2916math_isnan_impl(PyObject *module, double x)
2917/*[clinic end generated code: output=f537b4d6df878c3e input=935891e66083f46a]*/
2918{
2919 return PyBool_FromLong((long)Py_IS_NAN(x));
2920}
2921
2922
2923/*[clinic input]
2924math.isinf
2925
2926 x: double
2927 /
2928
2929Return True if x is a positive or negative infinity, and False otherwise.
2930[clinic start generated code]*/
2931
2932static PyObject *
2933math_isinf_impl(PyObject *module, double x)
2934/*[clinic end generated code: output=9f00cbec4de7b06b input=32630e4212cf961f]*/
2935{
2936 return PyBool_FromLong((long)Py_IS_INFINITY(x));
2937}
2938
2939
2940/*[clinic input]
2941math.isclose -> bool
2942
2943 a: double
2944 b: double
2945 *
2946 rel_tol: double = 1e-09
2947 maximum difference for being considered "close", relative to the
2948 magnitude of the input values
2949 abs_tol: double = 0.0
2950 maximum difference for being considered "close", regardless of the
2951 magnitude of the input values
2952
2953Determine whether two floating point numbers are close in value.
2954
2955Return True if a is close in value to b, and False otherwise.
2956
2957For the values to be considered close, the difference between them
2958must be smaller than at least one of the tolerances.
2959
2960-inf, inf and NaN behave similarly to the IEEE 754 Standard. That
2961is, NaN is not close to anything, even itself. inf and -inf are
2962only close to themselves.
2963[clinic start generated code]*/
2964
2965static int
2966math_isclose_impl(PyObject *module, double a, double b, double rel_tol,
2967 double abs_tol)
2968/*[clinic end generated code: output=b73070207511952d input=f28671871ea5bfba]*/
2969{
2970 double diff = 0.0;
2971
2972 /* sanity check on the inputs */
2973 if (rel_tol < 0.0 || abs_tol < 0.0 ) {
2974 PyErr_SetString(PyExc_ValueError,
2975 "tolerances must be non-negative");
2976 return -1;
2977 }
2978
2979 if ( a == b ) {
2980 /* short circuit exact equality -- needed to catch two infinities of
2981 the same sign. And perhaps speeds things up a bit sometimes.
2982 */
2983 return 1;
2984 }
2985
2986 /* This catches the case of two infinities of opposite sign, or
2987 one infinity and one finite number. Two infinities of opposite
2988 sign would otherwise have an infinite relative tolerance.
2989 Two infinities of the same sign are caught by the equality check
2990 above.
2991 */
2992
2993 if (Py_IS_INFINITY(a) || Py_IS_INFINITY(b)) {
2994 return 0;
2995 }
2996
2997 /* now do the regular computation
2998 this is essentially the "weak" test from the Boost library
2999 */
3000
3001 diff = fabs(b - a);
3002
3003 return (((diff <= fabs(rel_tol * b)) ||
3004 (diff <= fabs(rel_tol * a))) ||
3005 (diff <= abs_tol));
3006}
3007
3008static inline int
3009_check_long_mult_overflow(long a, long b) {
3010
3011 /* From Python2's int_mul code:
3012
3013 Integer overflow checking for * is painful: Python tried a couple ways, but
3014 they didn't work on all platforms, or failed in endcases (a product of
3015 -sys.maxint-1 has been a particular pain).
3016
3017 Here's another way:
3018
3019 The native long product x*y is either exactly right or *way* off, being
3020 just the last n bits of the true product, where n is the number of bits
3021 in a long (the delivered product is the true product plus i*2**n for
3022 some integer i).
3023
3024 The native double product (double)x * (double)y is subject to three
3025 rounding errors: on a sizeof(long)==8 box, each cast to double can lose
3026 info, and even on a sizeof(long)==4 box, the multiplication can lose info.
3027 But, unlike the native long product, it's not in *range* trouble: even
3028 if sizeof(long)==32 (256-bit longs), the product easily fits in the
3029 dynamic range of a double. So the leading 50 (or so) bits of the double
3030 product are correct.
3031
3032 We check these two ways against each other, and declare victory if they're
3033 approximately the same. Else, because the native long product is the only
3034 one that can lose catastrophic amounts of information, it's the native long
3035 product that must have overflowed.
3036
3037 */
3038
3039 long longprod = (long)((unsigned long)a * b);
3040 double doubleprod = (double)a * (double)b;
3041 double doubled_longprod = (double)longprod;
3042
3043 if (doubled_longprod == doubleprod) {
3044 return 0;
3045 }
3046
3047 const double diff = doubled_longprod - doubleprod;
3048 const double absdiff = diff >= 0.0 ? diff : -diff;
3049 const double absprod = doubleprod >= 0.0 ? doubleprod : -doubleprod;
3050
3051 if (32.0 * absdiff <= absprod) {
3052 return 0;
3053 }
3054
3055 return 1;
3056}
3057
3058/*[clinic input]
3059math.prod
3060
3061 iterable: object
3062 /
3063 *
3064 start: object(c_default="NULL") = 1
3065
3066Calculate the product of all the elements in the input iterable.
3067
3068The default start value for the product is 1.
3069
3070When the iterable is empty, return the start value. This function is
3071intended specifically for use with numeric values and may reject
3072non-numeric types.
3073[clinic start generated code]*/
3074
3075static PyObject *
3076math_prod_impl(PyObject *module, PyObject *iterable, PyObject *start)
3077/*[clinic end generated code: output=36153bedac74a198 input=4c5ab0682782ed54]*/
3078{
3079 PyObject *result = start;
3080 PyObject *temp, *item, *iter;
3081
3082 iter = PyObject_GetIter(iterable);
3083 if (iter == NULL) {
3084 return NULL;
3085 }
3086
3087 if (result == NULL) {
3088 result = _PyLong_GetOne();
3089 }
3090 Py_INCREF(result);
3091#ifndef SLOW_PROD
3092 /* Fast paths for integers keeping temporary products in C.
3093 * Assumes all inputs are the same type.
3094 * If the assumption fails, default to use PyObjects instead.
3095 */
3096 if (PyLong_CheckExact(result)) {
3097 int overflow;
3098 long i_result = PyLong_AsLongAndOverflow(result, &overflow);
3099 /* If this already overflowed, don't even enter the loop. */
3100 if (overflow == 0) {
3101 Py_DECREF(result);
3102 result = NULL;
3103 }
3104 /* Loop over all the items in the iterable until we finish, we overflow
3105 * or we found a non integer element */
3106 while (result == NULL) {
3107 item = PyIter_Next(iter);
3108 if (item == NULL) {
3109 Py_DECREF(iter);
3110 if (PyErr_Occurred()) {
3111 return NULL;
3112 }
3113 return PyLong_FromLong(i_result);
3114 }
3115 if (PyLong_CheckExact(item)) {
3116 long b = PyLong_AsLongAndOverflow(item, &overflow);
3117 if (overflow == 0 && !_check_long_mult_overflow(i_result, b)) {
3118 long x = i_result * b;
3119 i_result = x;
3120 Py_DECREF(item);
3121 continue;
3122 }
3123 }
3124 /* Either overflowed or is not an int.
3125 * Restore real objects and process normally */
3126 result = PyLong_FromLong(i_result);
3127 if (result == NULL) {
3128 Py_DECREF(item);
3129 Py_DECREF(iter);
3130 return NULL;
3131 }
3132 temp = PyNumber_Multiply(result, item);
3133 Py_DECREF(result);
3134 Py_DECREF(item);
3135 result = temp;
3136 if (result == NULL) {
3137 Py_DECREF(iter);
3138 return NULL;
3139 }
3140 }
3141 }
3142
3143 /* Fast paths for floats keeping temporary products in C.
3144 * Assumes all inputs are the same type.
3145 * If the assumption fails, default to use PyObjects instead.
3146 */
3147 if (PyFloat_CheckExact(result)) {
3148 double f_result = PyFloat_AS_DOUBLE(result);
3149 Py_DECREF(result);
3150 result = NULL;
3151 while(result == NULL) {
3152 item = PyIter_Next(iter);
3153 if (item == NULL) {
3154 Py_DECREF(iter);
3155 if (PyErr_Occurred()) {
3156 return NULL;
3157 }
3158 return PyFloat_FromDouble(f_result);
3159 }
3160 if (PyFloat_CheckExact(item)) {
3161 f_result *= PyFloat_AS_DOUBLE(item);
3162 Py_DECREF(item);
3163 continue;
3164 }
3165 if (PyLong_CheckExact(item)) {
3166 long value;
3167 int overflow;
3168 value = PyLong_AsLongAndOverflow(item, &overflow);
3169 if (!overflow) {
3170 f_result *= (double)value;
3171 Py_DECREF(item);
3172 continue;
3173 }
3174 }
3175 result = PyFloat_FromDouble(f_result);
3176 if (result == NULL) {
3177 Py_DECREF(item);
3178 Py_DECREF(iter);
3179 return NULL;
3180 }
3181 temp = PyNumber_Multiply(result, item);
3182 Py_DECREF(result);
3183 Py_DECREF(item);
3184 result = temp;
3185 if (result == NULL) {
3186 Py_DECREF(iter);
3187 return NULL;
3188 }
3189 }
3190 }
3191#endif
3192 /* Consume rest of the iterable (if any) that could not be handled
3193 * by specialized functions above.*/
3194 for(;;) {
3195 item = PyIter_Next(iter);
3196 if (item == NULL) {
3197 /* error, or end-of-sequence */
3198 if (PyErr_Occurred()) {
3199 Py_DECREF(result);
3200 result = NULL;
3201 }
3202 break;
3203 }
3204 temp = PyNumber_Multiply(result, item);
3205 Py_DECREF(result);
3206 Py_DECREF(item);
3207 result = temp;
3208 if (result == NULL)
3209 break;
3210 }
3211 Py_DECREF(iter);
3212 return result;
3213}
3214
3215
3216/*[clinic input]
3217math.perm
3218
3219 n: object
3220 k: object = None
3221 /
3222
3223Number of ways to choose k items from n items without repetition and with order.
3224
3225Evaluates to n! / (n - k)! when k <= n and evaluates
3226to zero when k > n.
3227
3228If k is not specified or is None, then k defaults to n
3229and the function returns n!.
3230
3231Raises TypeError if either of the arguments are not integers.
3232Raises ValueError if either of the arguments are negative.
3233[clinic start generated code]*/
3234
3235static PyObject *
3236math_perm_impl(PyObject *module, PyObject *n, PyObject *k)
3237/*[clinic end generated code: output=e021a25469653e23 input=5311c5a00f359b53]*/
3238{
3239 PyObject *result = NULL, *factor = NULL;
3240 int overflow, cmp;
3241 long long i, factors;
3242
3243 if (k == Py_None) {
3244 return math_factorial(module, n);
3245 }
3246 n = PyNumber_Index(n);
3247 if (n == NULL) {
3248 return NULL;
3249 }
3250 k = PyNumber_Index(k);
3251 if (k == NULL) {
3252 Py_DECREF(n);
3253 return NULL;
3254 }
3255
3256 if (Py_SIZE(n) < 0) {
3257 PyErr_SetString(PyExc_ValueError,
3258 "n must be a non-negative integer");
3259 goto error;
3260 }
3261 if (Py_SIZE(k) < 0) {
3262 PyErr_SetString(PyExc_ValueError,
3263 "k must be a non-negative integer");
3264 goto error;
3265 }
3266
3267 cmp = PyObject_RichCompareBool(n, k, Py_LT);
3268 if (cmp != 0) {
3269 if (cmp > 0) {
3270 result = PyLong_FromLong(0);
3271 goto done;
3272 }
3273 goto error;
3274 }
3275
3276 factors = PyLong_AsLongLongAndOverflow(k, &overflow);
3277 if (overflow > 0) {
3278 PyErr_Format(PyExc_OverflowError,
3279 "k must not exceed %lld",
3280 LLONG_MAX);
3281 goto error;
3282 }
3283 else if (factors == -1) {
3284 /* k is nonnegative, so a return value of -1 can only indicate error */
3285 goto error;
3286 }
3287
3288 if (factors == 0) {
3289 result = PyLong_FromLong(1);
3290 goto done;
3291 }
3292
3293 result = n;
3294 Py_INCREF(result);
3295 if (factors == 1) {
3296 goto done;
3297 }
3298
3299 factor = Py_NewRef(n);
3300 PyObject *one = _PyLong_GetOne(); // borrowed ref
3301 for (i = 1; i < factors; ++i) {
3302 Py_SETREF(factor, PyNumber_Subtract(factor, one));
3303 if (factor == NULL) {
3304 goto error;
3305 }
3306 Py_SETREF(result, PyNumber_Multiply(result, factor));
3307 if (result == NULL) {
3308 goto error;
3309 }
3310 }
3311 Py_DECREF(factor);
3312
3313done:
3314 Py_DECREF(n);
3315 Py_DECREF(k);
3316 return result;
3317
3318error:
3319 Py_XDECREF(factor);
3320 Py_XDECREF(result);
3321 Py_DECREF(n);
3322 Py_DECREF(k);
3323 return NULL;
3324}
3325
3326
3327/*[clinic input]
3328math.comb
3329
3330 n: object
3331 k: object
3332 /
3333
3334Number of ways to choose k items from n items without repetition and without order.
3335
3336Evaluates to n! / (k! * (n - k)!) when k <= n and evaluates
3337to zero when k > n.
3338
3339Also called the binomial coefficient because it is equivalent
3340to the coefficient of k-th term in polynomial expansion of the
3341expression (1 + x)**n.
3342
3343Raises TypeError if either of the arguments are not integers.
3344Raises ValueError if either of the arguments are negative.
3345
3346[clinic start generated code]*/
3347
3348static PyObject *
3349math_comb_impl(PyObject *module, PyObject *n, PyObject *k)
3350/*[clinic end generated code: output=bd2cec8d854f3493 input=9a05315af2518709]*/
3351{
3352 PyObject *result = NULL, *factor = NULL, *temp;
3353 int overflow, cmp;
3354 long long i, factors;
3355
3356 n = PyNumber_Index(n);
3357 if (n == NULL) {
3358 return NULL;
3359 }
3360 k = PyNumber_Index(k);
3361 if (k == NULL) {
3362 Py_DECREF(n);
3363 return NULL;
3364 }
3365
3366 if (Py_SIZE(n) < 0) {
3367 PyErr_SetString(PyExc_ValueError,
3368 "n must be a non-negative integer");
3369 goto error;
3370 }
3371 if (Py_SIZE(k) < 0) {
3372 PyErr_SetString(PyExc_ValueError,
3373 "k must be a non-negative integer");
3374 goto error;
3375 }
3376
3377 /* k = min(k, n - k) */
3378 temp = PyNumber_Subtract(n, k);
3379 if (temp == NULL) {
3380 goto error;
3381 }
3382 if (Py_SIZE(temp) < 0) {
3383 Py_DECREF(temp);
3384 result = PyLong_FromLong(0);
3385 goto done;
3386 }
3387 cmp = PyObject_RichCompareBool(temp, k, Py_LT);
3388 if (cmp > 0) {
3389 Py_SETREF(k, temp);
3390 }
3391 else {
3392 Py_DECREF(temp);
3393 if (cmp < 0) {
3394 goto error;
3395 }
3396 }
3397
3398 factors = PyLong_AsLongLongAndOverflow(k, &overflow);
3399 if (overflow > 0) {
3400 PyErr_Format(PyExc_OverflowError,
3401 "min(n - k, k) must not exceed %lld",
3402 LLONG_MAX);
3403 goto error;
3404 }
3405 if (factors == -1) {
3406 /* k is nonnegative, so a return value of -1 can only indicate error */
3407 goto error;
3408 }
3409
3410 if (factors == 0) {
3411 result = PyLong_FromLong(1);
3412 goto done;
3413 }
3414
3415 result = n;
3416 Py_INCREF(result);
3417 if (factors == 1) {
3418 goto done;
3419 }
3420
3421 factor = Py_NewRef(n);
3422 PyObject *one = _PyLong_GetOne(); // borrowed ref
3423 for (i = 1; i < factors; ++i) {
3424 Py_SETREF(factor, PyNumber_Subtract(factor, one));
3425 if (factor == NULL) {
3426 goto error;
3427 }
3428 Py_SETREF(result, PyNumber_Multiply(result, factor));
3429 if (result == NULL) {
3430 goto error;
3431 }
3432
3433 temp = PyLong_FromUnsignedLongLong((unsigned long long)i + 1);
3434 if (temp == NULL) {
3435 goto error;
3436 }
3437 Py_SETREF(result, PyNumber_FloorDivide(result, temp));
3438 Py_DECREF(temp);
3439 if (result == NULL) {
3440 goto error;
3441 }
3442 }
3443 Py_DECREF(factor);
3444
3445done:
3446 Py_DECREF(n);
3447 Py_DECREF(k);
3448 return result;
3449
3450error:
3451 Py_XDECREF(factor);
3452 Py_XDECREF(result);
3453 Py_DECREF(n);
3454 Py_DECREF(k);
3455 return NULL;
3456}
3457
3458
3459/*[clinic input]
3460math.nextafter
3461
3462 x: double
3463 y: double
3464 /
3465
3466Return the next floating-point value after x towards y.
3467[clinic start generated code]*/
3468
3469static PyObject *
3470math_nextafter_impl(PyObject *module, double x, double y)
3471/*[clinic end generated code: output=750c8266c1c540ce input=02b2d50cd1d9f9b6]*/
3472{
3473#if defined(_AIX)
3474 if (x == y) {
3475 /* On AIX 7.1, libm nextafter(-0.0, +0.0) returns -0.0.
3476 Bug fixed in bos.adt.libm 7.2.2.0 by APAR IV95512. */
3477 return PyFloat_FromDouble(y);
3478 }
3479 if (Py_IS_NAN(x)) {
3480 return PyFloat_FromDouble(x);
3481 }
3482 if (Py_IS_NAN(y)) {
3483 return PyFloat_FromDouble(y);
3484 }
3485#endif
3486 return PyFloat_FromDouble(nextafter(x, y));
3487}
3488
3489
3490/*[clinic input]
3491math.ulp -> double
3492
3493 x: double
3494 /
3495
3496Return the value of the least significant bit of the float x.
3497[clinic start generated code]*/
3498
3499static double
3500math_ulp_impl(PyObject *module, double x)
3501/*[clinic end generated code: output=f5207867a9384dd4 input=31f9bfbbe373fcaa]*/
3502{
3503 if (Py_IS_NAN(x)) {
3504 return x;
3505 }
3506 x = fabs(x);
3507 if (Py_IS_INFINITY(x)) {
3508 return x;
3509 }
3510 double inf = m_inf();
3511 double x2 = nextafter(x, inf);
3512 if (Py_IS_INFINITY(x2)) {
3513 /* special case: x is the largest positive representable float */
3514 x2 = nextafter(x, -inf);
3515 return x - x2;
3516 }
3517 return x2 - x;
3518}
3519
3520static int
3521math_exec(PyObject *module)
3522{
3523 if (PyModule_AddObject(module, "pi", PyFloat_FromDouble(Py_MATH_PI)) < 0) {
3524 return -1;
3525 }
3526 if (PyModule_AddObject(module, "e", PyFloat_FromDouble(Py_MATH_E)) < 0) {
3527 return -1;
3528 }
3529 // 2pi
3530 if (PyModule_AddObject(module, "tau", PyFloat_FromDouble(Py_MATH_TAU)) < 0) {
3531 return -1;
3532 }
3533 if (PyModule_AddObject(module, "inf", PyFloat_FromDouble(m_inf())) < 0) {
3534 return -1;
3535 }
3536#if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN)
3537 if (PyModule_AddObject(module, "nan", PyFloat_FromDouble(m_nan())) < 0) {
3538 return -1;
3539 }
3540#endif
3541 return 0;
3542}
3543
3544static PyMethodDef math_methods[] = {
3545 {"acos", math_acos, METH_O, math_acos_doc},
3546 {"acosh", math_acosh, METH_O, math_acosh_doc},
3547 {"asin", math_asin, METH_O, math_asin_doc},
3548 {"asinh", math_asinh, METH_O, math_asinh_doc},
3549 {"atan", math_atan, METH_O, math_atan_doc},
3550 {"atan2", (PyCFunction)(void(*)(void))math_atan2, METH_FASTCALL, math_atan2_doc},
3551 {"atanh", math_atanh, METH_O, math_atanh_doc},
3552 MATH_CEIL_METHODDEF
3553 {"copysign", (PyCFunction)(void(*)(void))math_copysign, METH_FASTCALL, math_copysign_doc},
3554 {"cos", math_cos, METH_O, math_cos_doc},
3555 {"cosh", math_cosh, METH_O, math_cosh_doc},
3556 MATH_DEGREES_METHODDEF
3557 MATH_DIST_METHODDEF
3558 {"erf", math_erf, METH_O, math_erf_doc},
3559 {"erfc", math_erfc, METH_O, math_erfc_doc},
3560 {"exp", math_exp, METH_O, math_exp_doc},
3561 {"expm1", math_expm1, METH_O, math_expm1_doc},
3562 {"fabs", math_fabs, METH_O, math_fabs_doc},
3563 MATH_FACTORIAL_METHODDEF
3564 MATH_FLOOR_METHODDEF
3565 MATH_FMOD_METHODDEF
3566 MATH_FREXP_METHODDEF
3567 MATH_FSUM_METHODDEF
3568 {"gamma", math_gamma, METH_O, math_gamma_doc},
3569 {"gcd", (PyCFunction)(void(*)(void))math_gcd, METH_FASTCALL, math_gcd_doc},
3570 {"hypot", (PyCFunction)(void(*)(void))math_hypot, METH_FASTCALL, math_hypot_doc},
3571 MATH_ISCLOSE_METHODDEF
3572 MATH_ISFINITE_METHODDEF
3573 MATH_ISINF_METHODDEF
3574 MATH_ISNAN_METHODDEF
3575 MATH_ISQRT_METHODDEF
3576 {"lcm", (PyCFunction)(void(*)(void))math_lcm, METH_FASTCALL, math_lcm_doc},
3577 MATH_LDEXP_METHODDEF
3578 {"lgamma", math_lgamma, METH_O, math_lgamma_doc},
3579 MATH_LOG_METHODDEF
3580 {"log1p", math_log1p, METH_O, math_log1p_doc},
3581 MATH_LOG10_METHODDEF
3582 MATH_LOG2_METHODDEF
3583 MATH_MODF_METHODDEF
3584 MATH_POW_METHODDEF
3585 MATH_RADIANS_METHODDEF
3586 {"remainder", (PyCFunction)(void(*)(void))math_remainder, METH_FASTCALL, math_remainder_doc},
3587 {"sin", math_sin, METH_O, math_sin_doc},
3588 {"sinh", math_sinh, METH_O, math_sinh_doc},
3589 {"sqrt", math_sqrt, METH_O, math_sqrt_doc},
3590 {"tan", math_tan, METH_O, math_tan_doc},
3591 {"tanh", math_tanh, METH_O, math_tanh_doc},
3592 MATH_TRUNC_METHODDEF
3593 MATH_PROD_METHODDEF
3594 MATH_PERM_METHODDEF
3595 MATH_COMB_METHODDEF
3596 MATH_NEXTAFTER_METHODDEF
3597 MATH_ULP_METHODDEF
3598 {NULL, NULL} /* sentinel */
3599};
3600
3601static PyModuleDef_Slot math_slots[] = {
3602 {Py_mod_exec, math_exec},
3603 {0, NULL}
3604};
3605
3606PyDoc_STRVAR(module_doc,
3607"This module provides access to the mathematical functions\n"
3608"defined by the C standard.");
3609
3610static struct PyModuleDef mathmodule = {
3611 PyModuleDef_HEAD_INIT,
3612 .m_name = "math",
3613 .m_doc = module_doc,
3614 .m_size = 0,
3615 .m_methods = math_methods,
3616 .m_slots = math_slots,
3617};
3618
3619PyMODINIT_FUNC
3620PyInit_math(void)
3621{
3622 return PyModuleDef_Init(&mathmodule);
3623}
3624