source: python/trunk/Modules/cmathmodule.c

Last change on this file was 391, checked in by dmik, 11 years ago

python: Merge vendor 2.7.6 to trunk.
  • Property svn:eol-style set to native
File size: 36.5 KB
RevLine 
[2]1/* Complex math module */
2
3/* much code borrowed from mathmodule.c */
4
5#include "Python.h"
[391]6#include "_math.h"
[2]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
[391]26 evaluation of exp, cos, cosh, sin, sinh, tan, and tanh to avoid unnecessary
[2]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
[391]35/*
[2]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 */
51static Py_complex c_asinh(Py_complex);
52static Py_complex c_atanh(Py_complex);
53static Py_complex c_cosh(Py_complex);
54static Py_complex c_sinh(Py_complex);
55static Py_complex c_sqrt(Py_complex);
56static Py_complex c_tanh(Py_complex);
57static 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
65enum special_types {
[391]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 */
[2]73};
74
75static enum special_types
76special_type(double d)
77{
[391]78 if (Py_IS_FINITE(d)) {
79 if (d != 0) {
80 if (copysign(1., d) == 1.)
81 return ST_POS;
82 else
83 return ST_NEG;
84 }
85 else {
86 if (copysign(1., d) == 1.)
87 return ST_PZERO;
88 else
89 return ST_NZERO;
90 }
91 }
92 if (Py_IS_NAN(d))
93 return ST_NAN;
94 if (copysign(1., d) == 1.)
95 return ST_PINF;
96 else
97 return ST_NINF;
[2]98}
99
[391]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 }
[2]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
124static Py_complex acos_special_values[7][7];
125
126static Py_complex
127c_acos(Py_complex z)
128{
[391]129 Py_complex s1, s2, r;
[2]130
[391]131 SPECIAL_VALUE(z, acos_special_values);
[2]132
[391]133 if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
134 /* avoid unnecessary overflow for large arguments */
135 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 if (z.real < 0.) {
139 r.imag = -copysign(log(hypot(z.real/2., z.imag/2.)) +
140 M_LN2*2., z.imag);
141 } else {
142 r.imag = copysign(log(hypot(z.real/2., z.imag/2.)) +
143 M_LN2*2., -z.imag);
144 }
145 } else {
146 s1.real = 1.-z.real;
147 s1.imag = -z.imag;
148 s1 = c_sqrt(s1);
149 s2.real = 1.+z.real;
150 s2.imag = z.imag;
151 s2 = c_sqrt(s2);
152 r.real = 2.*atan2(s1.real, s2.real);
153 r.imag = m_asinh(s2.real*s1.imag - s2.imag*s1.real);
154 }
155 errno = 0;
156 return r;
[2]157}
158
159PyDoc_STRVAR(c_acos_doc,
160"acos(x)\n"
161"\n"
162"Return the arc cosine of x.");
163
164
165static Py_complex acosh_special_values[7][7];
166
167static Py_complex
168c_acosh(Py_complex z)
169{
[391]170 Py_complex s1, s2, r;
[2]171
[391]172 SPECIAL_VALUE(z, acosh_special_values);
[2]173
[391]174 if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
175 /* avoid unnecessary overflow for large arguments */
176 r.real = log(hypot(z.real/2., z.imag/2.)) + M_LN2*2.;
177 r.imag = atan2(z.imag, z.real);
178 } else {
179 s1.real = z.real - 1.;
180 s1.imag = z.imag;
181 s1 = c_sqrt(s1);
182 s2.real = z.real + 1.;
183 s2.imag = z.imag;
184 s2 = c_sqrt(s2);
185 r.real = m_asinh(s1.real*s2.real + s1.imag*s2.imag);
186 r.imag = 2.*atan2(s1.imag, s2.real);
187 }
188 errno = 0;
189 return r;
[2]190}
191
192PyDoc_STRVAR(c_acosh_doc,
193"acosh(x)\n"
194"\n"
195"Return the hyperbolic arccosine of x.");
196
197
198static Py_complex
199c_asin(Py_complex z)
200{
[391]201 /* asin(z) = -i asinh(iz) */
202 Py_complex s, r;
203 s.real = -z.imag;
204 s.imag = z.real;
205 s = c_asinh(s);
206 r.real = s.imag;
207 r.imag = -s.real;
208 return r;
[2]209}
210
211PyDoc_STRVAR(c_asin_doc,
212"asin(x)\n"
213"\n"
214"Return the arc sine of x.");
215
216
217static Py_complex asinh_special_values[7][7];
218
219static Py_complex
220c_asinh(Py_complex z)
221{
[391]222 Py_complex s1, s2, r;
[2]223
[391]224 SPECIAL_VALUE(z, asinh_special_values);
[2]225
[391]226 if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
227 if (z.imag >= 0.) {
228 r.real = copysign(log(hypot(z.real/2., z.imag/2.)) +
229 M_LN2*2., z.real);
230 } else {
231 r.real = -copysign(log(hypot(z.real/2., z.imag/2.)) +
232 M_LN2*2., -z.real);
233 }
234 r.imag = atan2(z.imag, fabs(z.real));
235 } else {
236 s1.real = 1.+z.imag;
237 s1.imag = -z.real;
238 s1 = c_sqrt(s1);
239 s2.real = 1.-z.imag;
240 s2.imag = z.real;
241 s2 = c_sqrt(s2);
242 r.real = m_asinh(s1.real*s2.imag-s2.real*s1.imag);
243 r.imag = atan2(z.imag, s1.real*s2.real-s1.imag*s2.imag);
244 }
245 errno = 0;
246 return r;
[2]247}
248
249PyDoc_STRVAR(c_asinh_doc,
250"asinh(x)\n"
251"\n"
252"Return the hyperbolic arc sine of x.");
253
254
255static Py_complex
256c_atan(Py_complex z)
257{
[391]258 /* atan(z) = -i atanh(iz) */
259 Py_complex s, r;
260 s.real = -z.imag;
261 s.imag = z.real;
262 s = c_atanh(s);
263 r.real = s.imag;
264 r.imag = -s.real;
265 return r;
[2]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.). */
270static double
271c_atan2(Py_complex z)
272{
[391]273 if (Py_IS_NAN(z.real) || Py_IS_NAN(z.imag))
274 return Py_NAN;
275 if (Py_IS_INFINITY(z.imag)) {
276 if (Py_IS_INFINITY(z.real)) {
277 if (copysign(1., z.real) == 1.)
278 /* atan2(+-inf, +inf) == +-pi/4 */
279 return copysign(0.25*Py_MATH_PI, z.imag);
280 else
281 /* atan2(+-inf, -inf) == +-pi*3/4 */
282 return copysign(0.75*Py_MATH_PI, z.imag);
283 }
284 /* atan2(+-inf, x) == +-pi/2 for finite x */
285 return copysign(0.5*Py_MATH_PI, z.imag);
286 }
287 if (Py_IS_INFINITY(z.real) || z.imag == 0.) {
288 if (copysign(1., z.real) == 1.)
289 /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
290 return copysign(0., z.imag);
291 else
292 /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
293 return copysign(Py_MATH_PI, z.imag);
294 }
295 return atan2(z.imag, z.real);
[2]296}
297
298PyDoc_STRVAR(c_atan_doc,
299"atan(x)\n"
300"\n"
301"Return the arc tangent of x.");
302
303
304static Py_complex atanh_special_values[7][7];
305
306static Py_complex
307c_atanh(Py_complex z)
308{
[391]309 Py_complex r;
310 double ay, h;
[2]311
[391]312 SPECIAL_VALUE(z, atanh_special_values);
[2]313
[391]314 /* Reduce to case where z.real >= 0., using atanh(z) = -atanh(-z). */
315 if (z.real < 0.) {
316 return c_neg(c_atanh(c_neg(z)));
317 }
[2]318
[391]319 ay = fabs(z.imag);
320 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 h = hypot(z.real/2., z.imag/2.); /* safe from overflow */
327 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 r.imag = -copysign(Py_MATH_PI/2., -z.imag);
333 errno = 0;
334 } else if (z.real == 1. && ay < CM_SQRT_DBL_MIN) {
335 /* C99 standard says: atanh(1+/-0.) should be inf +/- 0i */
336 if (ay == 0.) {
337 r.real = INF;
338 r.imag = z.imag;
339 errno = EDOM;
340 } else {
341 r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.)));
342 r.imag = copysign(atan2(2., -ay)/2, z.imag);
343 errno = 0;
344 }
345 } else {
346 r.real = m_log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.;
347 r.imag = -atan2(-2.*z.imag, (1-z.real)*(1+z.real) - ay*ay)/2.;
348 errno = 0;
349 }
350 return r;
[2]351}
352
353PyDoc_STRVAR(c_atanh_doc,
354"atanh(x)\n"
355"\n"
356"Return the hyperbolic arc tangent of x.");
357
358
359static Py_complex
360c_cos(Py_complex z)
361{
[391]362 /* cos(z) = cosh(iz) */
363 Py_complex r;
364 r.real = -z.imag;
365 r.imag = z.real;
366 r = c_cosh(r);
367 return r;
[2]368}
369
370PyDoc_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 */
377static Py_complex cosh_special_values[7][7];
378
379static Py_complex
380c_cosh(Py_complex z)
381{
[391]382 Py_complex r;
383 double x_minus_one;
[2]384
[391]385 /* special treatment for cosh(+/-inf + iy) if y is not a NaN */
386 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
387 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) &&
388 (z.imag != 0.)) {
389 if (z.real > 0) {
390 r.real = copysign(INF, cos(z.imag));
391 r.imag = copysign(INF, sin(z.imag));
392 }
393 else {
394 r.real = copysign(INF, cos(z.imag));
395 r.imag = -copysign(INF, sin(z.imag));
396 }
397 }
398 else {
399 r = cosh_special_values[special_type(z.real)]
400 [special_type(z.imag)];
401 }
402 /* need to set errno = EDOM if y is +/- infinity and x is not
403 a NaN */
404 if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
405 errno = EDOM;
406 else
407 errno = 0;
408 return r;
409 }
[2]410
[391]411 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 x_minus_one = z.real - copysign(1., z.real);
415 r.real = cos(z.imag) * cosh(x_minus_one) * Py_MATH_E;
416 r.imag = sin(z.imag) * sinh(x_minus_one) * Py_MATH_E;
417 } else {
418 r.real = cos(z.imag) * cosh(z.real);
419 r.imag = sin(z.imag) * sinh(z.real);
420 }
421 /* detect overflow, and set errno accordingly */
422 if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
423 errno = ERANGE;
424 else
425 errno = 0;
426 return r;
[2]427}
428
429PyDoc_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 */
437static Py_complex exp_special_values[7][7];
438
439static Py_complex
440c_exp(Py_complex z)
441{
[391]442 Py_complex r;
443 double l;
[2]444
[391]445 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
446 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
447 && (z.imag != 0.)) {
448 if (z.real > 0) {
449 r.real = copysign(INF, cos(z.imag));
450 r.imag = copysign(INF, sin(z.imag));
451 }
452 else {
453 r.real = copysign(0., cos(z.imag));
454 r.imag = copysign(0., sin(z.imag));
455 }
456 }
457 else {
458 r = exp_special_values[special_type(z.real)]
459 [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 if (Py_IS_INFINITY(z.imag) &&
464 (Py_IS_FINITE(z.real) ||
465 (Py_IS_INFINITY(z.real) && z.real > 0)))
466 errno = EDOM;
467 else
468 errno = 0;
469 return r;
470 }
[2]471
[391]472 if (z.real > CM_LOG_LARGE_DOUBLE) {
473 l = exp(z.real-1.);
474 r.real = l*cos(z.imag)*Py_MATH_E;
475 r.imag = l*sin(z.imag)*Py_MATH_E;
476 } else {
477 l = exp(z.real);
478 r.real = l*cos(z.imag);
479 r.imag = l*sin(z.imag);
480 }
481 /* detect overflow, and set errno accordingly */
482 if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
483 errno = ERANGE;
484 else
485 errno = 0;
486 return r;
[2]487}
488
489PyDoc_STRVAR(c_exp_doc,
490"exp(x)\n"
491"\n"
492"Return the exponential value e**x.");
493
494
495static Py_complex log_special_values[7][7];
496
497static Py_complex
498c_log(Py_complex z)
499{
[391]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:
[2]504
[391]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.
[2]510
[391]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.
[2]514
[391]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.
[2]519
[391]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.
[2]525
[391]526 */
[2]527
[391]528 Py_complex r;
529 double ax, ay, am, an, h;
[2]530
[391]531 SPECIAL_VALUE(z, log_special_values);
[2]532
[391]533 ax = fabs(z.real);
534 ay = fabs(z.imag);
[2]535
[391]536 if (ax > CM_LARGE_DOUBLE || ay > CM_LARGE_DOUBLE) {
537 r.real = log(hypot(ax/2., ay/2.)) + M_LN2;
538 } else if (ax < DBL_MIN && ay < DBL_MIN) {
539 if (ax > 0. || ay > 0.) {
540 /* catch cases where hypot(ax, ay) is subnormal */
541 r.real = log(hypot(ldexp(ax, DBL_MANT_DIG),
542 ldexp(ay, DBL_MANT_DIG))) - DBL_MANT_DIG*M_LN2;
543 }
544 else {
545 /* log(+/-0. +/- 0i) */
546 r.real = -INF;
547 r.imag = atan2(z.imag, z.real);
548 errno = EDOM;
549 return r;
550 }
551 } else {
552 h = hypot(ax, ay);
553 if (0.71 <= h && h <= 1.73) {
554 am = ax > ay ? ax : ay; /* max(ax, ay) */
555 an = ax > ay ? ay : ax; /* min(ax, ay) */
556 r.real = m_log1p((am-1)*(am+1)+an*an)/2.;
557 } else {
558 r.real = log(h);
559 }
560 }
561 r.imag = atan2(z.imag, z.real);
562 errno = 0;
563 return r;
[2]564}
565
566
567static Py_complex
568c_log10(Py_complex z)
569{
[391]570 Py_complex r;
571 int errno_save;
[2]572
[391]573 r = c_log(z);
574 errno_save = errno; /* just in case the divisions affect errno */
575 r.real = r.real / M_LN10;
576 r.imag = r.imag / M_LN10;
577 errno = errno_save;
578 return r;
[2]579}
580
581PyDoc_STRVAR(c_log10_doc,
582"log10(x)\n"
583"\n"
584"Return the base-10 logarithm of x.");
585
586
587static Py_complex
588c_sin(Py_complex z)
589{
[391]590 /* sin(z) = -i sin(iz) */
591 Py_complex s, r;
592 s.real = -z.imag;
593 s.imag = z.real;
594 s = c_sinh(s);
595 r.real = s.imag;
596 r.imag = -s.real;
597 return r;
[2]598}
599
600PyDoc_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 */
607static Py_complex sinh_special_values[7][7];
608
609static Py_complex
610c_sinh(Py_complex z)
611{
[391]612 Py_complex r;
613 double x_minus_one;
[2]614
[391]615 /* special treatment for sinh(+/-inf + iy) if y is finite and
616 nonzero */
617 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
618 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
619 && (z.imag != 0.)) {
620 if (z.real > 0) {
621 r.real = copysign(INF, cos(z.imag));
622 r.imag = copysign(INF, sin(z.imag));
623 }
624 else {
625 r.real = -copysign(INF, cos(z.imag));
626 r.imag = copysign(INF, sin(z.imag));
627 }
628 }
629 else {
630 r = sinh_special_values[special_type(z.real)]
631 [special_type(z.imag)];
632 }
633 /* need to set errno = EDOM if y is +/- infinity and x is not
634 a NaN */
635 if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
636 errno = EDOM;
637 else
638 errno = 0;
639 return r;
640 }
[2]641
[391]642 if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
643 x_minus_one = z.real - copysign(1., z.real);
644 r.real = cos(z.imag) * sinh(x_minus_one) * Py_MATH_E;
645 r.imag = sin(z.imag) * cosh(x_minus_one) * Py_MATH_E;
646 } else {
647 r.real = cos(z.imag) * sinh(z.real);
648 r.imag = sin(z.imag) * cosh(z.real);
649 }
650 /* detect overflow, and set errno accordingly */
651 if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
652 errno = ERANGE;
653 else
654 errno = 0;
655 return r;
[2]656}
657
658PyDoc_STRVAR(c_sinh_doc,
659"sinh(x)\n"
660"\n"
661"Return the hyperbolic sine of x.");
662
663
664static Py_complex sqrt_special_values[7][7];
665
666static Py_complex
667c_sqrt(Py_complex z)
668{
[391]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
[2]673
[391]674 s = sqrt((x + hypot(x, y))/2)
[2]675
[391]676 and the imaginary part is
[2]677
[391]678 d = (y/2)/s
[2]679
[391]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:
[2]683
[391]684 s = 2*sqrt(x/8 + hypot(x/8, y/8))
[2]685
[391]686 This costs us two extra multiplications/divisions, but avoids the
687 overhead of checking for x and y large.
[2]688
[391]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 */
[2]694
695
[391]696 Py_complex r;
697 double s,d;
698 double ax, ay;
[2]699
[391]700 SPECIAL_VALUE(z, sqrt_special_values);
[2]701
[391]702 if (z.real == 0. && z.imag == 0.) {
703 r.real = 0.;
704 r.imag = z.imag;
705 return r;
706 }
[2]707
[391]708 ax = fabs(z.real);
709 ay = fabs(z.imag);
[2]710
[391]711 if (ax < DBL_MIN && ay < DBL_MIN && (ax > 0. || ay > 0.)) {
712 /* here we catch cases where hypot(ax, ay) is subnormal */
713 ax = ldexp(ax, CM_SCALE_UP);
714 s = ldexp(sqrt(ax + hypot(ax, ldexp(ay, CM_SCALE_UP))),
715 CM_SCALE_DOWN);
716 } else {
717 ax /= 8.;
718 s = 2.*sqrt(ax + hypot(ax, ay/8.));
719 }
720 d = ay/(2.*s);
[2]721
[391]722 if (z.real >= 0.) {
723 r.real = s;
724 r.imag = copysign(d, z.imag);
725 } else {
726 r.real = d;
727 r.imag = copysign(s, z.imag);
728 }
729 errno = 0;
730 return r;
[2]731}
732
733PyDoc_STRVAR(c_sqrt_doc,
734"sqrt(x)\n"
735"\n"
736"Return the square root of x.");
737
738
739static Py_complex
740c_tan(Py_complex z)
741{
[391]742 /* tan(z) = -i tanh(iz) */
743 Py_complex s, r;
744 s.real = -z.imag;
745 s.imag = z.real;
746 s = c_tanh(s);
747 r.real = s.imag;
748 r.imag = -s.real;
749 return r;
[2]750}
751
752PyDoc_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 */
759static Py_complex tanh_special_values[7][7];
760
761static Py_complex
762c_tanh(Py_complex z)
763{
[391]764 /* Formula:
[2]765
[391]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)
[2]768
[391]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).
[2]773
[391]774 */
[2]775
[391]776 Py_complex r;
777 double tx, ty, cx, txty, denom;
[2]778
[391]779 /* special treatment for tanh(+/-inf + iy) if y is finite and
780 nonzero */
781 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
782 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
783 && (z.imag != 0.)) {
784 if (z.real > 0) {
785 r.real = 1.0;
786 r.imag = copysign(0.,
787 2.*sin(z.imag)*cos(z.imag));
788 }
789 else {
790 r.real = -1.0;
791 r.imag = copysign(0.,
792 2.*sin(z.imag)*cos(z.imag));
793 }
794 }
795 else {
796 r = tanh_special_values[special_type(z.real)]
797 [special_type(z.imag)];
798 }
799 /* need to set errno = EDOM if z.imag is +/-infinity and
800 z.real is finite */
801 if (Py_IS_INFINITY(z.imag) && Py_IS_FINITE(z.real))
802 errno = EDOM;
803 else
804 errno = 0;
805 return r;
806 }
[2]807
[391]808 /* danger of overflow in 2.*z.imag !*/
809 if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
810 r.real = copysign(1., z.real);
811 r.imag = 4.*sin(z.imag)*cos(z.imag)*exp(-2.*fabs(z.real));
812 } else {
813 tx = tanh(z.real);
814 ty = tan(z.imag);
815 cx = 1./cosh(z.real);
816 txty = tx*ty;
817 denom = 1. + txty*txty;
818 r.real = tx*(1.+ty*ty)/denom;
819 r.imag = ((ty/denom)*cx)*cx;
820 }
821 errno = 0;
822 return r;
[2]823}
824
825PyDoc_STRVAR(c_tanh_doc,
826"tanh(x)\n"
827"\n"
828"Return the hyperbolic tangent of x.");
829
830
831static PyObject *
832cmath_log(PyObject *self, PyObject *args)
833{
[391]834 Py_complex x;
835 Py_complex y;
[2]836
[391]837 if (!PyArg_ParseTuple(args, "D|D", &x, &y))
838 return NULL;
[2]839
[391]840 errno = 0;
841 PyFPE_START_PROTECT("complex function", return 0)
842 x = c_log(x);
843 if (PyTuple_GET_SIZE(args) == 2) {
844 y = c_log(y);
845 x = c_quot(x, y);
846 }
847 PyFPE_END_PROTECT(x)
848 if (errno != 0)
849 return math_error();
850 return PyComplex_FromCComplex(x);
[2]851}
852
853PyDoc_STRVAR(cmath_log_doc,
854"log(x[, base]) -> the logarithm of x to the given base.\n\
855If 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
860static PyObject *
861math_error(void)
862{
[391]863 if (errno == EDOM)
864 PyErr_SetString(PyExc_ValueError, "math domain error");
865 else if (errno == ERANGE)
866 PyErr_SetString(PyExc_OverflowError, "math range error");
867 else /* Unexpected math error */
868 PyErr_SetFromErrno(PyExc_ValueError);
869 return NULL;
[2]870}
871
872static PyObject *
873math_1(PyObject *args, Py_complex (*func)(Py_complex))
874{
[391]875 Py_complex x,r ;
876 if (!PyArg_ParseTuple(args, "D", &x))
877 return NULL;
878 errno = 0;
879 PyFPE_START_PROTECT("complex function", return 0);
880 r = (*func)(x);
881 PyFPE_END_PROTECT(r);
882 if (errno == EDOM) {
883 PyErr_SetString(PyExc_ValueError, "math domain error");
884 return NULL;
885 }
886 else if (errno == ERANGE) {
887 PyErr_SetString(PyExc_OverflowError, "math range error");
888 return NULL;
889 }
890 else {
891 return PyComplex_FromCComplex(r);
892 }
[2]893}
894
895#define FUNC1(stubname, func) \
[391]896 static PyObject * stubname(PyObject *self, PyObject *args) { \
897 return math_1(args, func); \
898 }
[2]899
900FUNC1(cmath_acos, c_acos)
901FUNC1(cmath_acosh, c_acosh)
902FUNC1(cmath_asin, c_asin)
903FUNC1(cmath_asinh, c_asinh)
904FUNC1(cmath_atan, c_atan)
905FUNC1(cmath_atanh, c_atanh)
906FUNC1(cmath_cos, c_cos)
907FUNC1(cmath_cosh, c_cosh)
908FUNC1(cmath_exp, c_exp)
909FUNC1(cmath_log10, c_log10)
910FUNC1(cmath_sin, c_sin)
911FUNC1(cmath_sinh, c_sinh)
912FUNC1(cmath_sqrt, c_sqrt)
913FUNC1(cmath_tan, c_tan)
914FUNC1(cmath_tanh, c_tanh)
915
916static PyObject *
917cmath_phase(PyObject *self, PyObject *args)
918{
[391]919 Py_complex z;
920 double phi;
921 if (!PyArg_ParseTuple(args, "D:phase", &z))
922 return NULL;
923 errno = 0;
924 PyFPE_START_PROTECT("arg function", return 0)
925 phi = c_atan2(z);
926 PyFPE_END_PROTECT(phi)
927 if (errno != 0)
928 return math_error();
929 else
930 return PyFloat_FromDouble(phi);
[2]931}
932
933PyDoc_STRVAR(cmath_phase_doc,
934"phase(z) -> float\n\n\
935Return argument, also known as the phase angle, of a complex.");
936
937static PyObject *
938cmath_polar(PyObject *self, PyObject *args)
939{
[391]940 Py_complex z;
941 double r, phi;
942 if (!PyArg_ParseTuple(args, "D:polar", &z))
943 return NULL;
944 PyFPE_START_PROTECT("polar function", return 0)
945 phi = c_atan2(z); /* should not cause any exception */
946 r = c_abs(z); /* sets errno to ERANGE on overflow; otherwise 0 */
947 PyFPE_END_PROTECT(r)
948 if (errno != 0)
949 return math_error();
950 else
951 return Py_BuildValue("dd", r, phi);
[2]952}
953
954PyDoc_STRVAR(cmath_polar_doc,
955"polar(z) -> r: float, phi: float\n\n\
956Convert a complex from rectangular coordinates to polar coordinates. r is\n\
957the 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
970static Py_complex rect_special_values[7][7];
971
972static PyObject *
973cmath_rect(PyObject *self, PyObject *args)
974{
[391]975 Py_complex z;
976 double r, phi;
977 if (!PyArg_ParseTuple(args, "dd:rect", &r, &phi))
978 return NULL;
979 errno = 0;
980 PyFPE_START_PROTECT("rect function", return 0)
[2]981
[391]982 /* deal with special values */
983 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 if (Py_IS_INFINITY(r) && (Py_IS_FINITE(phi)
988 && (phi != 0.))) {
989 if (r > 0) {
990 z.real = copysign(INF, cos(phi));
991 z.imag = copysign(INF, sin(phi));
992 }
993 else {
994 z.real = -copysign(INF, cos(phi));
995 z.imag = -copysign(INF, sin(phi));
996 }
997 }
998 else {
999 z = rect_special_values[special_type(r)]
1000 [special_type(phi)];
1001 }
1002 /* need to set errno = EDOM if r is a nonzero number and phi
1003 is infinite */
1004 if (r != 0. && !Py_IS_NAN(r) && Py_IS_INFINITY(phi))
1005 errno = EDOM;
1006 else
1007 errno = 0;
1008 }
1009 else if (phi == 0.0) {
1010 /* Workaround for buggy results with phi=-0.0 on OS X 10.8. See
1011 bugs.python.org/issue18513. */
1012 z.real = r;
1013 z.imag = r * phi;
1014 errno = 0;
1015 }
1016 else {
1017 z.real = r * cos(phi);
1018 z.imag = r * sin(phi);
1019 errno = 0;
1020 }
[2]1021
[391]1022 PyFPE_END_PROTECT(z)
1023 if (errno != 0)
1024 return math_error();
1025 else
1026 return PyComplex_FromCComplex(z);
[2]1027}
1028
1029PyDoc_STRVAR(cmath_rect_doc,
1030"rect(r, phi) -> z: complex\n\n\
1031Convert from polar coordinates to rectangular coordinates.");
1032
1033static PyObject *
1034cmath_isnan(PyObject *self, PyObject *args)
1035{
[391]1036 Py_complex z;
1037 if (!PyArg_ParseTuple(args, "D:isnan", &z))
1038 return NULL;
1039 return PyBool_FromLong(Py_IS_NAN(z.real) || Py_IS_NAN(z.imag));
[2]1040}
1041
1042PyDoc_STRVAR(cmath_isnan_doc,
1043"isnan(z) -> bool\n\
1044Checks if the real or imaginary part of z not a number (NaN)");
1045
1046static PyObject *
1047cmath_isinf(PyObject *self, PyObject *args)
1048{
[391]1049 Py_complex z;
1050 if (!PyArg_ParseTuple(args, "D:isnan", &z))
1051 return NULL;
1052 return PyBool_FromLong(Py_IS_INFINITY(z.real) ||
1053 Py_IS_INFINITY(z.imag));
[2]1054}
1055
1056PyDoc_STRVAR(cmath_isinf_doc,
1057"isinf(z) -> bool\n\
1058Checks if the real or imaginary part of z is infinite.");
1059
1060
1061PyDoc_STRVAR(module_doc,
1062"This module is always available. It provides access to mathematical\n"
1063"functions for complex numbers.");
1064
1065static PyMethodDef cmath_methods[] = {
[391]1066 {"acos", cmath_acos, METH_VARARGS, c_acos_doc},
1067 {"acosh", cmath_acosh, METH_VARARGS, c_acosh_doc},
1068 {"asin", cmath_asin, METH_VARARGS, c_asin_doc},
1069 {"asinh", cmath_asinh, METH_VARARGS, c_asinh_doc},
1070 {"atan", cmath_atan, METH_VARARGS, c_atan_doc},
1071 {"atanh", cmath_atanh, METH_VARARGS, c_atanh_doc},
1072 {"cos", cmath_cos, METH_VARARGS, c_cos_doc},
1073 {"cosh", cmath_cosh, METH_VARARGS, c_cosh_doc},
1074 {"exp", cmath_exp, METH_VARARGS, c_exp_doc},
1075 {"isinf", cmath_isinf, METH_VARARGS, cmath_isinf_doc},
1076 {"isnan", cmath_isnan, METH_VARARGS, cmath_isnan_doc},
1077 {"log", cmath_log, METH_VARARGS, cmath_log_doc},
1078 {"log10", cmath_log10, METH_VARARGS, c_log10_doc},
1079 {"phase", cmath_phase, METH_VARARGS, cmath_phase_doc},
1080 {"polar", cmath_polar, METH_VARARGS, cmath_polar_doc},
1081 {"rect", cmath_rect, METH_VARARGS, cmath_rect_doc},
1082 {"sin", cmath_sin, METH_VARARGS, c_sin_doc},
1083 {"sinh", cmath_sinh, METH_VARARGS, c_sinh_doc},
1084 {"sqrt", cmath_sqrt, METH_VARARGS, c_sqrt_doc},
1085 {"tan", cmath_tan, METH_VARARGS, c_tan_doc},
1086 {"tanh", cmath_tanh, METH_VARARGS, c_tanh_doc},
1087 {NULL, NULL} /* sentinel */
[2]1088};
1089
1090PyMODINIT_FUNC
1091initcmath(void)
1092{
[391]1093 PyObject *m;
[2]1094
[391]1095 m = Py_InitModule3("cmath", cmath_methods, module_doc);
1096 if (m == NULL)
1097 return;
[2]1098
[391]1099 PyModule_AddObject(m, "pi",
1100 PyFloat_FromDouble(Py_MATH_PI));
1101 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
[2]1102
[391]1103 /* initialize special value tables */
[2]1104
1105#define INIT_SPECIAL_VALUES(NAME, BODY) { Py_complex* p = (Py_complex*)NAME; BODY }
1106#define C(REAL, IMAG) p->real = REAL; p->imag = IMAG; ++p;
1107
[391]1108 INIT_SPECIAL_VALUES(acos_special_values, {
1109 C(P34,INF) C(P,INF) C(P,INF) C(P,-INF) C(P,-INF) C(P34,-INF) C(N,INF)
1110 C(P12,INF) C(U,U) C(U,U) C(U,U) C(U,U) C(P12,-INF) C(N,N)
1111 C(P12,INF) C(U,U) C(P12,0.) C(P12,-0.) C(U,U) C(P12,-INF) C(P12,N)
1112 C(P12,INF) C(U,U) C(P12,0.) C(P12,-0.) C(U,U) C(P12,-INF) C(P12,N)
1113 C(P12,INF) C(U,U) C(U,U) C(U,U) C(U,U) C(P12,-INF) C(N,N)
1114 C(P14,INF) C(0.,INF) C(0.,INF) C(0.,-INF) C(0.,-INF) C(P14,-INF) C(N,INF)
1115 C(N,INF) C(N,N) C(N,N) C(N,N) C(N,N) C(N,-INF) C(N,N)
1116 })
[2]1117
[391]1118 INIT_SPECIAL_VALUES(acosh_special_values, {
1119 C(INF,-P34) C(INF,-P) C(INF,-P) C(INF,P) C(INF,P) C(INF,P34) C(INF,N)
1120 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1121 C(INF,-P12) C(U,U) C(0.,-P12) C(0.,P12) C(U,U) C(INF,P12) C(N,N)
1122 C(INF,-P12) C(U,U) C(0.,-P12) C(0.,P12) C(U,U) C(INF,P12) C(N,N)
1123 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1124 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
1125 C(INF,N) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,N) C(N,N)
1126 })
[2]1127
[391]1128 INIT_SPECIAL_VALUES(asinh_special_values, {
1129 C(-INF,-P14) C(-INF,-0.) C(-INF,-0.) C(-INF,0.) C(-INF,0.) C(-INF,P14) C(-INF,N)
1130 C(-INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(-INF,P12) C(N,N)
1131 C(-INF,-P12) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(-INF,P12) C(N,N)
1132 C(INF,-P12) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(INF,P12) C(N,N)
1133 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1134 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
1135 C(INF,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(INF,N) C(N,N)
1136 })
[2]1137
[391]1138 INIT_SPECIAL_VALUES(atanh_special_values, {
1139 C(-0.,-P12) C(-0.,-P12) C(-0.,-P12) C(-0.,P12) C(-0.,P12) C(-0.,P12) C(-0.,N)
1140 C(-0.,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(-0.,P12) C(N,N)
1141 C(-0.,-P12) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(-0.,P12) C(-0.,N)
1142 C(0.,-P12) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,P12) C(0.,N)
1143 C(0.,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(0.,P12) C(N,N)
1144 C(0.,-P12) C(0.,-P12) C(0.,-P12) C(0.,P12) C(0.,P12) C(0.,P12) C(0.,N)
1145 C(0.,-P12) C(N,N) C(N,N) C(N,N) C(N,N) C(0.,P12) C(N,N)
1146 })
[2]1147
[391]1148 INIT_SPECIAL_VALUES(cosh_special_values, {
1149 C(INF,N) C(U,U) C(INF,0.) C(INF,-0.) C(U,U) C(INF,N) C(INF,N)
1150 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1151 C(N,0.) C(U,U) C(1.,0.) C(1.,-0.) C(U,U) C(N,0.) C(N,0.)
1152 C(N,0.) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(N,0.) C(N,0.)
1153 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1154 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
1155 C(N,N) C(N,N) C(N,0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1156 })
[2]1157
[391]1158 INIT_SPECIAL_VALUES(exp_special_values, {
1159 C(0.,0.) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,0.) C(0.,0.)
1160 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1161 C(N,N) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(N,N) C(N,N)
1162 C(N,N) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(N,N) C(N,N)
1163 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1164 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
1165 C(N,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1166 })
[2]1167
[391]1168 INIT_SPECIAL_VALUES(log_special_values, {
1169 C(INF,-P34) C(INF,-P) C(INF,-P) C(INF,P) C(INF,P) C(INF,P34) C(INF,N)
1170 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1171 C(INF,-P12) C(U,U) C(-INF,-P) C(-INF,P) C(U,U) C(INF,P12) C(N,N)
1172 C(INF,-P12) C(U,U) C(-INF,-0.) C(-INF,0.) C(U,U) C(INF,P12) C(N,N)
1173 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1174 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
1175 C(INF,N) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,N) C(N,N)
1176 })
[2]1177
[391]1178 INIT_SPECIAL_VALUES(sinh_special_values, {
1179 C(INF,N) C(U,U) C(-INF,-0.) C(-INF,0.) C(U,U) C(INF,N) C(INF,N)
1180 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1181 C(0.,N) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(0.,N) C(0.,N)
1182 C(0.,N) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,N) C(0.,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 })
[2]1187
[391]1188 INIT_SPECIAL_VALUES(sqrt_special_values, {
1189 C(INF,-INF) C(0.,-INF) C(0.,-INF) C(0.,INF) C(0.,INF) C(INF,INF) C(N,INF)
1190 C(INF,-INF) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,INF) C(N,N)
1191 C(INF,-INF) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(INF,INF) C(N,N)
1192 C(INF,-INF) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(INF,INF) C(N,N)
1193 C(INF,-INF) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,INF) C(N,N)
1194 C(INF,-INF) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,INF) C(INF,N)
1195 C(INF,-INF) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,INF) C(N,N)
1196 })
[2]1197
[391]1198 INIT_SPECIAL_VALUES(tanh_special_values, {
1199 C(-1.,0.) C(U,U) C(-1.,-0.) C(-1.,0.) C(U,U) C(-1.,0.) C(-1.,0.)
1200 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1201 C(N,N) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(N,N) C(N,N)
1202 C(N,N) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(N,N) C(N,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(1.,0.) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(1.,0.) C(1.,0.)
1205 C(N,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1206 })
[2]1207
[391]1208 INIT_SPECIAL_VALUES(rect_special_values, {
1209 C(INF,N) C(U,U) C(-INF,0.) C(-INF,-0.) C(U,U) C(INF,N) C(INF,N)
1210 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1211 C(0.,0.) C(U,U) C(-0.,0.) C(-0.,-0.) C(U,U) C(0.,0.) C(0.,0.)
1212 C(0.,0.) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,0.) C(0.,0.)
1213 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1214 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
1215 C(N,N) C(N,N) C(N,0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1216 })
[2]1217}
Note: See TracBrowser for help on using the repository browser.