LCOV - code coverage report
Current view: | top level - Modules - mathmodule.c (source / functions) |  |  | Hit | Total | Coverage
Test: | CPython 3.12 LCOV report [commit 5e6661bce9] |  | Lines: | 1207 | 1345 | 89.7 %
Date: | 2023-03-20 08:15:36 |  | Functions: | 93 | 93 | 100.0 %
|  |  | Branches: | 753 | 918 | 82.0 %

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

Generated by: LCOV version 1.14