blob: 552cb78c88d19f02ebc23ff7e3e3741c36b14508 [file] [log] [blame]
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001/* Math module -- standard C math library functions, pi and e */
2
Christian Heimes53876d92008-04-19 00:31:39 +00003/* Here are some comments from Tim Peters, extracted from the
4 discussion attached to http://bugs.python.org/issue1640. They
5 describe the general aims of the math module with respect to
6 special values, IEEE-754 floating-point exceptions, and Python
7 exceptions.
8
9These are the "spirit of 754" rules:
10
111. If the mathematical result is a real number, but of magnitude too
12large to approximate by a machine float, overflow is signaled and the
13result is an infinity (with the appropriate sign).
14
152. If the mathematical result is a real number, but of magnitude too
16small to approximate by a machine float, underflow is signaled and the
17result is a zero (with the appropriate sign).
18
193. At a singularity (a value x such that the limit of f(y) as y
20approaches x exists and is an infinity), "divide by zero" is signaled
21and the result is an infinity (with the appropriate sign). This is
22complicated a little by that the left-side and right-side limits may
23not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0
24from the positive or negative directions. In that specific case, the
25sign of the zero determines the result of 1/0.
26
274. At a point where a function has no defined result in the extended
28reals (i.e., the reals plus an infinity or two), invalid operation is
29signaled and a NaN is returned.
30
31And these are what Python has historically /tried/ to do (but not
32always successfully, as platform libm behavior varies a lot):
33
34For #1, raise OverflowError.
35
36For #2, return a zero (with the appropriate sign if that happens by
37accident ;-)).
38
39For #3 and #4, raise ValueError. It may have made sense to raise
40Python's ZeroDivisionError in #3, but historically that's only been
41raised for division by zero and mod by zero.
42
43*/
44
45/*
46 In general, on an IEEE-754 platform the aim is to follow the C99
47 standard, including Annex 'F', whenever possible. Where the
48 standard recommends raising the 'divide-by-zero' or 'invalid'
49 floating-point exceptions, Python should raise a ValueError. Where
50 the standard recommends raising 'overflow', Python should raise an
51 OverflowError. In all other circumstances a value should be
52 returned.
53 */
54
Barry Warsaw8b43b191996-12-09 22:32:36 +000055#include "Python.h"
Mark Dickinson664b5112009-12-16 20:23:42 +000056#include "_math.h"
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000057
Serhiy Storchakac9ea9332017-01-19 18:13:09 +020058#include "clinic/mathmodule.c.h"
59
60/*[clinic input]
61module math
62[clinic start generated code]*/
63/*[clinic end generated code: output=da39a3ee5e6b4b0d input=76bc7002685dd942]*/
64
65
Mark Dickinson12c4bdb2009-09-28 19:21:11 +000066/*
67 sin(pi*x), giving accurate results for all finite x (especially x
68 integral or close to an integer). This is here for use in the
69 reflection formula for the gamma function. It conforms to IEEE
70 754-2008 for finite arguments, but not for infinities or nans.
71*/
Tim Petersa40c7932001-09-05 22:36:56 +000072
Mark Dickinson12c4bdb2009-09-28 19:21:11 +000073static const double pi = 3.141592653589793238462643383279502884197;
Mark Dickinson9c91eb82010-07-07 16:17:31 +000074static const double logpi = 1.144729885849400174143427351353058711647;
Louie Lu7a264642017-03-31 01:05:10 +080075#if !defined(HAVE_ERF) || !defined(HAVE_ERFC)
76static const double sqrtpi = 1.772453850905516027298167483341145182798;
77#endif /* !defined(HAVE_ERF) || !defined(HAVE_ERFC) */
Mark Dickinson12c4bdb2009-09-28 19:21:11 +000078
Raymond Hettingercfd735e2019-01-29 20:39:53 -080079
80/* Version of PyFloat_AsDouble() with in-line fast paths
81 for exact floats and integers. Gives a substantial
82 speed improvement for extracting float arguments.
83*/
84
85#define ASSIGN_DOUBLE(target_var, obj, error_label) \
86 if (PyFloat_CheckExact(obj)) { \
87 target_var = PyFloat_AS_DOUBLE(obj); \
88 } \
89 else if (PyLong_CheckExact(obj)) { \
90 target_var = PyLong_AsDouble(obj); \
91 if (target_var == -1.0 && PyErr_Occurred()) { \
92 goto error_label; \
93 } \
94 } \
95 else { \
96 target_var = PyFloat_AsDouble(obj); \
97 if (target_var == -1.0 && PyErr_Occurred()) { \
98 goto error_label; \
99 } \
100 }
101
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000102static double
Dima Pasechnikf57cd822019-02-26 06:36:11 +0000103m_sinpi(double x)
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000104{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000105 double y, r;
106 int n;
107 /* this function should only ever be called for finite arguments */
108 assert(Py_IS_FINITE(x));
109 y = fmod(fabs(x), 2.0);
110 n = (int)round(2.0*y);
111 assert(0 <= n && n <= 4);
112 switch (n) {
113 case 0:
114 r = sin(pi*y);
115 break;
116 case 1:
117 r = cos(pi*(y-0.5));
118 break;
119 case 2:
120 /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give
121 -0.0 instead of 0.0 when y == 1.0. */
122 r = sin(pi*(1.0-y));
123 break;
124 case 3:
125 r = -cos(pi*(y-1.5));
126 break;
127 case 4:
128 r = sin(pi*(y-2.0));
129 break;
130 default:
Barry Warsawb2e57942017-09-14 18:13:16 -0700131 Py_UNREACHABLE();
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000132 }
133 return copysign(1.0, x)*r;
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000134}
135
136/* Implementation of the real gamma function. In extensive but non-exhaustive
137 random tests, this function proved accurate to within <= 10 ulps across the
138 entire float domain. Note that accuracy may depend on the quality of the
139 system math functions, the pow function in particular. Special cases
140 follow C99 annex F. The parameters and method are tailored to platforms
141 whose double format is the IEEE 754 binary64 format.
142
143 Method: for x > 0.0 we use the Lanczos approximation with parameters N=13
144 and g=6.024680040776729583740234375; these parameters are amongst those
145 used by the Boost library. Following Boost (again), we re-express the
146 Lanczos sum as a rational function, and compute it that way. The
147 coefficients below were computed independently using MPFR, and have been
148 double-checked against the coefficients in the Boost source code.
149
150 For x < 0.0 we use the reflection formula.
151
152 There's one minor tweak that deserves explanation: Lanczos' formula for
153 Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x
154 values, x+g-0.5 can be represented exactly. However, in cases where it
155 can't be represented exactly the small error in x+g-0.5 can be magnified
156 significantly by the pow and exp calls, especially for large x. A cheap
157 correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error
158 involved in the computation of x+g-0.5 (that is, e = computed value of
159 x+g-0.5 - exact value of x+g-0.5). Here's the proof:
160
161 Correction factor
162 -----------------
163 Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754
164 double, and e is tiny. Then:
165
166 pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e)
167 = pow(y, x-0.5)/exp(y) * C,
168
169 where the correction_factor C is given by
170
171 C = pow(1-e/y, x-0.5) * exp(e)
172
173 Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so:
174
175 C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y
176
177 But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and
178
179 pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y),
180
181 Note that for accuracy, when computing r*C it's better to do
182
183 r + e*g/y*r;
184
185 than
186
187 r * (1 + e*g/y);
188
189 since the addition in the latter throws away most of the bits of
190 information in e*g/y.
191*/
192
193#define LANCZOS_N 13
194static const double lanczos_g = 6.024680040776729583740234375;
195static const double lanczos_g_minus_half = 5.524680040776729583740234375;
196static const double lanczos_num_coeffs[LANCZOS_N] = {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000197 23531376880.410759688572007674451636754734846804940,
198 42919803642.649098768957899047001988850926355848959,
199 35711959237.355668049440185451547166705960488635843,
200 17921034426.037209699919755754458931112671403265390,
201 6039542586.3520280050642916443072979210699388420708,
202 1439720407.3117216736632230727949123939715485786772,
203 248874557.86205415651146038641322942321632125127801,
204 31426415.585400194380614231628318205362874684987640,
205 2876370.6289353724412254090516208496135991145378768,
206 186056.26539522349504029498971604569928220784236328,
207 8071.6720023658162106380029022722506138218516325024,
208 210.82427775157934587250973392071336271166969580291,
209 2.5066282746310002701649081771338373386264310793408
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000210};
211
212/* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */
213static const double lanczos_den_coeffs[LANCZOS_N] = {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000214 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
215 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000216
217/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */
218#define NGAMMA_INTEGRAL 23
219static const double gamma_integral[NGAMMA_INTEGRAL] = {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000220 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
221 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
222 1307674368000.0, 20922789888000.0, 355687428096000.0,
223 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
224 51090942171709440000.0, 1124000727777607680000.0,
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000225};
226
227/* Lanczos' sum L_g(x), for positive x */
228
229static double
230lanczos_sum(double x)
231{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000232 double num = 0.0, den = 0.0;
233 int i;
234 assert(x > 0.0);
235 /* evaluate the rational function lanczos_sum(x). For large
236 x, the obvious algorithm risks overflow, so we instead
237 rescale the denominator and numerator of the rational
238 function by x**(1-LANCZOS_N) and treat this as a
239 rational function in 1/x. This also reduces the error for
240 larger x values. The choice of cutoff point (5.0 below) is
241 somewhat arbitrary; in tests, smaller cutoff values than
242 this resulted in lower accuracy. */
243 if (x < 5.0) {
244 for (i = LANCZOS_N; --i >= 0; ) {
245 num = num * x + lanczos_num_coeffs[i];
246 den = den * x + lanczos_den_coeffs[i];
247 }
248 }
249 else {
250 for (i = 0; i < LANCZOS_N; i++) {
251 num = num / x + lanczos_num_coeffs[i];
252 den = den / x + lanczos_den_coeffs[i];
253 }
254 }
255 return num/den;
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000256}
257
Mark Dickinsona5d0c7c2015-01-11 11:55:29 +0000258/* Constant for +infinity, generated in the same way as float('inf'). */
259
260static double
261m_inf(void)
262{
263#ifndef PY_NO_SHORT_FLOAT_REPR
264 return _Py_dg_infinity(0);
265#else
266 return Py_HUGE_VAL;
267#endif
268}
269
270/* Constant nan value, generated in the same way as float('nan'). */
271/* We don't currently assume that Py_NAN is defined everywhere. */
272
273#if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN)
274
275static double
276m_nan(void)
277{
278#ifndef PY_NO_SHORT_FLOAT_REPR
279 return _Py_dg_stdnan(0);
280#else
281 return Py_NAN;
282#endif
283}
284
285#endif
286
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000287static double
288m_tgamma(double x)
289{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000290 double absx, r, y, z, sqrtpow;
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000291
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000292 /* special cases */
293 if (!Py_IS_FINITE(x)) {
294 if (Py_IS_NAN(x) || x > 0.0)
295 return x; /* tgamma(nan) = nan, tgamma(inf) = inf */
296 else {
297 errno = EDOM;
298 return Py_NAN; /* tgamma(-inf) = nan, invalid */
299 }
300 }
301 if (x == 0.0) {
302 errno = EDOM;
Mark Dickinson50203a62011-09-25 15:26:43 +0100303 /* tgamma(+-0.0) = +-inf, divide-by-zero */
304 return copysign(Py_HUGE_VAL, x);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000305 }
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000306
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000307 /* integer arguments */
308 if (x == floor(x)) {
309 if (x < 0.0) {
310 errno = EDOM; /* tgamma(n) = nan, invalid for */
311 return Py_NAN; /* negative integers n */
312 }
313 if (x <= NGAMMA_INTEGRAL)
314 return gamma_integral[(int)x - 1];
315 }
316 absx = fabs(x);
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000317
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000318 /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */
319 if (absx < 1e-20) {
320 r = 1.0/x;
321 if (Py_IS_INFINITY(r))
322 errno = ERANGE;
323 return r;
324 }
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000325
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000326 /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for
327 x > 200, and underflows to +-0.0 for x < -200, not a negative
328 integer. */
329 if (absx > 200.0) {
330 if (x < 0.0) {
Dima Pasechnikf57cd822019-02-26 06:36:11 +0000331 return 0.0/m_sinpi(x);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000332 }
333 else {
334 errno = ERANGE;
335 return Py_HUGE_VAL;
336 }
337 }
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000338
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000339 y = absx + lanczos_g_minus_half;
340 /* compute error in sum */
341 if (absx > lanczos_g_minus_half) {
342 /* note: the correction can be foiled by an optimizing
343 compiler that (incorrectly) thinks that an expression like
344 a + b - a - b can be optimized to 0.0. This shouldn't
345 happen in a standards-conforming compiler. */
346 double q = y - absx;
347 z = q - lanczos_g_minus_half;
348 }
349 else {
350 double q = y - lanczos_g_minus_half;
351 z = q - absx;
352 }
353 z = z * lanczos_g / y;
354 if (x < 0.0) {
Dima Pasechnikf57cd822019-02-26 06:36:11 +0000355 r = -pi / m_sinpi(absx) / absx * exp(y) / lanczos_sum(absx);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000356 r -= z * r;
357 if (absx < 140.0) {
358 r /= pow(y, absx - 0.5);
359 }
360 else {
361 sqrtpow = pow(y, absx / 2.0 - 0.25);
362 r /= sqrtpow;
363 r /= sqrtpow;
364 }
365 }
366 else {
367 r = lanczos_sum(absx) / exp(y);
368 r += z * r;
369 if (absx < 140.0) {
370 r *= pow(y, absx - 0.5);
371 }
372 else {
373 sqrtpow = pow(y, absx / 2.0 - 0.25);
374 r *= sqrtpow;
375 r *= sqrtpow;
376 }
377 }
378 if (Py_IS_INFINITY(r))
379 errno = ERANGE;
380 return r;
Guido van Rossum8832b621991-12-16 15:44:24 +0000381}
382
Christian Heimes53876d92008-04-19 00:31:39 +0000383/*
Mark Dickinson05d2e082009-12-11 20:17:17 +0000384 lgamma: natural log of the absolute value of the Gamma function.
385 For large arguments, Lanczos' formula works extremely well here.
386*/
387
388static double
389m_lgamma(double x)
390{
Serhiy Storchaka97553fd2017-03-11 23:37:16 +0200391 double r;
Serhiy Storchaka97553fd2017-03-11 23:37:16 +0200392 double absx;
Mark Dickinson05d2e082009-12-11 20:17:17 +0000393
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000394 /* special cases */
395 if (!Py_IS_FINITE(x)) {
396 if (Py_IS_NAN(x))
397 return x; /* lgamma(nan) = nan */
398 else
399 return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */
400 }
Mark Dickinson05d2e082009-12-11 20:17:17 +0000401
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000402 /* integer arguments */
403 if (x == floor(x) && x <= 2.0) {
404 if (x <= 0.0) {
405 errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */
406 return Py_HUGE_VAL; /* integers n <= 0 */
407 }
408 else {
409 return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */
410 }
411 }
Mark Dickinson05d2e082009-12-11 20:17:17 +0000412
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000413 absx = fabs(x);
414 /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */
415 if (absx < 1e-20)
416 return -log(absx);
Mark Dickinson05d2e082009-12-11 20:17:17 +0000417
Mark Dickinson9c91eb82010-07-07 16:17:31 +0000418 /* Lanczos' formula. We could save a fraction of a ulp in accuracy by
419 having a second set of numerator coefficients for lanczos_sum that
420 absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g
421 subtraction below; it's probably not worth it. */
422 r = log(lanczos_sum(absx)) - lanczos_g;
423 r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1);
424 if (x < 0.0)
425 /* Use reflection formula to get value for negative x. */
Dima Pasechnikf57cd822019-02-26 06:36:11 +0000426 r = logpi - log(fabs(m_sinpi(absx))) - log(absx) - r;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000427 if (Py_IS_INFINITY(r))
428 errno = ERANGE;
429 return r;
Mark Dickinson05d2e082009-12-11 20:17:17 +0000430}
431
Serhiy Storchaka97553fd2017-03-11 23:37:16 +0200432#if !defined(HAVE_ERF) || !defined(HAVE_ERFC)
433
Mark Dickinson45f992a2009-12-19 11:20:49 +0000434/*
435 Implementations of the error function erf(x) and the complementary error
436 function erfc(x).
437
Brett Cannon45adb312016-01-15 09:38:24 -0800438 Method: we use a series approximation for erf for small x, and a continued
439 fraction approximation for erfc(x) for larger x;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000440 combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x),
441 this gives us erf(x) and erfc(x) for all x.
442
443 The series expansion used is:
444
445 erf(x) = x*exp(-x*x)/sqrt(pi) * [
446 2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...]
447
448 The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k).
449 This series converges well for smallish x, but slowly for larger x.
450
451 The continued fraction expansion used is:
452
453 erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - )
454 3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...]
455
456 after the first term, the general term has the form:
457
458 k*(k-0.5)/(2*k+0.5 + x**2 - ...).
459
460 This expansion converges fast for larger x, but convergence becomes
461 infinitely slow as x approaches 0.0. The (somewhat naive) continued
462 fraction evaluation algorithm used below also risks overflow for large x;
463 but for large x, erfc(x) == 0.0 to within machine precision. (For
464 example, erfc(30.0) is approximately 2.56e-393).
465
466 Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and
467 continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) <
468 ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the
469 numbers of terms to use for the relevant expansions. */
470
471#define ERF_SERIES_CUTOFF 1.5
472#define ERF_SERIES_TERMS 25
473#define ERFC_CONTFRAC_CUTOFF 30.0
474#define ERFC_CONTFRAC_TERMS 50
475
476/*
477 Error function, via power series.
478
479 Given a finite float x, return an approximation to erf(x).
480 Converges reasonably fast for small x.
481*/
482
483static double
484m_erf_series(double x)
485{
Mark Dickinsonbcdf9da2010-06-13 10:52:38 +0000486 double x2, acc, fk, result;
487 int i, saved_errno;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000488
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000489 x2 = x * x;
490 acc = 0.0;
491 fk = (double)ERF_SERIES_TERMS + 0.5;
492 for (i = 0; i < ERF_SERIES_TERMS; i++) {
493 acc = 2.0 + x2 * acc / fk;
494 fk -= 1.0;
495 }
Mark Dickinsonbcdf9da2010-06-13 10:52:38 +0000496 /* Make sure the exp call doesn't affect errno;
497 see m_erfc_contfrac for more. */
498 saved_errno = errno;
499 result = acc * x * exp(-x2) / sqrtpi;
500 errno = saved_errno;
501 return result;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000502}
503
504/*
505 Complementary error function, via continued fraction expansion.
506
507 Given a positive float x, return an approximation to erfc(x). Converges
508 reasonably fast for x large (say, x > 2.0), and should be safe from
509 overflow if x and nterms are not too large. On an IEEE 754 machine, with x
510 <= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller
511 than the smallest representable nonzero float. */
512
513static double
514m_erfc_contfrac(double x)
515{
Mark Dickinsonbcdf9da2010-06-13 10:52:38 +0000516 double x2, a, da, p, p_last, q, q_last, b, result;
517 int i, saved_errno;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000518
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000519 if (x >= ERFC_CONTFRAC_CUTOFF)
520 return 0.0;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000521
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000522 x2 = x*x;
523 a = 0.0;
524 da = 0.5;
525 p = 1.0; p_last = 0.0;
526 q = da + x2; q_last = 1.0;
527 for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) {
528 double temp;
529 a += da;
530 da += 2.0;
531 b = da + x2;
532 temp = p; p = b*p - a*p_last; p_last = temp;
533 temp = q; q = b*q - a*q_last; q_last = temp;
534 }
Mark Dickinsonbcdf9da2010-06-13 10:52:38 +0000535 /* Issue #8986: On some platforms, exp sets errno on underflow to zero;
536 save the current errno value so that we can restore it later. */
537 saved_errno = errno;
538 result = p / q * x * exp(-x2) / sqrtpi;
539 errno = saved_errno;
540 return result;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000541}
542
Serhiy Storchaka97553fd2017-03-11 23:37:16 +0200543#endif /* !defined(HAVE_ERF) || !defined(HAVE_ERFC) */
544
Mark Dickinson45f992a2009-12-19 11:20:49 +0000545/* Error function erf(x), for general x */
546
547static double
548m_erf(double x)
549{
Serhiy Storchaka97553fd2017-03-11 23:37:16 +0200550#ifdef HAVE_ERF
551 return erf(x);
552#else
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000553 double absx, cf;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000554
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000555 if (Py_IS_NAN(x))
556 return x;
557 absx = fabs(x);
558 if (absx < ERF_SERIES_CUTOFF)
559 return m_erf_series(x);
560 else {
561 cf = m_erfc_contfrac(absx);
562 return x > 0.0 ? 1.0 - cf : cf - 1.0;
563 }
Serhiy Storchaka97553fd2017-03-11 23:37:16 +0200564#endif
Mark Dickinson45f992a2009-12-19 11:20:49 +0000565}
566
567/* Complementary error function erfc(x), for general x. */
568
569static double
570m_erfc(double x)
571{
Serhiy Storchaka97553fd2017-03-11 23:37:16 +0200572#ifdef HAVE_ERFC
573 return erfc(x);
574#else
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000575 double absx, cf;
Mark Dickinson45f992a2009-12-19 11:20:49 +0000576
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000577 if (Py_IS_NAN(x))
578 return x;
579 absx = fabs(x);
580 if (absx < ERF_SERIES_CUTOFF)
581 return 1.0 - m_erf_series(x);
582 else {
583 cf = m_erfc_contfrac(absx);
584 return x > 0.0 ? cf : 2.0 - cf;
585 }
Serhiy Storchaka97553fd2017-03-11 23:37:16 +0200586#endif
Mark Dickinson45f992a2009-12-19 11:20:49 +0000587}
Mark Dickinson05d2e082009-12-11 20:17:17 +0000588
589/*
Christian Heimese57950f2008-04-21 13:08:03 +0000590 wrapper for atan2 that deals directly with special cases before
591 delegating to the platform libm for the remaining cases. This
592 is necessary to get consistent behaviour across platforms.
593 Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
594 always follow C99.
595*/
596
597static double
598m_atan2(double y, double x)
599{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000600 if (Py_IS_NAN(x) || Py_IS_NAN(y))
601 return Py_NAN;
602 if (Py_IS_INFINITY(y)) {
603 if (Py_IS_INFINITY(x)) {
604 if (copysign(1., x) == 1.)
605 /* atan2(+-inf, +inf) == +-pi/4 */
606 return copysign(0.25*Py_MATH_PI, y);
607 else
608 /* atan2(+-inf, -inf) == +-pi*3/4 */
609 return copysign(0.75*Py_MATH_PI, y);
610 }
611 /* atan2(+-inf, x) == +-pi/2 for finite x */
612 return copysign(0.5*Py_MATH_PI, y);
613 }
614 if (Py_IS_INFINITY(x) || y == 0.) {
615 if (copysign(1., x) == 1.)
616 /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
617 return copysign(0., y);
618 else
619 /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
620 return copysign(Py_MATH_PI, y);
621 }
622 return atan2(y, x);
Christian Heimese57950f2008-04-21 13:08:03 +0000623}
624
Mark Dickinsona0ce3752017-04-05 18:34:27 +0100625
626/* IEEE 754-style remainder operation: x - n*y where n*y is the nearest
627 multiple of y to x, taking n even in the case of a tie. Assuming an IEEE 754
628 binary floating-point format, the result is always exact. */
629
630static double
631m_remainder(double x, double y)
632{
633 /* Deal with most common case first. */
634 if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) {
635 double absx, absy, c, m, r;
636
637 if (y == 0.0) {
638 return Py_NAN;
639 }
640
641 absx = fabs(x);
642 absy = fabs(y);
643 m = fmod(absx, absy);
644
645 /*
646 Warning: some subtlety here. What we *want* to know at this point is
647 whether the remainder m is less than, equal to, or greater than half
648 of absy. However, we can't do that comparison directly because we
649 can't be sure that 0.5*absy is representable (the mutiplication
650 might incur precision loss due to underflow). So instead we compare
651 m with the complement c = absy - m: m < 0.5*absy if and only if m <
652 c, and so on. The catch is that absy - m might also not be
653 representable, but it turns out that it doesn't matter:
654
655 - if m > 0.5*absy then absy - m is exactly representable, by
656 Sterbenz's lemma, so m > c
657 - if m == 0.5*absy then again absy - m is exactly representable
658 and m == c
659 - if m < 0.5*absy then either (i) 0.5*absy is exactly representable,
660 in which case 0.5*absy < absy - m, so 0.5*absy <= c and hence m <
661 c, or (ii) absy is tiny, either subnormal or in the lowest normal
662 binade. Then absy - m is exactly representable and again m < c.
663 */
664
665 c = absy - m;
666 if (m < c) {
667 r = m;
668 }
669 else if (m > c) {
670 r = -c;
671 }
672 else {
673 /*
674 Here absx is exactly halfway between two multiples of absy,
675 and we need to choose the even multiple. x now has the form
676
677 absx = n * absy + m
678
679 for some integer n (recalling that m = 0.5*absy at this point).
680 If n is even we want to return m; if n is odd, we need to
681 return -m.
682
683 So
684
685 0.5 * (absx - m) = (n/2) * absy
686
687 and now reducing modulo absy gives us:
688
689 | m, if n is odd
690 fmod(0.5 * (absx - m), absy) = |
691 | 0, if n is even
692
693 Now m - 2.0 * fmod(...) gives the desired result: m
694 if n is even, -m if m is odd.
695
696 Note that all steps in fmod(0.5 * (absx - m), absy)
697 will be computed exactly, with no rounding error
698 introduced.
699 */
700 assert(m == c);
701 r = m - 2.0 * fmod(0.5 * (absx - m), absy);
702 }
703 return copysign(1.0, x) * r;
704 }
705
706 /* Special values. */
707 if (Py_IS_NAN(x)) {
708 return x;
709 }
710 if (Py_IS_NAN(y)) {
711 return y;
712 }
713 if (Py_IS_INFINITY(x)) {
714 return Py_NAN;
715 }
716 assert(Py_IS_INFINITY(y));
717 return x;
718}
719
720
Christian Heimese57950f2008-04-21 13:08:03 +0000721/*
Mark Dickinsone675f082008-12-11 21:56:00 +0000722 Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
723 log(-ve), log(NaN). Here are wrappers for log and log10 that deal with
724 special values directly, passing positive non-special values through to
725 the system log/log10.
726 */
727
728static double
729m_log(double x)
730{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000731 if (Py_IS_FINITE(x)) {
732 if (x > 0.0)
733 return log(x);
734 errno = EDOM;
735 if (x == 0.0)
736 return -Py_HUGE_VAL; /* log(0) = -inf */
737 else
738 return Py_NAN; /* log(-ve) = nan */
739 }
740 else if (Py_IS_NAN(x))
741 return x; /* log(nan) = nan */
742 else if (x > 0.0)
743 return x; /* log(inf) = inf */
744 else {
745 errno = EDOM;
746 return Py_NAN; /* log(-inf) = nan */
747 }
Mark Dickinsone675f082008-12-11 21:56:00 +0000748}
749
Victor Stinnerfa0e3d52011-05-09 01:01:09 +0200750/*
751 log2: log to base 2.
752
753 Uses an algorithm that should:
Mark Dickinson83b8c0b2011-05-09 08:40:20 +0100754
Victor Stinnerfa0e3d52011-05-09 01:01:09 +0200755 (a) produce exact results for powers of 2, and
Mark Dickinson83b8c0b2011-05-09 08:40:20 +0100756 (b) give a monotonic log2 (for positive finite floats),
757 assuming that the system log is monotonic.
Victor Stinnerfa0e3d52011-05-09 01:01:09 +0200758*/
759
760static double
761m_log2(double x)
762{
763 if (!Py_IS_FINITE(x)) {
764 if (Py_IS_NAN(x))
765 return x; /* log2(nan) = nan */
766 else if (x > 0.0)
767 return x; /* log2(+inf) = +inf */
768 else {
769 errno = EDOM;
770 return Py_NAN; /* log2(-inf) = nan, invalid-operation */
771 }
772 }
773
774 if (x > 0.0) {
Victor Stinner8f9f8d62011-05-09 12:45:41 +0200775#ifdef HAVE_LOG2
776 return log2(x);
777#else
Victor Stinnerfa0e3d52011-05-09 01:01:09 +0200778 double m;
779 int e;
780 m = frexp(x, &e);
781 /* We want log2(m * 2**e) == log(m) / log(2) + e. Care is needed when
782 * x is just greater than 1.0: in that case e is 1, log(m) is negative,
783 * and we get significant cancellation error from the addition of
784 * log(m) / log(2) to e. The slight rewrite of the expression below
785 * avoids this problem.
786 */
787 if (x >= 1.0) {
788 return log(2.0 * m) / log(2.0) + (e - 1);
789 }
790 else {
791 return log(m) / log(2.0) + e;
792 }
Victor Stinner8f9f8d62011-05-09 12:45:41 +0200793#endif
Victor Stinnerfa0e3d52011-05-09 01:01:09 +0200794 }
795 else if (x == 0.0) {
796 errno = EDOM;
797 return -Py_HUGE_VAL; /* log2(0) = -inf, divide-by-zero */
798 }
799 else {
800 errno = EDOM;
Mark Dickinson23442582011-05-09 08:05:00 +0100801 return Py_NAN; /* log2(-inf) = nan, invalid-operation */
Victor Stinnerfa0e3d52011-05-09 01:01:09 +0200802 }
803}
804
Mark Dickinsone675f082008-12-11 21:56:00 +0000805static double
806m_log10(double x)
807{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000808 if (Py_IS_FINITE(x)) {
809 if (x > 0.0)
810 return log10(x);
811 errno = EDOM;
812 if (x == 0.0)
813 return -Py_HUGE_VAL; /* log10(0) = -inf */
814 else
815 return Py_NAN; /* log10(-ve) = nan */
816 }
817 else if (Py_IS_NAN(x))
818 return x; /* log10(nan) = nan */
819 else if (x > 0.0)
820 return x; /* log10(inf) = inf */
821 else {
822 errno = EDOM;
823 return Py_NAN; /* log10(-inf) = nan */
824 }
Mark Dickinsone675f082008-12-11 21:56:00 +0000825}
826
827
Serhiy Storchakac9ea9332017-01-19 18:13:09 +0200828/*[clinic input]
829math.gcd
Serhiy Storchaka48e47aa2015-05-13 00:19:51 +0300830
Serhiy Storchakac9ea9332017-01-19 18:13:09 +0200831 x as a: object
832 y as b: object
833 /
834
835greatest common divisor of x and y
836[clinic start generated code]*/
837
838static PyObject *
839math_gcd_impl(PyObject *module, PyObject *a, PyObject *b)
840/*[clinic end generated code: output=7b2e0c151bd7a5d8 input=c2691e57fb2a98fa]*/
841{
842 PyObject *g;
Serhiy Storchaka48e47aa2015-05-13 00:19:51 +0300843
844 a = PyNumber_Index(a);
845 if (a == NULL)
846 return NULL;
847 b = PyNumber_Index(b);
848 if (b == NULL) {
849 Py_DECREF(a);
850 return NULL;
851 }
852 g = _PyLong_GCD(a, b);
853 Py_DECREF(a);
854 Py_DECREF(b);
855 return g;
856}
857
Serhiy Storchaka48e47aa2015-05-13 00:19:51 +0300858
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000859/* Call is_error when errno != 0, and where x is the result libm
860 * returned. is_error will usually set up an exception and return
861 * true (1), but may return false (0) without setting up an exception.
862 */
863static int
864is_error(double x)
865{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000866 int result = 1; /* presumption of guilt */
867 assert(errno); /* non-zero errno is a precondition for calling */
868 if (errno == EDOM)
869 PyErr_SetString(PyExc_ValueError, "math domain error");
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000870
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000871 else if (errno == ERANGE) {
872 /* ANSI C generally requires libm functions to set ERANGE
873 * on overflow, but also generally *allows* them to set
874 * ERANGE on underflow too. There's no consistency about
875 * the latter across platforms.
876 * Alas, C99 never requires that errno be set.
877 * Here we suppress the underflow errors (libm functions
878 * should return a zero on underflow, and +- HUGE_VAL on
879 * overflow, so testing the result for zero suffices to
880 * distinguish the cases).
881 *
882 * On some platforms (Ubuntu/ia64) it seems that errno can be
883 * set to ERANGE for subnormal results that do *not* underflow
884 * to zero. So to be safe, we'll ignore ERANGE whenever the
885 * function result is less than one in absolute value.
886 */
887 if (fabs(x) < 1.0)
888 result = 0;
889 else
890 PyErr_SetString(PyExc_OverflowError,
891 "math range error");
892 }
893 else
894 /* Unexpected math error */
895 PyErr_SetFromErrno(PyExc_ValueError);
896 return result;
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000897}
898
Mark Dickinsone675f082008-12-11 21:56:00 +0000899/*
Christian Heimes53876d92008-04-19 00:31:39 +0000900 math_1 is used to wrap a libm function f that takes a double
Serhiy Storchakac9ea9332017-01-19 18:13:09 +0200901 argument and returns a double.
Christian Heimes53876d92008-04-19 00:31:39 +0000902
903 The error reporting follows these rules, which are designed to do
904 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
905 platforms.
906
907 - a NaN result from non-NaN inputs causes ValueError to be raised
908 - an infinite result from finite inputs causes OverflowError to be
909 raised if can_overflow is 1, or raises ValueError if can_overflow
910 is 0.
911 - if the result is finite and errno == EDOM then ValueError is
912 raised
913 - if the result is finite and nonzero and errno == ERANGE then
914 OverflowError is raised
915
916 The last rule is used to catch overflow on platforms which follow
917 C89 but for which HUGE_VAL is not an infinity.
918
919 For the majority of one-argument functions these rules are enough
920 to ensure that Python's functions behave as specified in 'Annex F'
921 of the C99 standard, with the 'invalid' and 'divide-by-zero'
922 floating-point exceptions mapping to Python's ValueError and the
923 'overflow' floating-point exception mapping to OverflowError.
924 math_1 only works for functions that don't have singularities *and*
925 the possibility of overflow; fortunately, that covers everything we
926 care about right now.
927*/
928
Barry Warsaw8b43b191996-12-09 22:32:36 +0000929static PyObject *
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000930math_1_to_whatever(PyObject *arg, double (*func) (double),
Christian Heimes53876d92008-04-19 00:31:39 +0000931 PyObject *(*from_double_func) (double),
932 int can_overflow)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000933{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000934 double x, r;
935 x = PyFloat_AsDouble(arg);
936 if (x == -1.0 && PyErr_Occurred())
937 return NULL;
938 errno = 0;
939 PyFPE_START_PROTECT("in math_1", return 0);
940 r = (*func)(x);
941 PyFPE_END_PROTECT(r);
942 if (Py_IS_NAN(r) && !Py_IS_NAN(x)) {
943 PyErr_SetString(PyExc_ValueError,
944 "math domain error"); /* invalid arg */
945 return NULL;
946 }
947 if (Py_IS_INFINITY(r) && Py_IS_FINITE(x)) {
Benjamin Peterson2354a752012-03-13 16:13:09 -0500948 if (can_overflow)
949 PyErr_SetString(PyExc_OverflowError,
950 "math range error"); /* overflow */
951 else
952 PyErr_SetString(PyExc_ValueError,
953 "math domain error"); /* singularity */
954 return NULL;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000955 }
956 if (Py_IS_FINITE(r) && errno && is_error(r))
957 /* this branch unnecessary on most platforms */
958 return NULL;
Mark Dickinsonde429622008-05-01 00:19:23 +0000959
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000960 return (*from_double_func)(r);
Christian Heimes53876d92008-04-19 00:31:39 +0000961}
962
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000963/* variant of math_1, to be used when the function being wrapped is known to
964 set errno properly (that is, errno = EDOM for invalid or divide-by-zero,
965 errno = ERANGE for overflow). */
966
967static PyObject *
968math_1a(PyObject *arg, double (*func) (double))
969{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +0000970 double x, r;
971 x = PyFloat_AsDouble(arg);
972 if (x == -1.0 && PyErr_Occurred())
973 return NULL;
974 errno = 0;
975 PyFPE_START_PROTECT("in math_1a", return 0);
976 r = (*func)(x);
977 PyFPE_END_PROTECT(r);
978 if (errno && is_error(r))
979 return NULL;
980 return PyFloat_FromDouble(r);
Mark Dickinson12c4bdb2009-09-28 19:21:11 +0000981}
982
Christian Heimes53876d92008-04-19 00:31:39 +0000983/*
984 math_2 is used to wrap a libm function f that takes two double
985 arguments and returns a double.
986
987 The error reporting follows these rules, which are designed to do
988 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
989 platforms.
990
991 - a NaN result from non-NaN inputs causes ValueError to be raised
992 - an infinite result from finite inputs causes OverflowError to be
993 raised.
994 - if the result is finite and errno == EDOM then ValueError is
995 raised
996 - if the result is finite and nonzero and errno == ERANGE then
997 OverflowError is raised
998
999 The last rule is used to catch overflow on platforms which follow
1000 C89 but for which HUGE_VAL is not an infinity.
1001
1002 For most two-argument functions (copysign, fmod, hypot, atan2)
1003 these rules are enough to ensure that Python's functions behave as
1004 specified in 'Annex F' of the C99 standard, with the 'invalid' and
1005 'divide-by-zero' floating-point exceptions mapping to Python's
1006 ValueError and the 'overflow' floating-point exception mapping to
1007 OverflowError.
1008*/
1009
1010static PyObject *
1011math_1(PyObject *arg, double (*func) (double), int can_overflow)
1012{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001013 return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow);
Jeffrey Yasskinc2155832008-01-05 20:03:11 +00001014}
1015
1016static PyObject *
Christian Heimes53876d92008-04-19 00:31:39 +00001017math_1_to_int(PyObject *arg, double (*func) (double), int can_overflow)
Jeffrey Yasskinc2155832008-01-05 20:03:11 +00001018{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001019 return math_1_to_whatever(arg, func, PyLong_FromDouble, can_overflow);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001020}
1021
Barry Warsaw8b43b191996-12-09 22:32:36 +00001022static PyObject *
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02001023math_2(PyObject *const *args, Py_ssize_t nargs,
1024 double (*func) (double, double), const char *funcname)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001025{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001026 double x, y, r;
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02001027 if (!_PyArg_CheckPositional(funcname, nargs, 2, 2))
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001028 return NULL;
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02001029 x = PyFloat_AsDouble(args[0]);
1030 y = PyFloat_AsDouble(args[1]);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001031 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1032 return NULL;
1033 errno = 0;
1034 PyFPE_START_PROTECT("in math_2", return 0);
1035 r = (*func)(x, y);
1036 PyFPE_END_PROTECT(r);
1037 if (Py_IS_NAN(r)) {
1038 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1039 errno = EDOM;
1040 else
1041 errno = 0;
1042 }
1043 else if (Py_IS_INFINITY(r)) {
1044 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
1045 errno = ERANGE;
1046 else
1047 errno = 0;
1048 }
1049 if (errno && is_error(r))
1050 return NULL;
1051 else
1052 return PyFloat_FromDouble(r);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001053}
1054
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001055#define FUNC1(funcname, func, can_overflow, docstring) \
1056 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
1057 return math_1(args, func, can_overflow); \
1058 }\
1059 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001060
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001061#define FUNC1A(funcname, func, docstring) \
1062 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
1063 return math_1a(args, func); \
1064 }\
1065 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Mark Dickinson12c4bdb2009-09-28 19:21:11 +00001066
Fred Drake40c48682000-07-03 18:11:56 +00001067#define FUNC2(funcname, func, docstring) \
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02001068 static PyObject * math_##funcname(PyObject *self, PyObject *const *args, Py_ssize_t nargs) { \
1069 return math_2(args, nargs, func, #funcname); \
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001070 }\
1071 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001072
Christian Heimes53876d92008-04-19 00:31:39 +00001073FUNC1(acos, acos, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001074 "acos($module, x, /)\n--\n\n"
1075 "Return the arc cosine (measured in radians) of x.")
Mark Dickinsonf3718592009-12-21 15:27:41 +00001076FUNC1(acosh, m_acosh, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001077 "acosh($module, x, /)\n--\n\n"
1078 "Return the inverse hyperbolic cosine of x.")
Christian Heimes53876d92008-04-19 00:31:39 +00001079FUNC1(asin, asin, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001080 "asin($module, x, /)\n--\n\n"
1081 "Return the arc sine (measured in radians) of x.")
Mark Dickinsonf3718592009-12-21 15:27:41 +00001082FUNC1(asinh, m_asinh, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001083 "asinh($module, x, /)\n--\n\n"
1084 "Return the inverse hyperbolic sine of x.")
Christian Heimes53876d92008-04-19 00:31:39 +00001085FUNC1(atan, atan, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001086 "atan($module, x, /)\n--\n\n"
1087 "Return the arc tangent (measured in radians) of x.")
Christian Heimese57950f2008-04-21 13:08:03 +00001088FUNC2(atan2, m_atan2,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001089 "atan2($module, y, x, /)\n--\n\n"
1090 "Return the arc tangent (measured in radians) of y/x.\n\n"
Tim Petersfe71f812001-08-07 22:10:00 +00001091 "Unlike atan(y/x), the signs of both x and y are considered.")
Mark Dickinsonf3718592009-12-21 15:27:41 +00001092FUNC1(atanh, m_atanh, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001093 "atanh($module, x, /)\n--\n\n"
1094 "Return the inverse hyperbolic tangent of x.")
Guido van Rossum13e05de2007-08-23 22:56:55 +00001095
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001096/*[clinic input]
1097math.ceil
1098
1099 x as number: object
1100 /
1101
1102Return the ceiling of x as an Integral.
1103
1104This is the smallest integer >= x.
1105[clinic start generated code]*/
1106
1107static PyObject *
1108math_ceil(PyObject *module, PyObject *number)
1109/*[clinic end generated code: output=6c3b8a78bc201c67 input=2725352806399cab]*/
1110{
Benjamin Petersonce798522012-01-22 11:24:29 -05001111 _Py_IDENTIFIER(__ceil__);
Mark Dickinson6d02d9c2010-07-02 16:05:15 +00001112 PyObject *method, *result;
Guido van Rossum13e05de2007-08-23 22:56:55 +00001113
Benjamin Petersonce798522012-01-22 11:24:29 -05001114 method = _PyObject_LookupSpecial(number, &PyId___ceil__);
Benjamin Petersonf751bc92010-07-02 13:46:42 +00001115 if (method == NULL) {
1116 if (PyErr_Occurred())
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001117 return NULL;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001118 return math_1_to_int(number, ceil, 0);
Benjamin Petersonf751bc92010-07-02 13:46:42 +00001119 }
Victor Stinnerf17c3de2016-12-06 18:46:19 +01001120 result = _PyObject_CallNoArg(method);
Mark Dickinson6d02d9c2010-07-02 16:05:15 +00001121 Py_DECREF(method);
1122 return result;
Guido van Rossum13e05de2007-08-23 22:56:55 +00001123}
1124
Christian Heimes072c0f12008-01-03 23:01:04 +00001125FUNC2(copysign, copysign,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001126 "copysign($module, x, y, /)\n--\n\n"
1127 "Return a float with the magnitude (absolute value) of x but the sign of y.\n\n"
1128 "On platforms that support signed zeros, copysign(1.0, -0.0)\n"
1129 "returns -1.0.\n")
Christian Heimes53876d92008-04-19 00:31:39 +00001130FUNC1(cos, cos, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001131 "cos($module, x, /)\n--\n\n"
1132 "Return the cosine of x (measured in radians).")
Christian Heimes53876d92008-04-19 00:31:39 +00001133FUNC1(cosh, cosh, 1,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001134 "cosh($module, x, /)\n--\n\n"
1135 "Return the hyperbolic cosine of x.")
Mark Dickinson45f992a2009-12-19 11:20:49 +00001136FUNC1A(erf, m_erf,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001137 "erf($module, x, /)\n--\n\n"
1138 "Error function at x.")
Mark Dickinson45f992a2009-12-19 11:20:49 +00001139FUNC1A(erfc, m_erfc,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001140 "erfc($module, x, /)\n--\n\n"
1141 "Complementary error function at x.")
Christian Heimes53876d92008-04-19 00:31:39 +00001142FUNC1(exp, exp, 1,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001143 "exp($module, x, /)\n--\n\n"
1144 "Return e raised to the power of x.")
Mark Dickinson664b5112009-12-16 20:23:42 +00001145FUNC1(expm1, m_expm1, 1,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001146 "expm1($module, x, /)\n--\n\n"
1147 "Return exp(x)-1.\n\n"
Mark Dickinson664b5112009-12-16 20:23:42 +00001148 "This function avoids the loss of precision involved in the direct "
1149 "evaluation of exp(x)-1 for small x.")
Christian Heimes53876d92008-04-19 00:31:39 +00001150FUNC1(fabs, fabs, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001151 "fabs($module, x, /)\n--\n\n"
1152 "Return the absolute value of the float x.")
Guido van Rossum13e05de2007-08-23 22:56:55 +00001153
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001154/*[clinic input]
1155math.floor
1156
1157 x as number: object
1158 /
1159
1160Return the floor of x as an Integral.
1161
1162This is the largest integer <= x.
1163[clinic start generated code]*/
1164
1165static PyObject *
1166math_floor(PyObject *module, PyObject *number)
1167/*[clinic end generated code: output=c6a65c4884884b8a input=63af6b5d7ebcc3d6]*/
1168{
Benjamin Petersonce798522012-01-22 11:24:29 -05001169 _Py_IDENTIFIER(__floor__);
Benjamin Petersonb0125892010-07-02 13:35:17 +00001170 PyObject *method, *result;
Guido van Rossum13e05de2007-08-23 22:56:55 +00001171
Benjamin Petersonce798522012-01-22 11:24:29 -05001172 method = _PyObject_LookupSpecial(number, &PyId___floor__);
Benjamin Peterson8bb9cde2010-07-01 15:16:55 +00001173 if (method == NULL) {
1174 if (PyErr_Occurred())
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001175 return NULL;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001176 return math_1_to_int(number, floor, 0);
Benjamin Peterson8bb9cde2010-07-01 15:16:55 +00001177 }
Victor Stinnerf17c3de2016-12-06 18:46:19 +01001178 result = _PyObject_CallNoArg(method);
Benjamin Petersonb0125892010-07-02 13:35:17 +00001179 Py_DECREF(method);
1180 return result;
Guido van Rossum13e05de2007-08-23 22:56:55 +00001181}
1182
Mark Dickinson12c4bdb2009-09-28 19:21:11 +00001183FUNC1A(gamma, m_tgamma,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001184 "gamma($module, x, /)\n--\n\n"
1185 "Gamma function at x.")
Mark Dickinson05d2e082009-12-11 20:17:17 +00001186FUNC1A(lgamma, m_lgamma,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001187 "lgamma($module, x, /)\n--\n\n"
1188 "Natural logarithm of absolute value of Gamma function at x.")
Mark Dickinsonbe64d952010-07-07 16:21:29 +00001189FUNC1(log1p, m_log1p, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001190 "log1p($module, x, /)\n--\n\n"
1191 "Return the natural logarithm of 1+x (base e).\n\n"
Benjamin Petersona0dfa822009-11-13 02:25:08 +00001192 "The result is computed in a way which is accurate for x near zero.")
Mark Dickinsona0ce3752017-04-05 18:34:27 +01001193FUNC2(remainder, m_remainder,
1194 "remainder($module, x, y, /)\n--\n\n"
1195 "Difference between x and the closest integer multiple of y.\n\n"
1196 "Return x - n*y where n*y is the closest integer multiple of y.\n"
1197 "In the case where x is exactly halfway between two multiples of\n"
1198 "y, the nearest even value of n is used. The result is always exact.")
Christian Heimes53876d92008-04-19 00:31:39 +00001199FUNC1(sin, sin, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001200 "sin($module, x, /)\n--\n\n"
1201 "Return the sine of x (measured in radians).")
Christian Heimes53876d92008-04-19 00:31:39 +00001202FUNC1(sinh, sinh, 1,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001203 "sinh($module, x, /)\n--\n\n"
1204 "Return the hyperbolic sine of x.")
Christian Heimes53876d92008-04-19 00:31:39 +00001205FUNC1(sqrt, sqrt, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001206 "sqrt($module, x, /)\n--\n\n"
1207 "Return the square root of x.")
Christian Heimes53876d92008-04-19 00:31:39 +00001208FUNC1(tan, tan, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001209 "tan($module, x, /)\n--\n\n"
1210 "Return the tangent of x (measured in radians).")
Christian Heimes53876d92008-04-19 00:31:39 +00001211FUNC1(tanh, tanh, 0,
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001212 "tanh($module, x, /)\n--\n\n"
1213 "Return the hyperbolic tangent of x.")
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001214
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001215/* Precision summation function as msum() by Raymond Hettinger in
1216 <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
1217 enhanced with the exact partials sum and roundoff from Mark
1218 Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
1219 See those links for more details, proofs and other references.
1220
1221 Note 1: IEEE 754R floating point semantics are assumed,
1222 but the current implementation does not re-establish special
1223 value semantics across iterations (i.e. handling -Inf + Inf).
1224
1225 Note 2: No provision is made for intermediate overflow handling;
Georg Brandlf78e02b2008-06-10 17:40:04 +00001226 therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001227 sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
1228 overflow of the first partial sum.
1229
Benjamin Petersonfea6a942008-07-02 16:11:42 +00001230 Note 3: The intermediate values lo, yr, and hi are declared volatile so
1231 aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
Georg Brandlf78e02b2008-06-10 17:40:04 +00001232 Also, the volatile declaration forces the values to be stored in memory as
1233 regular doubles instead of extended long precision (80-bit) values. This
Benjamin Petersonfea6a942008-07-02 16:11:42 +00001234 prevents double rounding because any addition or subtraction of two doubles
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001235 can be resolved exactly into double-sized hi and lo values. As long as the
Georg Brandlf78e02b2008-06-10 17:40:04 +00001236 hi value gets forced into a double before yr and lo are computed, the extra
1237 bits in downstream extended precision operations (x87 for example) will be
1238 exactly zero and therefore can be losslessly stored back into a double,
1239 thereby preventing double rounding.
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001240
1241 Note 4: A similar implementation is in Modules/cmathmodule.c.
1242 Be sure to update both when making changes.
1243
Serhiy Storchakaa60c2fe2015-03-12 21:56:08 +02001244 Note 5: The signature of math.fsum() differs from builtins.sum()
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001245 because the start argument doesn't make sense in the context of
1246 accurate summation. Since the partials table is collapsed before
1247 returning a result, sum(seq2, start=sum(seq1)) may not equal the
1248 accurate result returned by sum(itertools.chain(seq1, seq2)).
1249*/
1250
1251#define NUM_PARTIALS 32 /* initial partials array size, on stack */
1252
1253/* Extend the partials array p[] by doubling its size. */
1254static int /* non-zero on error */
Mark Dickinsonaa7633a2008-08-01 08:16:13 +00001255_fsum_realloc(double **p_ptr, Py_ssize_t n,
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001256 double *ps, Py_ssize_t *m_ptr)
1257{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001258 void *v = NULL;
1259 Py_ssize_t m = *m_ptr;
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001260
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001261 m += m; /* double */
Victor Stinner049e5092014-08-17 22:20:00 +02001262 if (n < m && (size_t)m < ((size_t)PY_SSIZE_T_MAX / sizeof(double))) {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001263 double *p = *p_ptr;
1264 if (p == ps) {
1265 v = PyMem_Malloc(sizeof(double) * m);
1266 if (v != NULL)
1267 memcpy(v, ps, sizeof(double) * n);
1268 }
1269 else
1270 v = PyMem_Realloc(p, sizeof(double) * m);
1271 }
1272 if (v == NULL) { /* size overflow or no memory */
1273 PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
1274 return 1;
1275 }
1276 *p_ptr = (double*) v;
1277 *m_ptr = m;
1278 return 0;
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001279}
1280
1281/* Full precision summation of a sequence of floats.
1282
1283 def msum(iterable):
1284 partials = [] # sorted, non-overlapping partial sums
1285 for x in iterable:
Mark Dickinsonfdb0acc2010-06-25 20:22:24 +00001286 i = 0
1287 for y in partials:
1288 if abs(x) < abs(y):
1289 x, y = y, x
1290 hi = x + y
1291 lo = y - (hi - x)
1292 if lo:
1293 partials[i] = lo
1294 i += 1
1295 x = hi
1296 partials[i:] = [x]
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001297 return sum_exact(partials)
1298
1299 Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo
1300 are exactly equal to x+y. The inner loop applies hi/lo summation to each
1301 partial so that the list of partial sums remains exact.
1302
1303 Sum_exact() adds the partial sums exactly and correctly rounds the final
1304 result (using the round-half-to-even rule). The items in partials remain
1305 non-zero, non-special, non-overlapping and strictly increasing in
1306 magnitude, but possibly not all having the same sign.
1307
1308 Depends on IEEE 754 arithmetic guarantees and half-even rounding.
1309*/
1310
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001311/*[clinic input]
1312math.fsum
1313
1314 seq: object
1315 /
1316
1317Return an accurate floating point sum of values in the iterable seq.
1318
1319Assumes IEEE-754 floating point arithmetic.
1320[clinic start generated code]*/
1321
1322static PyObject *
1323math_fsum(PyObject *module, PyObject *seq)
1324/*[clinic end generated code: output=ba5c672b87fe34fc input=c51b7d8caf6f6e82]*/
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001325{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001326 PyObject *item, *iter, *sum = NULL;
1327 Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
1328 double x, y, t, ps[NUM_PARTIALS], *p = ps;
1329 double xsave, special_sum = 0.0, inf_sum = 0.0;
1330 volatile double hi, yr, lo;
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001331
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001332 iter = PyObject_GetIter(seq);
1333 if (iter == NULL)
1334 return NULL;
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001335
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001336 PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001337
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001338 for(;;) { /* for x in iterable */
1339 assert(0 <= n && n <= m);
1340 assert((m == NUM_PARTIALS && p == ps) ||
1341 (m > NUM_PARTIALS && p != NULL));
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001342
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001343 item = PyIter_Next(iter);
1344 if (item == NULL) {
1345 if (PyErr_Occurred())
1346 goto _fsum_error;
1347 break;
1348 }
Raymond Hettingercfd735e2019-01-29 20:39:53 -08001349 ASSIGN_DOUBLE(x, item, error_with_item);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001350 Py_DECREF(item);
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001351
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001352 xsave = x;
1353 for (i = j = 0; j < n; j++) { /* for y in partials */
1354 y = p[j];
1355 if (fabs(x) < fabs(y)) {
1356 t = x; x = y; y = t;
1357 }
1358 hi = x + y;
1359 yr = hi - x;
1360 lo = y - yr;
1361 if (lo != 0.0)
1362 p[i++] = lo;
1363 x = hi;
1364 }
Mark Dickinsonaa7633a2008-08-01 08:16:13 +00001365
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001366 n = i; /* ps[i:] = [x] */
1367 if (x != 0.0) {
1368 if (! Py_IS_FINITE(x)) {
1369 /* a nonfinite x could arise either as
1370 a result of intermediate overflow, or
1371 as a result of a nan or inf in the
1372 summands */
1373 if (Py_IS_FINITE(xsave)) {
1374 PyErr_SetString(PyExc_OverflowError,
1375 "intermediate overflow in fsum");
1376 goto _fsum_error;
1377 }
1378 if (Py_IS_INFINITY(xsave))
1379 inf_sum += xsave;
1380 special_sum += xsave;
1381 /* reset partials */
1382 n = 0;
1383 }
1384 else if (n >= m && _fsum_realloc(&p, n, ps, &m))
1385 goto _fsum_error;
1386 else
1387 p[n++] = x;
1388 }
1389 }
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001390
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001391 if (special_sum != 0.0) {
1392 if (Py_IS_NAN(inf_sum))
1393 PyErr_SetString(PyExc_ValueError,
1394 "-inf + inf in fsum");
1395 else
1396 sum = PyFloat_FromDouble(special_sum);
1397 goto _fsum_error;
1398 }
Mark Dickinsonaa7633a2008-08-01 08:16:13 +00001399
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001400 hi = 0.0;
1401 if (n > 0) {
1402 hi = p[--n];
1403 /* sum_exact(ps, hi) from the top, stop when the sum becomes
1404 inexact. */
1405 while (n > 0) {
1406 x = hi;
1407 y = p[--n];
1408 assert(fabs(y) < fabs(x));
1409 hi = x + y;
1410 yr = hi - x;
1411 lo = y - yr;
1412 if (lo != 0.0)
1413 break;
1414 }
1415 /* Make half-even rounding work across multiple partials.
1416 Needed so that sum([1e-16, 1, 1e16]) will round-up the last
1417 digit to two instead of down to zero (the 1e-16 makes the 1
1418 slightly closer to two). With a potential 1 ULP rounding
1419 error fixed-up, math.fsum() can guarantee commutativity. */
1420 if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
1421 (lo > 0.0 && p[n-1] > 0.0))) {
1422 y = lo * 2.0;
1423 x = hi + y;
1424 yr = x - hi;
1425 if (y == yr)
1426 hi = x;
1427 }
1428 }
1429 sum = PyFloat_FromDouble(hi);
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001430
Raymond Hettingercfd735e2019-01-29 20:39:53 -08001431 _fsum_error:
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001432 PyFPE_END_PROTECT(hi)
1433 Py_DECREF(iter);
1434 if (p != ps)
1435 PyMem_Free(p);
1436 return sum;
Raymond Hettingercfd735e2019-01-29 20:39:53 -08001437
1438 error_with_item:
1439 Py_DECREF(item);
1440 goto _fsum_error;
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001441}
1442
1443#undef NUM_PARTIALS
1444
Benjamin Peterson2b7411d2008-05-26 17:36:47 +00001445
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001446/* Return the smallest integer k such that n < 2**k, or 0 if n == 0.
1447 * Equivalent to floor(lg(x))+1. Also equivalent to: bitwidth_of_type -
1448 * count_leading_zero_bits(x)
1449 */
1450
1451/* XXX: This routine does more or less the same thing as
1452 * bits_in_digit() in Objects/longobject.c. Someday it would be nice to
1453 * consolidate them. On BSD, there's a library function called fls()
1454 * that we could use, and GCC provides __builtin_clz().
1455 */
1456
1457static unsigned long
1458bit_length(unsigned long n)
1459{
1460 unsigned long len = 0;
1461 while (n != 0) {
1462 ++len;
1463 n >>= 1;
1464 }
1465 return len;
1466}
1467
1468static unsigned long
1469count_set_bits(unsigned long n)
1470{
1471 unsigned long count = 0;
1472 while (n != 0) {
1473 ++count;
1474 n &= n - 1; /* clear least significant bit */
1475 }
1476 return count;
1477}
1478
Mark Dickinson73934b92019-05-18 12:29:50 +01001479/* Integer square root
1480
1481Given a nonnegative integer `n`, we want to compute the largest integer
1482`a` for which `a * a <= n`, or equivalently the integer part of the exact
1483square root of `n`.
1484
1485We use an adaptive-precision pure-integer version of Newton's iteration. Given
1486a positive integer `n`, the algorithm produces at each iteration an integer
1487approximation `a` to the square root of `n >> s` for some even integer `s`,
1488with `s` decreasing as the iterations progress. On the final iteration, `s` is
1489zero and we have an approximation to the square root of `n` itself.
1490
1491At every step, the approximation `a` is strictly within 1.0 of the true square
1492root, so we have
1493
1494 (a - 1)**2 < (n >> s) < (a + 1)**2
1495
1496After the final iteration, a check-and-correct step is needed to determine
1497whether `a` or `a - 1` gives the desired integer square root of `n`.
1498
1499The algorithm is remarkable in its simplicity. There's no need for a
1500per-iteration check-and-correct step, and termination is straightforward: the
1501number of iterations is known in advance (it's exactly `floor(log2(log2(n)))`
1502for `n > 1`). The only tricky part of the correctness proof is in establishing
1503that the bound `(a - 1)**2 < (n >> s) < (a + 1)**2` is maintained from one
1504iteration to the next. A sketch of the proof of this is given below.
1505
1506In addition to the proof sketch, a formal, computer-verified proof
1507of correctness (using Lean) of an equivalent recursive algorithm can be found
1508here:
1509
1510 https://github.com/mdickinson/snippets/blob/master/proofs/isqrt/src/isqrt.lean
1511
1512
1513Here's Python code equivalent to the C implementation below:
1514
1515 def isqrt(n):
1516 """
1517 Return the integer part of the square root of the input.
1518 """
1519 n = operator.index(n)
1520
1521 if n < 0:
1522 raise ValueError("isqrt() argument must be nonnegative")
1523 if n == 0:
1524 return 0
1525
1526 c = (n.bit_length() - 1) // 2
1527 a = 1
1528 d = 0
1529 for s in reversed(range(c.bit_length())):
Mark Dickinson2dfeaa92019-06-16 17:53:21 +01001530 # Loop invariant: (a-1)**2 < (n >> 2*(c - d)) < (a+1)**2
Mark Dickinson73934b92019-05-18 12:29:50 +01001531 e = d
1532 d = c >> s
1533 a = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a
Mark Dickinson73934b92019-05-18 12:29:50 +01001534
1535 return a - (a*a > n)
1536
1537
1538Sketch of proof of correctness
1539------------------------------
1540
1541The delicate part of the correctness proof is showing that the loop invariant
1542is preserved from one iteration to the next. That is, just before the line
1543
1544 a = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a
1545
1546is executed in the above code, we know that
1547
1548 (1) (a - 1)**2 < (n >> 2*(c - e)) < (a + 1)**2.
1549
1550(since `e` is always the value of `d` from the previous iteration). We must
1551prove that after that line is executed, we have
1552
1553 (a - 1)**2 < (n >> 2*(c - d)) < (a + 1)**2
1554
Min ho Kimf7d72e42019-07-06 07:39:32 +10001555To facilitate the proof, we make some changes of notation. Write `m` for
Mark Dickinson73934b92019-05-18 12:29:50 +01001556`n >> 2*(c-d)`, and write `b` for the new value of `a`, so
1557
1558 b = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a
1559
1560or equivalently:
1561
1562 (2) b = (a << d - e - 1) + (m >> d - e + 1) // a
1563
1564Then we can rewrite (1) as:
1565
1566 (3) (a - 1)**2 < (m >> 2*(d - e)) < (a + 1)**2
1567
1568and we must show that (b - 1)**2 < m < (b + 1)**2.
1569
1570From this point on, we switch to mathematical notation, so `/` means exact
1571division rather than integer division and `^` is used for exponentiation. We
1572use the `√` symbol for the exact square root. In (3), we can remove the
1573implicit floor operation to give:
1574
1575 (4) (a - 1)^2 < m / 4^(d - e) < (a + 1)^2
1576
1577Taking square roots throughout (4), scaling by `2^(d-e)`, and rearranging gives
1578
1579 (5) 0 <= | 2^(d-e)a - √m | < 2^(d-e)
1580
1581Squaring and dividing through by `2^(d-e+1) a` gives
1582
1583 (6) 0 <= 2^(d-e-1) a + m / (2^(d-e+1) a) - √m < 2^(d-e-1) / a
1584
1585We'll show below that `2^(d-e-1) <= a`. Given that, we can replace the
1586right-hand side of (6) with `1`, and now replacing the central
1587term `m / (2^(d-e+1) a)` with its floor in (6) gives
1588
1589 (7) -1 < 2^(d-e-1) a + m // 2^(d-e+1) a - √m < 1
1590
1591Or equivalently, from (2):
1592
1593 (7) -1 < b - √m < 1
1594
1595and rearranging gives that `(b-1)^2 < m < (b+1)^2`, which is what we needed
1596to prove.
1597
1598We're not quite done: we still have to prove the inequality `2^(d - e - 1) <=
1599a` that was used to get line (7) above. From the definition of `c`, we have
1600`4^c <= n`, which implies
1601
1602 (8) 4^d <= m
1603
1604also, since `e == d >> 1`, `d` is at most `2e + 1`, from which it follows
1605that `2d - 2e - 1 <= d` and hence that
1606
1607 (9) 4^(2d - 2e - 1) <= m
1608
1609Dividing both sides by `4^(d - e)` gives
1610
1611 (10) 4^(d - e - 1) <= m / 4^(d - e)
1612
1613But we know from (4) that `m / 4^(d-e) < (a + 1)^2`, hence
1614
1615 (11) 4^(d - e - 1) < (a + 1)^2
1616
1617Now taking square roots of both sides and observing that both `2^(d-e-1)` and
1618`a` are integers gives `2^(d - e - 1) <= a`, which is what we needed. This
1619completes the proof sketch.
1620
1621*/
1622
Mark Dickinson5c08ce92019-05-19 17:51:56 +01001623
1624/* Approximate square root of a large 64-bit integer.
1625
1626 Given `n` satisfying `2**62 <= n < 2**64`, return `a`
1627 satisfying `(a - 1)**2 < n < (a + 1)**2`. */
1628
1629static uint64_t
1630_approximate_isqrt(uint64_t n)
1631{
1632 uint32_t u = 1U + (n >> 62);
1633 u = (u << 1) + (n >> 59) / u;
1634 u = (u << 3) + (n >> 53) / u;
1635 u = (u << 7) + (n >> 41) / u;
1636 return (u << 15) + (n >> 17) / u;
1637}
1638
Mark Dickinson73934b92019-05-18 12:29:50 +01001639/*[clinic input]
1640math.isqrt
1641
1642 n: object
1643 /
1644
1645Return the integer part of the square root of the input.
1646[clinic start generated code]*/
1647
1648static PyObject *
1649math_isqrt(PyObject *module, PyObject *n)
1650/*[clinic end generated code: output=35a6f7f980beab26 input=5b6e7ae4fa6c43d6]*/
1651{
Mark Dickinson5c08ce92019-05-19 17:51:56 +01001652 int a_too_large, c_bit_length;
Mark Dickinson73934b92019-05-18 12:29:50 +01001653 size_t c, d;
Mark Dickinson5c08ce92019-05-19 17:51:56 +01001654 uint64_t m, u;
Mark Dickinson73934b92019-05-18 12:29:50 +01001655 PyObject *a = NULL, *b;
1656
1657 n = PyNumber_Index(n);
1658 if (n == NULL) {
1659 return NULL;
1660 }
1661
1662 if (_PyLong_Sign(n) < 0) {
1663 PyErr_SetString(
1664 PyExc_ValueError,
1665 "isqrt() argument must be nonnegative");
1666 goto error;
1667 }
1668 if (_PyLong_Sign(n) == 0) {
1669 Py_DECREF(n);
1670 return PyLong_FromLong(0);
1671 }
1672
Mark Dickinson5c08ce92019-05-19 17:51:56 +01001673 /* c = (n.bit_length() - 1) // 2 */
Mark Dickinson73934b92019-05-18 12:29:50 +01001674 c = _PyLong_NumBits(n);
1675 if (c == (size_t)(-1)) {
1676 goto error;
1677 }
1678 c = (c - 1U) / 2U;
1679
Mark Dickinson5c08ce92019-05-19 17:51:56 +01001680 /* Fast path: if c <= 31 then n < 2**64 and we can compute directly with a
1681 fast, almost branch-free algorithm. In the final correction, we use `u*u
1682 - 1 >= m` instead of the simpler `u*u > m` in order to get the correct
1683 result in the corner case where `u=2**32`. */
1684 if (c <= 31U) {
1685 m = (uint64_t)PyLong_AsUnsignedLongLong(n);
1686 Py_DECREF(n);
1687 if (m == (uint64_t)(-1) && PyErr_Occurred()) {
1688 return NULL;
1689 }
1690 u = _approximate_isqrt(m << (62U - 2U*c)) >> (31U - c);
1691 u -= u * u - 1U >= m;
1692 return PyLong_FromUnsignedLongLong((unsigned long long)u);
Mark Dickinson73934b92019-05-18 12:29:50 +01001693 }
1694
Mark Dickinson5c08ce92019-05-19 17:51:56 +01001695 /* Slow path: n >= 2**64. We perform the first five iterations in C integer
1696 arithmetic, then switch to using Python long integers. */
1697
1698 /* From n >= 2**64 it follows that c.bit_length() >= 6. */
1699 c_bit_length = 6;
1700 while ((c >> c_bit_length) > 0U) {
1701 ++c_bit_length;
1702 }
1703
1704 /* Initialise d and a. */
1705 d = c >> (c_bit_length - 5);
1706 b = _PyLong_Rshift(n, 2U*c - 62U);
1707 if (b == NULL) {
1708 goto error;
1709 }
1710 m = (uint64_t)PyLong_AsUnsignedLongLong(b);
1711 Py_DECREF(b);
1712 if (m == (uint64_t)(-1) && PyErr_Occurred()) {
1713 goto error;
1714 }
1715 u = _approximate_isqrt(m) >> (31U - d);
1716 a = PyLong_FromUnsignedLongLong((unsigned long long)u);
Mark Dickinson73934b92019-05-18 12:29:50 +01001717 if (a == NULL) {
1718 goto error;
1719 }
Mark Dickinson5c08ce92019-05-19 17:51:56 +01001720
1721 for (int s = c_bit_length - 6; s >= 0; --s) {
Serhiy Storchakaa5119e72019-05-19 14:14:38 +03001722 PyObject *q;
Mark Dickinson73934b92019-05-18 12:29:50 +01001723 size_t e = d;
1724
1725 d = c >> s;
1726
1727 /* q = (n >> 2*c - e - d + 1) // a */
Serhiy Storchakaa5119e72019-05-19 14:14:38 +03001728 q = _PyLong_Rshift(n, 2U*c - d - e + 1U);
Mark Dickinson73934b92019-05-18 12:29:50 +01001729 if (q == NULL) {
1730 goto error;
1731 }
1732 Py_SETREF(q, PyNumber_FloorDivide(q, a));
1733 if (q == NULL) {
1734 goto error;
1735 }
1736
1737 /* a = (a << d - 1 - e) + q */
Serhiy Storchakaa5119e72019-05-19 14:14:38 +03001738 Py_SETREF(a, _PyLong_Lshift(a, d - 1U - e));
Mark Dickinson73934b92019-05-18 12:29:50 +01001739 if (a == NULL) {
1740 Py_DECREF(q);
1741 goto error;
1742 }
1743 Py_SETREF(a, PyNumber_Add(a, q));
1744 Py_DECREF(q);
1745 if (a == NULL) {
1746 goto error;
1747 }
1748 }
1749
1750 /* The correct result is either a or a - 1. Figure out which, and
1751 decrement a if necessary. */
1752
1753 /* a_too_large = n < a * a */
1754 b = PyNumber_Multiply(a, a);
1755 if (b == NULL) {
1756 goto error;
1757 }
1758 a_too_large = PyObject_RichCompareBool(n, b, Py_LT);
1759 Py_DECREF(b);
1760 if (a_too_large == -1) {
1761 goto error;
1762 }
1763
1764 if (a_too_large) {
1765 Py_SETREF(a, PyNumber_Subtract(a, _PyLong_One));
1766 }
1767 Py_DECREF(n);
1768 return a;
1769
1770 error:
1771 Py_XDECREF(a);
1772 Py_DECREF(n);
1773 return NULL;
1774}
1775
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001776/* Divide-and-conquer factorial algorithm
1777 *
Raymond Hettinger15f44ab2016-08-30 10:47:49 -07001778 * Based on the formula and pseudo-code provided at:
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001779 * http://www.luschny.de/math/factorial/binarysplitfact.html
1780 *
1781 * Faster algorithms exist, but they're more complicated and depend on
Ezio Melotti9527afd2010-07-08 15:03:02 +00001782 * a fast prime factorization algorithm.
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001783 *
1784 * Notes on the algorithm
1785 * ----------------------
1786 *
1787 * factorial(n) is written in the form 2**k * m, with m odd. k and m are
1788 * computed separately, and then combined using a left shift.
1789 *
1790 * The function factorial_odd_part computes the odd part m (i.e., the greatest
1791 * odd divisor) of factorial(n), using the formula:
1792 *
1793 * factorial_odd_part(n) =
1794 *
1795 * product_{i >= 0} product_{0 < j <= n / 2**i, j odd} j
1796 *
1797 * Example: factorial_odd_part(20) =
1798 *
1799 * (1) *
1800 * (1) *
1801 * (1 * 3 * 5) *
1802 * (1 * 3 * 5 * 7 * 9)
1803 * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1804 *
1805 * Here i goes from large to small: the first term corresponds to i=4 (any
1806 * larger i gives an empty product), and the last term corresponds to i=0.
1807 * Each term can be computed from the last by multiplying by the extra odd
1808 * numbers required: e.g., to get from the penultimate term to the last one,
1809 * we multiply by (11 * 13 * 15 * 17 * 19).
1810 *
1811 * To see a hint of why this formula works, here are the same numbers as above
1812 * but with the even parts (i.e., the appropriate powers of 2) included. For
1813 * each subterm in the product for i, we multiply that subterm by 2**i:
1814 *
1815 * factorial(20) =
1816 *
1817 * (16) *
1818 * (8) *
1819 * (4 * 12 * 20) *
1820 * (2 * 6 * 10 * 14 * 18) *
1821 * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1822 *
1823 * The factorial_partial_product function computes the product of all odd j in
1824 * range(start, stop) for given start and stop. It's used to compute the
1825 * partial products like (11 * 13 * 15 * 17 * 19) in the example above. It
1826 * operates recursively, repeatedly splitting the range into two roughly equal
1827 * pieces until the subranges are small enough to be computed using only C
1828 * integer arithmetic.
1829 *
1830 * The two-valuation k (i.e., the exponent of the largest power of 2 dividing
1831 * the factorial) is computed independently in the main math_factorial
1832 * function. By standard results, its value is:
1833 *
1834 * two_valuation = n//2 + n//4 + n//8 + ....
1835 *
1836 * It can be shown (e.g., by complete induction on n) that two_valuation is
1837 * equal to n - count_set_bits(n), where count_set_bits(n) gives the number of
1838 * '1'-bits in the binary expansion of n.
1839 */
1840
1841/* factorial_partial_product: Compute product(range(start, stop, 2)) using
1842 * divide and conquer. Assumes start and stop are odd and stop > start.
1843 * max_bits must be >= bit_length(stop - 2). */
1844
1845static PyObject *
1846factorial_partial_product(unsigned long start, unsigned long stop,
1847 unsigned long max_bits)
1848{
1849 unsigned long midpoint, num_operands;
1850 PyObject *left = NULL, *right = NULL, *result = NULL;
1851
1852 /* If the return value will fit an unsigned long, then we can
1853 * multiply in a tight, fast loop where each multiply is O(1).
1854 * Compute an upper bound on the number of bits required to store
1855 * the answer.
1856 *
1857 * Storing some integer z requires floor(lg(z))+1 bits, which is
1858 * conveniently the value returned by bit_length(z). The
1859 * product x*y will require at most
1860 * bit_length(x) + bit_length(y) bits to store, based
1861 * on the idea that lg product = lg x + lg y.
1862 *
1863 * We know that stop - 2 is the largest number to be multiplied. From
1864 * there, we have: bit_length(answer) <= num_operands *
1865 * bit_length(stop - 2)
1866 */
1867
1868 num_operands = (stop - start) / 2;
1869 /* The "num_operands <= 8 * SIZEOF_LONG" check guards against the
1870 * unlikely case of an overflow in num_operands * max_bits. */
1871 if (num_operands <= 8 * SIZEOF_LONG &&
1872 num_operands * max_bits <= 8 * SIZEOF_LONG) {
1873 unsigned long j, total;
1874 for (total = start, j = start + 2; j < stop; j += 2)
1875 total *= j;
1876 return PyLong_FromUnsignedLong(total);
1877 }
1878
1879 /* find midpoint of range(start, stop), rounded up to next odd number. */
1880 midpoint = (start + num_operands) | 1;
1881 left = factorial_partial_product(start, midpoint,
1882 bit_length(midpoint - 2));
1883 if (left == NULL)
1884 goto error;
1885 right = factorial_partial_product(midpoint, stop, max_bits);
1886 if (right == NULL)
1887 goto error;
1888 result = PyNumber_Multiply(left, right);
1889
1890 error:
1891 Py_XDECREF(left);
1892 Py_XDECREF(right);
1893 return result;
1894}
1895
1896/* factorial_odd_part: compute the odd part of factorial(n). */
1897
1898static PyObject *
1899factorial_odd_part(unsigned long n)
1900{
1901 long i;
1902 unsigned long v, lower, upper;
1903 PyObject *partial, *tmp, *inner, *outer;
1904
1905 inner = PyLong_FromLong(1);
1906 if (inner == NULL)
1907 return NULL;
1908 outer = inner;
1909 Py_INCREF(outer);
1910
1911 upper = 3;
1912 for (i = bit_length(n) - 2; i >= 0; i--) {
1913 v = n >> i;
1914 if (v <= 2)
1915 continue;
1916 lower = upper;
1917 /* (v + 1) | 1 = least odd integer strictly larger than n / 2**i */
1918 upper = (v + 1) | 1;
1919 /* Here inner is the product of all odd integers j in the range (0,
1920 n/2**(i+1)]. The factorial_partial_product call below gives the
1921 product of all odd integers j in the range (n/2**(i+1), n/2**i]. */
1922 partial = factorial_partial_product(lower, upper, bit_length(upper-2));
1923 /* inner *= partial */
1924 if (partial == NULL)
1925 goto error;
1926 tmp = PyNumber_Multiply(inner, partial);
1927 Py_DECREF(partial);
1928 if (tmp == NULL)
1929 goto error;
1930 Py_DECREF(inner);
1931 inner = tmp;
1932 /* Now inner is the product of all odd integers j in the range (0,
1933 n/2**i], giving the inner product in the formula above. */
1934
1935 /* outer *= inner; */
1936 tmp = PyNumber_Multiply(outer, inner);
1937 if (tmp == NULL)
1938 goto error;
1939 Py_DECREF(outer);
1940 outer = tmp;
1941 }
Mark Dickinson76464492012-10-25 10:46:28 +01001942 Py_DECREF(inner);
1943 return outer;
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001944
1945 error:
1946 Py_DECREF(outer);
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001947 Py_DECREF(inner);
Mark Dickinson76464492012-10-25 10:46:28 +01001948 return NULL;
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001949}
1950
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001951
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001952/* Lookup table for small factorial values */
1953
1954static const unsigned long SmallFactorials[] = {
1955 1, 1, 2, 6, 24, 120, 720, 5040, 40320,
1956 362880, 3628800, 39916800, 479001600,
1957#if SIZEOF_LONG >= 8
1958 6227020800, 87178291200, 1307674368000,
1959 20922789888000, 355687428096000, 6402373705728000,
1960 121645100408832000, 2432902008176640000
1961#endif
1962};
1963
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001964/*[clinic input]
1965math.factorial
1966
1967 x as arg: object
1968 /
1969
1970Find x!.
1971
1972Raise a ValueError if x is negative or non-integral.
1973[clinic start generated code]*/
1974
Barry Warsaw8b43b191996-12-09 22:32:36 +00001975static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02001976math_factorial(PyObject *module, PyObject *arg)
1977/*[clinic end generated code: output=6686f26fae00e9ca input=6d1c8105c0d91fb4]*/
Georg Brandlc28e1fa2008-06-10 19:20:26 +00001978{
Serhiy Storchakaa5119e72019-05-19 14:14:38 +03001979 long x, two_valuation;
Mark Dickinson5990d282014-04-10 09:29:39 -04001980 int overflow;
Serhiy Storchakaa5119e72019-05-19 14:14:38 +03001981 PyObject *result, *odd_part, *pyint_form;
Georg Brandlc28e1fa2008-06-10 19:20:26 +00001982
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001983 if (PyFloat_Check(arg)) {
Serhiy Storchaka231aad32019-06-17 16:57:27 +03001984 if (PyErr_WarnEx(PyExc_DeprecationWarning,
1985 "Using factorial() with floats is deprecated",
1986 1) < 0)
1987 {
1988 return NULL;
1989 }
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001990 PyObject *lx;
1991 double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);
1992 if (!(Py_IS_FINITE(dx) && dx == floor(dx))) {
1993 PyErr_SetString(PyExc_ValueError,
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00001994 "factorial() only accepts integral values");
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00001995 return NULL;
1996 }
1997 lx = PyLong_FromDouble(dx);
1998 if (lx == NULL)
1999 return NULL;
Mark Dickinson5990d282014-04-10 09:29:39 -04002000 x = PyLong_AsLongAndOverflow(lx, &overflow);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002001 Py_DECREF(lx);
2002 }
Pablo Galindoe9ba3702018-09-03 22:20:06 +01002003 else {
2004 pyint_form = PyNumber_Index(arg);
2005 if (pyint_form == NULL) {
2006 return NULL;
2007 }
2008 x = PyLong_AsLongAndOverflow(pyint_form, &overflow);
2009 Py_DECREF(pyint_form);
2010 }
Georg Brandlc28e1fa2008-06-10 19:20:26 +00002011
Mark Dickinson5990d282014-04-10 09:29:39 -04002012 if (x == -1 && PyErr_Occurred()) {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002013 return NULL;
Mark Dickinson5990d282014-04-10 09:29:39 -04002014 }
2015 else if (overflow == 1) {
2016 PyErr_Format(PyExc_OverflowError,
2017 "factorial() argument should not exceed %ld",
2018 LONG_MAX);
2019 return NULL;
2020 }
2021 else if (overflow == -1 || x < 0) {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002022 PyErr_SetString(PyExc_ValueError,
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00002023 "factorial() not defined for negative values");
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002024 return NULL;
2025 }
Georg Brandlc28e1fa2008-06-10 19:20:26 +00002026
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00002027 /* use lookup table if x is small */
Victor Stinner63941882011-09-29 00:42:28 +02002028 if (x < (long)Py_ARRAY_LENGTH(SmallFactorials))
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00002029 return PyLong_FromUnsignedLong(SmallFactorials[x]);
2030
2031 /* else express in the form odd_part * 2**two_valuation, and compute as
2032 odd_part << two_valuation. */
2033 odd_part = factorial_odd_part(x);
2034 if (odd_part == NULL)
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002035 return NULL;
Serhiy Storchakaa5119e72019-05-19 14:14:38 +03002036 two_valuation = x - count_set_bits(x);
2037 result = _PyLong_Lshift(odd_part, two_valuation);
Mark Dickinson4c8a9a22010-05-15 17:02:38 +00002038 Py_DECREF(odd_part);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002039 return result;
Georg Brandlc28e1fa2008-06-10 19:20:26 +00002040}
2041
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002042
2043/*[clinic input]
2044math.trunc
2045
2046 x: object
2047 /
2048
2049Truncates the Real x to the nearest Integral toward 0.
2050
2051Uses the __trunc__ magic method.
2052[clinic start generated code]*/
Georg Brandlc28e1fa2008-06-10 19:20:26 +00002053
2054static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002055math_trunc(PyObject *module, PyObject *x)
2056/*[clinic end generated code: output=34b9697b707e1031 input=2168b34e0a09134d]*/
Christian Heimes400adb02008-02-01 08:12:03 +00002057{
Benjamin Petersonce798522012-01-22 11:24:29 -05002058 _Py_IDENTIFIER(__trunc__);
Benjamin Petersonb0125892010-07-02 13:35:17 +00002059 PyObject *trunc, *result;
Christian Heimes400adb02008-02-01 08:12:03 +00002060
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002061 if (Py_TYPE(x)->tp_dict == NULL) {
2062 if (PyType_Ready(Py_TYPE(x)) < 0)
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002063 return NULL;
2064 }
Christian Heimes400adb02008-02-01 08:12:03 +00002065
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002066 trunc = _PyObject_LookupSpecial(x, &PyId___trunc__);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002067 if (trunc == NULL) {
Benjamin Peterson8bb9cde2010-07-01 15:16:55 +00002068 if (!PyErr_Occurred())
2069 PyErr_Format(PyExc_TypeError,
2070 "type %.100s doesn't define __trunc__ method",
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002071 Py_TYPE(x)->tp_name);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002072 return NULL;
2073 }
Victor Stinnerf17c3de2016-12-06 18:46:19 +01002074 result = _PyObject_CallNoArg(trunc);
Benjamin Petersonb0125892010-07-02 13:35:17 +00002075 Py_DECREF(trunc);
2076 return result;
Christian Heimes400adb02008-02-01 08:12:03 +00002077}
2078
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002079
2080/*[clinic input]
2081math.frexp
2082
2083 x: double
2084 /
2085
2086Return the mantissa and exponent of x, as pair (m, e).
2087
2088m is a float and e is an int, such that x = m * 2.**e.
2089If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0.
2090[clinic start generated code]*/
Christian Heimes400adb02008-02-01 08:12:03 +00002091
2092static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002093math_frexp_impl(PyObject *module, double x)
2094/*[clinic end generated code: output=03e30d252a15ad4a input=96251c9e208bc6e9]*/
Guido van Rossumd18ad581991-10-24 14:57:21 +00002095{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002096 int i;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002097 /* deal with special cases directly, to sidestep platform
2098 differences */
2099 if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
2100 i = 0;
2101 }
2102 else {
2103 PyFPE_START_PROTECT("in math_frexp", return 0);
2104 x = frexp(x, &i);
2105 PyFPE_END_PROTECT(x);
2106 }
2107 return Py_BuildValue("(di)", x, i);
Guido van Rossumd18ad581991-10-24 14:57:21 +00002108}
2109
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002110
2111/*[clinic input]
2112math.ldexp
2113
2114 x: double
2115 i: object
2116 /
2117
2118Return x * (2**i).
2119
2120This is essentially the inverse of frexp().
2121[clinic start generated code]*/
Guido van Rossumc6e22901998-12-04 19:26:43 +00002122
Barry Warsaw8b43b191996-12-09 22:32:36 +00002123static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002124math_ldexp_impl(PyObject *module, double x, PyObject *i)
2125/*[clinic end generated code: output=b6892f3c2df9cc6a input=17d5970c1a40a8c1]*/
Guido van Rossumd18ad581991-10-24 14:57:21 +00002126{
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002127 double r;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002128 long exp;
2129 int overflow;
Alexandre Vassalotti6461e102008-05-15 22:09:29 +00002130
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002131 if (PyLong_Check(i)) {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002132 /* on overflow, replace exponent with either LONG_MAX
2133 or LONG_MIN, depending on the sign. */
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002134 exp = PyLong_AsLongAndOverflow(i, &overflow);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002135 if (exp == -1 && PyErr_Occurred())
2136 return NULL;
2137 if (overflow)
2138 exp = overflow < 0 ? LONG_MIN : LONG_MAX;
2139 }
2140 else {
2141 PyErr_SetString(PyExc_TypeError,
Serhiy Storchaka95949422013-08-27 19:40:23 +03002142 "Expected an int as second argument to ldexp.");
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002143 return NULL;
2144 }
Alexandre Vassalotti6461e102008-05-15 22:09:29 +00002145
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002146 if (x == 0. || !Py_IS_FINITE(x)) {
2147 /* NaNs, zeros and infinities are returned unchanged */
2148 r = x;
2149 errno = 0;
2150 } else if (exp > INT_MAX) {
2151 /* overflow */
2152 r = copysign(Py_HUGE_VAL, x);
2153 errno = ERANGE;
2154 } else if (exp < INT_MIN) {
2155 /* underflow to +-0 */
2156 r = copysign(0., x);
2157 errno = 0;
2158 } else {
2159 errno = 0;
2160 PyFPE_START_PROTECT("in math_ldexp", return 0);
2161 r = ldexp(x, (int)exp);
2162 PyFPE_END_PROTECT(r);
2163 if (Py_IS_INFINITY(r))
2164 errno = ERANGE;
2165 }
Alexandre Vassalotti6461e102008-05-15 22:09:29 +00002166
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002167 if (errno && is_error(r))
2168 return NULL;
2169 return PyFloat_FromDouble(r);
Guido van Rossumd18ad581991-10-24 14:57:21 +00002170}
2171
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002172
2173/*[clinic input]
2174math.modf
2175
2176 x: double
2177 /
2178
2179Return the fractional and integer parts of x.
2180
2181Both results carry the sign of x and are floats.
2182[clinic start generated code]*/
Guido van Rossumc6e22901998-12-04 19:26:43 +00002183
Barry Warsaw8b43b191996-12-09 22:32:36 +00002184static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002185math_modf_impl(PyObject *module, double x)
2186/*[clinic end generated code: output=90cee0260014c3c0 input=b4cfb6786afd9035]*/
Guido van Rossumd18ad581991-10-24 14:57:21 +00002187{
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002188 double y;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002189 /* some platforms don't do the right thing for NaNs and
2190 infinities, so we take care of special cases directly. */
2191 if (!Py_IS_FINITE(x)) {
2192 if (Py_IS_INFINITY(x))
2193 return Py_BuildValue("(dd)", copysign(0., x), x);
2194 else if (Py_IS_NAN(x))
2195 return Py_BuildValue("(dd)", x, x);
2196 }
Christian Heimesa342c012008-04-20 21:01:16 +00002197
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002198 errno = 0;
2199 PyFPE_START_PROTECT("in math_modf", return 0);
2200 x = modf(x, &y);
2201 PyFPE_END_PROTECT(x);
2202 return Py_BuildValue("(dd)", x, y);
Guido van Rossumd18ad581991-10-24 14:57:21 +00002203}
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00002204
Guido van Rossumc6e22901998-12-04 19:26:43 +00002205
Serhiy Storchaka95949422013-08-27 19:40:23 +03002206/* A decent logarithm is easy to compute even for huge ints, but libm can't
Tim Peters78526162001-09-05 00:53:45 +00002207 do that by itself -- loghelper can. func is log or log10, and name is
Serhiy Storchaka95949422013-08-27 19:40:23 +03002208 "log" or "log10". Note that overflow of the result isn't possible: an int
Mark Dickinson6ecd9e52010-01-02 15:33:56 +00002209 can contain no more than INT_MAX * SHIFT bits, so has value certainly less
2210 than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
Tim Peters78526162001-09-05 00:53:45 +00002211 small enough to fit in an IEEE single. log and log10 are even smaller.
Serhiy Storchaka95949422013-08-27 19:40:23 +03002212 However, intermediate overflow is possible for an int if the number of bits
2213 in that int is larger than PY_SSIZE_T_MAX. */
Tim Peters78526162001-09-05 00:53:45 +00002214
2215static PyObject*
Serhiy Storchakaef1585e2015-12-25 20:01:53 +02002216loghelper(PyObject* arg, double (*func)(double), const char *funcname)
Tim Peters78526162001-09-05 00:53:45 +00002217{
Serhiy Storchaka95949422013-08-27 19:40:23 +03002218 /* If it is int, do it ourselves. */
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002219 if (PyLong_Check(arg)) {
Mark Dickinsonc6037172010-09-29 19:06:36 +00002220 double x, result;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002221 Py_ssize_t e;
Mark Dickinsonc6037172010-09-29 19:06:36 +00002222
2223 /* Negative or zero inputs give a ValueError. */
2224 if (Py_SIZE(arg) <= 0) {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002225 PyErr_SetString(PyExc_ValueError,
2226 "math domain error");
2227 return NULL;
2228 }
Mark Dickinsonfa41e602010-09-28 07:22:27 +00002229
Mark Dickinsonc6037172010-09-29 19:06:36 +00002230 x = PyLong_AsDouble(arg);
2231 if (x == -1.0 && PyErr_Occurred()) {
2232 if (!PyErr_ExceptionMatches(PyExc_OverflowError))
2233 return NULL;
2234 /* Here the conversion to double overflowed, but it's possible
2235 to compute the log anyway. Clear the exception and continue. */
2236 PyErr_Clear();
2237 x = _PyLong_Frexp((PyLongObject *)arg, &e);
2238 if (x == -1.0 && PyErr_Occurred())
2239 return NULL;
2240 /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */
2241 result = func(x) + func(2.0) * e;
2242 }
2243 else
2244 /* Successfully converted x to a double. */
2245 result = func(x);
2246 return PyFloat_FromDouble(result);
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002247 }
Tim Peters78526162001-09-05 00:53:45 +00002248
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002249 /* Else let libm handle it by itself. */
2250 return math_1(arg, func, 0);
Tim Peters78526162001-09-05 00:53:45 +00002251}
2252
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002253
2254/*[clinic input]
2255math.log
2256
2257 x: object
2258 [
2259 base: object(c_default="NULL") = math.e
2260 ]
2261 /
2262
2263Return the logarithm of x to the given base.
2264
2265If the base not specified, returns the natural logarithm (base e) of x.
2266[clinic start generated code]*/
2267
Tim Peters78526162001-09-05 00:53:45 +00002268static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002269math_log_impl(PyObject *module, PyObject *x, int group_right_1,
2270 PyObject *base)
2271/*[clinic end generated code: output=7b5a39e526b73fc9 input=0f62d5726cbfebbd]*/
Tim Peters78526162001-09-05 00:53:45 +00002272{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002273 PyObject *num, *den;
2274 PyObject *ans;
Raymond Hettinger866964c2002-12-14 19:51:34 +00002275
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002276 num = loghelper(x, m_log, "log");
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002277 if (num == NULL || base == NULL)
2278 return num;
Raymond Hettinger866964c2002-12-14 19:51:34 +00002279
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002280 den = loghelper(base, m_log, "log");
2281 if (den == NULL) {
2282 Py_DECREF(num);
2283 return NULL;
2284 }
Raymond Hettinger866964c2002-12-14 19:51:34 +00002285
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002286 ans = PyNumber_TrueDivide(num, den);
2287 Py_DECREF(num);
2288 Py_DECREF(den);
2289 return ans;
Tim Peters78526162001-09-05 00:53:45 +00002290}
2291
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002292
2293/*[clinic input]
2294math.log2
2295
2296 x: object
2297 /
2298
2299Return the base 2 logarithm of x.
2300[clinic start generated code]*/
Tim Peters78526162001-09-05 00:53:45 +00002301
2302static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002303math_log2(PyObject *module, PyObject *x)
2304/*[clinic end generated code: output=5425899a4d5d6acb input=08321262bae4f39b]*/
Victor Stinnerfa0e3d52011-05-09 01:01:09 +02002305{
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002306 return loghelper(x, m_log2, "log2");
Victor Stinnerfa0e3d52011-05-09 01:01:09 +02002307}
2308
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002309
2310/*[clinic input]
2311math.log10
2312
2313 x: object
2314 /
2315
2316Return the base 10 logarithm of x.
2317[clinic start generated code]*/
Victor Stinnerfa0e3d52011-05-09 01:01:09 +02002318
2319static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002320math_log10(PyObject *module, PyObject *x)
2321/*[clinic end generated code: output=be72a64617df9c6f input=b2469d02c6469e53]*/
Tim Peters78526162001-09-05 00:53:45 +00002322{
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002323 return loghelper(x, m_log10, "log10");
Tim Peters78526162001-09-05 00:53:45 +00002324}
2325
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002326
2327/*[clinic input]
2328math.fmod
2329
2330 x: double
2331 y: double
2332 /
2333
2334Return fmod(x, y), according to platform C.
2335
2336x % y may differ.
2337[clinic start generated code]*/
Tim Peters78526162001-09-05 00:53:45 +00002338
Christian Heimes53876d92008-04-19 00:31:39 +00002339static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002340math_fmod_impl(PyObject *module, double x, double y)
2341/*[clinic end generated code: output=7559d794343a27b5 input=4f84caa8cfc26a03]*/
Christian Heimes53876d92008-04-19 00:31:39 +00002342{
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002343 double r;
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002344 /* fmod(x, +/-Inf) returns x for finite x. */
2345 if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
2346 return PyFloat_FromDouble(x);
2347 errno = 0;
2348 PyFPE_START_PROTECT("in math_fmod", return 0);
2349 r = fmod(x, y);
2350 PyFPE_END_PROTECT(r);
2351 if (Py_IS_NAN(r)) {
2352 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
2353 errno = EDOM;
2354 else
2355 errno = 0;
2356 }
2357 if (errno && is_error(r))
2358 return NULL;
2359 else
2360 return PyFloat_FromDouble(r);
Christian Heimes53876d92008-04-19 00:31:39 +00002361}
2362
Raymond Hettinger13990742018-08-11 11:26:36 -07002363/*
Raymond Hettinger745c0f32018-08-31 11:22:13 -07002364Given an *n* length *vec* of values and a value *max*, compute:
Raymond Hettinger13990742018-08-11 11:26:36 -07002365
Raymond Hettingerc630e102018-08-11 18:39:05 -07002366 max * sqrt(sum((x / max) ** 2 for x in vec))
Raymond Hettinger13990742018-08-11 11:26:36 -07002367
Raymond Hettinger745c0f32018-08-31 11:22:13 -07002368The value of the *max* variable must be non-negative and
Raymond Hettinger216aaaa2018-11-09 01:06:02 -08002369equal to the absolute value of the largest magnitude
Raymond Hettinger745c0f32018-08-31 11:22:13 -07002370entry in the vector. If n==0, then *max* should be 0.0.
2371If an infinity is present in the vec, *max* should be INF.
Raymond Hettingerc630e102018-08-11 18:39:05 -07002372
2373The *found_nan* variable indicates whether some member of
2374the *vec* is a NaN.
Raymond Hettinger21786f52018-08-28 22:47:24 -07002375
2376To improve accuracy and to increase the number of cases where
2377vector_norm() is commutative, we use a variant of Neumaier
2378summation specialized to exploit that we always know that
2379|csum| >= |x|.
2380
2381The *csum* variable tracks the cumulative sum and *frac* tracks
2382the cumulative fractional errors at each step. Since this
2383variant assumes that |csum| >= |x| at each step, we establish
2384the precondition by starting the accumulation from 1.0 which
Raymond Hettinger745c0f32018-08-31 11:22:13 -07002385represents the largest possible value of (x/max)**2.
2386
2387After the loop is finished, the initial 1.0 is subtracted out
2388for a net zero effect on the final sum. Since *csum* will be
2389greater than 1.0, the subtraction of 1.0 will not cause
2390fractional digits to be dropped from *csum*.
Raymond Hettinger21786f52018-08-28 22:47:24 -07002391
Raymond Hettinger13990742018-08-11 11:26:36 -07002392*/
2393
2394static inline double
Raymond Hettingerc630e102018-08-11 18:39:05 -07002395vector_norm(Py_ssize_t n, double *vec, double max, int found_nan)
Raymond Hettinger13990742018-08-11 11:26:36 -07002396{
Raymond Hettinger745c0f32018-08-31 11:22:13 -07002397 double x, csum = 1.0, oldcsum, frac = 0.0;
Raymond Hettinger13990742018-08-11 11:26:36 -07002398 Py_ssize_t i;
2399
Raymond Hettingerc630e102018-08-11 18:39:05 -07002400 if (Py_IS_INFINITY(max)) {
2401 return max;
2402 }
2403 if (found_nan) {
2404 return Py_NAN;
2405 }
Raymond Hettingerf3267142018-09-02 13:34:21 -07002406 if (max == 0.0 || n <= 1) {
Raymond Hettinger745c0f32018-08-31 11:22:13 -07002407 return max;
Raymond Hettinger13990742018-08-11 11:26:36 -07002408 }
Raymond Hettinger745c0f32018-08-31 11:22:13 -07002409 for (i=0 ; i < n ; i++) {
Raymond Hettinger13990742018-08-11 11:26:36 -07002410 x = vec[i];
Raymond Hettinger745c0f32018-08-31 11:22:13 -07002411 assert(Py_IS_FINITE(x) && fabs(x) <= max);
Raymond Hettinger13990742018-08-11 11:26:36 -07002412 x /= max;
Raymond Hettinger21786f52018-08-28 22:47:24 -07002413 x = x*x;
Raymond Hettinger13990742018-08-11 11:26:36 -07002414 oldcsum = csum;
2415 csum += x;
Raymond Hettinger745c0f32018-08-31 11:22:13 -07002416 assert(csum >= x);
Raymond Hettinger21786f52018-08-28 22:47:24 -07002417 frac += (oldcsum - csum) + x;
Raymond Hettinger13990742018-08-11 11:26:36 -07002418 }
Raymond Hettinger745c0f32018-08-31 11:22:13 -07002419 return max * sqrt(csum - 1.0 + frac);
Raymond Hettinger13990742018-08-11 11:26:36 -07002420}
2421
Raymond Hettingerc630e102018-08-11 18:39:05 -07002422#define NUM_STACK_ELEMS 16
2423
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002424/*[clinic input]
2425math.dist
2426
Ammar Askarcb08a712019-01-12 01:23:41 -05002427 p: object(subclass_of='&PyTuple_Type')
2428 q: object(subclass_of='&PyTuple_Type')
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002429 /
2430
2431Return the Euclidean distance between two points p and q.
2432
2433The points should be specified as tuples of coordinates.
2434Both tuples must be the same size.
2435
2436Roughly equivalent to:
2437 sqrt(sum((px - qx) ** 2.0 for px, qx in zip(p, q)))
2438[clinic start generated code]*/
2439
2440static PyObject *
2441math_dist_impl(PyObject *module, PyObject *p, PyObject *q)
Ammar Askarcb08a712019-01-12 01:23:41 -05002442/*[clinic end generated code: output=56bd9538d06bbcfe input=937122eaa5f19272]*/
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002443{
2444 PyObject *item;
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002445 double max = 0.0;
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002446 double x, px, qx, result;
2447 Py_ssize_t i, m, n;
2448 int found_nan = 0;
Raymond Hettingerc630e102018-08-11 18:39:05 -07002449 double diffs_on_stack[NUM_STACK_ELEMS];
2450 double *diffs = diffs_on_stack;
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002451
2452 m = PyTuple_GET_SIZE(p);
2453 n = PyTuple_GET_SIZE(q);
2454 if (m != n) {
2455 PyErr_SetString(PyExc_ValueError,
2456 "both points must have the same number of dimensions");
2457 return NULL;
2458
2459 }
Raymond Hettingerc630e102018-08-11 18:39:05 -07002460 if (n > NUM_STACK_ELEMS) {
2461 diffs = (double *) PyObject_Malloc(n * sizeof(double));
2462 if (diffs == NULL) {
Zackery Spytz4c49da02018-12-07 03:11:30 -07002463 return PyErr_NoMemory();
Raymond Hettingerc630e102018-08-11 18:39:05 -07002464 }
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002465 }
2466 for (i=0 ; i<n ; i++) {
2467 item = PyTuple_GET_ITEM(p, i);
Raymond Hettingercfd735e2019-01-29 20:39:53 -08002468 ASSIGN_DOUBLE(px, item, error_exit);
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002469 item = PyTuple_GET_ITEM(q, i);
Raymond Hettingercfd735e2019-01-29 20:39:53 -08002470 ASSIGN_DOUBLE(qx, item, error_exit);
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002471 x = fabs(px - qx);
2472 diffs[i] = x;
2473 found_nan |= Py_IS_NAN(x);
2474 if (x > max) {
2475 max = x;
2476 }
2477 }
Raymond Hettingerc630e102018-08-11 18:39:05 -07002478 result = vector_norm(n, diffs, max, found_nan);
2479 if (diffs != diffs_on_stack) {
2480 PyObject_Free(diffs);
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002481 }
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002482 return PyFloat_FromDouble(result);
Raymond Hettingerc630e102018-08-11 18:39:05 -07002483
2484 error_exit:
2485 if (diffs != diffs_on_stack) {
2486 PyObject_Free(diffs);
2487 }
2488 return NULL;
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07002489}
2490
Raymond Hettingerc6dabe32018-07-28 07:48:04 -07002491/* AC: cannot convert yet, waiting for *args support */
Christian Heimes53876d92008-04-19 00:31:39 +00002492static PyObject *
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02002493math_hypot(PyObject *self, PyObject *const *args, Py_ssize_t nargs)
Christian Heimes53876d92008-04-19 00:31:39 +00002494{
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02002495 Py_ssize_t i;
Raymond Hettingerc6dabe32018-07-28 07:48:04 -07002496 PyObject *item;
Raymond Hettingerc6dabe32018-07-28 07:48:04 -07002497 double max = 0.0;
Raymond Hettingerc6dabe32018-07-28 07:48:04 -07002498 double x, result;
2499 int found_nan = 0;
Raymond Hettingerc630e102018-08-11 18:39:05 -07002500 double coord_on_stack[NUM_STACK_ELEMS];
2501 double *coordinates = coord_on_stack;
Raymond Hettingerc6dabe32018-07-28 07:48:04 -07002502
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02002503 if (nargs > NUM_STACK_ELEMS) {
2504 coordinates = (double *) PyObject_Malloc(nargs * sizeof(double));
Zackery Spytz4c49da02018-12-07 03:11:30 -07002505 if (coordinates == NULL) {
2506 return PyErr_NoMemory();
2507 }
Raymond Hettingerc630e102018-08-11 18:39:05 -07002508 }
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02002509 for (i = 0; i < nargs; i++) {
2510 item = args[i];
Raymond Hettingercfd735e2019-01-29 20:39:53 -08002511 ASSIGN_DOUBLE(x, item, error_exit);
Raymond Hettingerc6dabe32018-07-28 07:48:04 -07002512 x = fabs(x);
2513 coordinates[i] = x;
2514 found_nan |= Py_IS_NAN(x);
2515 if (x > max) {
2516 max = x;
2517 }
2518 }
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02002519 result = vector_norm(nargs, coordinates, max, found_nan);
Raymond Hettingerc630e102018-08-11 18:39:05 -07002520 if (coordinates != coord_on_stack) {
2521 PyObject_Free(coordinates);
Raymond Hettingerc6dabe32018-07-28 07:48:04 -07002522 }
Raymond Hettingerc6dabe32018-07-28 07:48:04 -07002523 return PyFloat_FromDouble(result);
Raymond Hettingerc630e102018-08-11 18:39:05 -07002524
2525 error_exit:
2526 if (coordinates != coord_on_stack) {
2527 PyObject_Free(coordinates);
2528 }
2529 return NULL;
Christian Heimes53876d92008-04-19 00:31:39 +00002530}
2531
Raymond Hettingerc630e102018-08-11 18:39:05 -07002532#undef NUM_STACK_ELEMS
2533
Raymond Hettingerc6dabe32018-07-28 07:48:04 -07002534PyDoc_STRVAR(math_hypot_doc,
2535 "hypot(*coordinates) -> value\n\n\
2536Multidimensional Euclidean distance from the origin to a point.\n\
2537\n\
2538Roughly equivalent to:\n\
2539 sqrt(sum(x**2 for x in coordinates))\n\
2540\n\
2541For a two dimensional point (x, y), gives the hypotenuse\n\
2542using the Pythagorean theorem: sqrt(x*x + y*y).\n\
2543\n\
2544For example, the hypotenuse of a 3/4/5 right triangle is:\n\
2545\n\
2546 >>> hypot(3.0, 4.0)\n\
2547 5.0\n\
2548");
Christian Heimes53876d92008-04-19 00:31:39 +00002549
2550/* pow can't use math_2, but needs its own wrapper: the problem is
2551 that an infinite result can arise either as a result of overflow
2552 (in which case OverflowError should be raised) or as a result of
2553 e.g. 0.**-5. (for which ValueError needs to be raised.)
2554*/
2555
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002556/*[clinic input]
2557math.pow
Christian Heimes53876d92008-04-19 00:31:39 +00002558
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002559 x: double
2560 y: double
2561 /
2562
2563Return x**y (x to the power of y).
2564[clinic start generated code]*/
2565
2566static PyObject *
2567math_pow_impl(PyObject *module, double x, double y)
2568/*[clinic end generated code: output=fff93e65abccd6b0 input=c26f1f6075088bfd]*/
2569{
2570 double r;
2571 int odd_y;
Christian Heimesa342c012008-04-20 21:01:16 +00002572
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002573 /* deal directly with IEEE specials, to cope with problems on various
2574 platforms whose semantics don't exactly match C99 */
2575 r = 0.; /* silence compiler warning */
2576 if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
2577 errno = 0;
2578 if (Py_IS_NAN(x))
2579 r = y == 0. ? 1. : x; /* NaN**0 = 1 */
2580 else if (Py_IS_NAN(y))
2581 r = x == 1. ? 1. : y; /* 1**NaN = 1 */
2582 else if (Py_IS_INFINITY(x)) {
2583 odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
2584 if (y > 0.)
2585 r = odd_y ? x : fabs(x);
2586 else if (y == 0.)
2587 r = 1.;
2588 else /* y < 0. */
2589 r = odd_y ? copysign(0., x) : 0.;
2590 }
2591 else if (Py_IS_INFINITY(y)) {
2592 if (fabs(x) == 1.0)
2593 r = 1.;
2594 else if (y > 0. && fabs(x) > 1.0)
2595 r = y;
2596 else if (y < 0. && fabs(x) < 1.0) {
2597 r = -y; /* result is +inf */
2598 if (x == 0.) /* 0**-inf: divide-by-zero */
2599 errno = EDOM;
2600 }
2601 else
2602 r = 0.;
2603 }
2604 }
2605 else {
2606 /* let libm handle finite**finite */
2607 errno = 0;
2608 PyFPE_START_PROTECT("in math_pow", return 0);
2609 r = pow(x, y);
2610 PyFPE_END_PROTECT(r);
2611 /* a NaN result should arise only from (-ve)**(finite
2612 non-integer); in this case we want to raise ValueError. */
2613 if (!Py_IS_FINITE(r)) {
2614 if (Py_IS_NAN(r)) {
2615 errno = EDOM;
2616 }
2617 /*
2618 an infinite result here arises either from:
2619 (A) (+/-0.)**negative (-> divide-by-zero)
2620 (B) overflow of x**y with x and y finite
2621 */
2622 else if (Py_IS_INFINITY(r)) {
2623 if (x == 0.)
2624 errno = EDOM;
2625 else
2626 errno = ERANGE;
2627 }
2628 }
2629 }
Christian Heimes53876d92008-04-19 00:31:39 +00002630
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002631 if (errno && is_error(r))
2632 return NULL;
2633 else
2634 return PyFloat_FromDouble(r);
Christian Heimes53876d92008-04-19 00:31:39 +00002635}
2636
Christian Heimes53876d92008-04-19 00:31:39 +00002637
Christian Heimes072c0f12008-01-03 23:01:04 +00002638static const double degToRad = Py_MATH_PI / 180.0;
2639static const double radToDeg = 180.0 / Py_MATH_PI;
Raymond Hettingerd6f22672002-05-13 03:56:10 +00002640
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002641/*[clinic input]
2642math.degrees
2643
2644 x: double
2645 /
2646
2647Convert angle x from radians to degrees.
2648[clinic start generated code]*/
2649
Raymond Hettingerd6f22672002-05-13 03:56:10 +00002650static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002651math_degrees_impl(PyObject *module, double x)
2652/*[clinic end generated code: output=7fea78b294acd12f input=81e016555d6e3660]*/
Raymond Hettingerd6f22672002-05-13 03:56:10 +00002653{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002654 return PyFloat_FromDouble(x * radToDeg);
Raymond Hettingerd6f22672002-05-13 03:56:10 +00002655}
2656
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002657
2658/*[clinic input]
2659math.radians
2660
2661 x: double
2662 /
2663
2664Convert angle x from degrees to radians.
2665[clinic start generated code]*/
Raymond Hettingerd6f22672002-05-13 03:56:10 +00002666
2667static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002668math_radians_impl(PyObject *module, double x)
2669/*[clinic end generated code: output=34daa47caf9b1590 input=91626fc489fe3d63]*/
Raymond Hettingerd6f22672002-05-13 03:56:10 +00002670{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002671 return PyFloat_FromDouble(x * degToRad);
Raymond Hettingerd6f22672002-05-13 03:56:10 +00002672}
2673
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002674
2675/*[clinic input]
2676math.isfinite
2677
2678 x: double
2679 /
2680
2681Return True if x is neither an infinity nor a NaN, and False otherwise.
2682[clinic start generated code]*/
Tim Peters78526162001-09-05 00:53:45 +00002683
Christian Heimes072c0f12008-01-03 23:01:04 +00002684static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002685math_isfinite_impl(PyObject *module, double x)
2686/*[clinic end generated code: output=8ba1f396440c9901 input=46967d254812e54a]*/
Mark Dickinson8e0c9962010-07-11 17:38:24 +00002687{
Mark Dickinson8e0c9962010-07-11 17:38:24 +00002688 return PyBool_FromLong((long)Py_IS_FINITE(x));
2689}
2690
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002691
2692/*[clinic input]
2693math.isnan
2694
2695 x: double
2696 /
2697
2698Return True if x is a NaN (not a number), and False otherwise.
2699[clinic start generated code]*/
Mark Dickinson8e0c9962010-07-11 17:38:24 +00002700
2701static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002702math_isnan_impl(PyObject *module, double x)
2703/*[clinic end generated code: output=f537b4d6df878c3e input=935891e66083f46a]*/
Christian Heimes072c0f12008-01-03 23:01:04 +00002704{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002705 return PyBool_FromLong((long)Py_IS_NAN(x));
Christian Heimes072c0f12008-01-03 23:01:04 +00002706}
2707
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002708
2709/*[clinic input]
2710math.isinf
2711
2712 x: double
2713 /
2714
2715Return True if x is a positive or negative infinity, and False otherwise.
2716[clinic start generated code]*/
Christian Heimes072c0f12008-01-03 23:01:04 +00002717
2718static PyObject *
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002719math_isinf_impl(PyObject *module, double x)
2720/*[clinic end generated code: output=9f00cbec4de7b06b input=32630e4212cf961f]*/
Christian Heimes072c0f12008-01-03 23:01:04 +00002721{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00002722 return PyBool_FromLong((long)Py_IS_INFINITY(x));
Christian Heimes072c0f12008-01-03 23:01:04 +00002723}
2724
Christian Heimes072c0f12008-01-03 23:01:04 +00002725
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002726/*[clinic input]
2727math.isclose -> bool
2728
2729 a: double
2730 b: double
2731 *
2732 rel_tol: double = 1e-09
2733 maximum difference for being considered "close", relative to the
2734 magnitude of the input values
2735 abs_tol: double = 0.0
2736 maximum difference for being considered "close", regardless of the
2737 magnitude of the input values
2738
2739Determine whether two floating point numbers are close in value.
2740
2741Return True if a is close in value to b, and False otherwise.
2742
2743For the values to be considered close, the difference between them
2744must be smaller than at least one of the tolerances.
2745
2746-inf, inf and NaN behave similarly to the IEEE 754 Standard. That
2747is, NaN is not close to anything, even itself. inf and -inf are
2748only close to themselves.
2749[clinic start generated code]*/
2750
2751static int
2752math_isclose_impl(PyObject *module, double a, double b, double rel_tol,
2753 double abs_tol)
2754/*[clinic end generated code: output=b73070207511952d input=f28671871ea5bfba]*/
Tal Einatd5519ed2015-05-31 22:05:00 +03002755{
Tal Einatd5519ed2015-05-31 22:05:00 +03002756 double diff = 0.0;
Tal Einatd5519ed2015-05-31 22:05:00 +03002757
2758 /* sanity check on the inputs */
2759 if (rel_tol < 0.0 || abs_tol < 0.0 ) {
2760 PyErr_SetString(PyExc_ValueError,
2761 "tolerances must be non-negative");
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002762 return -1;
Tal Einatd5519ed2015-05-31 22:05:00 +03002763 }
2764
2765 if ( a == b ) {
2766 /* short circuit exact equality -- needed to catch two infinities of
2767 the same sign. And perhaps speeds things up a bit sometimes.
2768 */
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002769 return 1;
Tal Einatd5519ed2015-05-31 22:05:00 +03002770 }
2771
2772 /* This catches the case of two infinities of opposite sign, or
2773 one infinity and one finite number. Two infinities of opposite
2774 sign would otherwise have an infinite relative tolerance.
2775 Two infinities of the same sign are caught by the equality check
2776 above.
2777 */
2778
2779 if (Py_IS_INFINITY(a) || Py_IS_INFINITY(b)) {
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002780 return 0;
Tal Einatd5519ed2015-05-31 22:05:00 +03002781 }
2782
2783 /* now do the regular computation
2784 this is essentially the "weak" test from the Boost library
2785 */
2786
2787 diff = fabs(b - a);
2788
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02002789 return (((diff <= fabs(rel_tol * b)) ||
2790 (diff <= fabs(rel_tol * a))) ||
2791 (diff <= abs_tol));
Tal Einatd5519ed2015-05-31 22:05:00 +03002792}
2793
Pablo Galindo04114112019-03-09 19:18:08 +00002794static inline int
2795_check_long_mult_overflow(long a, long b) {
2796
2797 /* From Python2's int_mul code:
2798
2799 Integer overflow checking for * is painful: Python tried a couple ways, but
2800 they didn't work on all platforms, or failed in endcases (a product of
2801 -sys.maxint-1 has been a particular pain).
2802
2803 Here's another way:
2804
2805 The native long product x*y is either exactly right or *way* off, being
2806 just the last n bits of the true product, where n is the number of bits
2807 in a long (the delivered product is the true product plus i*2**n for
2808 some integer i).
2809
2810 The native double product (double)x * (double)y is subject to three
2811 rounding errors: on a sizeof(long)==8 box, each cast to double can lose
2812 info, and even on a sizeof(long)==4 box, the multiplication can lose info.
2813 But, unlike the native long product, it's not in *range* trouble: even
2814 if sizeof(long)==32 (256-bit longs), the product easily fits in the
2815 dynamic range of a double. So the leading 50 (or so) bits of the double
2816 product are correct.
2817
2818 We check these two ways against each other, and declare victory if they're
2819 approximately the same. Else, because the native long product is the only
2820 one that can lose catastrophic amounts of information, it's the native long
2821 product that must have overflowed.
2822
2823 */
2824
2825 long longprod = (long)((unsigned long)a * b);
2826 double doubleprod = (double)a * (double)b;
2827 double doubled_longprod = (double)longprod;
2828
2829 if (doubled_longprod == doubleprod) {
2830 return 0;
2831 }
2832
2833 const double diff = doubled_longprod - doubleprod;
2834 const double absdiff = diff >= 0.0 ? diff : -diff;
2835 const double absprod = doubleprod >= 0.0 ? doubleprod : -doubleprod;
2836
2837 if (32.0 * absdiff <= absprod) {
2838 return 0;
2839 }
2840
2841 return 1;
2842}
Tal Einatd5519ed2015-05-31 22:05:00 +03002843
Pablo Galindobc098512019-02-07 07:04:02 +00002844/*[clinic input]
2845math.prod
2846
2847 iterable: object
2848 /
2849 *
2850 start: object(c_default="NULL") = 1
2851
2852Calculate the product of all the elements in the input iterable.
2853
2854The default start value for the product is 1.
2855
2856When the iterable is empty, return the start value. This function is
2857intended specifically for use with numeric values and may reject
2858non-numeric types.
2859[clinic start generated code]*/
2860
2861static PyObject *
2862math_prod_impl(PyObject *module, PyObject *iterable, PyObject *start)
2863/*[clinic end generated code: output=36153bedac74a198 input=4c5ab0682782ed54]*/
2864{
2865 PyObject *result = start;
2866 PyObject *temp, *item, *iter;
2867
2868 iter = PyObject_GetIter(iterable);
2869 if (iter == NULL) {
2870 return NULL;
2871 }
2872
2873 if (result == NULL) {
2874 result = PyLong_FromLong(1);
2875 if (result == NULL) {
2876 Py_DECREF(iter);
2877 return NULL;
2878 }
2879 } else {
2880 Py_INCREF(result);
2881 }
2882#ifndef SLOW_PROD
2883 /* Fast paths for integers keeping temporary products in C.
2884 * Assumes all inputs are the same type.
2885 * If the assumption fails, default to use PyObjects instead.
2886 */
2887 if (PyLong_CheckExact(result)) {
2888 int overflow;
2889 long i_result = PyLong_AsLongAndOverflow(result, &overflow);
2890 /* If this already overflowed, don't even enter the loop. */
2891 if (overflow == 0) {
2892 Py_DECREF(result);
2893 result = NULL;
2894 }
2895 /* Loop over all the items in the iterable until we finish, we overflow
2896 * or we found a non integer element */
2897 while(result == NULL) {
2898 item = PyIter_Next(iter);
2899 if (item == NULL) {
2900 Py_DECREF(iter);
2901 if (PyErr_Occurred()) {
2902 return NULL;
2903 }
2904 return PyLong_FromLong(i_result);
2905 }
2906 if (PyLong_CheckExact(item)) {
2907 long b = PyLong_AsLongAndOverflow(item, &overflow);
Pablo Galindo04114112019-03-09 19:18:08 +00002908 if (overflow == 0 && !_check_long_mult_overflow(i_result, b)) {
2909 long x = i_result * b;
Pablo Galindobc098512019-02-07 07:04:02 +00002910 i_result = x;
2911 Py_DECREF(item);
2912 continue;
2913 }
2914 }
2915 /* Either overflowed or is not an int.
2916 * Restore real objects and process normally */
2917 result = PyLong_FromLong(i_result);
2918 if (result == NULL) {
2919 Py_DECREF(item);
2920 Py_DECREF(iter);
2921 return NULL;
2922 }
2923 temp = PyNumber_Multiply(result, item);
2924 Py_DECREF(result);
2925 Py_DECREF(item);
2926 result = temp;
2927 if (result == NULL) {
2928 Py_DECREF(iter);
2929 return NULL;
2930 }
2931 }
2932 }
2933
2934 /* Fast paths for floats keeping temporary products in C.
2935 * Assumes all inputs are the same type.
2936 * If the assumption fails, default to use PyObjects instead.
2937 */
2938 if (PyFloat_CheckExact(result)) {
2939 double f_result = PyFloat_AS_DOUBLE(result);
2940 Py_DECREF(result);
2941 result = NULL;
2942 while(result == NULL) {
2943 item = PyIter_Next(iter);
2944 if (item == NULL) {
2945 Py_DECREF(iter);
2946 if (PyErr_Occurred()) {
2947 return NULL;
2948 }
2949 return PyFloat_FromDouble(f_result);
2950 }
2951 if (PyFloat_CheckExact(item)) {
2952 f_result *= PyFloat_AS_DOUBLE(item);
2953 Py_DECREF(item);
2954 continue;
2955 }
2956 if (PyLong_CheckExact(item)) {
2957 long value;
2958 int overflow;
2959 value = PyLong_AsLongAndOverflow(item, &overflow);
2960 if (!overflow) {
2961 f_result *= (double)value;
2962 Py_DECREF(item);
2963 continue;
2964 }
2965 }
2966 result = PyFloat_FromDouble(f_result);
2967 if (result == NULL) {
2968 Py_DECREF(item);
2969 Py_DECREF(iter);
2970 return NULL;
2971 }
2972 temp = PyNumber_Multiply(result, item);
2973 Py_DECREF(result);
2974 Py_DECREF(item);
2975 result = temp;
2976 if (result == NULL) {
2977 Py_DECREF(iter);
2978 return NULL;
2979 }
2980 }
2981 }
2982#endif
2983 /* Consume rest of the iterable (if any) that could not be handled
2984 * by specialized functions above.*/
2985 for(;;) {
2986 item = PyIter_Next(iter);
2987 if (item == NULL) {
2988 /* error, or end-of-sequence */
2989 if (PyErr_Occurred()) {
2990 Py_DECREF(result);
2991 result = NULL;
2992 }
2993 break;
2994 }
2995 temp = PyNumber_Multiply(result, item);
2996 Py_DECREF(result);
2997 Py_DECREF(item);
2998 result = temp;
2999 if (result == NULL)
3000 break;
3001 }
3002 Py_DECREF(iter);
3003 return result;
3004}
3005
3006
Yash Aggarwal4a686502019-06-01 12:51:27 +05303007/*[clinic input]
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003008math.perm
3009
3010 n: object
Raymond Hettingere119b3d2019-06-08 08:58:11 -07003011 k: object = None
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003012 /
3013
3014Number of ways to choose k items from n items without repetition and with order.
3015
Raymond Hettinger963eb0f2019-06-04 01:23:06 -07003016Evaluates to n! / (n - k)! when k <= n and evaluates
3017to zero when k > n.
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003018
Raymond Hettingere119b3d2019-06-08 08:58:11 -07003019If k is not specified or is None, then k defaults to n
3020and the function returns n!.
3021
Raymond Hettinger963eb0f2019-06-04 01:23:06 -07003022Raises TypeError if either of the arguments are not integers.
3023Raises ValueError if either of the arguments are negative.
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003024[clinic start generated code]*/
3025
3026static PyObject *
3027math_perm_impl(PyObject *module, PyObject *n, PyObject *k)
Raymond Hettingere119b3d2019-06-08 08:58:11 -07003028/*[clinic end generated code: output=e021a25469653e23 input=5311c5a00f359b53]*/
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003029{
3030 PyObject *result = NULL, *factor = NULL;
3031 int overflow, cmp;
3032 long long i, factors;
3033
Raymond Hettingere119b3d2019-06-08 08:58:11 -07003034 if (k == Py_None) {
3035 return math_factorial(module, n);
3036 }
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003037 n = PyNumber_Index(n);
3038 if (n == NULL) {
3039 return NULL;
3040 }
3041 if (!PyLong_CheckExact(n)) {
3042 Py_SETREF(n, _PyLong_Copy((PyLongObject *)n));
3043 if (n == NULL) {
3044 return NULL;
3045 }
3046 }
3047 k = PyNumber_Index(k);
3048 if (k == NULL) {
3049 Py_DECREF(n);
3050 return NULL;
3051 }
3052 if (!PyLong_CheckExact(k)) {
3053 Py_SETREF(k, _PyLong_Copy((PyLongObject *)k));
3054 if (k == NULL) {
3055 Py_DECREF(n);
3056 return NULL;
3057 }
3058 }
3059
3060 if (Py_SIZE(n) < 0) {
3061 PyErr_SetString(PyExc_ValueError,
3062 "n must be a non-negative integer");
3063 goto error;
3064 }
Mark Dickinson45e04112019-06-16 11:06:06 +01003065 if (Py_SIZE(k) < 0) {
3066 PyErr_SetString(PyExc_ValueError,
3067 "k must be a non-negative integer");
3068 goto error;
3069 }
3070
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003071 cmp = PyObject_RichCompareBool(n, k, Py_LT);
3072 if (cmp != 0) {
3073 if (cmp > 0) {
Raymond Hettinger963eb0f2019-06-04 01:23:06 -07003074 result = PyLong_FromLong(0);
3075 goto done;
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003076 }
3077 goto error;
3078 }
3079
3080 factors = PyLong_AsLongLongAndOverflow(k, &overflow);
3081 if (overflow > 0) {
3082 PyErr_Format(PyExc_OverflowError,
3083 "k must not exceed %lld",
3084 LLONG_MAX);
3085 goto error;
3086 }
Mark Dickinson45e04112019-06-16 11:06:06 +01003087 else if (factors == -1) {
3088 /* k is nonnegative, so a return value of -1 can only indicate error */
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003089 goto error;
3090 }
3091
3092 if (factors == 0) {
3093 result = PyLong_FromLong(1);
3094 goto done;
3095 }
3096
3097 result = n;
3098 Py_INCREF(result);
3099 if (factors == 1) {
3100 goto done;
3101 }
3102
3103 factor = n;
3104 Py_INCREF(factor);
3105 for (i = 1; i < factors; ++i) {
3106 Py_SETREF(factor, PyNumber_Subtract(factor, _PyLong_One));
3107 if (factor == NULL) {
3108 goto error;
3109 }
3110 Py_SETREF(result, PyNumber_Multiply(result, factor));
3111 if (result == NULL) {
3112 goto error;
3113 }
3114 }
3115 Py_DECREF(factor);
3116
3117done:
3118 Py_DECREF(n);
3119 Py_DECREF(k);
3120 return result;
3121
3122error:
3123 Py_XDECREF(factor);
3124 Py_XDECREF(result);
3125 Py_DECREF(n);
3126 Py_DECREF(k);
3127 return NULL;
3128}
3129
3130
3131/*[clinic input]
Yash Aggarwal4a686502019-06-01 12:51:27 +05303132math.comb
3133
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003134 n: object
3135 k: object
3136 /
Yash Aggarwal4a686502019-06-01 12:51:27 +05303137
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003138Number of ways to choose k items from n items without repetition and without order.
Yash Aggarwal4a686502019-06-01 12:51:27 +05303139
Raymond Hettinger963eb0f2019-06-04 01:23:06 -07003140Evaluates to n! / (k! * (n - k)!) when k <= n and evaluates
3141to zero when k > n.
Yash Aggarwal4a686502019-06-01 12:51:27 +05303142
Raymond Hettinger963eb0f2019-06-04 01:23:06 -07003143Also called the binomial coefficient because it is equivalent
3144to the coefficient of k-th term in polynomial expansion of the
3145expression (1 + x)**n.
3146
3147Raises TypeError if either of the arguments are not integers.
3148Raises ValueError if either of the arguments are negative.
Yash Aggarwal4a686502019-06-01 12:51:27 +05303149
3150[clinic start generated code]*/
3151
3152static PyObject *
3153math_comb_impl(PyObject *module, PyObject *n, PyObject *k)
Raymond Hettinger963eb0f2019-06-04 01:23:06 -07003154/*[clinic end generated code: output=bd2cec8d854f3493 input=9a05315af2518709]*/
Yash Aggarwal4a686502019-06-01 12:51:27 +05303155{
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003156 PyObject *result = NULL, *factor = NULL, *temp;
Yash Aggarwal4a686502019-06-01 12:51:27 +05303157 int overflow, cmp;
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003158 long long i, factors;
Yash Aggarwal4a686502019-06-01 12:51:27 +05303159
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003160 n = PyNumber_Index(n);
3161 if (n == NULL) {
3162 return NULL;
Yash Aggarwal4a686502019-06-01 12:51:27 +05303163 }
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003164 if (!PyLong_CheckExact(n)) {
3165 Py_SETREF(n, _PyLong_Copy((PyLongObject *)n));
3166 if (n == NULL) {
3167 return NULL;
3168 }
3169 }
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003170 k = PyNumber_Index(k);
3171 if (k == NULL) {
3172 Py_DECREF(n);
3173 return NULL;
Yash Aggarwal4a686502019-06-01 12:51:27 +05303174 }
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003175 if (!PyLong_CheckExact(k)) {
3176 Py_SETREF(k, _PyLong_Copy((PyLongObject *)k));
3177 if (k == NULL) {
3178 Py_DECREF(n);
3179 return NULL;
3180 }
3181 }
Yash Aggarwal4a686502019-06-01 12:51:27 +05303182
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003183 if (Py_SIZE(n) < 0) {
3184 PyErr_SetString(PyExc_ValueError,
3185 "n must be a non-negative integer");
3186 goto error;
Yash Aggarwal4a686502019-06-01 12:51:27 +05303187 }
Mark Dickinson45e04112019-06-16 11:06:06 +01003188 if (Py_SIZE(k) < 0) {
3189 PyErr_SetString(PyExc_ValueError,
3190 "k must be a non-negative integer");
3191 goto error;
3192 }
3193
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003194 /* k = min(k, n - k) */
3195 temp = PyNumber_Subtract(n, k);
3196 if (temp == NULL) {
3197 goto error;
Yash Aggarwal4a686502019-06-01 12:51:27 +05303198 }
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003199 if (Py_SIZE(temp) < 0) {
3200 Py_DECREF(temp);
Raymond Hettinger963eb0f2019-06-04 01:23:06 -07003201 result = PyLong_FromLong(0);
3202 goto done;
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003203 }
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003204 cmp = PyObject_RichCompareBool(temp, k, Py_LT);
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003205 if (cmp > 0) {
3206 Py_SETREF(k, temp);
Yash Aggarwal4a686502019-06-01 12:51:27 +05303207 }
3208 else {
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003209 Py_DECREF(temp);
3210 if (cmp < 0) {
3211 goto error;
3212 }
Yash Aggarwal4a686502019-06-01 12:51:27 +05303213 }
3214
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003215 factors = PyLong_AsLongLongAndOverflow(k, &overflow);
3216 if (overflow > 0) {
Yash Aggarwal4a686502019-06-01 12:51:27 +05303217 PyErr_Format(PyExc_OverflowError,
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003218 "min(n - k, k) must not exceed %lld",
Yash Aggarwal4a686502019-06-01 12:51:27 +05303219 LLONG_MAX);
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003220 goto error;
Yash Aggarwal4a686502019-06-01 12:51:27 +05303221 }
Mark Dickinson45e04112019-06-16 11:06:06 +01003222 if (factors == -1) {
3223 /* k is nonnegative, so a return value of -1 can only indicate error */
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003224 goto error;
Yash Aggarwal4a686502019-06-01 12:51:27 +05303225 }
3226
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003227 if (factors == 0) {
3228 result = PyLong_FromLong(1);
3229 goto done;
Yash Aggarwal4a686502019-06-01 12:51:27 +05303230 }
3231
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003232 result = n;
3233 Py_INCREF(result);
3234 if (factors == 1) {
3235 goto done;
Yash Aggarwal4a686502019-06-01 12:51:27 +05303236 }
3237
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003238 factor = n;
3239 Py_INCREF(factor);
3240 for (i = 1; i < factors; ++i) {
3241 Py_SETREF(factor, PyNumber_Subtract(factor, _PyLong_One));
3242 if (factor == NULL) {
3243 goto error;
3244 }
3245 Py_SETREF(result, PyNumber_Multiply(result, factor));
3246 if (result == NULL) {
3247 goto error;
3248 }
Yash Aggarwal4a686502019-06-01 12:51:27 +05303249
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003250 temp = PyLong_FromUnsignedLongLong((unsigned long long)i + 1);
3251 if (temp == NULL) {
3252 goto error;
3253 }
3254 Py_SETREF(result, PyNumber_FloorDivide(result, temp));
3255 Py_DECREF(temp);
3256 if (result == NULL) {
3257 goto error;
3258 }
3259 }
3260 Py_DECREF(factor);
Yash Aggarwal4a686502019-06-01 12:51:27 +05303261
Serhiy Storchaka2b843ac2019-06-01 22:09:02 +03003262done:
3263 Py_DECREF(n);
3264 Py_DECREF(k);
3265 return result;
3266
3267error:
3268 Py_XDECREF(factor);
3269 Py_XDECREF(result);
3270 Py_DECREF(n);
3271 Py_DECREF(k);
Yash Aggarwal4a686502019-06-01 12:51:27 +05303272 return NULL;
3273}
3274
3275
Barry Warsaw8b43b191996-12-09 22:32:36 +00003276static PyMethodDef math_methods[] = {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003277 {"acos", math_acos, METH_O, math_acos_doc},
3278 {"acosh", math_acosh, METH_O, math_acosh_doc},
3279 {"asin", math_asin, METH_O, math_asin_doc},
3280 {"asinh", math_asinh, METH_O, math_asinh_doc},
3281 {"atan", math_atan, METH_O, math_atan_doc},
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02003282 {"atan2", (PyCFunction)(void(*)(void))math_atan2, METH_FASTCALL, math_atan2_doc},
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003283 {"atanh", math_atanh, METH_O, math_atanh_doc},
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02003284 MATH_CEIL_METHODDEF
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02003285 {"copysign", (PyCFunction)(void(*)(void))math_copysign, METH_FASTCALL, math_copysign_doc},
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003286 {"cos", math_cos, METH_O, math_cos_doc},
3287 {"cosh", math_cosh, METH_O, math_cosh_doc},
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02003288 MATH_DEGREES_METHODDEF
Raymond Hettinger9c18b1a2018-07-31 00:45:49 -07003289 MATH_DIST_METHODDEF
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003290 {"erf", math_erf, METH_O, math_erf_doc},
3291 {"erfc", math_erfc, METH_O, math_erfc_doc},
3292 {"exp", math_exp, METH_O, math_exp_doc},
3293 {"expm1", math_expm1, METH_O, math_expm1_doc},
3294 {"fabs", math_fabs, METH_O, math_fabs_doc},
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02003295 MATH_FACTORIAL_METHODDEF
3296 MATH_FLOOR_METHODDEF
3297 MATH_FMOD_METHODDEF
3298 MATH_FREXP_METHODDEF
3299 MATH_FSUM_METHODDEF
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003300 {"gamma", math_gamma, METH_O, math_gamma_doc},
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02003301 MATH_GCD_METHODDEF
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02003302 {"hypot", (PyCFunction)(void(*)(void))math_hypot, METH_FASTCALL, math_hypot_doc},
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02003303 MATH_ISCLOSE_METHODDEF
3304 MATH_ISFINITE_METHODDEF
3305 MATH_ISINF_METHODDEF
3306 MATH_ISNAN_METHODDEF
Mark Dickinson73934b92019-05-18 12:29:50 +01003307 MATH_ISQRT_METHODDEF
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02003308 MATH_LDEXP_METHODDEF
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003309 {"lgamma", math_lgamma, METH_O, math_lgamma_doc},
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02003310 MATH_LOG_METHODDEF
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003311 {"log1p", math_log1p, METH_O, math_log1p_doc},
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02003312 MATH_LOG10_METHODDEF
3313 MATH_LOG2_METHODDEF
3314 MATH_MODF_METHODDEF
3315 MATH_POW_METHODDEF
3316 MATH_RADIANS_METHODDEF
Serhiy Storchakad0d3e992019-01-12 08:26:34 +02003317 {"remainder", (PyCFunction)(void(*)(void))math_remainder, METH_FASTCALL, math_remainder_doc},
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003318 {"sin", math_sin, METH_O, math_sin_doc},
3319 {"sinh", math_sinh, METH_O, math_sinh_doc},
3320 {"sqrt", math_sqrt, METH_O, math_sqrt_doc},
3321 {"tan", math_tan, METH_O, math_tan_doc},
3322 {"tanh", math_tanh, METH_O, math_tanh_doc},
Serhiy Storchakac9ea9332017-01-19 18:13:09 +02003323 MATH_TRUNC_METHODDEF
Pablo Galindobc098512019-02-07 07:04:02 +00003324 MATH_PROD_METHODDEF
Serhiy Storchaka5ae299a2019-06-02 11:16:49 +03003325 MATH_PERM_METHODDEF
Yash Aggarwal4a686502019-06-01 12:51:27 +05303326 MATH_COMB_METHODDEF
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003327 {NULL, NULL} /* sentinel */
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00003328};
3329
Guido van Rossumc6e22901998-12-04 19:26:43 +00003330
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +00003331PyDoc_STRVAR(module_doc,
Ned Batchelder6faad352019-05-17 05:59:14 -04003332"This module provides access to the mathematical functions\n"
3333"defined by the C standard.");
Guido van Rossumc6e22901998-12-04 19:26:43 +00003334
Martin v. Löwis1a214512008-06-11 05:26:20 +00003335
3336static struct PyModuleDef mathmodule = {
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003337 PyModuleDef_HEAD_INIT,
3338 "math",
3339 module_doc,
3340 -1,
3341 math_methods,
3342 NULL,
3343 NULL,
3344 NULL,
3345 NULL
Martin v. Löwis1a214512008-06-11 05:26:20 +00003346};
3347
Mark Hammondfe51c6d2002-08-02 02:27:13 +00003348PyMODINIT_FUNC
Martin v. Löwis1a214512008-06-11 05:26:20 +00003349PyInit_math(void)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00003350{
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003351 PyObject *m;
Tim Petersfe71f812001-08-07 22:10:00 +00003352
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003353 m = PyModule_Create(&mathmodule);
3354 if (m == NULL)
3355 goto finally;
Barry Warsawfc93f751996-12-17 00:47:03 +00003356
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003357 PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
3358 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
Guido van Rossum0a891d72016-08-15 09:12:52 -07003359 PyModule_AddObject(m, "tau", PyFloat_FromDouble(Py_MATH_TAU)); /* 2pi */
Mark Dickinsona5d0c7c2015-01-11 11:55:29 +00003360 PyModule_AddObject(m, "inf", PyFloat_FromDouble(m_inf()));
3361#if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN)
3362 PyModule_AddObject(m, "nan", PyFloat_FromDouble(m_nan()));
3363#endif
Barry Warsawfc93f751996-12-17 00:47:03 +00003364
Mark Dickinsona5d0c7c2015-01-11 11:55:29 +00003365 finally:
Antoine Pitrouf95a1b32010-05-09 15:52:27 +00003366 return m;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00003367}