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