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