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