Home | History | Annotate | Download | only in Modules
      1 /* Complex math module */
      2 
      3 /* much code borrowed from mathmodule.c */
      4 
      5 #include "Python.h"
      6 #include "_math.h"
      7 /* we need DBL_MAX, DBL_MIN, DBL_EPSILON, DBL_MANT_DIG and FLT_RADIX from
      8    float.h.  We assume that FLT_RADIX is either 2 or 16. */
      9 #include <float.h>
     10 
     11 #include "clinic/cmathmodule.c.h"
     12 /*[clinic input]
     13 module cmath
     14 [clinic start generated code]*/
     15 /*[clinic end generated code: output=da39a3ee5e6b4b0d input=308d6839f4a46333]*/
     16 
     17 /*[python input]
     18 class Py_complex_protected_converter(Py_complex_converter):
     19     def modify(self):
     20         return 'errno = 0; PyFPE_START_PROTECT("complex function", goto exit);'
     21 
     22 
     23 class Py_complex_protected_return_converter(CReturnConverter):
     24     type = "Py_complex"
     25 
     26     def render(self, function, data):
     27         self.declare(data)
     28         data.return_conversion.append("""
     29 PyFPE_END_PROTECT(_return_value);
     30 if (errno == EDOM) {
     31     PyErr_SetString(PyExc_ValueError, "math domain error");
     32     goto exit;
     33 }
     34 else if (errno == ERANGE) {
     35     PyErr_SetString(PyExc_OverflowError, "math range error");
     36     goto exit;
     37 }
     38 else {
     39     return_value = PyComplex_FromCComplex(_return_value);
     40 }
     41 """.strip())
     42 [python start generated code]*/
     43 /*[python end generated code: output=da39a3ee5e6b4b0d input=345daa075b1028e7]*/
     44 
     45 #if (FLT_RADIX != 2 && FLT_RADIX != 16)
     46 #error "Modules/cmathmodule.c expects FLT_RADIX to be 2 or 16"
     47 #endif
     48 
     49 #ifndef M_LN2
     50 #define M_LN2 (0.6931471805599453094) /* natural log of 2 */
     51 #endif
     52 
     53 #ifndef M_LN10
     54 #define M_LN10 (2.302585092994045684) /* natural log of 10 */
     55 #endif
     56 
     57 /*
     58    CM_LARGE_DOUBLE is used to avoid spurious overflow in the sqrt, log,
     59    inverse trig and inverse hyperbolic trig functions.  Its log is used in the
     60    evaluation of exp, cos, cosh, sin, sinh, tan, and tanh to avoid unnecessary
     61    overflow.
     62  */
     63 
     64 #define CM_LARGE_DOUBLE (DBL_MAX/4.)
     65 #define CM_SQRT_LARGE_DOUBLE (sqrt(CM_LARGE_DOUBLE))
     66 #define CM_LOG_LARGE_DOUBLE (log(CM_LARGE_DOUBLE))
     67 #define CM_SQRT_DBL_MIN (sqrt(DBL_MIN))
     68 
     69 /*
     70    CM_SCALE_UP is an odd integer chosen such that multiplication by
     71    2**CM_SCALE_UP is sufficient to turn a subnormal into a normal.
     72    CM_SCALE_DOWN is (-(CM_SCALE_UP+1)/2).  These scalings are used to compute
     73    square roots accurately when the real and imaginary parts of the argument
     74    are subnormal.
     75 */
     76 
     77 #if FLT_RADIX==2
     78 #define CM_SCALE_UP (2*(DBL_MANT_DIG/2) + 1)
     79 #elif FLT_RADIX==16
     80 #define CM_SCALE_UP (4*DBL_MANT_DIG+1)
     81 #endif
     82 #define CM_SCALE_DOWN (-(CM_SCALE_UP+1)/2)
     83 
     84 /* Constants cmath.inf, cmath.infj, cmath.nan, cmath.nanj.
     85    cmath.nan and cmath.nanj are defined only when either
     86    PY_NO_SHORT_FLOAT_REPR is *not* defined (which should be
     87    the most common situation on machines using an IEEE 754
     88    representation), or Py_NAN is defined. */
     89 
     90 static double
     91 m_inf(void)
     92 {
     93 #ifndef PY_NO_SHORT_FLOAT_REPR
     94     return _Py_dg_infinity(0);
     95 #else
     96     return Py_HUGE_VAL;
     97 #endif
     98 }
     99 
    100 static Py_complex
    101 c_infj(void)
    102 {
    103     Py_complex r;
    104     r.real = 0.0;
    105     r.imag = m_inf();
    106     return r;
    107 }
    108 
    109 #if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN)
    110 
    111 static double
    112 m_nan(void)
    113 {
    114 #ifndef PY_NO_SHORT_FLOAT_REPR
    115     return _Py_dg_stdnan(0);
    116 #else
    117     return Py_NAN;
    118 #endif
    119 }
    120 
    121 static Py_complex
    122 c_nanj(void)
    123 {
    124     Py_complex r;
    125     r.real = 0.0;
    126     r.imag = m_nan();
    127     return r;
    128 }
    129 
    130 #endif
    131 
    132 /* forward declarations */
    133 static Py_complex cmath_asinh_impl(PyObject *, Py_complex);
    134 static Py_complex cmath_atanh_impl(PyObject *, Py_complex);
    135 static Py_complex cmath_cosh_impl(PyObject *, Py_complex);
    136 static Py_complex cmath_sinh_impl(PyObject *, Py_complex);
    137 static Py_complex cmath_sqrt_impl(PyObject *, Py_complex);
    138 static Py_complex cmath_tanh_impl(PyObject *, Py_complex);
    139 static PyObject * math_error(void);
    140 
    141 /* Code to deal with special values (infinities, NaNs, etc.). */
    142 
    143 /* special_type takes a double and returns an integer code indicating
    144    the type of the double as follows:
    145 */
    146 
    147 enum special_types {
    148     ST_NINF,            /* 0, negative infinity */
    149     ST_NEG,             /* 1, negative finite number (nonzero) */
    150     ST_NZERO,           /* 2, -0. */
    151     ST_PZERO,           /* 3, +0. */
    152     ST_POS,             /* 4, positive finite number (nonzero) */
    153     ST_PINF,            /* 5, positive infinity */
    154     ST_NAN              /* 6, Not a Number */
    155 };
    156 
    157 static enum special_types
    158 special_type(double d)
    159 {
    160     if (Py_IS_FINITE(d)) {
    161         if (d != 0) {
    162             if (copysign(1., d) == 1.)
    163                 return ST_POS;
    164             else
    165                 return ST_NEG;
    166         }
    167         else {
    168             if (copysign(1., d) == 1.)
    169                 return ST_PZERO;
    170             else
    171                 return ST_NZERO;
    172         }
    173     }
    174     if (Py_IS_NAN(d))
    175         return ST_NAN;
    176     if (copysign(1., d) == 1.)
    177         return ST_PINF;
    178     else
    179         return ST_NINF;
    180 }
    181 
    182 #define SPECIAL_VALUE(z, table)                                         \
    183     if (!Py_IS_FINITE((z).real) || !Py_IS_FINITE((z).imag)) {           \
    184         errno = 0;                                              \
    185         return table[special_type((z).real)]                            \
    186                     [special_type((z).imag)];                           \
    187     }
    188 
    189 #define P Py_MATH_PI
    190 #define P14 0.25*Py_MATH_PI
    191 #define P12 0.5*Py_MATH_PI
    192 #define P34 0.75*Py_MATH_PI
    193 #define INF Py_HUGE_VAL
    194 #define N Py_NAN
    195 #define U -9.5426319407711027e33 /* unlikely value, used as placeholder */
    196 
    197 /* First, the C functions that do the real work.  Each of the c_*
    198    functions computes and returns the C99 Annex G recommended result
    199    and also sets errno as follows: errno = 0 if no floating-point
    200    exception is associated with the result; errno = EDOM if C99 Annex
    201    G recommends raising divide-by-zero or invalid for this result; and
    202    errno = ERANGE where the overflow floating-point signal should be
    203    raised.
    204 */
    205 
    206 static Py_complex acos_special_values[7][7];
    207 
    208 /*[clinic input]
    209 cmath.acos -> Py_complex_protected
    210 
    211     z: Py_complex_protected
    212     /
    213 
    214 Return the arc cosine of z.
    215 [clinic start generated code]*/
    216 
    217 static Py_complex
    218 cmath_acos_impl(PyObject *module, Py_complex z)
    219 /*[clinic end generated code: output=40bd42853fd460ae input=bd6cbd78ae851927]*/
    220 {
    221     Py_complex s1, s2, r;
    222 
    223     SPECIAL_VALUE(z, acos_special_values);
    224 
    225     if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
    226         /* avoid unnecessary overflow for large arguments */
    227         r.real = atan2(fabs(z.imag), z.real);
    228         /* split into cases to make sure that the branch cut has the
    229            correct continuity on systems with unsigned zeros */
    230         if (z.real < 0.) {
    231             r.imag = -copysign(log(hypot(z.real/2., z.imag/2.)) +
    232                                M_LN2*2., z.imag);
    233         } else {
    234             r.imag = copysign(log(hypot(z.real/2., z.imag/2.)) +
    235                               M_LN2*2., -z.imag);
    236         }
    237     } else {
    238         s1.real = 1.-z.real;
    239         s1.imag = -z.imag;
    240         s1 = cmath_sqrt_impl(module, s1);
    241         s2.real = 1.+z.real;
    242         s2.imag = z.imag;
    243         s2 = cmath_sqrt_impl(module, s2);
    244         r.real = 2.*atan2(s1.real, s2.real);
    245         r.imag = m_asinh(s2.real*s1.imag - s2.imag*s1.real);
    246     }
    247     errno = 0;
    248     return r;
    249 }
    250 
    251 
    252 static Py_complex acosh_special_values[7][7];
    253 
    254 /*[clinic input]
    255 cmath.acosh = cmath.acos
    256 
    257 Return the inverse hyperbolic cosine of z.
    258 [clinic start generated code]*/
    259 
    260 static Py_complex
    261 cmath_acosh_impl(PyObject *module, Py_complex z)
    262 /*[clinic end generated code: output=3e2454d4fcf404ca input=3f61bee7d703e53c]*/
    263 {
    264     Py_complex s1, s2, r;
    265 
    266     SPECIAL_VALUE(z, acosh_special_values);
    267 
    268     if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
    269         /* avoid unnecessary overflow for large arguments */
    270         r.real = log(hypot(z.real/2., z.imag/2.)) + M_LN2*2.;
    271         r.imag = atan2(z.imag, z.real);
    272     } else {
    273         s1.real = z.real - 1.;
    274         s1.imag = z.imag;
    275         s1 = cmath_sqrt_impl(module, s1);
    276         s2.real = z.real + 1.;
    277         s2.imag = z.imag;
    278         s2 = cmath_sqrt_impl(module, s2);
    279         r.real = m_asinh(s1.real*s2.real + s1.imag*s2.imag);
    280         r.imag = 2.*atan2(s1.imag, s2.real);
    281     }
    282     errno = 0;
    283     return r;
    284 }
    285 
    286 /*[clinic input]
    287 cmath.asin = cmath.acos
    288 
    289 Return the arc sine of z.
    290 [clinic start generated code]*/
    291 
    292 static Py_complex
    293 cmath_asin_impl(PyObject *module, Py_complex z)
    294 /*[clinic end generated code: output=3b264cd1b16bf4e1 input=be0bf0cfdd5239c5]*/
    295 {
    296     /* asin(z) = -i asinh(iz) */
    297     Py_complex s, r;
    298     s.real = -z.imag;
    299     s.imag = z.real;
    300     s = cmath_asinh_impl(module, s);
    301     r.real = s.imag;
    302     r.imag = -s.real;
    303     return r;
    304 }
    305 
    306 
    307 static Py_complex asinh_special_values[7][7];
    308 
    309 /*[clinic input]
    310 cmath.asinh = cmath.acos
    311 
    312 Return the inverse hyperbolic sine of z.
    313 [clinic start generated code]*/
    314 
    315 static Py_complex
    316 cmath_asinh_impl(PyObject *module, Py_complex z)
    317 /*[clinic end generated code: output=733d8107841a7599 input=5c09448fcfc89a79]*/
    318 {
    319     Py_complex s1, s2, r;
    320 
    321     SPECIAL_VALUE(z, asinh_special_values);
    322 
    323     if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
    324         if (z.imag >= 0.) {
    325             r.real = copysign(log(hypot(z.real/2., z.imag/2.)) +
    326                               M_LN2*2., z.real);
    327         } else {
    328             r.real = -copysign(log(hypot(z.real/2., z.imag/2.)) +
    329                                M_LN2*2., -z.real);
    330         }
    331         r.imag = atan2(z.imag, fabs(z.real));
    332     } else {
    333         s1.real = 1.+z.imag;
    334         s1.imag = -z.real;
    335         s1 = cmath_sqrt_impl(module, s1);
    336         s2.real = 1.-z.imag;
    337         s2.imag = z.real;
    338         s2 = cmath_sqrt_impl(module, s2);
    339         r.real = m_asinh(s1.real*s2.imag-s2.real*s1.imag);
    340         r.imag = atan2(z.imag, s1.real*s2.real-s1.imag*s2.imag);
    341     }
    342     errno = 0;
    343     return r;
    344 }
    345 
    346 
    347 /*[clinic input]
    348 cmath.atan = cmath.acos
    349 
    350 Return the arc tangent of z.
    351 [clinic start generated code]*/
    352 
    353 static Py_complex
    354 cmath_atan_impl(PyObject *module, Py_complex z)
    355 /*[clinic end generated code: output=b6bfc497058acba4 input=3b21ff7d5eac632a]*/
    356 {
    357     /* atan(z) = -i atanh(iz) */
    358     Py_complex s, r;
    359     s.real = -z.imag;
    360     s.imag = z.real;
    361     s = cmath_atanh_impl(module, s);
    362     r.real = s.imag;
    363     r.imag = -s.real;
    364     return r;
    365 }
    366 
    367 /* Windows screws up atan2 for inf and nan, and alpha Tru64 5.1 doesn't follow
    368    C99 for atan2(0., 0.). */
    369 static double
    370 c_atan2(Py_complex z)
    371 {
    372     if (Py_IS_NAN(z.real) || Py_IS_NAN(z.imag))
    373         return Py_NAN;
    374     if (Py_IS_INFINITY(z.imag)) {
    375         if (Py_IS_INFINITY(z.real)) {
    376             if (copysign(1., z.real) == 1.)
    377                 /* atan2(+-inf, +inf) == +-pi/4 */
    378                 return copysign(0.25*Py_MATH_PI, z.imag);
    379             else
    380                 /* atan2(+-inf, -inf) == +-pi*3/4 */
    381                 return copysign(0.75*Py_MATH_PI, z.imag);
    382         }
    383         /* atan2(+-inf, x) == +-pi/2 for finite x */
    384         return copysign(0.5*Py_MATH_PI, z.imag);
    385     }
    386     if (Py_IS_INFINITY(z.real) || z.imag == 0.) {
    387         if (copysign(1., z.real) == 1.)
    388             /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
    389             return copysign(0., z.imag);
    390         else
    391             /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
    392             return copysign(Py_MATH_PI, z.imag);
    393     }
    394     return atan2(z.imag, z.real);
    395 }
    396 
    397 
    398 static Py_complex atanh_special_values[7][7];
    399 
    400 /*[clinic input]
    401 cmath.atanh = cmath.acos
    402 
    403 Return the inverse hyperbolic tangent of z.
    404 [clinic start generated code]*/
    405 
    406 static Py_complex
    407 cmath_atanh_impl(PyObject *module, Py_complex z)
    408 /*[clinic end generated code: output=e83355f93a989c9e input=2b3fdb82fb34487b]*/
    409 {
    410     Py_complex r;
    411     double ay, h;
    412 
    413     SPECIAL_VALUE(z, atanh_special_values);
    414 
    415     /* Reduce to case where z.real >= 0., using atanh(z) = -atanh(-z). */
    416     if (z.real < 0.) {
    417         return _Py_c_neg(cmath_atanh_impl(module, _Py_c_neg(z)));
    418     }
    419 
    420     ay = fabs(z.imag);
    421     if (z.real > CM_SQRT_LARGE_DOUBLE || ay > CM_SQRT_LARGE_DOUBLE) {
    422         /*
    423            if abs(z) is large then we use the approximation
    424            atanh(z) ~ 1/z +/- i*pi/2 (+/- depending on the sign
    425            of z.imag)
    426         */
    427         h = hypot(z.real/2., z.imag/2.);  /* safe from overflow */
    428         r.real = z.real/4./h/h;
    429         /* the two negations in the next line cancel each other out
    430            except when working with unsigned zeros: they're there to
    431            ensure that the branch cut has the correct continuity on
    432            systems that don't support signed zeros */
    433         r.imag = -copysign(Py_MATH_PI/2., -z.imag);
    434         errno = 0;
    435     } else if (z.real == 1. && ay < CM_SQRT_DBL_MIN) {
    436         /* C99 standard says:  atanh(1+/-0.) should be inf +/- 0i */
    437         if (ay == 0.) {
    438             r.real = INF;
    439             r.imag = z.imag;
    440             errno = EDOM;
    441         } else {
    442             r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.)));
    443             r.imag = copysign(atan2(2., -ay)/2, z.imag);
    444             errno = 0;
    445         }
    446     } else {
    447         r.real = m_log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.;
    448         r.imag = -atan2(-2.*z.imag, (1-z.real)*(1+z.real) - ay*ay)/2.;
    449         errno = 0;
    450     }
    451     return r;
    452 }
    453 
    454 
    455 /*[clinic input]
    456 cmath.cos = cmath.acos
    457 
    458 Return the cosine of z.
    459 [clinic start generated code]*/
    460 
    461 static Py_complex
    462 cmath_cos_impl(PyObject *module, Py_complex z)
    463 /*[clinic end generated code: output=fd64918d5b3186db input=6022e39b77127ac7]*/
    464 {
    465     /* cos(z) = cosh(iz) */
    466     Py_complex r;
    467     r.real = -z.imag;
    468     r.imag = z.real;
    469     r = cmath_cosh_impl(module, r);
    470     return r;
    471 }
    472 
    473 
    474 /* cosh(infinity + i*y) needs to be dealt with specially */
    475 static Py_complex cosh_special_values[7][7];
    476 
    477 /*[clinic input]
    478 cmath.cosh = cmath.acos
    479 
    480 Return the hyperbolic cosine of z.
    481 [clinic start generated code]*/
    482 
    483 static Py_complex
    484 cmath_cosh_impl(PyObject *module, Py_complex z)
    485 /*[clinic end generated code: output=2e969047da601bdb input=d6b66339e9cc332b]*/
    486 {
    487     Py_complex r;
    488     double x_minus_one;
    489 
    490     /* special treatment for cosh(+/-inf + iy) if y is not a NaN */
    491     if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
    492         if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) &&
    493             (z.imag != 0.)) {
    494             if (z.real > 0) {
    495                 r.real = copysign(INF, cos(z.imag));
    496                 r.imag = copysign(INF, sin(z.imag));
    497             }
    498             else {
    499                 r.real = copysign(INF, cos(z.imag));
    500                 r.imag = -copysign(INF, sin(z.imag));
    501             }
    502         }
    503         else {
    504             r = cosh_special_values[special_type(z.real)]
    505                                    [special_type(z.imag)];
    506         }
    507         /* need to set errno = EDOM if y is +/- infinity and x is not
    508            a NaN */
    509         if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
    510             errno = EDOM;
    511         else
    512             errno = 0;
    513         return r;
    514     }
    515 
    516     if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
    517         /* deal correctly with cases where cosh(z.real) overflows but
    518            cosh(z) does not. */
    519         x_minus_one = z.real - copysign(1., z.real);
    520         r.real = cos(z.imag) * cosh(x_minus_one) * Py_MATH_E;
    521         r.imag = sin(z.imag) * sinh(x_minus_one) * Py_MATH_E;
    522     } else {
    523         r.real = cos(z.imag) * cosh(z.real);
    524         r.imag = sin(z.imag) * sinh(z.real);
    525     }
    526     /* detect overflow, and set errno accordingly */
    527     if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
    528         errno = ERANGE;
    529     else
    530         errno = 0;
    531     return r;
    532 }
    533 
    534 
    535 /* exp(infinity + i*y) and exp(-infinity + i*y) need special treatment for
    536    finite y */
    537 static Py_complex exp_special_values[7][7];
    538 
    539 /*[clinic input]
    540 cmath.exp = cmath.acos
    541 
    542 Return the exponential value e**z.
    543 [clinic start generated code]*/
    544 
    545 static Py_complex
    546 cmath_exp_impl(PyObject *module, Py_complex z)
    547 /*[clinic end generated code: output=edcec61fb9dfda6c input=8b9e6cf8a92174c3]*/
    548 {
    549     Py_complex r;
    550     double l;
    551 
    552     if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
    553         if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
    554             && (z.imag != 0.)) {
    555             if (z.real > 0) {
    556                 r.real = copysign(INF, cos(z.imag));
    557                 r.imag = copysign(INF, sin(z.imag));
    558             }
    559             else {
    560                 r.real = copysign(0., cos(z.imag));
    561                 r.imag = copysign(0., sin(z.imag));
    562             }
    563         }
    564         else {
    565             r = exp_special_values[special_type(z.real)]
    566                                   [special_type(z.imag)];
    567         }
    568         /* need to set errno = EDOM if y is +/- infinity and x is not
    569            a NaN and not -infinity */
    570         if (Py_IS_INFINITY(z.imag) &&
    571             (Py_IS_FINITE(z.real) ||
    572              (Py_IS_INFINITY(z.real) && z.real > 0)))
    573             errno = EDOM;
    574         else
    575             errno = 0;
    576         return r;
    577     }
    578 
    579     if (z.real > CM_LOG_LARGE_DOUBLE) {
    580         l = exp(z.real-1.);
    581         r.real = l*cos(z.imag)*Py_MATH_E;
    582         r.imag = l*sin(z.imag)*Py_MATH_E;
    583     } else {
    584         l = exp(z.real);
    585         r.real = l*cos(z.imag);
    586         r.imag = l*sin(z.imag);
    587     }
    588     /* detect overflow, and set errno accordingly */
    589     if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
    590         errno = ERANGE;
    591     else
    592         errno = 0;
    593     return r;
    594 }
    595 
    596 static Py_complex log_special_values[7][7];
    597 
    598 static Py_complex
    599 c_log(Py_complex z)
    600 {
    601     /*
    602        The usual formula for the real part is log(hypot(z.real, z.imag)).
    603        There are four situations where this formula is potentially
    604        problematic:
    605 
    606        (1) the absolute value of z is subnormal.  Then hypot is subnormal,
    607        so has fewer than the usual number of bits of accuracy, hence may
    608        have large relative error.  This then gives a large absolute error
    609        in the log.  This can be solved by rescaling z by a suitable power
    610        of 2.
    611 
    612        (2) the absolute value of z is greater than DBL_MAX (e.g. when both
    613        z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX)
    614        Again, rescaling solves this.
    615 
    616        (3) the absolute value of z is close to 1.  In this case it's
    617        difficult to achieve good accuracy, at least in part because a
    618        change of 1ulp in the real or imaginary part of z can result in a
    619        change of billions of ulps in the correctly rounded answer.
    620 
    621        (4) z = 0.  The simplest thing to do here is to call the
    622        floating-point log with an argument of 0, and let its behaviour
    623        (returning -infinity, signaling a floating-point exception, setting
    624        errno, or whatever) determine that of c_log.  So the usual formula
    625        is fine here.
    626 
    627      */
    628 
    629     Py_complex r;
    630     double ax, ay, am, an, h;
    631 
    632     SPECIAL_VALUE(z, log_special_values);
    633 
    634     ax = fabs(z.real);
    635     ay = fabs(z.imag);
    636 
    637     if (ax > CM_LARGE_DOUBLE || ay > CM_LARGE_DOUBLE) {
    638         r.real = log(hypot(ax/2., ay/2.)) + M_LN2;
    639     } else if (ax < DBL_MIN && ay < DBL_MIN) {
    640         if (ax > 0. || ay > 0.) {
    641             /* catch cases where hypot(ax, ay) is subnormal */
    642             r.real = log(hypot(ldexp(ax, DBL_MANT_DIG),
    643                      ldexp(ay, DBL_MANT_DIG))) - DBL_MANT_DIG*M_LN2;
    644         }
    645         else {
    646             /* log(+/-0. +/- 0i) */
    647             r.real = -INF;
    648             r.imag = atan2(z.imag, z.real);
    649             errno = EDOM;
    650             return r;
    651         }
    652     } else {
    653         h = hypot(ax, ay);
    654         if (0.71 <= h && h <= 1.73) {
    655             am = ax > ay ? ax : ay;  /* max(ax, ay) */
    656             an = ax > ay ? ay : ax;  /* min(ax, ay) */
    657             r.real = m_log1p((am-1)*(am+1)+an*an)/2.;
    658         } else {
    659             r.real = log(h);
    660         }
    661     }
    662     r.imag = atan2(z.imag, z.real);
    663     errno = 0;
    664     return r;
    665 }
    666 
    667 
    668 /*[clinic input]
    669 cmath.log10 = cmath.acos
    670 
    671 Return the base-10 logarithm of z.
    672 [clinic start generated code]*/
    673 
    674 static Py_complex
    675 cmath_log10_impl(PyObject *module, Py_complex z)
    676 /*[clinic end generated code: output=2922779a7c38cbe1 input=cff5644f73c1519c]*/
    677 {
    678     Py_complex r;
    679     int errno_save;
    680 
    681     r = c_log(z);
    682     errno_save = errno; /* just in case the divisions affect errno */
    683     r.real = r.real / M_LN10;
    684     r.imag = r.imag / M_LN10;
    685     errno = errno_save;
    686     return r;
    687 }
    688 
    689 
    690 /*[clinic input]
    691 cmath.sin = cmath.acos
    692 
    693 Return the sine of z.
    694 [clinic start generated code]*/
    695 
    696 static Py_complex
    697 cmath_sin_impl(PyObject *module, Py_complex z)
    698 /*[clinic end generated code: output=980370d2ff0bb5aa input=2d3519842a8b4b85]*/
    699 {
    700     /* sin(z) = -i sin(iz) */
    701     Py_complex s, r;
    702     s.real = -z.imag;
    703     s.imag = z.real;
    704     s = cmath_sinh_impl(module, s);
    705     r.real = s.imag;
    706     r.imag = -s.real;
    707     return r;
    708 }
    709 
    710 
    711 /* sinh(infinity + i*y) needs to be dealt with specially */
    712 static Py_complex sinh_special_values[7][7];
    713 
    714 /*[clinic input]
    715 cmath.sinh = cmath.acos
    716 
    717 Return the hyperbolic sine of z.
    718 [clinic start generated code]*/
    719 
    720 static Py_complex
    721 cmath_sinh_impl(PyObject *module, Py_complex z)
    722 /*[clinic end generated code: output=38b0a6cce26f3536 input=d2d3fc8c1ddfd2dd]*/
    723 {
    724     Py_complex r;
    725     double x_minus_one;
    726 
    727     /* special treatment for sinh(+/-inf + iy) if y is finite and
    728        nonzero */
    729     if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
    730         if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
    731             && (z.imag != 0.)) {
    732             if (z.real > 0) {
    733                 r.real = copysign(INF, cos(z.imag));
    734                 r.imag = copysign(INF, sin(z.imag));
    735             }
    736             else {
    737                 r.real = -copysign(INF, cos(z.imag));
    738                 r.imag = copysign(INF, sin(z.imag));
    739             }
    740         }
    741         else {
    742             r = sinh_special_values[special_type(z.real)]
    743                                    [special_type(z.imag)];
    744         }
    745         /* need to set errno = EDOM if y is +/- infinity and x is not
    746            a NaN */
    747         if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
    748             errno = EDOM;
    749         else
    750             errno = 0;
    751         return r;
    752     }
    753 
    754     if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
    755         x_minus_one = z.real - copysign(1., z.real);
    756         r.real = cos(z.imag) * sinh(x_minus_one) * Py_MATH_E;
    757         r.imag = sin(z.imag) * cosh(x_minus_one) * Py_MATH_E;
    758     } else {
    759         r.real = cos(z.imag) * sinh(z.real);
    760         r.imag = sin(z.imag) * cosh(z.real);
    761     }
    762     /* detect overflow, and set errno accordingly */
    763     if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
    764         errno = ERANGE;
    765     else
    766         errno = 0;
    767     return r;
    768 }
    769 
    770 
    771 static Py_complex sqrt_special_values[7][7];
    772 
    773 /*[clinic input]
    774 cmath.sqrt = cmath.acos
    775 
    776 Return the square root of z.
    777 [clinic start generated code]*/
    778 
    779 static Py_complex
    780 cmath_sqrt_impl(PyObject *module, Py_complex z)
    781 /*[clinic end generated code: output=b6507b3029c339fc input=7088b166fc9a58c7]*/
    782 {
    783     /*
    784        Method: use symmetries to reduce to the case when x = z.real and y
    785        = z.imag are nonnegative.  Then the real part of the result is
    786        given by
    787 
    788          s = sqrt((x + hypot(x, y))/2)
    789 
    790        and the imaginary part is
    791 
    792          d = (y/2)/s
    793 
    794        If either x or y is very large then there's a risk of overflow in
    795        computation of the expression x + hypot(x, y).  We can avoid this
    796        by rewriting the formula for s as:
    797 
    798          s = 2*sqrt(x/8 + hypot(x/8, y/8))
    799 
    800        This costs us two extra multiplications/divisions, but avoids the
    801        overhead of checking for x and y large.
    802 
    803        If both x and y are subnormal then hypot(x, y) may also be
    804        subnormal, so will lack full precision.  We solve this by rescaling
    805        x and y by a sufficiently large power of 2 to ensure that x and y
    806        are normal.
    807     */
    808 
    809 
    810     Py_complex r;
    811     double s,d;
    812     double ax, ay;
    813 
    814     SPECIAL_VALUE(z, sqrt_special_values);
    815 
    816     if (z.real == 0. && z.imag == 0.) {
    817         r.real = 0.;
    818         r.imag = z.imag;
    819         return r;
    820     }
    821 
    822     ax = fabs(z.real);
    823     ay = fabs(z.imag);
    824 
    825     if (ax < DBL_MIN && ay < DBL_MIN && (ax > 0. || ay > 0.)) {
    826         /* here we catch cases where hypot(ax, ay) is subnormal */
    827         ax = ldexp(ax, CM_SCALE_UP);
    828         s = ldexp(sqrt(ax + hypot(ax, ldexp(ay, CM_SCALE_UP))),
    829                   CM_SCALE_DOWN);
    830     } else {
    831         ax /= 8.;
    832         s = 2.*sqrt(ax + hypot(ax, ay/8.));
    833     }
    834     d = ay/(2.*s);
    835 
    836     if (z.real >= 0.) {
    837         r.real = s;
    838         r.imag = copysign(d, z.imag);
    839     } else {
    840         r.real = d;
    841         r.imag = copysign(s, z.imag);
    842     }
    843     errno = 0;
    844     return r;
    845 }
    846 
    847 
    848 /*[clinic input]
    849 cmath.tan = cmath.acos
    850 
    851 Return the tangent of z.
    852 [clinic start generated code]*/
    853 
    854 static Py_complex
    855 cmath_tan_impl(PyObject *module, Py_complex z)
    856 /*[clinic end generated code: output=7c5f13158a72eb13 input=fc167e528767888e]*/
    857 {
    858     /* tan(z) = -i tanh(iz) */
    859     Py_complex s, r;
    860     s.real = -z.imag;
    861     s.imag = z.real;
    862     s = cmath_tanh_impl(module, s);
    863     r.real = s.imag;
    864     r.imag = -s.real;
    865     return r;
    866 }
    867 
    868 
    869 /* tanh(infinity + i*y) needs to be dealt with specially */
    870 static Py_complex tanh_special_values[7][7];
    871 
    872 /*[clinic input]
    873 cmath.tanh = cmath.acos
    874 
    875 Return the hyperbolic tangent of z.
    876 [clinic start generated code]*/
    877 
    878 static Py_complex
    879 cmath_tanh_impl(PyObject *module, Py_complex z)
    880 /*[clinic end generated code: output=36d547ef7aca116c input=22f67f9dc6d29685]*/
    881 {
    882     /* Formula:
    883 
    884        tanh(x+iy) = (tanh(x)(1+tan(y)^2) + i tan(y)(1-tanh(x))^2) /
    885        (1+tan(y)^2 tanh(x)^2)
    886 
    887        To avoid excessive roundoff error, 1-tanh(x)^2 is better computed
    888        as 1/cosh(x)^2.  When abs(x) is large, we approximate 1-tanh(x)^2
    889        by 4 exp(-2*x) instead, to avoid possible overflow in the
    890        computation of cosh(x).
    891 
    892     */
    893 
    894     Py_complex r;
    895     double tx, ty, cx, txty, denom;
    896 
    897     /* special treatment for tanh(+/-inf + iy) if y is finite and
    898        nonzero */
    899     if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
    900         if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
    901             && (z.imag != 0.)) {
    902             if (z.real > 0) {
    903                 r.real = 1.0;
    904                 r.imag = copysign(0.,
    905                                   2.*sin(z.imag)*cos(z.imag));
    906             }
    907             else {
    908                 r.real = -1.0;
    909                 r.imag = copysign(0.,
    910                                   2.*sin(z.imag)*cos(z.imag));
    911             }
    912         }
    913         else {
    914             r = tanh_special_values[special_type(z.real)]
    915                                    [special_type(z.imag)];
    916         }
    917         /* need to set errno = EDOM if z.imag is +/-infinity and
    918            z.real is finite */
    919         if (Py_IS_INFINITY(z.imag) && Py_IS_FINITE(z.real))
    920             errno = EDOM;
    921         else
    922             errno = 0;
    923         return r;
    924     }
    925 
    926     /* danger of overflow in 2.*z.imag !*/
    927     if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
    928         r.real = copysign(1., z.real);
    929         r.imag = 4.*sin(z.imag)*cos(z.imag)*exp(-2.*fabs(z.real));
    930     } else {
    931         tx = tanh(z.real);
    932         ty = tan(z.imag);
    933         cx = 1./cosh(z.real);
    934         txty = tx*ty;
    935         denom = 1. + txty*txty;
    936         r.real = tx*(1.+ty*ty)/denom;
    937         r.imag = ((ty/denom)*cx)*cx;
    938     }
    939     errno = 0;
    940     return r;
    941 }
    942 
    943 
    944 /*[clinic input]
    945 cmath.log
    946 
    947     x: Py_complex
    948     y_obj: object = NULL
    949     /
    950 
    951 The logarithm of z to the given base.
    952 
    953 If the base not specified, returns the natural logarithm (base e) of z.
    954 [clinic start generated code]*/
    955 
    956 static PyObject *
    957 cmath_log_impl(PyObject *module, Py_complex x, PyObject *y_obj)
    958 /*[clinic end generated code: output=4effdb7d258e0d94 input=ee0e823a7c6e68ea]*/
    959 {
    960     Py_complex y;
    961 
    962     errno = 0;
    963     PyFPE_START_PROTECT("complex function", return 0)
    964     x = c_log(x);
    965     if (y_obj != NULL) {
    966         y = PyComplex_AsCComplex(y_obj);
    967         if (PyErr_Occurred()) {
    968             return NULL;
    969         }
    970         y = c_log(y);
    971         x = _Py_c_quot(x, y);
    972     }
    973     PyFPE_END_PROTECT(x)
    974     if (errno != 0)
    975         return math_error();
    976     return PyComplex_FromCComplex(x);
    977 }
    978 
    979 
    980 /* And now the glue to make them available from Python: */
    981 
    982 static PyObject *
    983 math_error(void)
    984 {
    985     if (errno == EDOM)
    986         PyErr_SetString(PyExc_ValueError, "math domain error");
    987     else if (errno == ERANGE)
    988         PyErr_SetString(PyExc_OverflowError, "math range error");
    989     else    /* Unexpected math error */
    990         PyErr_SetFromErrno(PyExc_ValueError);
    991     return NULL;
    992 }
    993 
    994 
    995 /*[clinic input]
    996 cmath.phase
    997 
    998     z: Py_complex
    999     /
   1000 
   1001 Return argument, also known as the phase angle, of a complex.
   1002 [clinic start generated code]*/
   1003 
   1004 static PyObject *
   1005 cmath_phase_impl(PyObject *module, Py_complex z)
   1006 /*[clinic end generated code: output=50725086a7bfd253 input=5cf75228ba94b69d]*/
   1007 {
   1008     double phi;
   1009 
   1010     errno = 0;
   1011     PyFPE_START_PROTECT("arg function", return 0)
   1012     phi = c_atan2(z);
   1013     PyFPE_END_PROTECT(phi)
   1014     if (errno != 0)
   1015         return math_error();
   1016     else
   1017         return PyFloat_FromDouble(phi);
   1018 }
   1019 
   1020 /*[clinic input]
   1021 cmath.polar
   1022 
   1023     z: Py_complex
   1024     /
   1025 
   1026 Convert a complex from rectangular coordinates to polar coordinates.
   1027 
   1028 r is the distance from 0 and phi the phase angle.
   1029 [clinic start generated code]*/
   1030 
   1031 static PyObject *
   1032 cmath_polar_impl(PyObject *module, Py_complex z)
   1033 /*[clinic end generated code: output=d0a8147c41dbb654 input=26c353574fd1a861]*/
   1034 {
   1035     double r, phi;
   1036 
   1037     errno = 0;
   1038     PyFPE_START_PROTECT("polar function", return 0)
   1039     phi = c_atan2(z); /* should not cause any exception */
   1040     r = _Py_c_abs(z); /* sets errno to ERANGE on overflow */
   1041     PyFPE_END_PROTECT(r)
   1042     if (errno != 0)
   1043         return math_error();
   1044     else
   1045         return Py_BuildValue("dd", r, phi);
   1046 }
   1047 
   1048 /*
   1049   rect() isn't covered by the C99 standard, but it's not too hard to
   1050   figure out 'spirit of C99' rules for special value handing:
   1051 
   1052     rect(x, t) should behave like exp(log(x) + it) for positive-signed x
   1053     rect(x, t) should behave like -exp(log(-x) + it) for negative-signed x
   1054     rect(nan, t) should behave like exp(nan + it), except that rect(nan, 0)
   1055       gives nan +- i0 with the sign of the imaginary part unspecified.
   1056 
   1057 */
   1058 
   1059 static Py_complex rect_special_values[7][7];
   1060 
   1061 /*[clinic input]
   1062 cmath.rect
   1063 
   1064     r: double
   1065     phi: double
   1066     /
   1067 
   1068 Convert from polar coordinates to rectangular coordinates.
   1069 [clinic start generated code]*/
   1070 
   1071 static PyObject *
   1072 cmath_rect_impl(PyObject *module, double r, double phi)
   1073 /*[clinic end generated code: output=385a0690925df2d5 input=24c5646d147efd69]*/
   1074 {
   1075     Py_complex z;
   1076     errno = 0;
   1077     PyFPE_START_PROTECT("rect function", return 0)
   1078 
   1079     /* deal with special values */
   1080     if (!Py_IS_FINITE(r) || !Py_IS_FINITE(phi)) {
   1081         /* if r is +/-infinity and phi is finite but nonzero then
   1082            result is (+-INF +-INF i), but we need to compute cos(phi)
   1083            and sin(phi) to figure out the signs. */
   1084         if (Py_IS_INFINITY(r) && (Py_IS_FINITE(phi)
   1085                                   && (phi != 0.))) {
   1086             if (r > 0) {
   1087                 z.real = copysign(INF, cos(phi));
   1088                 z.imag = copysign(INF, sin(phi));
   1089             }
   1090             else {
   1091                 z.real = -copysign(INF, cos(phi));
   1092                 z.imag = -copysign(INF, sin(phi));
   1093             }
   1094         }
   1095         else {
   1096             z = rect_special_values[special_type(r)]
   1097                                    [special_type(phi)];
   1098         }
   1099         /* need to set errno = EDOM if r is a nonzero number and phi
   1100            is infinite */
   1101         if (r != 0. && !Py_IS_NAN(r) && Py_IS_INFINITY(phi))
   1102             errno = EDOM;
   1103         else
   1104             errno = 0;
   1105     }
   1106     else if (phi == 0.0) {
   1107         /* Workaround for buggy results with phi=-0.0 on OS X 10.8.  See
   1108            bugs.python.org/issue18513. */
   1109         z.real = r;
   1110         z.imag = r * phi;
   1111         errno = 0;
   1112     }
   1113     else {
   1114         z.real = r * cos(phi);
   1115         z.imag = r * sin(phi);
   1116         errno = 0;
   1117     }
   1118 
   1119     PyFPE_END_PROTECT(z)
   1120     if (errno != 0)
   1121         return math_error();
   1122     else
   1123         return PyComplex_FromCComplex(z);
   1124 }
   1125 
   1126 /*[clinic input]
   1127 cmath.isfinite = cmath.polar
   1128 
   1129 Return True if both the real and imaginary parts of z are finite, else False.
   1130 [clinic start generated code]*/
   1131 
   1132 static PyObject *
   1133 cmath_isfinite_impl(PyObject *module, Py_complex z)
   1134 /*[clinic end generated code: output=ac76611e2c774a36 input=848e7ee701895815]*/
   1135 {
   1136     return PyBool_FromLong(Py_IS_FINITE(z.real) && Py_IS_FINITE(z.imag));
   1137 }
   1138 
   1139 /*[clinic input]
   1140 cmath.isnan = cmath.polar
   1141 
   1142 Checks if the real or imaginary part of z not a number (NaN).
   1143 [clinic start generated code]*/
   1144 
   1145 static PyObject *
   1146 cmath_isnan_impl(PyObject *module, Py_complex z)
   1147 /*[clinic end generated code: output=e7abf6e0b28beab7 input=71799f5d284c9baf]*/
   1148 {
   1149     return PyBool_FromLong(Py_IS_NAN(z.real) || Py_IS_NAN(z.imag));
   1150 }
   1151 
   1152 /*[clinic input]
   1153 cmath.isinf = cmath.polar
   1154 
   1155 Checks if the real or imaginary part of z is infinite.
   1156 [clinic start generated code]*/
   1157 
   1158 static PyObject *
   1159 cmath_isinf_impl(PyObject *module, Py_complex z)
   1160 /*[clinic end generated code: output=502a75a79c773469 input=363df155c7181329]*/
   1161 {
   1162     return PyBool_FromLong(Py_IS_INFINITY(z.real) ||
   1163                            Py_IS_INFINITY(z.imag));
   1164 }
   1165 
   1166 /*[clinic input]
   1167 cmath.isclose -> bool
   1168 
   1169     a: Py_complex
   1170     b: Py_complex
   1171     *
   1172     rel_tol: double = 1e-09
   1173         maximum difference for being considered "close", relative to the
   1174         magnitude of the input values
   1175     abs_tol: double = 0.0
   1176         maximum difference for being considered "close", regardless of the
   1177         magnitude of the input values
   1178 
   1179 Determine whether two complex numbers are close in value.
   1180 
   1181 Return True if a is close in value to b, and False otherwise.
   1182 
   1183 For the values to be considered close, the difference between them must be
   1184 smaller than at least one of the tolerances.
   1185 
   1186 -inf, inf and NaN behave similarly to the IEEE 754 Standard. That is, NaN is
   1187 not close to anything, even itself. inf and -inf are only close to themselves.
   1188 [clinic start generated code]*/
   1189 
   1190 static int
   1191 cmath_isclose_impl(PyObject *module, Py_complex a, Py_complex b,
   1192                    double rel_tol, double abs_tol)
   1193 /*[clinic end generated code: output=8a2486cc6e0014d1 input=df9636d7de1d4ac3]*/
   1194 {
   1195     double diff;
   1196 
   1197     /* sanity check on the inputs */
   1198     if (rel_tol < 0.0 || abs_tol < 0.0 ) {
   1199         PyErr_SetString(PyExc_ValueError,
   1200                         "tolerances must be non-negative");
   1201         return -1;
   1202     }
   1203 
   1204     if ( (a.real == b.real) && (a.imag == b.imag) ) {
   1205         /* short circuit exact equality -- needed to catch two infinities of
   1206            the same sign. And perhaps speeds things up a bit sometimes.
   1207         */
   1208         return 1;
   1209     }
   1210 
   1211     /* This catches the case of two infinities of opposite sign, or
   1212        one infinity and one finite number. Two infinities of opposite
   1213        sign would otherwise have an infinite relative tolerance.
   1214        Two infinities of the same sign are caught by the equality check
   1215        above.
   1216     */
   1217 
   1218     if (Py_IS_INFINITY(a.real) || Py_IS_INFINITY(a.imag) ||
   1219         Py_IS_INFINITY(b.real) || Py_IS_INFINITY(b.imag)) {
   1220         return 0;
   1221     }
   1222 
   1223     /* now do the regular computation
   1224        this is essentially the "weak" test from the Boost library
   1225     */
   1226 
   1227     diff = _Py_c_abs(_Py_c_diff(a, b));
   1228 
   1229     return (((diff <= rel_tol * _Py_c_abs(b)) ||
   1230              (diff <= rel_tol * _Py_c_abs(a))) ||
   1231             (diff <= abs_tol));
   1232 }
   1233 
   1234 PyDoc_STRVAR(module_doc,
   1235 "This module is always available. It provides access to mathematical\n"
   1236 "functions for complex numbers.");
   1237 
   1238 static PyMethodDef cmath_methods[] = {
   1239     CMATH_ACOS_METHODDEF
   1240     CMATH_ACOSH_METHODDEF
   1241     CMATH_ASIN_METHODDEF
   1242     CMATH_ASINH_METHODDEF
   1243     CMATH_ATAN_METHODDEF
   1244     CMATH_ATANH_METHODDEF
   1245     CMATH_COS_METHODDEF
   1246     CMATH_COSH_METHODDEF
   1247     CMATH_EXP_METHODDEF
   1248     CMATH_ISCLOSE_METHODDEF
   1249     CMATH_ISFINITE_METHODDEF
   1250     CMATH_ISINF_METHODDEF
   1251     CMATH_ISNAN_METHODDEF
   1252     CMATH_LOG_METHODDEF
   1253     CMATH_LOG10_METHODDEF
   1254     CMATH_PHASE_METHODDEF
   1255     CMATH_POLAR_METHODDEF
   1256     CMATH_RECT_METHODDEF
   1257     CMATH_SIN_METHODDEF
   1258     CMATH_SINH_METHODDEF
   1259     CMATH_SQRT_METHODDEF
   1260     CMATH_TAN_METHODDEF
   1261     CMATH_TANH_METHODDEF
   1262     {NULL, NULL}  /* sentinel */
   1263 };
   1264 
   1265 
   1266 static struct PyModuleDef cmathmodule = {
   1267     PyModuleDef_HEAD_INIT,
   1268     "cmath",
   1269     module_doc,
   1270     -1,
   1271     cmath_methods,
   1272     NULL,
   1273     NULL,
   1274     NULL,
   1275     NULL
   1276 };
   1277 
   1278 PyMODINIT_FUNC
   1279 PyInit_cmath(void)
   1280 {
   1281     PyObject *m;
   1282 
   1283     m = PyModule_Create(&cmathmodule);
   1284     if (m == NULL)
   1285         return NULL;
   1286 
   1287     PyModule_AddObject(m, "pi",
   1288                        PyFloat_FromDouble(Py_MATH_PI));
   1289     PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
   1290     PyModule_AddObject(m, "tau", PyFloat_FromDouble(Py_MATH_TAU)); /* 2pi */
   1291     PyModule_AddObject(m, "inf", PyFloat_FromDouble(m_inf()));
   1292     PyModule_AddObject(m, "infj", PyComplex_FromCComplex(c_infj()));
   1293 #if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN)
   1294     PyModule_AddObject(m, "nan", PyFloat_FromDouble(m_nan()));
   1295     PyModule_AddObject(m, "nanj", PyComplex_FromCComplex(c_nanj()));
   1296 #endif
   1297 
   1298     /* initialize special value tables */
   1299 
   1300 #define INIT_SPECIAL_VALUES(NAME, BODY) { Py_complex* p = (Py_complex*)NAME; BODY }
   1301 #define C(REAL, IMAG) p->real = REAL; p->imag = IMAG; ++p;
   1302 
   1303     INIT_SPECIAL_VALUES(acos_special_values, {
   1304       C(P34,INF) C(P,INF)  C(P,INF)  C(P,-INF)  C(P,-INF)  C(P34,-INF) C(N,INF)
   1305       C(P12,INF) C(U,U)    C(U,U)    C(U,U)     C(U,U)     C(P12,-INF) C(N,N)
   1306       C(P12,INF) C(U,U)    C(P12,0.) C(P12,-0.) C(U,U)     C(P12,-INF) C(P12,N)
   1307       C(P12,INF) C(U,U)    C(P12,0.) C(P12,-0.) C(U,U)     C(P12,-INF) C(P12,N)
   1308       C(P12,INF) C(U,U)    C(U,U)    C(U,U)     C(U,U)     C(P12,-INF) C(N,N)
   1309       C(P14,INF) C(0.,INF) C(0.,INF) C(0.,-INF) C(0.,-INF) C(P14,-INF) C(N,INF)
   1310       C(N,INF)   C(N,N)    C(N,N)    C(N,N)     C(N,N)     C(N,-INF)   C(N,N)
   1311     })
   1312 
   1313     INIT_SPECIAL_VALUES(acosh_special_values, {
   1314       C(INF,-P34) C(INF,-P)  C(INF,-P)  C(INF,P)  C(INF,P)  C(INF,P34) C(INF,N)
   1315       C(INF,-P12) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,P12) C(N,N)
   1316       C(INF,-P12) C(U,U)     C(0.,-P12) C(0.,P12) C(U,U)    C(INF,P12) C(N,N)
   1317       C(INF,-P12) C(U,U)     C(0.,-P12) C(0.,P12) C(U,U)    C(INF,P12) C(N,N)
   1318       C(INF,-P12) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,P12) C(N,N)
   1319       C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
   1320       C(INF,N)    C(N,N)     C(N,N)     C(N,N)    C(N,N)    C(INF,N)   C(N,N)
   1321     })
   1322 
   1323     INIT_SPECIAL_VALUES(asinh_special_values, {
   1324       C(-INF,-P14) C(-INF,-0.) C(-INF,-0.) C(-INF,0.) C(-INF,0.) C(-INF,P14) C(-INF,N)
   1325       C(-INF,-P12) C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(-INF,P12) C(N,N)
   1326       C(-INF,-P12) C(U,U)      C(-0.,-0.)  C(-0.,0.)  C(U,U)     C(-INF,P12) C(N,N)
   1327       C(INF,-P12)  C(U,U)      C(0.,-0.)   C(0.,0.)   C(U,U)     C(INF,P12)  C(N,N)
   1328       C(INF,-P12)  C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(INF,P12)  C(N,N)
   1329       C(INF,-P14)  C(INF,-0.)  C(INF,-0.)  C(INF,0.)  C(INF,0.)  C(INF,P14)  C(INF,N)
   1330       C(INF,N)     C(N,N)      C(N,-0.)    C(N,0.)    C(N,N)     C(INF,N)    C(N,N)
   1331     })
   1332 
   1333     INIT_SPECIAL_VALUES(atanh_special_values, {
   1334       C(-0.,-P12) C(-0.,-P12) C(-0.,-P12) C(-0.,P12) C(-0.,P12) C(-0.,P12) C(-0.,N)
   1335       C(-0.,-P12) C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(-0.,P12) C(N,N)
   1336       C(-0.,-P12) C(U,U)      C(-0.,-0.)  C(-0.,0.)  C(U,U)     C(-0.,P12) C(-0.,N)
   1337       C(0.,-P12)  C(U,U)      C(0.,-0.)   C(0.,0.)   C(U,U)     C(0.,P12)  C(0.,N)
   1338       C(0.,-P12)  C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(0.,P12)  C(N,N)
   1339       C(0.,-P12)  C(0.,-P12)  C(0.,-P12)  C(0.,P12)  C(0.,P12)  C(0.,P12)  C(0.,N)
   1340       C(0.,-P12)  C(N,N)      C(N,N)      C(N,N)     C(N,N)     C(0.,P12)  C(N,N)
   1341     })
   1342 
   1343     INIT_SPECIAL_VALUES(cosh_special_values, {
   1344       C(INF,N) C(U,U) C(INF,0.)  C(INF,-0.) C(U,U) C(INF,N) C(INF,N)
   1345       C(N,N)   C(U,U) C(U,U)     C(U,U)     C(U,U) C(N,N)   C(N,N)
   1346       C(N,0.)  C(U,U) C(1.,0.)   C(1.,-0.)  C(U,U) C(N,0.)  C(N,0.)
   1347       C(N,0.)  C(U,U) C(1.,-0.)  C(1.,0.)   C(U,U) C(N,0.)  C(N,0.)
   1348       C(N,N)   C(U,U) C(U,U)     C(U,U)     C(U,U) C(N,N)   C(N,N)
   1349       C(INF,N) C(U,U) C(INF,-0.) C(INF,0.)  C(U,U) C(INF,N) C(INF,N)
   1350       C(N,N)   C(N,N) C(N,0.)    C(N,0.)    C(N,N) C(N,N)   C(N,N)
   1351     })
   1352 
   1353     INIT_SPECIAL_VALUES(exp_special_values, {
   1354       C(0.,0.) C(U,U) C(0.,-0.)  C(0.,0.)  C(U,U) C(0.,0.) C(0.,0.)
   1355       C(N,N)   C(U,U) C(U,U)     C(U,U)    C(U,U) C(N,N)   C(N,N)
   1356       C(N,N)   C(U,U) C(1.,-0.)  C(1.,0.)  C(U,U) C(N,N)   C(N,N)
   1357       C(N,N)   C(U,U) C(1.,-0.)  C(1.,0.)  C(U,U) C(N,N)   C(N,N)
   1358       C(N,N)   C(U,U) C(U,U)     C(U,U)    C(U,U) C(N,N)   C(N,N)
   1359       C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
   1360       C(N,N)   C(N,N) C(N,-0.)   C(N,0.)   C(N,N) C(N,N)   C(N,N)
   1361     })
   1362 
   1363     INIT_SPECIAL_VALUES(log_special_values, {
   1364       C(INF,-P34) C(INF,-P)  C(INF,-P)   C(INF,P)   C(INF,P)  C(INF,P34)  C(INF,N)
   1365       C(INF,-P12) C(U,U)     C(U,U)      C(U,U)     C(U,U)    C(INF,P12)  C(N,N)
   1366       C(INF,-P12) C(U,U)     C(-INF,-P)  C(-INF,P)  C(U,U)    C(INF,P12)  C(N,N)
   1367       C(INF,-P12) C(U,U)     C(-INF,-0.) C(-INF,0.) C(U,U)    C(INF,P12)  C(N,N)
   1368       C(INF,-P12) C(U,U)     C(U,U)      C(U,U)     C(U,U)    C(INF,P12)  C(N,N)
   1369       C(INF,-P14) C(INF,-0.) C(INF,-0.)  C(INF,0.)  C(INF,0.) C(INF,P14)  C(INF,N)
   1370       C(INF,N)    C(N,N)     C(N,N)      C(N,N)     C(N,N)    C(INF,N)    C(N,N)
   1371     })
   1372 
   1373     INIT_SPECIAL_VALUES(sinh_special_values, {
   1374       C(INF,N) C(U,U) C(-INF,-0.) C(-INF,0.) C(U,U) C(INF,N) C(INF,N)
   1375       C(N,N)   C(U,U) C(U,U)      C(U,U)     C(U,U) C(N,N)   C(N,N)
   1376       C(0.,N)  C(U,U) C(-0.,-0.)  C(-0.,0.)  C(U,U) C(0.,N)  C(0.,N)
   1377       C(0.,N)  C(U,U) C(0.,-0.)   C(0.,0.)   C(U,U) C(0.,N)  C(0.,N)
   1378       C(N,N)   C(U,U) C(U,U)      C(U,U)     C(U,U) C(N,N)   C(N,N)
   1379       C(INF,N) C(U,U) C(INF,-0.)  C(INF,0.)  C(U,U) C(INF,N) C(INF,N)
   1380       C(N,N)   C(N,N) C(N,-0.)    C(N,0.)    C(N,N) C(N,N)   C(N,N)
   1381     })
   1382 
   1383     INIT_SPECIAL_VALUES(sqrt_special_values, {
   1384       C(INF,-INF) C(0.,-INF) C(0.,-INF) C(0.,INF) C(0.,INF) C(INF,INF) C(N,INF)
   1385       C(INF,-INF) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,INF) C(N,N)
   1386       C(INF,-INF) C(U,U)     C(0.,-0.)  C(0.,0.)  C(U,U)    C(INF,INF) C(N,N)
   1387       C(INF,-INF) C(U,U)     C(0.,-0.)  C(0.,0.)  C(U,U)    C(INF,INF) C(N,N)
   1388       C(INF,-INF) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,INF) C(N,N)
   1389       C(INF,-INF) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,INF) C(INF,N)
   1390       C(INF,-INF) C(N,N)     C(N,N)     C(N,N)    C(N,N)    C(INF,INF) C(N,N)
   1391     })
   1392 
   1393     INIT_SPECIAL_VALUES(tanh_special_values, {
   1394       C(-1.,0.) C(U,U) C(-1.,-0.) C(-1.,0.) C(U,U) C(-1.,0.) C(-1.,0.)
   1395       C(N,N)    C(U,U) C(U,U)     C(U,U)    C(U,U) C(N,N)    C(N,N)
   1396       C(N,N)    C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(N,N)    C(N,N)
   1397       C(N,N)    C(U,U) C(0.,-0.)  C(0.,0.)  C(U,U) C(N,N)    C(N,N)
   1398       C(N,N)    C(U,U) C(U,U)     C(U,U)    C(U,U) C(N,N)    C(N,N)
   1399       C(1.,0.)  C(U,U) C(1.,-0.)  C(1.,0.)  C(U,U) C(1.,0.)  C(1.,0.)
   1400       C(N,N)    C(N,N) C(N,-0.)   C(N,0.)   C(N,N) C(N,N)    C(N,N)
   1401     })
   1402 
   1403     INIT_SPECIAL_VALUES(rect_special_values, {
   1404       C(INF,N) C(U,U) C(-INF,0.) C(-INF,-0.) C(U,U) C(INF,N) C(INF,N)
   1405       C(N,N)   C(U,U) C(U,U)     C(U,U)      C(U,U) C(N,N)   C(N,N)
   1406       C(0.,0.) C(U,U) C(-0.,0.)  C(-0.,-0.)  C(U,U) C(0.,0.) C(0.,0.)
   1407       C(0.,0.) C(U,U) C(0.,-0.)  C(0.,0.)    C(U,U) C(0.,0.) C(0.,0.)
   1408       C(N,N)   C(U,U) C(U,U)     C(U,U)      C(U,U) C(N,N)   C(N,N)
   1409       C(INF,N) C(U,U) C(INF,-0.) C(INF,0.)   C(U,U) C(INF,N) C(INF,N)
   1410       C(N,N)   C(N,N) C(N,0.)    C(N,0.)     C(N,N) C(N,N)   C(N,N)
   1411     })
   1412     return m;
   1413 }
   1414