Add better error reporting for MemoryErrors caused by str->float conversions.
[python.git] / Modules / cmathmodule.c
blob788a19dad5c4182abaa092fd8b58552a1781e529
1 /* Complex math module */
3 /* much code borrowed from mathmodule.c */
5 #include "Python.h"
6 #include "_math.h"
7 /* we need DBL_MAX, DBL_MIN, DBL_EPSILON, DBL_MANT_DIG and FLT_RADIX from
8 float.h. We assume that FLT_RADIX is either 2 or 16. */
9 #include <float.h>
11 #if (FLT_RADIX != 2 && FLT_RADIX != 16)
12 #error "Modules/cmathmodule.c expects FLT_RADIX to be 2 or 16"
13 #endif
15 #ifndef M_LN2
16 #define M_LN2 (0.6931471805599453094) /* natural log of 2 */
17 #endif
19 #ifndef M_LN10
20 #define M_LN10 (2.302585092994045684) /* natural log of 10 */
21 #endif
24 CM_LARGE_DOUBLE is used to avoid spurious overflow in the sqrt, log,
25 inverse trig and inverse hyperbolic trig functions. Its log is used in the
26 evaluation of exp, cos, cosh, sin, sinh, tan, and tanh to avoid unecessary
27 overflow.
30 #define CM_LARGE_DOUBLE (DBL_MAX/4.)
31 #define CM_SQRT_LARGE_DOUBLE (sqrt(CM_LARGE_DOUBLE))
32 #define CM_LOG_LARGE_DOUBLE (log(CM_LARGE_DOUBLE))
33 #define CM_SQRT_DBL_MIN (sqrt(DBL_MIN))
35 /*
36 CM_SCALE_UP is an odd integer chosen such that multiplication by
37 2**CM_SCALE_UP is sufficient to turn a subnormal into a normal.
38 CM_SCALE_DOWN is (-(CM_SCALE_UP+1)/2). These scalings are used to compute
39 square roots accurately when the real and imaginary parts of the argument
40 are subnormal.
43 #if FLT_RADIX==2
44 #define CM_SCALE_UP (2*(DBL_MANT_DIG/2) + 1)
45 #elif FLT_RADIX==16
46 #define CM_SCALE_UP (4*DBL_MANT_DIG+1)
47 #endif
48 #define CM_SCALE_DOWN (-(CM_SCALE_UP+1)/2)
50 /* forward declarations */
51 static Py_complex c_asinh(Py_complex);
52 static Py_complex c_atanh(Py_complex);
53 static Py_complex c_cosh(Py_complex);
54 static Py_complex c_sinh(Py_complex);
55 static Py_complex c_sqrt(Py_complex);
56 static Py_complex c_tanh(Py_complex);
57 static PyObject * math_error(void);
59 /* Code to deal with special values (infinities, NaNs, etc.). */
61 /* special_type takes a double and returns an integer code indicating
62 the type of the double as follows:
65 enum special_types {
66 ST_NINF, /* 0, negative infinity */
67 ST_NEG, /* 1, negative finite number (nonzero) */
68 ST_NZERO, /* 2, -0. */
69 ST_PZERO, /* 3, +0. */
70 ST_POS, /* 4, positive finite number (nonzero) */
71 ST_PINF, /* 5, positive infinity */
72 ST_NAN /* 6, Not a Number */
75 static enum special_types
76 special_type(double d)
78 if (Py_IS_FINITE(d)) {
79 if (d != 0) {
80 if (copysign(1., d) == 1.)
81 return ST_POS;
82 else
83 return ST_NEG;
85 else {
86 if (copysign(1., d) == 1.)
87 return ST_PZERO;
88 else
89 return ST_NZERO;
92 if (Py_IS_NAN(d))
93 return ST_NAN;
94 if (copysign(1., d) == 1.)
95 return ST_PINF;
96 else
97 return ST_NINF;
100 #define SPECIAL_VALUE(z, table) \
101 if (!Py_IS_FINITE((z).real) || !Py_IS_FINITE((z).imag)) { \
102 errno = 0; \
103 return table[special_type((z).real)] \
104 [special_type((z).imag)]; \
107 #define P Py_MATH_PI
108 #define P14 0.25*Py_MATH_PI
109 #define P12 0.5*Py_MATH_PI
110 #define P34 0.75*Py_MATH_PI
111 #define INF Py_HUGE_VAL
112 #define N Py_NAN
113 #define U -9.5426319407711027e33 /* unlikely value, used as placeholder */
115 /* First, the C functions that do the real work. Each of the c_*
116 functions computes and returns the C99 Annex G recommended result
117 and also sets errno as follows: errno = 0 if no floating-point
118 exception is associated with the result; errno = EDOM if C99 Annex
119 G recommends raising divide-by-zero or invalid for this result; and
120 errno = ERANGE where the overflow floating-point signal should be
121 raised.
124 static Py_complex acos_special_values[7][7];
126 static Py_complex
127 c_acos(Py_complex z)
129 Py_complex s1, s2, r;
131 SPECIAL_VALUE(z, acos_special_values);
133 if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
134 /* avoid unnecessary overflow for large arguments */
135 r.real = atan2(fabs(z.imag), z.real);
136 /* split into cases to make sure that the branch cut has the
137 correct continuity on systems with unsigned zeros */
138 if (z.real < 0.) {
139 r.imag = -copysign(log(hypot(z.real/2., z.imag/2.)) +
140 M_LN2*2., z.imag);
141 } else {
142 r.imag = copysign(log(hypot(z.real/2., z.imag/2.)) +
143 M_LN2*2., -z.imag);
145 } else {
146 s1.real = 1.-z.real;
147 s1.imag = -z.imag;
148 s1 = c_sqrt(s1);
149 s2.real = 1.+z.real;
150 s2.imag = z.imag;
151 s2 = c_sqrt(s2);
152 r.real = 2.*atan2(s1.real, s2.real);
153 r.imag = m_asinh(s2.real*s1.imag - s2.imag*s1.real);
155 errno = 0;
156 return r;
159 PyDoc_STRVAR(c_acos_doc,
160 "acos(x)\n"
161 "\n"
162 "Return the arc cosine of x.");
165 static Py_complex acosh_special_values[7][7];
167 static Py_complex
168 c_acosh(Py_complex z)
170 Py_complex s1, s2, r;
172 SPECIAL_VALUE(z, acosh_special_values);
174 if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
175 /* avoid unnecessary overflow for large arguments */
176 r.real = log(hypot(z.real/2., z.imag/2.)) + M_LN2*2.;
177 r.imag = atan2(z.imag, z.real);
178 } else {
179 s1.real = z.real - 1.;
180 s1.imag = z.imag;
181 s1 = c_sqrt(s1);
182 s2.real = z.real + 1.;
183 s2.imag = z.imag;
184 s2 = c_sqrt(s2);
185 r.real = m_asinh(s1.real*s2.real + s1.imag*s2.imag);
186 r.imag = 2.*atan2(s1.imag, s2.real);
188 errno = 0;
189 return r;
192 PyDoc_STRVAR(c_acosh_doc,
193 "acosh(x)\n"
194 "\n"
195 "Return the hyperbolic arccosine of x.");
198 static Py_complex
199 c_asin(Py_complex z)
201 /* asin(z) = -i asinh(iz) */
202 Py_complex s, r;
203 s.real = -z.imag;
204 s.imag = z.real;
205 s = c_asinh(s);
206 r.real = s.imag;
207 r.imag = -s.real;
208 return r;
211 PyDoc_STRVAR(c_asin_doc,
212 "asin(x)\n"
213 "\n"
214 "Return the arc sine of x.");
217 static Py_complex asinh_special_values[7][7];
219 static Py_complex
220 c_asinh(Py_complex z)
222 Py_complex s1, s2, r;
224 SPECIAL_VALUE(z, asinh_special_values);
226 if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
227 if (z.imag >= 0.) {
228 r.real = copysign(log(hypot(z.real/2., z.imag/2.)) +
229 M_LN2*2., z.real);
230 } else {
231 r.real = -copysign(log(hypot(z.real/2., z.imag/2.)) +
232 M_LN2*2., -z.real);
234 r.imag = atan2(z.imag, fabs(z.real));
235 } else {
236 s1.real = 1.+z.imag;
237 s1.imag = -z.real;
238 s1 = c_sqrt(s1);
239 s2.real = 1.-z.imag;
240 s2.imag = z.real;
241 s2 = c_sqrt(s2);
242 r.real = m_asinh(s1.real*s2.imag-s2.real*s1.imag);
243 r.imag = atan2(z.imag, s1.real*s2.real-s1.imag*s2.imag);
245 errno = 0;
246 return r;
249 PyDoc_STRVAR(c_asinh_doc,
250 "asinh(x)\n"
251 "\n"
252 "Return the hyperbolic arc sine of x.");
255 static Py_complex
256 c_atan(Py_complex z)
258 /* atan(z) = -i atanh(iz) */
259 Py_complex s, r;
260 s.real = -z.imag;
261 s.imag = z.real;
262 s = c_atanh(s);
263 r.real = s.imag;
264 r.imag = -s.real;
265 return r;
268 /* Windows screws up atan2 for inf and nan, and alpha Tru64 5.1 doesn't follow
269 C99 for atan2(0., 0.). */
270 static double
271 c_atan2(Py_complex z)
273 if (Py_IS_NAN(z.real) || Py_IS_NAN(z.imag))
274 return Py_NAN;
275 if (Py_IS_INFINITY(z.imag)) {
276 if (Py_IS_INFINITY(z.real)) {
277 if (copysign(1., z.real) == 1.)
278 /* atan2(+-inf, +inf) == +-pi/4 */
279 return copysign(0.25*Py_MATH_PI, z.imag);
280 else
281 /* atan2(+-inf, -inf) == +-pi*3/4 */
282 return copysign(0.75*Py_MATH_PI, z.imag);
284 /* atan2(+-inf, x) == +-pi/2 for finite x */
285 return copysign(0.5*Py_MATH_PI, z.imag);
287 if (Py_IS_INFINITY(z.real) || z.imag == 0.) {
288 if (copysign(1., z.real) == 1.)
289 /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
290 return copysign(0., z.imag);
291 else
292 /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
293 return copysign(Py_MATH_PI, z.imag);
295 return atan2(z.imag, z.real);
298 PyDoc_STRVAR(c_atan_doc,
299 "atan(x)\n"
300 "\n"
301 "Return the arc tangent of x.");
304 static Py_complex atanh_special_values[7][7];
306 static Py_complex
307 c_atanh(Py_complex z)
309 Py_complex r;
310 double ay, h;
312 SPECIAL_VALUE(z, atanh_special_values);
314 /* Reduce to case where z.real >= 0., using atanh(z) = -atanh(-z). */
315 if (z.real < 0.) {
316 return c_neg(c_atanh(c_neg(z)));
319 ay = fabs(z.imag);
320 if (z.real > CM_SQRT_LARGE_DOUBLE || ay > CM_SQRT_LARGE_DOUBLE) {
322 if abs(z) is large then we use the approximation
323 atanh(z) ~ 1/z +/- i*pi/2 (+/- depending on the sign
324 of z.imag)
326 h = hypot(z.real/2., z.imag/2.); /* safe from overflow */
327 r.real = z.real/4./h/h;
328 /* the two negations in the next line cancel each other out
329 except when working with unsigned zeros: they're there to
330 ensure that the branch cut has the correct continuity on
331 systems that don't support signed zeros */
332 r.imag = -copysign(Py_MATH_PI/2., -z.imag);
333 errno = 0;
334 } else if (z.real == 1. && ay < CM_SQRT_DBL_MIN) {
335 /* C99 standard says: atanh(1+/-0.) should be inf +/- 0i */
336 if (ay == 0.) {
337 r.real = INF;
338 r.imag = z.imag;
339 errno = EDOM;
340 } else {
341 r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.)));
342 r.imag = copysign(atan2(2., -ay)/2, z.imag);
343 errno = 0;
345 } else {
346 r.real = m_log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.;
347 r.imag = -atan2(-2.*z.imag, (1-z.real)*(1+z.real) - ay*ay)/2.;
348 errno = 0;
350 return r;
353 PyDoc_STRVAR(c_atanh_doc,
354 "atanh(x)\n"
355 "\n"
356 "Return the hyperbolic arc tangent of x.");
359 static Py_complex
360 c_cos(Py_complex z)
362 /* cos(z) = cosh(iz) */
363 Py_complex r;
364 r.real = -z.imag;
365 r.imag = z.real;
366 r = c_cosh(r);
367 return r;
370 PyDoc_STRVAR(c_cos_doc,
371 "cos(x)\n"
372 "\n"
373 "Return the cosine of x.");
376 /* cosh(infinity + i*y) needs to be dealt with specially */
377 static Py_complex cosh_special_values[7][7];
379 static Py_complex
380 c_cosh(Py_complex z)
382 Py_complex r;
383 double x_minus_one;
385 /* special treatment for cosh(+/-inf + iy) if y is not a NaN */
386 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
387 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) &&
388 (z.imag != 0.)) {
389 if (z.real > 0) {
390 r.real = copysign(INF, cos(z.imag));
391 r.imag = copysign(INF, sin(z.imag));
393 else {
394 r.real = copysign(INF, cos(z.imag));
395 r.imag = -copysign(INF, sin(z.imag));
398 else {
399 r = cosh_special_values[special_type(z.real)]
400 [special_type(z.imag)];
402 /* need to set errno = EDOM if y is +/- infinity and x is not
403 a NaN */
404 if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
405 errno = EDOM;
406 else
407 errno = 0;
408 return r;
411 if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
412 /* deal correctly with cases where cosh(z.real) overflows but
413 cosh(z) does not. */
414 x_minus_one = z.real - copysign(1., z.real);
415 r.real = cos(z.imag) * cosh(x_minus_one) * Py_MATH_E;
416 r.imag = sin(z.imag) * sinh(x_minus_one) * Py_MATH_E;
417 } else {
418 r.real = cos(z.imag) * cosh(z.real);
419 r.imag = sin(z.imag) * sinh(z.real);
421 /* detect overflow, and set errno accordingly */
422 if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
423 errno = ERANGE;
424 else
425 errno = 0;
426 return r;
429 PyDoc_STRVAR(c_cosh_doc,
430 "cosh(x)\n"
431 "\n"
432 "Return the hyperbolic cosine of x.");
435 /* exp(infinity + i*y) and exp(-infinity + i*y) need special treatment for
436 finite y */
437 static Py_complex exp_special_values[7][7];
439 static Py_complex
440 c_exp(Py_complex z)
442 Py_complex r;
443 double l;
445 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
446 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
447 && (z.imag != 0.)) {
448 if (z.real > 0) {
449 r.real = copysign(INF, cos(z.imag));
450 r.imag = copysign(INF, sin(z.imag));
452 else {
453 r.real = copysign(0., cos(z.imag));
454 r.imag = copysign(0., sin(z.imag));
457 else {
458 r = exp_special_values[special_type(z.real)]
459 [special_type(z.imag)];
461 /* need to set errno = EDOM if y is +/- infinity and x is not
462 a NaN and not -infinity */
463 if (Py_IS_INFINITY(z.imag) &&
464 (Py_IS_FINITE(z.real) ||
465 (Py_IS_INFINITY(z.real) && z.real > 0)))
466 errno = EDOM;
467 else
468 errno = 0;
469 return r;
472 if (z.real > CM_LOG_LARGE_DOUBLE) {
473 l = exp(z.real-1.);
474 r.real = l*cos(z.imag)*Py_MATH_E;
475 r.imag = l*sin(z.imag)*Py_MATH_E;
476 } else {
477 l = exp(z.real);
478 r.real = l*cos(z.imag);
479 r.imag = l*sin(z.imag);
481 /* detect overflow, and set errno accordingly */
482 if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
483 errno = ERANGE;
484 else
485 errno = 0;
486 return r;
489 PyDoc_STRVAR(c_exp_doc,
490 "exp(x)\n"
491 "\n"
492 "Return the exponential value e**x.");
495 static Py_complex log_special_values[7][7];
497 static Py_complex
498 c_log(Py_complex z)
501 The usual formula for the real part is log(hypot(z.real, z.imag)).
502 There are four situations where this formula is potentially
503 problematic:
505 (1) the absolute value of z is subnormal. Then hypot is subnormal,
506 so has fewer than the usual number of bits of accuracy, hence may
507 have large relative error. This then gives a large absolute error
508 in the log. This can be solved by rescaling z by a suitable power
509 of 2.
511 (2) the absolute value of z is greater than DBL_MAX (e.g. when both
512 z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX)
513 Again, rescaling solves this.
515 (3) the absolute value of z is close to 1. In this case it's
516 difficult to achieve good accuracy, at least in part because a
517 change of 1ulp in the real or imaginary part of z can result in a
518 change of billions of ulps in the correctly rounded answer.
520 (4) z = 0. The simplest thing to do here is to call the
521 floating-point log with an argument of 0, and let its behaviour
522 (returning -infinity, signaling a floating-point exception, setting
523 errno, or whatever) determine that of c_log. So the usual formula
524 is fine here.
528 Py_complex r;
529 double ax, ay, am, an, h;
531 SPECIAL_VALUE(z, log_special_values);
533 ax = fabs(z.real);
534 ay = fabs(z.imag);
536 if (ax > CM_LARGE_DOUBLE || ay > CM_LARGE_DOUBLE) {
537 r.real = log(hypot(ax/2., ay/2.)) + M_LN2;
538 } else if (ax < DBL_MIN && ay < DBL_MIN) {
539 if (ax > 0. || ay > 0.) {
540 /* catch cases where hypot(ax, ay) is subnormal */
541 r.real = log(hypot(ldexp(ax, DBL_MANT_DIG),
542 ldexp(ay, DBL_MANT_DIG))) - DBL_MANT_DIG*M_LN2;
544 else {
545 /* log(+/-0. +/- 0i) */
546 r.real = -INF;
547 r.imag = atan2(z.imag, z.real);
548 errno = EDOM;
549 return r;
551 } else {
552 h = hypot(ax, ay);
553 if (0.71 <= h && h <= 1.73) {
554 am = ax > ay ? ax : ay; /* max(ax, ay) */
555 an = ax > ay ? ay : ax; /* min(ax, ay) */
556 r.real = m_log1p((am-1)*(am+1)+an*an)/2.;
557 } else {
558 r.real = log(h);
561 r.imag = atan2(z.imag, z.real);
562 errno = 0;
563 return r;
567 static Py_complex
568 c_log10(Py_complex z)
570 Py_complex r;
571 int errno_save;
573 r = c_log(z);
574 errno_save = errno; /* just in case the divisions affect errno */
575 r.real = r.real / M_LN10;
576 r.imag = r.imag / M_LN10;
577 errno = errno_save;
578 return r;
581 PyDoc_STRVAR(c_log10_doc,
582 "log10(x)\n"
583 "\n"
584 "Return the base-10 logarithm of x.");
587 static Py_complex
588 c_sin(Py_complex z)
590 /* sin(z) = -i sin(iz) */
591 Py_complex s, r;
592 s.real = -z.imag;
593 s.imag = z.real;
594 s = c_sinh(s);
595 r.real = s.imag;
596 r.imag = -s.real;
597 return r;
600 PyDoc_STRVAR(c_sin_doc,
601 "sin(x)\n"
602 "\n"
603 "Return the sine of x.");
606 /* sinh(infinity + i*y) needs to be dealt with specially */
607 static Py_complex sinh_special_values[7][7];
609 static Py_complex
610 c_sinh(Py_complex z)
612 Py_complex r;
613 double x_minus_one;
615 /* special treatment for sinh(+/-inf + iy) if y is finite and
616 nonzero */
617 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
618 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
619 && (z.imag != 0.)) {
620 if (z.real > 0) {
621 r.real = copysign(INF, cos(z.imag));
622 r.imag = copysign(INF, sin(z.imag));
624 else {
625 r.real = -copysign(INF, cos(z.imag));
626 r.imag = copysign(INF, sin(z.imag));
629 else {
630 r = sinh_special_values[special_type(z.real)]
631 [special_type(z.imag)];
633 /* need to set errno = EDOM if y is +/- infinity and x is not
634 a NaN */
635 if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
636 errno = EDOM;
637 else
638 errno = 0;
639 return r;
642 if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
643 x_minus_one = z.real - copysign(1., z.real);
644 r.real = cos(z.imag) * sinh(x_minus_one) * Py_MATH_E;
645 r.imag = sin(z.imag) * cosh(x_minus_one) * Py_MATH_E;
646 } else {
647 r.real = cos(z.imag) * sinh(z.real);
648 r.imag = sin(z.imag) * cosh(z.real);
650 /* detect overflow, and set errno accordingly */
651 if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
652 errno = ERANGE;
653 else
654 errno = 0;
655 return r;
658 PyDoc_STRVAR(c_sinh_doc,
659 "sinh(x)\n"
660 "\n"
661 "Return the hyperbolic sine of x.");
664 static Py_complex sqrt_special_values[7][7];
666 static Py_complex
667 c_sqrt(Py_complex z)
670 Method: use symmetries to reduce to the case when x = z.real and y
671 = z.imag are nonnegative. Then the real part of the result is
672 given by
674 s = sqrt((x + hypot(x, y))/2)
676 and the imaginary part is
678 d = (y/2)/s
680 If either x or y is very large then there's a risk of overflow in
681 computation of the expression x + hypot(x, y). We can avoid this
682 by rewriting the formula for s as:
684 s = 2*sqrt(x/8 + hypot(x/8, y/8))
686 This costs us two extra multiplications/divisions, but avoids the
687 overhead of checking for x and y large.
689 If both x and y are subnormal then hypot(x, y) may also be
690 subnormal, so will lack full precision. We solve this by rescaling
691 x and y by a sufficiently large power of 2 to ensure that x and y
692 are normal.
696 Py_complex r;
697 double s,d;
698 double ax, ay;
700 SPECIAL_VALUE(z, sqrt_special_values);
702 if (z.real == 0. && z.imag == 0.) {
703 r.real = 0.;
704 r.imag = z.imag;
705 return r;
708 ax = fabs(z.real);
709 ay = fabs(z.imag);
711 if (ax < DBL_MIN && ay < DBL_MIN && (ax > 0. || ay > 0.)) {
712 /* here we catch cases where hypot(ax, ay) is subnormal */
713 ax = ldexp(ax, CM_SCALE_UP);
714 s = ldexp(sqrt(ax + hypot(ax, ldexp(ay, CM_SCALE_UP))),
715 CM_SCALE_DOWN);
716 } else {
717 ax /= 8.;
718 s = 2.*sqrt(ax + hypot(ax, ay/8.));
720 d = ay/(2.*s);
722 if (z.real >= 0.) {
723 r.real = s;
724 r.imag = copysign(d, z.imag);
725 } else {
726 r.real = d;
727 r.imag = copysign(s, z.imag);
729 errno = 0;
730 return r;
733 PyDoc_STRVAR(c_sqrt_doc,
734 "sqrt(x)\n"
735 "\n"
736 "Return the square root of x.");
739 static Py_complex
740 c_tan(Py_complex z)
742 /* tan(z) = -i tanh(iz) */
743 Py_complex s, r;
744 s.real = -z.imag;
745 s.imag = z.real;
746 s = c_tanh(s);
747 r.real = s.imag;
748 r.imag = -s.real;
749 return r;
752 PyDoc_STRVAR(c_tan_doc,
753 "tan(x)\n"
754 "\n"
755 "Return the tangent of x.");
758 /* tanh(infinity + i*y) needs to be dealt with specially */
759 static Py_complex tanh_special_values[7][7];
761 static Py_complex
762 c_tanh(Py_complex z)
764 /* Formula:
766 tanh(x+iy) = (tanh(x)(1+tan(y)^2) + i tan(y)(1-tanh(x))^2) /
767 (1+tan(y)^2 tanh(x)^2)
769 To avoid excessive roundoff error, 1-tanh(x)^2 is better computed
770 as 1/cosh(x)^2. When abs(x) is large, we approximate 1-tanh(x)^2
771 by 4 exp(-2*x) instead, to avoid possible overflow in the
772 computation of cosh(x).
776 Py_complex r;
777 double tx, ty, cx, txty, denom;
779 /* special treatment for tanh(+/-inf + iy) if y is finite and
780 nonzero */
781 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
782 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
783 && (z.imag != 0.)) {
784 if (z.real > 0) {
785 r.real = 1.0;
786 r.imag = copysign(0.,
787 2.*sin(z.imag)*cos(z.imag));
789 else {
790 r.real = -1.0;
791 r.imag = copysign(0.,
792 2.*sin(z.imag)*cos(z.imag));
795 else {
796 r = tanh_special_values[special_type(z.real)]
797 [special_type(z.imag)];
799 /* need to set errno = EDOM if z.imag is +/-infinity and
800 z.real is finite */
801 if (Py_IS_INFINITY(z.imag) && Py_IS_FINITE(z.real))
802 errno = EDOM;
803 else
804 errno = 0;
805 return r;
808 /* danger of overflow in 2.*z.imag !*/
809 if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
810 r.real = copysign(1., z.real);
811 r.imag = 4.*sin(z.imag)*cos(z.imag)*exp(-2.*fabs(z.real));
812 } else {
813 tx = tanh(z.real);
814 ty = tan(z.imag);
815 cx = 1./cosh(z.real);
816 txty = tx*ty;
817 denom = 1. + txty*txty;
818 r.real = tx*(1.+ty*ty)/denom;
819 r.imag = ((ty/denom)*cx)*cx;
821 errno = 0;
822 return r;
825 PyDoc_STRVAR(c_tanh_doc,
826 "tanh(x)\n"
827 "\n"
828 "Return the hyperbolic tangent of x.");
831 static PyObject *
832 cmath_log(PyObject *self, PyObject *args)
834 Py_complex x;
835 Py_complex y;
837 if (!PyArg_ParseTuple(args, "D|D", &x, &y))
838 return NULL;
840 errno = 0;
841 PyFPE_START_PROTECT("complex function", return 0)
842 x = c_log(x);
843 if (PyTuple_GET_SIZE(args) == 2) {
844 y = c_log(y);
845 x = c_quot(x, y);
847 PyFPE_END_PROTECT(x)
848 if (errno != 0)
849 return math_error();
850 return PyComplex_FromCComplex(x);
853 PyDoc_STRVAR(cmath_log_doc,
854 "log(x[, base]) -> the logarithm of x to the given base.\n\
855 If the base not specified, returns the natural logarithm (base e) of x.");
858 /* And now the glue to make them available from Python: */
860 static PyObject *
861 math_error(void)
863 if (errno == EDOM)
864 PyErr_SetString(PyExc_ValueError, "math domain error");
865 else if (errno == ERANGE)
866 PyErr_SetString(PyExc_OverflowError, "math range error");
867 else /* Unexpected math error */
868 PyErr_SetFromErrno(PyExc_ValueError);
869 return NULL;
872 static PyObject *
873 math_1(PyObject *args, Py_complex (*func)(Py_complex))
875 Py_complex x,r ;
876 if (!PyArg_ParseTuple(args, "D", &x))
877 return NULL;
878 errno = 0;
879 PyFPE_START_PROTECT("complex function", return 0);
880 r = (*func)(x);
881 PyFPE_END_PROTECT(r);
882 if (errno == EDOM) {
883 PyErr_SetString(PyExc_ValueError, "math domain error");
884 return NULL;
886 else if (errno == ERANGE) {
887 PyErr_SetString(PyExc_OverflowError, "math range error");
888 return NULL;
890 else {
891 return PyComplex_FromCComplex(r);
895 #define FUNC1(stubname, func) \
896 static PyObject * stubname(PyObject *self, PyObject *args) { \
897 return math_1(args, func); \
900 FUNC1(cmath_acos, c_acos)
901 FUNC1(cmath_acosh, c_acosh)
902 FUNC1(cmath_asin, c_asin)
903 FUNC1(cmath_asinh, c_asinh)
904 FUNC1(cmath_atan, c_atan)
905 FUNC1(cmath_atanh, c_atanh)
906 FUNC1(cmath_cos, c_cos)
907 FUNC1(cmath_cosh, c_cosh)
908 FUNC1(cmath_exp, c_exp)
909 FUNC1(cmath_log10, c_log10)
910 FUNC1(cmath_sin, c_sin)
911 FUNC1(cmath_sinh, c_sinh)
912 FUNC1(cmath_sqrt, c_sqrt)
913 FUNC1(cmath_tan, c_tan)
914 FUNC1(cmath_tanh, c_tanh)
916 static PyObject *
917 cmath_phase(PyObject *self, PyObject *args)
919 Py_complex z;
920 double phi;
921 if (!PyArg_ParseTuple(args, "D:phase", &z))
922 return NULL;
923 errno = 0;
924 PyFPE_START_PROTECT("arg function", return 0)
925 phi = c_atan2(z);
926 PyFPE_END_PROTECT(phi)
927 if (errno != 0)
928 return math_error();
929 else
930 return PyFloat_FromDouble(phi);
933 PyDoc_STRVAR(cmath_phase_doc,
934 "phase(z) -> float\n\n\
935 Return argument, also known as the phase angle, of a complex.");
937 static PyObject *
938 cmath_polar(PyObject *self, PyObject *args)
940 Py_complex z;
941 double r, phi;
942 if (!PyArg_ParseTuple(args, "D:polar", &z))
943 return NULL;
944 PyFPE_START_PROTECT("polar function", return 0)
945 phi = c_atan2(z); /* should not cause any exception */
946 r = c_abs(z); /* sets errno to ERANGE on overflow; otherwise 0 */
947 PyFPE_END_PROTECT(r)
948 if (errno != 0)
949 return math_error();
950 else
951 return Py_BuildValue("dd", r, phi);
954 PyDoc_STRVAR(cmath_polar_doc,
955 "polar(z) -> r: float, phi: float\n\n\
956 Convert a complex from rectangular coordinates to polar coordinates. r is\n\
957 the distance from 0 and phi the phase angle.");
960 rect() isn't covered by the C99 standard, but it's not too hard to
961 figure out 'spirit of C99' rules for special value handing:
963 rect(x, t) should behave like exp(log(x) + it) for positive-signed x
964 rect(x, t) should behave like -exp(log(-x) + it) for negative-signed x
965 rect(nan, t) should behave like exp(nan + it), except that rect(nan, 0)
966 gives nan +- i0 with the sign of the imaginary part unspecified.
970 static Py_complex rect_special_values[7][7];
972 static PyObject *
973 cmath_rect(PyObject *self, PyObject *args)
975 Py_complex z;
976 double r, phi;
977 if (!PyArg_ParseTuple(args, "dd:rect", &r, &phi))
978 return NULL;
979 errno = 0;
980 PyFPE_START_PROTECT("rect function", return 0)
982 /* deal with special values */
983 if (!Py_IS_FINITE(r) || !Py_IS_FINITE(phi)) {
984 /* if r is +/-infinity and phi is finite but nonzero then
985 result is (+-INF +-INF i), but we need to compute cos(phi)
986 and sin(phi) to figure out the signs. */
987 if (Py_IS_INFINITY(r) && (Py_IS_FINITE(phi)
988 && (phi != 0.))) {
989 if (r > 0) {
990 z.real = copysign(INF, cos(phi));
991 z.imag = copysign(INF, sin(phi));
993 else {
994 z.real = -copysign(INF, cos(phi));
995 z.imag = -copysign(INF, sin(phi));
998 else {
999 z = rect_special_values[special_type(r)]
1000 [special_type(phi)];
1002 /* need to set errno = EDOM if r is a nonzero number and phi
1003 is infinite */
1004 if (r != 0. && !Py_IS_NAN(r) && Py_IS_INFINITY(phi))
1005 errno = EDOM;
1006 else
1007 errno = 0;
1009 else {
1010 z.real = r * cos(phi);
1011 z.imag = r * sin(phi);
1012 errno = 0;
1015 PyFPE_END_PROTECT(z)
1016 if (errno != 0)
1017 return math_error();
1018 else
1019 return PyComplex_FromCComplex(z);
1022 PyDoc_STRVAR(cmath_rect_doc,
1023 "rect(r, phi) -> z: complex\n\n\
1024 Convert from polar coordinates to rectangular coordinates.");
1026 static PyObject *
1027 cmath_isnan(PyObject *self, PyObject *args)
1029 Py_complex z;
1030 if (!PyArg_ParseTuple(args, "D:isnan", &z))
1031 return NULL;
1032 return PyBool_FromLong(Py_IS_NAN(z.real) || Py_IS_NAN(z.imag));
1035 PyDoc_STRVAR(cmath_isnan_doc,
1036 "isnan(z) -> bool\n\
1037 Checks if the real or imaginary part of z not a number (NaN)");
1039 static PyObject *
1040 cmath_isinf(PyObject *self, PyObject *args)
1042 Py_complex z;
1043 if (!PyArg_ParseTuple(args, "D:isnan", &z))
1044 return NULL;
1045 return PyBool_FromLong(Py_IS_INFINITY(z.real) ||
1046 Py_IS_INFINITY(z.imag));
1049 PyDoc_STRVAR(cmath_isinf_doc,
1050 "isinf(z) -> bool\n\
1051 Checks if the real or imaginary part of z is infinite.");
1054 PyDoc_STRVAR(module_doc,
1055 "This module is always available. It provides access to mathematical\n"
1056 "functions for complex numbers.");
1058 static PyMethodDef cmath_methods[] = {
1059 {"acos", cmath_acos, METH_VARARGS, c_acos_doc},
1060 {"acosh", cmath_acosh, METH_VARARGS, c_acosh_doc},
1061 {"asin", cmath_asin, METH_VARARGS, c_asin_doc},
1062 {"asinh", cmath_asinh, METH_VARARGS, c_asinh_doc},
1063 {"atan", cmath_atan, METH_VARARGS, c_atan_doc},
1064 {"atanh", cmath_atanh, METH_VARARGS, c_atanh_doc},
1065 {"cos", cmath_cos, METH_VARARGS, c_cos_doc},
1066 {"cosh", cmath_cosh, METH_VARARGS, c_cosh_doc},
1067 {"exp", cmath_exp, METH_VARARGS, c_exp_doc},
1068 {"isinf", cmath_isinf, METH_VARARGS, cmath_isinf_doc},
1069 {"isnan", cmath_isnan, METH_VARARGS, cmath_isnan_doc},
1070 {"log", cmath_log, METH_VARARGS, cmath_log_doc},
1071 {"log10", cmath_log10, METH_VARARGS, c_log10_doc},
1072 {"phase", cmath_phase, METH_VARARGS, cmath_phase_doc},
1073 {"polar", cmath_polar, METH_VARARGS, cmath_polar_doc},
1074 {"rect", cmath_rect, METH_VARARGS, cmath_rect_doc},
1075 {"sin", cmath_sin, METH_VARARGS, c_sin_doc},
1076 {"sinh", cmath_sinh, METH_VARARGS, c_sinh_doc},
1077 {"sqrt", cmath_sqrt, METH_VARARGS, c_sqrt_doc},
1078 {"tan", cmath_tan, METH_VARARGS, c_tan_doc},
1079 {"tanh", cmath_tanh, METH_VARARGS, c_tanh_doc},
1080 {NULL, NULL} /* sentinel */
1083 PyMODINIT_FUNC
1084 initcmath(void)
1086 PyObject *m;
1088 m = Py_InitModule3("cmath", cmath_methods, module_doc);
1089 if (m == NULL)
1090 return;
1092 PyModule_AddObject(m, "pi",
1093 PyFloat_FromDouble(Py_MATH_PI));
1094 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
1096 /* initialize special value tables */
1098 #define INIT_SPECIAL_VALUES(NAME, BODY) { Py_complex* p = (Py_complex*)NAME; BODY }
1099 #define C(REAL, IMAG) p->real = REAL; p->imag = IMAG; ++p;
1101 INIT_SPECIAL_VALUES(acos_special_values, {
1102 C(P34,INF) C(P,INF) C(P,INF) C(P,-INF) C(P,-INF) C(P34,-INF) C(N,INF)
1103 C(P12,INF) C(U,U) C(U,U) C(U,U) C(U,U) C(P12,-INF) C(N,N)
1104 C(P12,INF) C(U,U) C(P12,0.) C(P12,-0.) C(U,U) C(P12,-INF) C(P12,N)
1105 C(P12,INF) C(U,U) C(P12,0.) C(P12,-0.) C(U,U) C(P12,-INF) C(P12,N)
1106 C(P12,INF) C(U,U) C(U,U) C(U,U) C(U,U) C(P12,-INF) C(N,N)
1107 C(P14,INF) C(0.,INF) C(0.,INF) C(0.,-INF) C(0.,-INF) C(P14,-INF) C(N,INF)
1108 C(N,INF) C(N,N) C(N,N) C(N,N) C(N,N) C(N,-INF) C(N,N)
1111 INIT_SPECIAL_VALUES(acosh_special_values, {
1112 C(INF,-P34) C(INF,-P) C(INF,-P) C(INF,P) C(INF,P) C(INF,P34) C(INF,N)
1113 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1114 C(INF,-P12) C(U,U) C(0.,-P12) C(0.,P12) C(U,U) C(INF,P12) C(N,N)
1115 C(INF,-P12) C(U,U) C(0.,-P12) C(0.,P12) C(U,U) C(INF,P12) C(N,N)
1116 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1117 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
1118 C(INF,N) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,N) C(N,N)
1121 INIT_SPECIAL_VALUES(asinh_special_values, {
1122 C(-INF,-P14) C(-INF,-0.) C(-INF,-0.) C(-INF,0.) C(-INF,0.) C(-INF,P14) C(-INF,N)
1123 C(-INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(-INF,P12) C(N,N)
1124 C(-INF,-P12) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(-INF,P12) C(N,N)
1125 C(INF,-P12) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(INF,P12) C(N,N)
1126 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1127 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
1128 C(INF,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(INF,N) C(N,N)
1131 INIT_SPECIAL_VALUES(atanh_special_values, {
1132 C(-0.,-P12) C(-0.,-P12) C(-0.,-P12) C(-0.,P12) C(-0.,P12) C(-0.,P12) C(-0.,N)
1133 C(-0.,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(-0.,P12) C(N,N)
1134 C(-0.,-P12) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(-0.,P12) C(-0.,N)
1135 C(0.,-P12) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,P12) C(0.,N)
1136 C(0.,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(0.,P12) C(N,N)
1137 C(0.,-P12) C(0.,-P12) C(0.,-P12) C(0.,P12) C(0.,P12) C(0.,P12) C(0.,N)
1138 C(0.,-P12) C(N,N) C(N,N) C(N,N) C(N,N) C(0.,P12) C(N,N)
1141 INIT_SPECIAL_VALUES(cosh_special_values, {
1142 C(INF,N) C(U,U) C(INF,0.) C(INF,-0.) C(U,U) C(INF,N) C(INF,N)
1143 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1144 C(N,0.) C(U,U) C(1.,0.) C(1.,-0.) C(U,U) C(N,0.) C(N,0.)
1145 C(N,0.) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(N,0.) C(N,0.)
1146 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1147 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
1148 C(N,N) C(N,N) C(N,0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1151 INIT_SPECIAL_VALUES(exp_special_values, {
1152 C(0.,0.) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,0.) C(0.,0.)
1153 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1154 C(N,N) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(N,N) C(N,N)
1155 C(N,N) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(N,N) C(N,N)
1156 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1157 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
1158 C(N,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1161 INIT_SPECIAL_VALUES(log_special_values, {
1162 C(INF,-P34) C(INF,-P) C(INF,-P) C(INF,P) C(INF,P) C(INF,P34) C(INF,N)
1163 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1164 C(INF,-P12) C(U,U) C(-INF,-P) C(-INF,P) C(U,U) C(INF,P12) C(N,N)
1165 C(INF,-P12) C(U,U) C(-INF,-0.) C(-INF,0.) C(U,U) C(INF,P12) C(N,N)
1166 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1167 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
1168 C(INF,N) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,N) C(N,N)
1171 INIT_SPECIAL_VALUES(sinh_special_values, {
1172 C(INF,N) C(U,U) C(-INF,-0.) C(-INF,0.) C(U,U) C(INF,N) C(INF,N)
1173 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1174 C(0.,N) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(0.,N) C(0.,N)
1175 C(0.,N) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,N) C(0.,N)
1176 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1177 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
1178 C(N,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1181 INIT_SPECIAL_VALUES(sqrt_special_values, {
1182 C(INF,-INF) C(0.,-INF) C(0.,-INF) C(0.,INF) C(0.,INF) C(INF,INF) C(N,INF)
1183 C(INF,-INF) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,INF) C(N,N)
1184 C(INF,-INF) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(INF,INF) C(N,N)
1185 C(INF,-INF) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(INF,INF) C(N,N)
1186 C(INF,-INF) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,INF) C(N,N)
1187 C(INF,-INF) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,INF) C(INF,N)
1188 C(INF,-INF) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,INF) C(N,N)
1191 INIT_SPECIAL_VALUES(tanh_special_values, {
1192 C(-1.,0.) C(U,U) C(-1.,-0.) C(-1.,0.) C(U,U) C(-1.,0.) C(-1.,0.)
1193 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1194 C(N,N) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(N,N) C(N,N)
1195 C(N,N) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(N,N) C(N,N)
1196 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1197 C(1.,0.) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(1.,0.) C(1.,0.)
1198 C(N,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1201 INIT_SPECIAL_VALUES(rect_special_values, {
1202 C(INF,N) C(U,U) C(-INF,0.) C(-INF,-0.) C(U,U) C(INF,N) C(INF,N)
1203 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1204 C(0.,0.) C(U,U) C(-0.,0.) C(-0.,-0.) C(U,U) C(0.,0.) C(0.,0.)
1205 C(0.,0.) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,0.) C(0.,0.)
1206 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1207 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
1208 C(N,N) C(N,N) C(N,0.) C(N,0.) C(N,N) C(N,N) C(N,N)