source: python/vendor/Python-2.7.6/Modules/cmathmodule.c

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

python: Update vendor to 2.7.6.
  • Property svn:eol-style set to native
File size: 36.5 KB
Line 
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 */
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 {
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
75static enum special_types
76special_type(double d)
77{
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;
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
124static Py_complex acos_special_values[7][7];
125
126static Py_complex
127c_acos(Py_complex z)
128{
129 Py_complex s1, s2, r;
130
131 SPECIAL_VALUE(z, acos_special_values);
132
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;
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{
170 Py_complex s1, s2, r;
171
172 SPECIAL_VALUE(z, acosh_special_values);
173
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;
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{
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;
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{
222 Py_complex s1, s2, r;
223
224 SPECIAL_VALUE(z, asinh_special_values);
225
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;
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{
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;
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{
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);
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{
309 Py_complex r;
310 double ay, h;
311
312 SPECIAL_VALUE(z, atanh_special_values);
313
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 }
318
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;
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{
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;
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{
382 Py_complex r;
383 double x_minus_one;
384
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 }
410
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;
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{
442 Py_complex r;
443 double l;
444
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 }
471
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;
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{
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 SPECIAL_VALUE(z, log_special_values);
532
533 ax = fabs(z.real);
534 ay = fabs(z.imag);
535
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;
564}
565
566
567static Py_complex
568c_log10(Py_complex z)
569{
570 Py_complex r;
571 int errno_save;
572
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;
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{
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;
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{
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 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 }
641
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;
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{
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 SPECIAL_VALUE(z, sqrt_special_values);
701
702 if (z.real == 0. && z.imag == 0.) {
703 r.real = 0.;
704 r.imag = z.imag;
705 return r;
706 }
707
708 ax = fabs(z.real);
709 ay = fabs(z.imag);
710
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);
721
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;
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{
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;
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{
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 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 }
807
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;
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{
834 Py_complex x;
835 Py_complex y;
836
837 if (!PyArg_ParseTuple(args, "D|D", &x, &y))
838 return NULL;
839
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);
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{
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;
870}
871
872static PyObject *
873math_1(PyObject *args, Py_complex (*func)(Py_complex))
874{
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 }
893}
894
895#define FUNC1(stubname, func) \
896 static PyObject * stubname(PyObject *self, PyObject *args) { \
897 return math_1(args, func); \
898 }
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{
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);
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{
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);
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{
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)
981
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 }
1021
1022 PyFPE_END_PROTECT(z)
1023 if (errno != 0)
1024 return math_error();
1025 else
1026 return PyComplex_FromCComplex(z);
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{
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));
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{
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));
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[] = {
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 */
1088};
1089
1090PyMODINIT_FUNC
1091initcmath(void)
1092{
1093 PyObject *m;
1094
1095 m = Py_InitModule3("cmath", cmath_methods, module_doc);
1096 if (m == NULL)
1097 return;
1098
1099 PyModule_AddObject(m, "pi",
1100 PyFloat_FromDouble(Py_MATH_PI));
1101 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
1102
1103 /* initialize special value tables */
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
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 })
1117
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 })
1127
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 })
1137
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 })
1147
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 })
1157
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 })
1167
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 })
1177
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 })
1187
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 })
1197
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 })
1207
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 })
1217}
Note: See TracBrowser for help on using the repository browser.