• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /* Complex math module */
2 
3 /* much code borrowed from mathmodule.c */
4 
5 #include "Python.h"
6 #include "pycore_dtoa.h"
7 #include "_math.h"
8 /* we need DBL_MAX, DBL_MIN, DBL_EPSILON, DBL_MANT_DIG and FLT_RADIX from
9    float.h.  We assume that FLT_RADIX is either 2 or 16. */
10 #include <float.h>
11 
12 #include "clinic/cmathmodule.c.h"
13 /*[clinic input]
14 module cmath
15 [clinic start generated code]*/
16 /*[clinic end generated code: output=da39a3ee5e6b4b0d input=308d6839f4a46333]*/
17 
18 /*[python input]
19 class Py_complex_protected_converter(Py_complex_converter):
20     def modify(self):
21         return 'errno = 0;'
22 
23 
24 class Py_complex_protected_return_converter(CReturnConverter):
25     type = "Py_complex"
26 
27     def render(self, function, data):
28         self.declare(data)
29         data.return_conversion.append("""
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=8b27adb674c08321]*/
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
m_inf(void)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
c_infj(void)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
m_nan(void)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
c_nanj(void)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
special_type(double d)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
cmath_acos_impl(PyObject * module,Py_complex z)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
cmath_acosh_impl(PyObject * module,Py_complex z)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
cmath_asin_impl(PyObject * module,Py_complex z)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
cmath_asinh_impl(PyObject * module,Py_complex z)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
cmath_atan_impl(PyObject * module,Py_complex z)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
c_atan2(Py_complex z)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
cmath_atanh_impl(PyObject * module,Py_complex z)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
cmath_cos_impl(PyObject * module,Py_complex z)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
cmath_cosh_impl(PyObject * module,Py_complex z)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
cmath_exp_impl(PyObject * module,Py_complex z)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
c_log(Py_complex z)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
cmath_log10_impl(PyObject * module,Py_complex z)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
cmath_sin_impl(PyObject * module,Py_complex z)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
cmath_sinh_impl(PyObject * module,Py_complex z)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
cmath_sqrt_impl(PyObject * module,Py_complex z)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
cmath_tan_impl(PyObject * module,Py_complex z)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
cmath_tanh_impl(PyObject * module,Py_complex z)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     z as x: Py_complex
948     base as y_obj: object = NULL
949     /
950 
951 log(z[, base]) -> 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 *
cmath_log_impl(PyObject * module,Py_complex x,PyObject * y_obj)957 cmath_log_impl(PyObject *module, Py_complex x, PyObject *y_obj)
958 /*[clinic end generated code: output=4effdb7d258e0d94 input=230ed3a71ecd000a]*/
959 {
960     Py_complex y;
961 
962     errno = 0;
963     x = c_log(x);
964     if (y_obj != NULL) {
965         y = PyComplex_AsCComplex(y_obj);
966         if (PyErr_Occurred()) {
967             return NULL;
968         }
969         y = c_log(y);
970         x = _Py_c_quot(x, y);
971     }
972     if (errno != 0)
973         return math_error();
974     return PyComplex_FromCComplex(x);
975 }
976 
977 
978 /* And now the glue to make them available from Python: */
979 
980 static PyObject *
math_error(void)981 math_error(void)
982 {
983     if (errno == EDOM)
984         PyErr_SetString(PyExc_ValueError, "math domain error");
985     else if (errno == ERANGE)
986         PyErr_SetString(PyExc_OverflowError, "math range error");
987     else    /* Unexpected math error */
988         PyErr_SetFromErrno(PyExc_ValueError);
989     return NULL;
990 }
991 
992 
993 /*[clinic input]
994 cmath.phase
995 
996     z: Py_complex
997     /
998 
999 Return argument, also known as the phase angle, of a complex.
1000 [clinic start generated code]*/
1001 
1002 static PyObject *
cmath_phase_impl(PyObject * module,Py_complex z)1003 cmath_phase_impl(PyObject *module, Py_complex z)
1004 /*[clinic end generated code: output=50725086a7bfd253 input=5cf75228ba94b69d]*/
1005 {
1006     double phi;
1007 
1008     errno = 0;
1009     phi = c_atan2(z);
1010     if (errno != 0)
1011         return math_error();
1012     else
1013         return PyFloat_FromDouble(phi);
1014 }
1015 
1016 /*[clinic input]
1017 cmath.polar
1018 
1019     z: Py_complex
1020     /
1021 
1022 Convert a complex from rectangular coordinates to polar coordinates.
1023 
1024 r is the distance from 0 and phi the phase angle.
1025 [clinic start generated code]*/
1026 
1027 static PyObject *
cmath_polar_impl(PyObject * module,Py_complex z)1028 cmath_polar_impl(PyObject *module, Py_complex z)
1029 /*[clinic end generated code: output=d0a8147c41dbb654 input=26c353574fd1a861]*/
1030 {
1031     double r, phi;
1032 
1033     errno = 0;
1034     phi = c_atan2(z); /* should not cause any exception */
1035     r = _Py_c_abs(z); /* sets errno to ERANGE on overflow */
1036     if (errno != 0)
1037         return math_error();
1038     else
1039         return Py_BuildValue("dd", r, phi);
1040 }
1041 
1042 /*
1043   rect() isn't covered by the C99 standard, but it's not too hard to
1044   figure out 'spirit of C99' rules for special value handing:
1045 
1046     rect(x, t) should behave like exp(log(x) + it) for positive-signed x
1047     rect(x, t) should behave like -exp(log(-x) + it) for negative-signed x
1048     rect(nan, t) should behave like exp(nan + it), except that rect(nan, 0)
1049       gives nan +- i0 with the sign of the imaginary part unspecified.
1050 
1051 */
1052 
1053 static Py_complex rect_special_values[7][7];
1054 
1055 /*[clinic input]
1056 cmath.rect
1057 
1058     r: double
1059     phi: double
1060     /
1061 
1062 Convert from polar coordinates to rectangular coordinates.
1063 [clinic start generated code]*/
1064 
1065 static PyObject *
cmath_rect_impl(PyObject * module,double r,double phi)1066 cmath_rect_impl(PyObject *module, double r, double phi)
1067 /*[clinic end generated code: output=385a0690925df2d5 input=24c5646d147efd69]*/
1068 {
1069     Py_complex z;
1070     errno = 0;
1071 
1072     /* deal with special values */
1073     if (!Py_IS_FINITE(r) || !Py_IS_FINITE(phi)) {
1074         /* if r is +/-infinity and phi is finite but nonzero then
1075            result is (+-INF +-INF i), but we need to compute cos(phi)
1076            and sin(phi) to figure out the signs. */
1077         if (Py_IS_INFINITY(r) && (Py_IS_FINITE(phi)
1078                                   && (phi != 0.))) {
1079             if (r > 0) {
1080                 z.real = copysign(INF, cos(phi));
1081                 z.imag = copysign(INF, sin(phi));
1082             }
1083             else {
1084                 z.real = -copysign(INF, cos(phi));
1085                 z.imag = -copysign(INF, sin(phi));
1086             }
1087         }
1088         else {
1089             z = rect_special_values[special_type(r)]
1090                                    [special_type(phi)];
1091         }
1092         /* need to set errno = EDOM if r is a nonzero number and phi
1093            is infinite */
1094         if (r != 0. && !Py_IS_NAN(r) && Py_IS_INFINITY(phi))
1095             errno = EDOM;
1096         else
1097             errno = 0;
1098     }
1099     else if (phi == 0.0) {
1100         /* Workaround for buggy results with phi=-0.0 on OS X 10.8.  See
1101            bugs.python.org/issue18513. */
1102         z.real = r;
1103         z.imag = r * phi;
1104         errno = 0;
1105     }
1106     else {
1107         z.real = r * cos(phi);
1108         z.imag = r * sin(phi);
1109         errno = 0;
1110     }
1111 
1112     if (errno != 0)
1113         return math_error();
1114     else
1115         return PyComplex_FromCComplex(z);
1116 }
1117 
1118 /*[clinic input]
1119 cmath.isfinite = cmath.polar
1120 
1121 Return True if both the real and imaginary parts of z are finite, else False.
1122 [clinic start generated code]*/
1123 
1124 static PyObject *
cmath_isfinite_impl(PyObject * module,Py_complex z)1125 cmath_isfinite_impl(PyObject *module, Py_complex z)
1126 /*[clinic end generated code: output=ac76611e2c774a36 input=848e7ee701895815]*/
1127 {
1128     return PyBool_FromLong(Py_IS_FINITE(z.real) && Py_IS_FINITE(z.imag));
1129 }
1130 
1131 /*[clinic input]
1132 cmath.isnan = cmath.polar
1133 
1134 Checks if the real or imaginary part of z not a number (NaN).
1135 [clinic start generated code]*/
1136 
1137 static PyObject *
cmath_isnan_impl(PyObject * module,Py_complex z)1138 cmath_isnan_impl(PyObject *module, Py_complex z)
1139 /*[clinic end generated code: output=e7abf6e0b28beab7 input=71799f5d284c9baf]*/
1140 {
1141     return PyBool_FromLong(Py_IS_NAN(z.real) || Py_IS_NAN(z.imag));
1142 }
1143 
1144 /*[clinic input]
1145 cmath.isinf = cmath.polar
1146 
1147 Checks if the real or imaginary part of z is infinite.
1148 [clinic start generated code]*/
1149 
1150 static PyObject *
cmath_isinf_impl(PyObject * module,Py_complex z)1151 cmath_isinf_impl(PyObject *module, Py_complex z)
1152 /*[clinic end generated code: output=502a75a79c773469 input=363df155c7181329]*/
1153 {
1154     return PyBool_FromLong(Py_IS_INFINITY(z.real) ||
1155                            Py_IS_INFINITY(z.imag));
1156 }
1157 
1158 /*[clinic input]
1159 cmath.isclose -> bool
1160 
1161     a: Py_complex
1162     b: Py_complex
1163     *
1164     rel_tol: double = 1e-09
1165         maximum difference for being considered "close", relative to the
1166         magnitude of the input values
1167     abs_tol: double = 0.0
1168         maximum difference for being considered "close", regardless of the
1169         magnitude of the input values
1170 
1171 Determine whether two complex numbers are close in value.
1172 
1173 Return True if a is close in value to b, and False otherwise.
1174 
1175 For the values to be considered close, the difference between them must be
1176 smaller than at least one of the tolerances.
1177 
1178 -inf, inf and NaN behave similarly to the IEEE 754 Standard. That is, NaN is
1179 not close to anything, even itself. inf and -inf are only close to themselves.
1180 [clinic start generated code]*/
1181 
1182 static int
cmath_isclose_impl(PyObject * module,Py_complex a,Py_complex b,double rel_tol,double abs_tol)1183 cmath_isclose_impl(PyObject *module, Py_complex a, Py_complex b,
1184                    double rel_tol, double abs_tol)
1185 /*[clinic end generated code: output=8a2486cc6e0014d1 input=df9636d7de1d4ac3]*/
1186 {
1187     double diff;
1188 
1189     /* sanity check on the inputs */
1190     if (rel_tol < 0.0 || abs_tol < 0.0 ) {
1191         PyErr_SetString(PyExc_ValueError,
1192                         "tolerances must be non-negative");
1193         return -1;
1194     }
1195 
1196     if ( (a.real == b.real) && (a.imag == b.imag) ) {
1197         /* short circuit exact equality -- needed to catch two infinities of
1198            the same sign. And perhaps speeds things up a bit sometimes.
1199         */
1200         return 1;
1201     }
1202 
1203     /* This catches the case of two infinities of opposite sign, or
1204        one infinity and one finite number. Two infinities of opposite
1205        sign would otherwise have an infinite relative tolerance.
1206        Two infinities of the same sign are caught by the equality check
1207        above.
1208     */
1209 
1210     if (Py_IS_INFINITY(a.real) || Py_IS_INFINITY(a.imag) ||
1211         Py_IS_INFINITY(b.real) || Py_IS_INFINITY(b.imag)) {
1212         return 0;
1213     }
1214 
1215     /* now do the regular computation
1216        this is essentially the "weak" test from the Boost library
1217     */
1218 
1219     diff = _Py_c_abs(_Py_c_diff(a, b));
1220 
1221     return (((diff <= rel_tol * _Py_c_abs(b)) ||
1222              (diff <= rel_tol * _Py_c_abs(a))) ||
1223             (diff <= abs_tol));
1224 }
1225 
1226 PyDoc_STRVAR(module_doc,
1227 "This module provides access to mathematical functions for complex\n"
1228 "numbers.");
1229 
1230 static PyMethodDef cmath_methods[] = {
1231     CMATH_ACOS_METHODDEF
1232     CMATH_ACOSH_METHODDEF
1233     CMATH_ASIN_METHODDEF
1234     CMATH_ASINH_METHODDEF
1235     CMATH_ATAN_METHODDEF
1236     CMATH_ATANH_METHODDEF
1237     CMATH_COS_METHODDEF
1238     CMATH_COSH_METHODDEF
1239     CMATH_EXP_METHODDEF
1240     CMATH_ISCLOSE_METHODDEF
1241     CMATH_ISFINITE_METHODDEF
1242     CMATH_ISINF_METHODDEF
1243     CMATH_ISNAN_METHODDEF
1244     CMATH_LOG_METHODDEF
1245     CMATH_LOG10_METHODDEF
1246     CMATH_PHASE_METHODDEF
1247     CMATH_POLAR_METHODDEF
1248     CMATH_RECT_METHODDEF
1249     CMATH_SIN_METHODDEF
1250     CMATH_SINH_METHODDEF
1251     CMATH_SQRT_METHODDEF
1252     CMATH_TAN_METHODDEF
1253     CMATH_TANH_METHODDEF
1254     {NULL, NULL}  /* sentinel */
1255 };
1256 
1257 static int
cmath_exec(PyObject * mod)1258 cmath_exec(PyObject *mod)
1259 {
1260     if (PyModule_AddObject(mod, "pi", PyFloat_FromDouble(Py_MATH_PI)) < 0) {
1261         return -1;
1262     }
1263     if (PyModule_AddObject(mod, "e", PyFloat_FromDouble(Py_MATH_E)) < 0) {
1264         return -1;
1265     }
1266     // 2pi
1267     if (PyModule_AddObject(mod, "tau", PyFloat_FromDouble(Py_MATH_TAU)) < 0) {
1268         return -1;
1269     }
1270     if (PyModule_AddObject(mod, "inf", PyFloat_FromDouble(m_inf())) < 0) {
1271         return -1;
1272     }
1273 
1274     if (PyModule_AddObject(mod, "infj",
1275                            PyComplex_FromCComplex(c_infj())) < 0) {
1276         return -1;
1277     }
1278 #if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN)
1279     if (PyModule_AddObject(mod, "nan", PyFloat_FromDouble(m_nan())) < 0) {
1280         return -1;
1281     }
1282     if (PyModule_AddObject(mod, "nanj",
1283                            PyComplex_FromCComplex(c_nanj())) < 0) {
1284         return -1;
1285     }
1286 #endif
1287 
1288     /* initialize special value tables */
1289 
1290 #define INIT_SPECIAL_VALUES(NAME, BODY) { Py_complex* p = (Py_complex*)NAME; BODY }
1291 #define C(REAL, IMAG) p->real = REAL; p->imag = IMAG; ++p;
1292 
1293     INIT_SPECIAL_VALUES(acos_special_values, {
1294       C(P34,INF) C(P,INF)  C(P,INF)  C(P,-INF)  C(P,-INF)  C(P34,-INF) C(N,INF)
1295       C(P12,INF) C(U,U)    C(U,U)    C(U,U)     C(U,U)     C(P12,-INF) C(N,N)
1296       C(P12,INF) C(U,U)    C(P12,0.) C(P12,-0.) C(U,U)     C(P12,-INF) C(P12,N)
1297       C(P12,INF) C(U,U)    C(P12,0.) C(P12,-0.) C(U,U)     C(P12,-INF) C(P12,N)
1298       C(P12,INF) C(U,U)    C(U,U)    C(U,U)     C(U,U)     C(P12,-INF) C(N,N)
1299       C(P14,INF) C(0.,INF) C(0.,INF) C(0.,-INF) C(0.,-INF) C(P14,-INF) C(N,INF)
1300       C(N,INF)   C(N,N)    C(N,N)    C(N,N)     C(N,N)     C(N,-INF)   C(N,N)
1301     })
1302 
1303     INIT_SPECIAL_VALUES(acosh_special_values, {
1304       C(INF,-P34) C(INF,-P)  C(INF,-P)  C(INF,P)  C(INF,P)  C(INF,P34) C(INF,N)
1305       C(INF,-P12) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,P12) C(N,N)
1306       C(INF,-P12) C(U,U)     C(0.,-P12) C(0.,P12) C(U,U)    C(INF,P12) C(N,N)
1307       C(INF,-P12) C(U,U)     C(0.,-P12) C(0.,P12) C(U,U)    C(INF,P12) C(N,N)
1308       C(INF,-P12) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,P12) C(N,N)
1309       C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
1310       C(INF,N)    C(N,N)     C(N,N)     C(N,N)    C(N,N)    C(INF,N)   C(N,N)
1311     })
1312 
1313     INIT_SPECIAL_VALUES(asinh_special_values, {
1314       C(-INF,-P14) C(-INF,-0.) C(-INF,-0.) C(-INF,0.) C(-INF,0.) C(-INF,P14) 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.,-0.)  C(-0.,0.)  C(U,U)     C(-INF,P12) C(N,N)
1317       C(INF,-P12)  C(U,U)      C(0.,-0.)   C(0.,0.)   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,-0.)    C(N,0.)    C(N,N)     C(INF,N)    C(N,N)
1321     })
1322 
1323     INIT_SPECIAL_VALUES(atanh_special_values, {
1324       C(-0.,-P12) C(-0.,-P12) C(-0.,-P12) C(-0.,P12) C(-0.,P12) C(-0.,P12) C(-0.,N)
1325       C(-0.,-P12) C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(-0.,P12) C(N,N)
1326       C(-0.,-P12) C(U,U)      C(-0.,-0.)  C(-0.,0.)  C(U,U)     C(-0.,P12) C(-0.,N)
1327       C(0.,-P12)  C(U,U)      C(0.,-0.)   C(0.,0.)   C(U,U)     C(0.,P12)  C(0.,N)
1328       C(0.,-P12)  C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(0.,P12)  C(N,N)
1329       C(0.,-P12)  C(0.,-P12)  C(0.,-P12)  C(0.,P12)  C(0.,P12)  C(0.,P12)  C(0.,N)
1330       C(0.,-P12)  C(N,N)      C(N,N)      C(N,N)     C(N,N)     C(0.,P12)  C(N,N)
1331     })
1332 
1333     INIT_SPECIAL_VALUES(cosh_special_values, {
1334       C(INF,N) C(U,U) C(INF,0.)  C(INF,-0.) C(U,U) C(INF,N) C(INF,N)
1335       C(N,N)   C(U,U) C(U,U)     C(U,U)     C(U,U) C(N,N)   C(N,N)
1336       C(N,0.)  C(U,U) C(1.,0.)   C(1.,-0.)  C(U,U) C(N,0.)  C(N,0.)
1337       C(N,0.)  C(U,U) C(1.,-0.)  C(1.,0.)   C(U,U) C(N,0.)  C(N,0.)
1338       C(N,N)   C(U,U) C(U,U)     C(U,U)     C(U,U) C(N,N)   C(N,N)
1339       C(INF,N) C(U,U) C(INF,-0.) C(INF,0.)  C(U,U) C(INF,N) C(INF,N)
1340       C(N,N)   C(N,N) C(N,0.)    C(N,0.)    C(N,N) C(N,N)   C(N,N)
1341     })
1342 
1343     INIT_SPECIAL_VALUES(exp_special_values, {
1344       C(0.,0.) C(U,U) C(0.,-0.)  C(0.,0.)  C(U,U) C(0.,0.) C(0.,0.)
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,N)   C(U,U) C(1.,-0.)  C(1.,0.)  C(U,U) C(N,N)   C(N,N)
1347       C(N,N)   C(U,U) C(1.,-0.)  C(1.,0.)  C(U,U) C(N,N)   C(N,N)
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(log_special_values, {
1354       C(INF,-P34) C(INF,-P)  C(INF,-P)   C(INF,P)   C(INF,P)  C(INF,P34)  C(INF,N)
1355       C(INF,-P12) C(U,U)     C(U,U)      C(U,U)     C(U,U)    C(INF,P12)  C(N,N)
1356       C(INF,-P12) C(U,U)     C(-INF,-P)  C(-INF,P)  C(U,U)    C(INF,P12)  C(N,N)
1357       C(INF,-P12) C(U,U)     C(-INF,-0.) C(-INF,0.) C(U,U)    C(INF,P12)  C(N,N)
1358       C(INF,-P12) C(U,U)     C(U,U)      C(U,U)     C(U,U)    C(INF,P12)  C(N,N)
1359       C(INF,-P14) C(INF,-0.) C(INF,-0.)  C(INF,0.)  C(INF,0.) C(INF,P14)  C(INF,N)
1360       C(INF,N)    C(N,N)     C(N,N)      C(N,N)     C(N,N)    C(INF,N)    C(N,N)
1361     })
1362 
1363     INIT_SPECIAL_VALUES(sinh_special_values, {
1364       C(INF,N) C(U,U) C(-INF,-0.) C(-INF,0.) C(U,U) C(INF,N) C(INF,N)
1365       C(N,N)   C(U,U) C(U,U)      C(U,U)     C(U,U) C(N,N)   C(N,N)
1366       C(0.,N)  C(U,U) C(-0.,-0.)  C(-0.,0.)  C(U,U) C(0.,N)  C(0.,N)
1367       C(0.,N)  C(U,U) C(0.,-0.)   C(0.,0.)   C(U,U) C(0.,N)  C(0.,N)
1368       C(N,N)   C(U,U) C(U,U)      C(U,U)     C(U,U) C(N,N)   C(N,N)
1369       C(INF,N) C(U,U) C(INF,-0.)  C(INF,0.)  C(U,U) C(INF,N) C(INF,N)
1370       C(N,N)   C(N,N) C(N,-0.)    C(N,0.)    C(N,N) C(N,N)   C(N,N)
1371     })
1372 
1373     INIT_SPECIAL_VALUES(sqrt_special_values, {
1374       C(INF,-INF) C(0.,-INF) C(0.,-INF) C(0.,INF) C(0.,INF) C(INF,INF) C(N,INF)
1375       C(INF,-INF) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,INF) C(N,N)
1376       C(INF,-INF) C(U,U)     C(0.,-0.)  C(0.,0.)  C(U,U)    C(INF,INF) C(N,N)
1377       C(INF,-INF) C(U,U)     C(0.,-0.)  C(0.,0.)  C(U,U)    C(INF,INF) C(N,N)
1378       C(INF,-INF) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,INF) C(N,N)
1379       C(INF,-INF) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,INF) C(INF,N)
1380       C(INF,-INF) C(N,N)     C(N,N)     C(N,N)    C(N,N)    C(INF,INF) C(N,N)
1381     })
1382 
1383     INIT_SPECIAL_VALUES(tanh_special_values, {
1384       C(-1.,0.) C(U,U) C(-1.,-0.) C(-1.,0.) C(U,U) C(-1.,0.) C(-1.,0.)
1385       C(N,N)    C(U,U) C(U,U)     C(U,U)    C(U,U) C(N,N)    C(N,N)
1386       C(N,N)    C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(N,N)    C(N,N)
1387       C(N,N)    C(U,U) C(0.,-0.)  C(0.,0.)  C(U,U) C(N,N)    C(N,N)
1388       C(N,N)    C(U,U) C(U,U)     C(U,U)    C(U,U) C(N,N)    C(N,N)
1389       C(1.,0.)  C(U,U) C(1.,-0.)  C(1.,0.)  C(U,U) C(1.,0.)  C(1.,0.)
1390       C(N,N)    C(N,N) C(N,-0.)   C(N,0.)   C(N,N) C(N,N)    C(N,N)
1391     })
1392 
1393     INIT_SPECIAL_VALUES(rect_special_values, {
1394       C(INF,N) C(U,U) C(-INF,0.) C(-INF,-0.) C(U,U) C(INF,N) C(INF,N)
1395       C(N,N)   C(U,U) C(U,U)     C(U,U)      C(U,U) C(N,N)   C(N,N)
1396       C(0.,0.) C(U,U) C(-0.,0.)  C(-0.,-0.)  C(U,U) C(0.,0.) C(0.,0.)
1397       C(0.,0.) C(U,U) C(0.,-0.)  C(0.,0.)    C(U,U) C(0.,0.) C(0.,0.)
1398       C(N,N)   C(U,U) C(U,U)     C(U,U)      C(U,U) C(N,N)   C(N,N)
1399       C(INF,N) C(U,U) C(INF,-0.) C(INF,0.)   C(U,U) C(INF,N) C(INF,N)
1400       C(N,N)   C(N,N) C(N,0.)    C(N,0.)     C(N,N) C(N,N)   C(N,N)
1401     })
1402     return 0;
1403 }
1404 
1405 static PyModuleDef_Slot cmath_slots[] = {
1406     {Py_mod_exec, cmath_exec},
1407     {0, NULL}
1408 };
1409 
1410 static struct PyModuleDef cmathmodule = {
1411     PyModuleDef_HEAD_INIT,
1412     .m_name = "cmath",
1413     .m_doc = module_doc,
1414     .m_size = 0,
1415     .m_methods = cmath_methods,
1416     .m_slots = cmath_slots
1417 };
1418 
1419 PyMODINIT_FUNC
PyInit_cmath(void)1420 PyInit_cmath(void)
1421 {
1422     return PyModuleDef_Init(&cmathmodule);
1423 }