ยปCore Development>Code coverage>Modules/cmathmodule.c

Python code coverage for Modules/cmathmodule.c

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