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

Last change on this file was 2, checked in by Yuri Dario, 15 years ago

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