Branch data Line data Source code
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]
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 : : /* 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 : :
97 : : static double
98 : 4 : m_inf(void)
99 : : {
100 : : #if _PY_SHORT_FLOAT_REPR == 1
101 : 4 : return _Py_dg_infinity(0);
102 : : #else
103 : : return Py_HUGE_VAL;
104 : : #endif
105 : : }
106 : :
107 : : static Py_complex
108 : 2 : c_infj(void)
109 : : {
110 : : Py_complex r;
111 : 2 : r.real = 0.0;
112 : 2 : r.imag = m_inf();
113 : 2 : return r;
114 : : }
115 : :
116 : : #if _PY_SHORT_FLOAT_REPR == 1
117 : :
118 : : static double
119 : 4 : m_nan(void)
120 : : {
121 : : #if _PY_SHORT_FLOAT_REPR == 1
122 : 4 : return _Py_dg_stdnan(0);
123 : : #else
124 : : return Py_NAN;
125 : : #endif
126 : : }
127 : :
128 : : static Py_complex
129 : 2 : c_nanj(void)
130 : : {
131 : : Py_complex r;
132 : 2 : r.real = 0.0;
133 : 2 : r.imag = m_nan();
134 : 2 : return r;
135 : : }
136 : :
137 : : #endif
138 : :
139 : : /* forward declarations */
140 : : static Py_complex cmath_asinh_impl(PyObject *, Py_complex);
141 : : static Py_complex cmath_atanh_impl(PyObject *, Py_complex);
142 : : static Py_complex cmath_cosh_impl(PyObject *, Py_complex);
143 : : static Py_complex cmath_sinh_impl(PyObject *, Py_complex);
144 : : static Py_complex cmath_sqrt_impl(PyObject *, Py_complex);
145 : : static Py_complex cmath_tanh_impl(PyObject *, Py_complex);
146 : : static 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 : :
154 : : enum 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 : :
164 : : static enum special_types
165 : 1058 : special_type(double d)
166 : : {
167 [ + + ]: 1058 : if (Py_IS_FINITE(d)) {
168 [ + + ]: 376 : if (d != 0) {
169 [ + + ]: 172 : if (copysign(1., d) == 1.)
170 : 86 : return ST_POS;
171 : : else
172 : 86 : return ST_NEG;
173 : : }
174 : : else {
175 [ + + ]: 204 : if (copysign(1., d) == 1.)
176 : 102 : return ST_PZERO;
177 : : else
178 : 102 : return ST_NZERO;
179 : : }
180 : : }
181 [ + + ]: 682 : if (Py_IS_NAN(d))
182 : 238 : return ST_NAN;
183 [ + + ]: 444 : if (copysign(1., d) == 1.)
184 : 222 : return ST_PINF;
185 : : else
186 : 222 : 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 : :
213 : : static Py_complex acos_special_values[7][7];
214 : :
215 : : /*[clinic input]
216 : : cmath.acos -> Py_complex_protected
217 : :
218 : : z: Py_complex_protected
219 : : /
220 : :
221 : : Return the arc cosine of z.
222 : : [clinic start generated code]*/
223 : :
224 : : static Py_complex
225 : 186 : cmath_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 [ + + + + ]: 186 : SPECIAL_VALUE(z, acos_special_values);
231 : :
232 [ + + + + ]: 153 : if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
233 : : /* avoid unnecessary overflow for large arguments */
234 : 20 : 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 [ + + ]: 20 : if (z.real < 0.) {
238 : 8 : r.imag = -copysign(log(hypot(z.real/2., z.imag/2.)) +
239 : : M_LN2*2., z.imag);
240 : : } else {
241 : 12 : r.imag = copysign(log(hypot(z.real/2., z.imag/2.)) +
242 : 12 : M_LN2*2., -z.imag);
243 : : }
244 : : } else {
245 : 133 : s1.real = 1.-z.real;
246 : 133 : s1.imag = -z.imag;
247 : 133 : s1 = cmath_sqrt_impl(module, s1);
248 : 133 : s2.real = 1.+z.real;
249 : 133 : s2.imag = z.imag;
250 : 133 : s2 = cmath_sqrt_impl(module, s2);
251 : 133 : r.real = 2.*atan2(s1.real, s2.real);
252 : 133 : r.imag = asinh(s2.real*s1.imag - s2.imag*s1.real);
253 : : }
254 : 153 : errno = 0;
255 : 153 : return r;
256 : : }
257 : :
258 : :
259 : : static Py_complex acosh_special_values[7][7];
260 : :
261 : : /*[clinic input]
262 : : cmath.acosh = cmath.acos
263 : :
264 : : Return the inverse hyperbolic cosine of z.
265 : : [clinic start generated code]*/
266 : :
267 : : static Py_complex
268 : 171 : cmath_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 [ + + + + ]: 171 : SPECIAL_VALUE(z, acosh_special_values);
274 : :
275 [ + + + + ]: 138 : if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
276 : : /* avoid unnecessary overflow for large arguments */
277 : 20 : r.real = log(hypot(z.real/2., z.imag/2.)) + M_LN2*2.;
278 : 20 : r.imag = atan2(z.imag, z.real);
279 : : } else {
280 : 118 : s1.real = z.real - 1.;
281 : 118 : s1.imag = z.imag;
282 : 118 : s1 = cmath_sqrt_impl(module, s1);
283 : 118 : s2.real = z.real + 1.;
284 : 118 : s2.imag = z.imag;
285 : 118 : s2 = cmath_sqrt_impl(module, s2);
286 : 118 : r.real = asinh(s1.real*s2.real + s1.imag*s2.imag);
287 : 118 : r.imag = 2.*atan2(s1.imag, s2.real);
288 : : }
289 : 138 : errno = 0;
290 : 138 : return r;
291 : : }
292 : :
293 : : /*[clinic input]
294 : : cmath.asin = cmath.acos
295 : :
296 : : Return the arc sine of z.
297 : : [clinic start generated code]*/
298 : :
299 : : static Py_complex
300 : 174 : cmath_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 : 174 : s.real = -z.imag;
306 : 174 : s.imag = z.real;
307 : 174 : s = cmath_asinh_impl(module, s);
308 : 174 : r.real = s.imag;
309 : 174 : r.imag = -s.real;
310 : 174 : return r;
311 : : }
312 : :
313 : :
314 : : static Py_complex asinh_special_values[7][7];
315 : :
316 : : /*[clinic input]
317 : : cmath.asinh = cmath.acos
318 : :
319 : : Return the inverse hyperbolic sine of z.
320 : : [clinic start generated code]*/
321 : :
322 : : static Py_complex
323 : 345 : cmath_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 [ + + + + ]: 345 : SPECIAL_VALUE(z, asinh_special_values);
329 : :
330 [ + + + + ]: 279 : if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
331 [ + + ]: 40 : if (z.imag >= 0.) {
332 : 24 : r.real = copysign(log(hypot(z.real/2., z.imag/2.)) +
333 : : M_LN2*2., z.real);
334 : : } else {
335 : 16 : r.real = -copysign(log(hypot(z.real/2., z.imag/2.)) +
336 : 16 : M_LN2*2., -z.real);
337 : : }
338 : 40 : r.imag = atan2(z.imag, fabs(z.real));
339 : : } else {
340 : 239 : s1.real = 1.+z.imag;
341 : 239 : s1.imag = -z.real;
342 : 239 : s1 = cmath_sqrt_impl(module, s1);
343 : 239 : s2.real = 1.-z.imag;
344 : 239 : s2.imag = z.real;
345 : 239 : s2 = cmath_sqrt_impl(module, s2);
346 : 239 : r.real = asinh(s1.real*s2.imag-s2.real*s1.imag);
347 : 239 : r.imag = atan2(z.imag, s1.real*s2.real-s1.imag*s2.imag);
348 : : }
349 : 279 : errno = 0;
350 : 279 : return r;
351 : : }
352 : :
353 : :
354 : : /*[clinic input]
355 : : cmath.atan = cmath.acos
356 : :
357 : : Return the arc tangent of z.
358 : : [clinic start generated code]*/
359 : :
360 : : static Py_complex
361 : 201 : cmath_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 : 201 : s.real = -z.imag;
367 : 201 : s.imag = z.real;
368 : 201 : s = cmath_atanh_impl(module, s);
369 : 201 : r.real = s.imag;
370 : 201 : r.imag = -s.real;
371 : 201 : 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.). */
376 : : static double
377 : 121 : c_atan2(Py_complex z)
378 : : {
379 [ + + + + ]: 121 : if (Py_IS_NAN(z.real) || Py_IS_NAN(z.imag))
380 : 39 : return Py_NAN;
381 [ + + + + ]: 82 : if (Py_IS_INFINITY(z.imag)) {
382 [ + + + + ]: 36 : if (Py_IS_INFINITY(z.real)) {
383 [ + + ]: 16 : if (copysign(1., z.real) == 1.)
384 : : /* atan2(+-inf, +inf) == +-pi/4 */
385 : 8 : return copysign(0.25*Py_MATH_PI, z.imag);
386 : : else
387 : : /* atan2(+-inf, -inf) == +-pi*3/4 */
388 : 8 : return copysign(0.75*Py_MATH_PI, z.imag);
389 : : }
390 : : /* atan2(+-inf, x) == +-pi/2 for finite x */
391 : 20 : return copysign(0.5*Py_MATH_PI, z.imag);
392 : : }
393 [ + + + + ]: 46 : if (Py_IS_INFINITY(z.real) || z.imag == 0.) {
394 [ + + ]: 37 : if (copysign(1., z.real) == 1.)
395 : : /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
396 : 20 : return copysign(0., z.imag);
397 : : else
398 : : /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
399 : 17 : return copysign(Py_MATH_PI, z.imag);
400 : : }
401 : 9 : return atan2(z.imag, z.real);
402 : : }
403 : :
404 : :
405 : : static Py_complex atanh_special_values[7][7];
406 : :
407 : : /*[clinic input]
408 : : cmath.atanh = cmath.acos
409 : :
410 : : Return the inverse hyperbolic tangent of z.
411 : : [clinic start generated code]*/
412 : :
413 : : static Py_complex
414 : 493 : cmath_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 [ + + + + ]: 493 : SPECIAL_VALUE(z, atanh_special_values);
421 : :
422 : : /* Reduce to case where z.real >= 0., using atanh(z) = -atanh(-z). */
423 [ + + ]: 427 : if (z.real < 0.) {
424 : 116 : return _Py_c_neg(cmath_atanh_impl(module, _Py_c_neg(z)));
425 : : }
426 : :
427 : 311 : ay = fabs(z.imag);
428 [ + + + + ]: 311 : 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 : 50 : h = hypot(z.real/2., z.imag/2.); /* safe from overflow */
435 : 50 : 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 : 50 : r.imag = -copysign(Py_MATH_PI/2., -z.imag);
441 : 50 : errno = 0;
442 [ + + + + ]: 261 : } else if (z.real == 1. && ay < CM_SQRT_DBL_MIN) {
443 : : /* C99 standard says: atanh(1+/-0.) should be inf +/- 0i */
444 [ + + ]: 16 : if (ay == 0.) {
445 : 8 : r.real = INF;
446 : 8 : r.imag = z.imag;
447 : 8 : errno = EDOM;
448 : : } else {
449 : 8 : r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.)));
450 : 8 : r.imag = copysign(atan2(2., -ay)/2, z.imag);
451 : 8 : errno = 0;
452 : : }
453 : : } else {
454 : 245 : r.real = m_log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.;
455 : 245 : r.imag = -atan2(-2.*z.imag, (1-z.real)*(1+z.real) - ay*ay)/2.;
456 : 245 : errno = 0;
457 : : }
458 : 311 : return r;
459 : : }
460 : :
461 : :
462 : : /*[clinic input]
463 : : cmath.cos = cmath.acos
464 : :
465 : : Return the cosine of z.
466 : : [clinic start generated code]*/
467 : :
468 : : static Py_complex
469 : 136 : cmath_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 : 136 : r.real = -z.imag;
475 : 136 : r.imag = z.real;
476 : 136 : r = cmath_cosh_impl(module, r);
477 : 136 : return r;
478 : : }
479 : :
480 : :
481 : : /* cosh(infinity + i*y) needs to be dealt with specially */
482 : : static Py_complex cosh_special_values[7][7];
483 : :
484 : : /*[clinic input]
485 : : cmath.cosh = cmath.acos
486 : :
487 : : Return the hyperbolic cosine of z.
488 : : [clinic start generated code]*/
489 : :
490 : : static Py_complex
491 : 275 : cmath_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 [ + + + + ]: 275 : if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
499 [ + + + + ]: 98 : if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) &&
500 [ + + ]: 48 : (z.imag != 0.)) {
501 [ + + ]: 40 : if (z.real > 0) {
502 : 20 : r.real = copysign(INF, cos(z.imag));
503 : 20 : r.imag = copysign(INF, sin(z.imag));
504 : : }
505 : : else {
506 : 20 : r.real = copysign(INF, cos(z.imag));
507 : 20 : r.imag = -copysign(INF, sin(z.imag));
508 : : }
509 : : }
510 : : else {
511 : 116 : r = cosh_special_values[special_type(z.real)]
512 : 58 : [special_type(z.imag)];
513 : : }
514 : : /* need to set errno = EDOM if y is +/- infinity and x is not
515 : : a NaN */
516 [ + + + + ]: 98 : if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
517 : 24 : errno = EDOM;
518 : : else
519 : 74 : errno = 0;
520 : 98 : return r;
521 : : }
522 : :
523 [ + + ]: 177 : 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 : 4 : x_minus_one = z.real - copysign(1., z.real);
527 : 4 : r.real = cos(z.imag) * cosh(x_minus_one) * Py_MATH_E;
528 : 4 : r.imag = sin(z.imag) * sinh(x_minus_one) * Py_MATH_E;
529 : : } else {
530 : 173 : r.real = cos(z.imag) * cosh(z.real);
531 : 173 : r.imag = sin(z.imag) * sinh(z.real);
532 : : }
533 : : /* detect overflow, and set errno accordingly */
534 [ - + - + ]: 177 : if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
535 : 0 : errno = ERANGE;
536 : : else
537 : 177 : errno = 0;
538 : 177 : return r;
539 : : }
540 : :
541 : :
542 : : /* exp(infinity + i*y) and exp(-infinity + i*y) need special treatment for
543 : : finite y */
544 : : static Py_complex exp_special_values[7][7];
545 : :
546 : : /*[clinic input]
547 : : cmath.exp = cmath.acos
548 : :
549 : : Return the exponential value e**z.
550 : : [clinic start generated code]*/
551 : :
552 : : static Py_complex
553 : 148 : cmath_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 [ + + + + ]: 148 : if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
560 [ + + + + ]: 49 : if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
561 [ + + ]: 24 : && (z.imag != 0.)) {
562 [ + + ]: 20 : if (z.real > 0) {
563 : 10 : r.real = copysign(INF, cos(z.imag));
564 : 10 : r.imag = copysign(INF, sin(z.imag));
565 : : }
566 : : else {
567 : 10 : r.real = copysign(0., cos(z.imag));
568 : 10 : r.imag = copysign(0., sin(z.imag));
569 : : }
570 : : }
571 : : else {
572 : 58 : r = exp_special_values[special_type(z.real)]
573 : 29 : [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 [ + + ]: 49 : if (Py_IS_INFINITY(z.imag) &&
578 [ + + ]: 14 : (Py_IS_FINITE(z.real) ||
579 [ + + + + ]: 6 : (Py_IS_INFINITY(z.real) && z.real > 0)))
580 : 10 : errno = EDOM;
581 : : else
582 : 39 : errno = 0;
583 : 49 : return r;
584 : : }
585 : :
586 [ + + ]: 99 : if (z.real > CM_LOG_LARGE_DOUBLE) {
587 : 8 : l = exp(z.real-1.);
588 : 8 : r.real = l*cos(z.imag)*Py_MATH_E;
589 : 8 : r.imag = l*sin(z.imag)*Py_MATH_E;
590 : : } else {
591 : 91 : l = exp(z.real);
592 : 91 : r.real = l*cos(z.imag);
593 : 91 : r.imag = l*sin(z.imag);
594 : : }
595 : : /* detect overflow, and set errno accordingly */
596 [ + + + + ]: 99 : if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
597 : 5 : errno = ERANGE;
598 : : else
599 : 94 : errno = 0;
600 : 99 : return r;
601 : : }
602 : :
603 : : static Py_complex log_special_values[7][7];
604 : :
605 : : static Py_complex
606 : 538 : c_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 [ + + + + ]: 538 : SPECIAL_VALUE(z, log_special_values);
640 : :
641 : 472 : ax = fabs(z.real);
642 : 472 : ay = fabs(z.imag);
643 : :
644 [ + + + + ]: 472 : if (ax > CM_LARGE_DOUBLE || ay > CM_LARGE_DOUBLE) {
645 : 40 : r.real = log(hypot(ax/2., ay/2.)) + M_LN2;
646 [ + + + + ]: 432 : } else if (ax < DBL_MIN && ay < DBL_MIN) {
647 [ + + + + ]: 36 : if (ax > 0. || ay > 0.) {
648 : : /* catch cases where hypot(ax, ay) is subnormal */
649 : 28 : r.real = log(hypot(ldexp(ax, DBL_MANT_DIG),
650 : 28 : ldexp(ay, DBL_MANT_DIG))) - DBL_MANT_DIG*M_LN2;
651 : : }
652 : : else {
653 : : /* log(+/-0. +/- 0i) */
654 : 8 : r.real = -INF;
655 : 8 : r.imag = atan2(z.imag, z.real);
656 : 8 : errno = EDOM;
657 : 8 : return r;
658 : : }
659 : : } else {
660 : 396 : h = hypot(ax, ay);
661 [ + + + + ]: 396 : if (0.71 <= h && h <= 1.73) {
662 [ + + ]: 65 : am = ax > ay ? ax : ay; /* max(ax, ay) */
663 [ + + ]: 65 : an = ax > ay ? ay : ax; /* min(ax, ay) */
664 : 65 : r.real = m_log1p((am-1)*(am+1)+an*an)/2.;
665 : : } else {
666 : 331 : r.real = log(h);
667 : : }
668 : : }
669 : 464 : r.imag = atan2(z.imag, z.real);
670 : 464 : errno = 0;
671 : 464 : return r;
672 : : }
673 : :
674 : :
675 : : /*[clinic input]
676 : : cmath.log10 = cmath.acos
677 : :
678 : : Return the base-10 logarithm of z.
679 : : [clinic start generated code]*/
680 : :
681 : : static Py_complex
682 : 182 : cmath_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 : 182 : r = c_log(z);
689 : 182 : errno_save = errno; /* just in case the divisions affect errno */
690 : 182 : r.real = r.real / M_LN10;
691 : 182 : r.imag = r.imag / M_LN10;
692 : 182 : errno = errno_save;
693 : 182 : return r;
694 : : }
695 : :
696 : :
697 : : /*[clinic input]
698 : : cmath.sin = cmath.acos
699 : :
700 : : Return the sine of z.
701 : : [clinic start generated code]*/
702 : :
703 : : static Py_complex
704 : 136 : cmath_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 : 136 : s.real = -z.imag;
710 : 136 : s.imag = z.real;
711 : 136 : s = cmath_sinh_impl(module, s);
712 : 136 : r.real = s.imag;
713 : 136 : r.imag = -s.real;
714 : 136 : return r;
715 : : }
716 : :
717 : :
718 : : /* sinh(infinity + i*y) needs to be dealt with specially */
719 : : static Py_complex sinh_special_values[7][7];
720 : :
721 : : /*[clinic input]
722 : : cmath.sinh = cmath.acos
723 : :
724 : : Return the hyperbolic sine of z.
725 : : [clinic start generated code]*/
726 : :
727 : : static Py_complex
728 : 276 : cmath_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 [ + + + + ]: 276 : if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
737 [ + + + + ]: 98 : if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
738 [ + + ]: 48 : && (z.imag != 0.)) {
739 [ + + ]: 40 : if (z.real > 0) {
740 : 20 : r.real = copysign(INF, cos(z.imag));
741 : 20 : r.imag = copysign(INF, sin(z.imag));
742 : : }
743 : : else {
744 : 20 : r.real = -copysign(INF, cos(z.imag));
745 : 20 : r.imag = copysign(INF, sin(z.imag));
746 : : }
747 : : }
748 : : else {
749 : 116 : r = sinh_special_values[special_type(z.real)]
750 : 58 : [special_type(z.imag)];
751 : : }
752 : : /* need to set errno = EDOM if y is +/- infinity and x is not
753 : : a NaN */
754 [ + + + + ]: 98 : if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
755 : 24 : errno = EDOM;
756 : : else
757 : 74 : errno = 0;
758 : 98 : return r;
759 : : }
760 : :
761 [ + + ]: 178 : if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
762 : 4 : x_minus_one = z.real - copysign(1., z.real);
763 : 4 : r.real = cos(z.imag) * sinh(x_minus_one) * Py_MATH_E;
764 : 4 : r.imag = sin(z.imag) * cosh(x_minus_one) * Py_MATH_E;
765 : : } else {
766 : 174 : r.real = cos(z.imag) * sinh(z.real);
767 : 174 : r.imag = sin(z.imag) * cosh(z.real);
768 : : }
769 : : /* detect overflow, and set errno accordingly */
770 [ - + - + ]: 178 : if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
771 : 0 : errno = ERANGE;
772 : : else
773 : 178 : errno = 0;
774 : 178 : return r;
775 : : }
776 : :
777 : :
778 : : static Py_complex sqrt_special_values[7][7];
779 : :
780 : : /*[clinic input]
781 : : cmath.sqrt = cmath.acos
782 : :
783 : : Return the square root of z.
784 : : [clinic start generated code]*/
785 : :
786 : : static Py_complex
787 : 1146 : cmath_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 [ + + + + ]: 1146 : SPECIAL_VALUE(z, sqrt_special_values);
822 : :
823 [ + + + + ]: 1113 : if (z.real == 0. && z.imag == 0.) {
824 : 29 : r.real = 0.;
825 : 29 : r.imag = z.imag;
826 : 29 : return r;
827 : : }
828 : :
829 : 1084 : ax = fabs(z.real);
830 : 1084 : ay = fabs(z.imag);
831 : :
832 [ + + + - : 1084 : if (ax < DBL_MIN && ay < DBL_MIN && (ax > 0. || ay > 0.)) {
+ + + - ]
833 : : /* here we catch cases where hypot(ax, ay) is subnormal */
834 : 16 : ax = ldexp(ax, CM_SCALE_UP);
835 : 16 : s = ldexp(sqrt(ax + hypot(ax, ldexp(ay, CM_SCALE_UP))),
836 : : CM_SCALE_DOWN);
837 : : } else {
838 : 1068 : ax /= 8.;
839 : 1068 : s = 2.*sqrt(ax + hypot(ax, ay/8.));
840 : : }
841 : 1084 : d = ay/(2.*s);
842 : :
843 [ + + ]: 1084 : if (z.real >= 0.) {
844 : 767 : r.real = s;
845 : 767 : r.imag = copysign(d, z.imag);
846 : : } else {
847 : 317 : r.real = d;
848 : 317 : r.imag = copysign(s, z.imag);
849 : : }
850 : 1084 : errno = 0;
851 : 1084 : return r;
852 : : }
853 : :
854 : :
855 : : /*[clinic input]
856 : : cmath.tan = cmath.acos
857 : :
858 : : Return the tangent of z.
859 : : [clinic start generated code]*/
860 : :
861 : : static Py_complex
862 : 139 : cmath_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 : 139 : s.real = -z.imag;
868 : 139 : s.imag = z.real;
869 : 139 : s = cmath_tanh_impl(module, s);
870 : 139 : r.real = s.imag;
871 : 139 : r.imag = -s.real;
872 : 139 : return r;
873 : : }
874 : :
875 : :
876 : : /* tanh(infinity + i*y) needs to be dealt with specially */
877 : : static Py_complex tanh_special_values[7][7];
878 : :
879 : : /*[clinic input]
880 : : cmath.tanh = cmath.acos
881 : :
882 : : Return the hyperbolic tangent of z.
883 : : [clinic start generated code]*/
884 : :
885 : : static Py_complex
886 : 281 : cmath_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 [ + + + + ]: 281 : if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
907 [ + + + + ]: 98 : if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
908 [ + + ]: 48 : && (z.imag != 0.)) {
909 [ + + ]: 40 : if (z.real > 0) {
910 : 20 : r.real = 1.0;
911 : 20 : r.imag = copysign(0.,
912 : 20 : 2.*sin(z.imag)*cos(z.imag));
913 : : }
914 : : else {
915 : 20 : r.real = -1.0;
916 : 20 : r.imag = copysign(0.,
917 : 20 : 2.*sin(z.imag)*cos(z.imag));
918 : : }
919 : : }
920 : : else {
921 : 116 : r = tanh_special_values[special_type(z.real)]
922 : 58 : [special_type(z.imag)];
923 : : }
924 : : /* need to set errno = EDOM if z.imag is +/-infinity and
925 : : z.real is finite */
926 [ + + + + ]: 98 : if (Py_IS_INFINITY(z.imag) && Py_IS_FINITE(z.real))
927 : 16 : errno = EDOM;
928 : : else
929 : 82 : errno = 0;
930 : 98 : return r;
931 : : }
932 : :
933 : : /* danger of overflow in 2.*z.imag !*/
934 [ + + ]: 183 : if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
935 : 10 : r.real = copysign(1., z.real);
936 : 10 : r.imag = 4.*sin(z.imag)*cos(z.imag)*exp(-2.*fabs(z.real));
937 : : } else {
938 : 173 : tx = tanh(z.real);
939 : 173 : ty = tan(z.imag);
940 : 173 : cx = 1./cosh(z.real);
941 : 173 : txty = tx*ty;
942 : 173 : denom = 1. + txty*txty;
943 : 173 : r.real = tx*(1.+ty*ty)/denom;
944 : 173 : r.imag = ((ty/denom)*cx)*cx;
945 : : }
946 : 183 : errno = 0;
947 : 183 : return r;
948 : : }
949 : :
950 : :
951 : : /*[clinic input]
952 : : cmath.log
953 : :
954 : : z as x: Py_complex
955 : : base as y_obj: object = NULL
956 : : /
957 : :
958 : : log(z[, base]) -> the logarithm of z to the given base.
959 : :
960 : : If the base not specified, returns the natural logarithm (base e) of z.
961 : : [clinic start generated code]*/
962 : :
963 : : static PyObject *
964 : 281 : cmath_log_impl(PyObject *module, Py_complex x, PyObject *y_obj)
965 : : /*[clinic end generated code: output=4effdb7d258e0d94 input=230ed3a71ecd000a]*/
966 : : {
967 : : Py_complex y;
968 : :
969 : 281 : errno = 0;
970 : 281 : x = c_log(x);
971 [ + + ]: 281 : if (y_obj != NULL) {
972 : 99 : y = PyComplex_AsCComplex(y_obj);
973 [ + + ]: 99 : if (PyErr_Occurred()) {
974 : 24 : return NULL;
975 : : }
976 : 75 : y = c_log(y);
977 : 75 : x = _Py_c_quot(x, y);
978 : : }
979 [ + + ]: 257 : if (errno != 0)
980 : 4 : return math_error();
981 : 253 : return PyComplex_FromCComplex(x);
982 : : }
983 : :
984 : :
985 : : /* And now the glue to make them available from Python: */
986 : :
987 : : static PyObject *
988 : 13 : math_error(void)
989 : : {
990 [ + + ]: 13 : if (errno == EDOM)
991 : 12 : PyErr_SetString(PyExc_ValueError, "math domain error");
992 [ + - ]: 1 : else if (errno == ERANGE)
993 : 1 : PyErr_SetString(PyExc_OverflowError, "math range error");
994 : : else /* Unexpected math error */
995 : 0 : PyErr_SetFromErrno(PyExc_ValueError);
996 : 13 : return NULL;
997 : : }
998 : :
999 : :
1000 : : /*[clinic input]
1001 : : cmath.phase
1002 : :
1003 : : z: Py_complex
1004 : : /
1005 : :
1006 : : Return argument, also known as the phase angle, of a complex.
1007 : : [clinic start generated code]*/
1008 : :
1009 : : static PyObject *
1010 : 43 : cmath_phase_impl(PyObject *module, Py_complex z)
1011 : : /*[clinic end generated code: output=50725086a7bfd253 input=5cf75228ba94b69d]*/
1012 : : {
1013 : : double phi;
1014 : :
1015 : 43 : errno = 0;
1016 : 43 : phi = c_atan2(z);
1017 [ - + ]: 43 : if (errno != 0)
1018 : 0 : return math_error();
1019 : : else
1020 : 43 : return PyFloat_FromDouble(phi);
1021 : : }
1022 : :
1023 : : /*[clinic input]
1024 : : cmath.polar
1025 : :
1026 : : z: Py_complex
1027 : : /
1028 : :
1029 : : Convert a complex from rectangular coordinates to polar coordinates.
1030 : :
1031 : : r is the distance from 0 and phi the phase angle.
1032 : : [clinic start generated code]*/
1033 : :
1034 : : static PyObject *
1035 : 78 : cmath_polar_impl(PyObject *module, Py_complex z)
1036 : : /*[clinic end generated code: output=d0a8147c41dbb654 input=26c353574fd1a861]*/
1037 : : {
1038 : : double r, phi;
1039 : :
1040 : 78 : errno = 0;
1041 : 78 : phi = c_atan2(z); /* should not cause any exception */
1042 : 78 : r = _Py_c_abs(z); /* sets errno to ERANGE on overflow */
1043 [ + + ]: 78 : if (errno != 0)
1044 : 1 : return math_error();
1045 : : else
1046 : 77 : 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 : :
1060 : : static Py_complex rect_special_values[7][7];
1061 : :
1062 : : /*[clinic input]
1063 : : cmath.rect
1064 : :
1065 : : r: double
1066 : : phi: double
1067 : : /
1068 : :
1069 : : Convert from polar coordinates to rectangular coordinates.
1070 : : [clinic start generated code]*/
1071 : :
1072 : : static PyObject *
1073 : 58 : cmath_rect_impl(PyObject *module, double r, double phi)
1074 : : /*[clinic end generated code: output=385a0690925df2d5 input=24c5646d147efd69]*/
1075 : : {
1076 : : Py_complex z;
1077 : 58 : errno = 0;
1078 : :
1079 : : /* deal with special values */
1080 [ + + + + ]: 58 : 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 [ + + + + ]: 49 : if (Py_IS_INFINITY(r) && (Py_IS_FINITE(phi)
1085 [ + + ]: 24 : && (phi != 0.))) {
1086 [ + + ]: 20 : if (r > 0) {
1087 : 10 : z.real = copysign(INF, cos(phi));
1088 : 10 : z.imag = copysign(INF, sin(phi));
1089 : : }
1090 : : else {
1091 : 10 : z.real = -copysign(INF, cos(phi));
1092 : 10 : z.imag = -copysign(INF, sin(phi));
1093 : : }
1094 : : }
1095 : : else {
1096 : 58 : z = rect_special_values[special_type(r)]
1097 : 29 : [special_type(phi)];
1098 : : }
1099 : : /* need to set errno = EDOM if r is a nonzero number and phi
1100 : : is infinite */
1101 [ + + + + : 49 : if (r != 0. && !Py_IS_NAN(r) && Py_IS_INFINITY(phi))
+ + ]
1102 : 8 : errno = EDOM;
1103 : : else
1104 : 41 : errno = 0;
1105 : : }
1106 [ + + ]: 9 : 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 : 6 : z.real = r;
1110 : 6 : z.imag = r * phi;
1111 : 6 : errno = 0;
1112 : : }
1113 : : else {
1114 : 3 : z.real = r * cos(phi);
1115 : 3 : z.imag = r * sin(phi);
1116 : 3 : errno = 0;
1117 : : }
1118 : :
1119 [ + + ]: 58 : if (errno != 0)
1120 : 8 : return math_error();
1121 : : else
1122 : 50 : return PyComplex_FromCComplex(z);
1123 : : }
1124 : :
1125 : : /*[clinic input]
1126 : : cmath.isfinite = cmath.polar
1127 : :
1128 : : Return True if both the real and imaginary parts of z are finite, else False.
1129 : : [clinic start generated code]*/
1130 : :
1131 : : static PyObject *
1132 : 49 : cmath_isfinite_impl(PyObject *module, Py_complex z)
1133 : : /*[clinic end generated code: output=ac76611e2c774a36 input=848e7ee701895815]*/
1134 : : {
1135 [ + + + + ]: 49 : return PyBool_FromLong(Py_IS_FINITE(z.real) && Py_IS_FINITE(z.imag));
1136 : : }
1137 : :
1138 : : /*[clinic input]
1139 : : cmath.isnan = cmath.polar
1140 : :
1141 : : Checks if the real or imaginary part of z not a number (NaN).
1142 : : [clinic start generated code]*/
1143 : :
1144 : : static PyObject *
1145 : 9 : cmath_isnan_impl(PyObject *module, Py_complex z)
1146 : : /*[clinic end generated code: output=e7abf6e0b28beab7 input=71799f5d284c9baf]*/
1147 : : {
1148 [ + + + + ]: 9 : return PyBool_FromLong(Py_IS_NAN(z.real) || Py_IS_NAN(z.imag));
1149 : : }
1150 : :
1151 : : /*[clinic input]
1152 : : cmath.isinf = cmath.polar
1153 : :
1154 : : Checks if the real or imaginary part of z is infinite.
1155 : : [clinic start generated code]*/
1156 : :
1157 : : static PyObject *
1158 : 9 : cmath_isinf_impl(PyObject *module, Py_complex z)
1159 : : /*[clinic end generated code: output=502a75a79c773469 input=363df155c7181329]*/
1160 : : {
1161 [ + + ]: 14 : return PyBool_FromLong(Py_IS_INFINITY(z.real) ||
1162 [ + + ]: 5 : Py_IS_INFINITY(z.imag));
1163 : : }
1164 : :
1165 : : /*[clinic input]
1166 : : cmath.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 : :
1178 : : Determine whether two complex numbers are close in value.
1179 : :
1180 : : Return True if a is close in value to b, and False otherwise.
1181 : :
1182 : : For the values to be considered close, the difference between them must be
1183 : : smaller than at least one of the tolerances.
1184 : :
1185 : : -inf, inf and NaN behave similarly to the IEEE 754 Standard. That is, NaN is
1186 : : not close to anything, even itself. inf and -inf are only close to themselves.
1187 : : [clinic start generated code]*/
1188 : :
1189 : : static int
1190 : 82 : cmath_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 [ + + + + ]: 82 : if (rel_tol < 0.0 || abs_tol < 0.0 ) {
1198 : 2 : PyErr_SetString(PyExc_ValueError,
1199 : : "tolerances must be non-negative");
1200 : 2 : return -1;
1201 : : }
1202 : :
1203 [ + + + + ]: 80 : 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 : 13 : 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 [ + + - + ]: 67 : if (Py_IS_INFINITY(a.real) || Py_IS_INFINITY(a.imag) ||
1218 [ + + - + ]: 63 : Py_IS_INFINITY(b.real) || Py_IS_INFINITY(b.imag)) {
1219 : 7 : return 0;
1220 : : }
1221 : :
1222 : : /* now do the regular computation
1223 : : this is essentially the "weak" test from the Boost library
1224 : : */
1225 : :
1226 : 60 : diff = _Py_c_abs(_Py_c_diff(a, b));
1227 : :
1228 : 60 : return (((diff <= rel_tol * _Py_c_abs(b)) ||
1229 [ + + + + : 60 : (diff <= rel_tol * _Py_c_abs(a))) ||
+ + ]
1230 : : (diff <= abs_tol));
1231 : : }
1232 : :
1233 : : PyDoc_STRVAR(module_doc,
1234 : : "This module provides access to mathematical functions for complex\n"
1235 : : "numbers.");
1236 : :
1237 : : static 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 : :
1264 : : static int
1265 : 2 : cmath_exec(PyObject *mod)
1266 : : {
1267 [ - + ]: 2 : if (PyModule_AddObject(mod, "pi", PyFloat_FromDouble(Py_MATH_PI)) < 0) {
1268 : 0 : return -1;
1269 : : }
1270 [ - + ]: 2 : if (PyModule_AddObject(mod, "e", PyFloat_FromDouble(Py_MATH_E)) < 0) {
1271 : 0 : return -1;
1272 : : }
1273 : : // 2pi
1274 [ - + ]: 2 : if (PyModule_AddObject(mod, "tau", PyFloat_FromDouble(Py_MATH_TAU)) < 0) {
1275 : 0 : return -1;
1276 : : }
1277 [ - + ]: 2 : if (PyModule_AddObject(mod, "inf", PyFloat_FromDouble(m_inf())) < 0) {
1278 : 0 : return -1;
1279 : : }
1280 : :
1281 [ - + ]: 2 : if (PyModule_AddObject(mod, "infj",
1282 : : PyComplex_FromCComplex(c_infj())) < 0) {
1283 : 0 : return -1;
1284 : : }
1285 : : #if _PY_SHORT_FLOAT_REPR == 1
1286 [ - + ]: 2 : if (PyModule_AddObject(mod, "nan", PyFloat_FromDouble(m_nan())) < 0) {
1287 : 0 : return -1;
1288 : : }
1289 [ - + ]: 2 : if (PyModule_AddObject(mod, "nanj",
1290 : : PyComplex_FromCComplex(c_nanj())) < 0) {
1291 : 0 : 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 : 2 : 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 : 2 : 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 : 2 : 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 : 2 : 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 : 2 : 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 : 2 : 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 : 2 : 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 : 2 : 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 : 2 : 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 : 2 : 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 : 2 : 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 : 2 : return 0;
1410 : : }
1411 : :
1412 : : static PyModuleDef_Slot cmath_slots[] = {
1413 : : {Py_mod_exec, cmath_exec},
1414 : : {0, NULL}
1415 : : };
1416 : :
1417 : : static 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 : :
1426 : : PyMODINIT_FUNC
1427 : 2 : PyInit_cmath(void)
1428 : : {
1429 : 2 : return PyModuleDef_Init(&cmathmodule);
1430 : : }
|
