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