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