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