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