?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
move sections
[python/dscho.git] / Modules / cmathmodule.c
blob825f7ad4dfc9db81cff7e47ad6e43814ac8ada20
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 unecessary
27 overflow.
28 */
29
30 #define CM_LARGE_DOUBLE (DBL_MAX/4.)
31 #define CM_SQRT_LARGE_DOUBLE (sqrt(CM_LARGE_DOUBLE))
32 #define CM_LOG_LARGE_DOUBLE (log(CM_LARGE_DOUBLE))
33 #define CM_SQRT_DBL_MIN (sqrt(DBL_MIN))
34
35 /*
36 CM_SCALE_UP is an odd integer chosen such that multiplication by
37 2**CM_SCALE_UP is sufficient to turn a subnormal into a normal.
38 CM_SCALE_DOWN is (-(CM_SCALE_UP+1)/2). These scalings are used to compute
39 square roots accurately when the real and imaginary parts of the argument
40 are subnormal.
41 */
42
43 #if FLT_RADIX==2
44 #define CM_SCALE_UP (2*(DBL_MANT_DIG/2) + 1)
45 #elif FLT_RADIX==16
46 #define CM_SCALE_UP (4*DBL_MANT_DIG+1)
47 #endif
48 #define CM_SCALE_DOWN (-(CM_SCALE_UP+1)/2)
49
50 /* forward declarations */
51 static Py_complex c_asinh(Py_complex);
52 static Py_complex c_atanh(Py_complex);
53 static Py_complex c_cosh(Py_complex);
54 static Py_complex c_sinh(Py_complex);
55 static Py_complex c_sqrt(Py_complex);
56 static Py_complex c_tanh(Py_complex);
57 static PyObject * math_error(void);
58
59 /* Code to deal with special values (infinities, NaNs, etc.). */
60
61 /* special_type takes a double and returns an integer code indicating
62 the type of the double as follows:
63 */
64
65 enum special_types {
66 ST_NINF, /* 0, negative infinity */
67 ST_NEG, /* 1, negative finite number (nonzero) */
68 ST_NZERO, /* 2, -0. */
69 ST_PZERO, /* 3, +0. */
70 ST_POS, /* 4, positive finite number (nonzero) */
71 ST_PINF, /* 5, positive infinity */
72 ST_NAN /* 6, Not a Number */
73 };
74
75 static enum special_types
76 special_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
124 static Py_complex acos_special_values[7][7];
125
126 static Py_complex
127 c_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
159 PyDoc_STRVAR(c_acos_doc,
160 "acos(x)\n"
161 "\n"
162 "Return the arc cosine of x.");
163
164
165 static Py_complex acosh_special_values[7][7];
166
167 static Py_complex
168 c_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
192 PyDoc_STRVAR(c_acosh_doc,
193 "acosh(x)\n"
194 "\n"
195 "Return the hyperbolic arccosine of x.");
196
197
198 static Py_complex
199 c_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
211 PyDoc_STRVAR(c_asin_doc,
212 "asin(x)\n"
213 "\n"
214 "Return the arc sine of x.");
215
216
217 static Py_complex asinh_special_values[7][7];
218
219 static Py_complex
220 c_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
249 PyDoc_STRVAR(c_asinh_doc,
250 "asinh(x)\n"
251 "\n"
252 "Return the hyperbolic arc sine of x.");
253
254
255 static Py_complex
256 c_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.). */
270 static double
271 c_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
298 PyDoc_STRVAR(c_atan_doc,
299 "atan(x)\n"
300 "\n"
301 "Return the arc tangent of x.");
302
303
304 static Py_complex atanh_special_values[7][7];
305
306 static Py_complex
307 c_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
353 PyDoc_STRVAR(c_atanh_doc,
354 "atanh(x)\n"
355 "\n"
356 "Return the hyperbolic arc tangent of x.");
357
358
359 static Py_complex
360 c_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
370 PyDoc_STRVAR(c_cos_doc,
371 "cos(x)\n"
372 "\n"
373 "Return the cosine of x.");
374
375
376 /* cosh(infinity + i*y) needs to be dealt with specially */
377 static Py_complex cosh_special_values[7][7];
378
379 static Py_complex
380 c_cosh(Py_complex z)
381 {
382 Py_complex r;
383 double x_minus_one;
384
385 /* special treatment for cosh(+/-inf + iy) if y is not a NaN */
386 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
429 PyDoc_STRVAR(c_cosh_doc,
430 "cosh(x)\n"
431 "\n"
432 "Return the hyperbolic cosine of x.");
433
434
435 /* exp(infinity + i*y) and exp(-infinity + i*y) need special treatment for
436 finite y */
437 static Py_complex exp_special_values[7][7];
438
439 static Py_complex
440 c_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
489 PyDoc_STRVAR(c_exp_doc,
490 "exp(x)\n"
491 "\n"
492 "Return the exponential value e**x.");
493
494
495 static Py_complex log_special_values[7][7];
496
497 static Py_complex
498 c_log(Py_complex z)
499 {
500 /*
501 The usual formula for the real part is log(hypot(z.real, z.imag)).
502 There are four situations where this formula is potentially
503 problematic:
504
505 (1) the absolute value of z is subnormal. Then hypot is subnormal,
506 so has fewer than the usual number of bits of accuracy, hence may
507 have large relative error. This then gives a large absolute error
508 in the log. This can be solved by rescaling z by a suitable power
509 of 2.
510
511 (2) the absolute value of z is greater than DBL_MAX (e.g. when both
512 z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX)
513 Again, rescaling solves this.
514
515 (3) the absolute value of z is close to 1. In this case it's
516 difficult to achieve good accuracy, at least in part because a
517 change of 1ulp in the real or imaginary part of z can result in a
518 change of billions of ulps in the correctly rounded answer.
519
520 (4) z = 0. The simplest thing to do here is to call the
521 floating-point log with an argument of 0, and let its behaviour
522 (returning -infinity, signaling a floating-point exception, setting
523 errno, or whatever) determine that of c_log. So the usual formula
524 is fine here.
525
526 */
527
528 Py_complex r;
529 double ax, ay, am, an, h;
530
531 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
567 static Py_complex
568 c_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
581 PyDoc_STRVAR(c_log10_doc,
582 "log10(x)\n"
583 "\n"
584 "Return the base-10 logarithm of x.");
585
586
587 static Py_complex
588 c_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
600 PyDoc_STRVAR(c_sin_doc,
601 "sin(x)\n"
602 "\n"
603 "Return the sine of x.");
604
605
606 /* sinh(infinity + i*y) needs to be dealt with specially */
607 static Py_complex sinh_special_values[7][7];
608
609 static Py_complex
610 c_sinh(Py_complex z)
611 {
612 Py_complex r;
613 double x_minus_one;
614
615 /* special treatment for sinh(+/-inf + iy) if y is finite and
616 nonzero */
617 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
658 PyDoc_STRVAR(c_sinh_doc,
659 "sinh(x)\n"
660 "\n"
661 "Return the hyperbolic sine of x.");
662
663
664 static Py_complex sqrt_special_values[7][7];
665
666 static Py_complex
667 c_sqrt(Py_complex z)
668 {
669 /*
670 Method: use symmetries to reduce to the case when x = z.real and y
671 = z.imag are nonnegative. Then the real part of the result is
672 given by
673
674 s = sqrt((x + hypot(x, y))/2)
675
676 and the imaginary part is
677
678 d = (y/2)/s
679
680 If either x or y is very large then there's a risk of overflow in
681 computation of the expression x + hypot(x, y). We can avoid this
682 by rewriting the formula for s as:
683
684 s = 2*sqrt(x/8 + hypot(x/8, y/8))
685
686 This costs us two extra multiplications/divisions, but avoids the
687 overhead of checking for x and y large.
688
689 If both x and y are subnormal then hypot(x, y) may also be
690 subnormal, so will lack full precision. We solve this by rescaling
691 x and y by a sufficiently large power of 2 to ensure that x and y
692 are normal.
693 */
694
695
696 Py_complex r;
697 double s,d;
698 double ax, ay;
699
700 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
733 PyDoc_STRVAR(c_sqrt_doc,
734 "sqrt(x)\n"
735 "\n"
736 "Return the square root of x.");
737
738
739 static Py_complex
740 c_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
752 PyDoc_STRVAR(c_tan_doc,
753 "tan(x)\n"
754 "\n"
755 "Return the tangent of x.");
756
757
758 /* tanh(infinity + i*y) needs to be dealt with specially */
759 static Py_complex tanh_special_values[7][7];
760
761 static Py_complex
762 c_tanh(Py_complex z)
763 {
764 /* Formula:
765
766 tanh(x+iy) = (tanh(x)(1+tan(y)^2) + i tan(y)(1-tanh(x))^2) /
767 (1+tan(y)^2 tanh(x)^2)
768
769 To avoid excessive roundoff error, 1-tanh(x)^2 is better computed
770 as 1/cosh(x)^2. When abs(x) is large, we approximate 1-tanh(x)^2
771 by 4 exp(-2*x) instead, to avoid possible overflow in the
772 computation of cosh(x).
773
774 */
775
776 Py_complex r;
777 double tx, ty, cx, txty, denom;
778
779 /* special treatment for tanh(+/-inf + iy) if y is finite and
780 nonzero */
781 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
825 PyDoc_STRVAR(c_tanh_doc,
826 "tanh(x)\n"
827 "\n"
828 "Return the hyperbolic tangent of x.");
829
830
831 static PyObject *
832 cmath_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
853 PyDoc_STRVAR(cmath_log_doc,
854 "log(x[, base]) -> the logarithm of x to the given base.\n\
855 If the base not specified, returns the natural logarithm (base e) of x.");
856
857
858 /* And now the glue to make them available from Python: */
859
860 static PyObject *
861 math_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
872 static PyObject *
873 math_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
900 FUNC1(cmath_acos, c_acos)
901 FUNC1(cmath_acosh, c_acosh)
902 FUNC1(cmath_asin, c_asin)
903 FUNC1(cmath_asinh, c_asinh)
904 FUNC1(cmath_atan, c_atan)
905 FUNC1(cmath_atanh, c_atanh)
906 FUNC1(cmath_cos, c_cos)
907 FUNC1(cmath_cosh, c_cosh)
908 FUNC1(cmath_exp, c_exp)
909 FUNC1(cmath_log10, c_log10)
910 FUNC1(cmath_sin, c_sin)
911 FUNC1(cmath_sinh, c_sinh)
912 FUNC1(cmath_sqrt, c_sqrt)
913 FUNC1(cmath_tan, c_tan)
914 FUNC1(cmath_tanh, c_tanh)
915
916 static PyObject *
917 cmath_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
933 PyDoc_STRVAR(cmath_phase_doc,
934 "phase(z) -> float\n\n\
935 Return argument, also known as the phase angle, of a complex.");
936
937 static PyObject *
938 cmath_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
954 PyDoc_STRVAR(cmath_polar_doc,
955 "polar(z) -> r: float, phi: float\n\n\
956 Convert a complex from rectangular coordinates to polar coordinates. r is\n\
957 the distance from 0 and phi the phase angle.");
958
959 /*
960 rect() isn't covered by the C99 standard, but it's not too hard to
961 figure out 'spirit of C99' rules for special value handing:
962
963 rect(x, t) should behave like exp(log(x) + it) for positive-signed x
964 rect(x, t) should behave like -exp(log(-x) + it) for negative-signed x
965 rect(nan, t) should behave like exp(nan + it), except that rect(nan, 0)
966 gives nan +- i0 with the sign of the imaginary part unspecified.
967
968 */
969
970 static Py_complex rect_special_values[7][7];
971
972 static PyObject *
973 cmath_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 {
1010 z.real = r * cos(phi);
1011 z.imag = r * sin(phi);
1012 errno = 0;
1013 }
1014
1015 PyFPE_END_PROTECT(z)
1016 if (errno != 0)
1017 return math_error();
1018 else
1019 return PyComplex_FromCComplex(z);
1020 }
1021
1022 PyDoc_STRVAR(cmath_rect_doc,
1023 "rect(r, phi) -> z: complex\n\n\
1024 Convert from polar coordinates to rectangular coordinates.");
1025
1026 static PyObject *
1027 cmath_isnan(PyObject *self, PyObject *args)
1028 {
1029 Py_complex z;
1030 if (!PyArg_ParseTuple(args, "D:isnan", &z))
1031 return NULL;
1032 return PyBool_FromLong(Py_IS_NAN(z.real) || Py_IS_NAN(z.imag));
1033 }
1034
1035 PyDoc_STRVAR(cmath_isnan_doc,
1036 "isnan(z) -> bool\n\
1037 Checks if the real or imaginary part of z not a number (NaN)");
1038
1039 static PyObject *
1040 cmath_isinf(PyObject *self, PyObject *args)
1041 {
1042 Py_complex z;
1043 if (!PyArg_ParseTuple(args, "D:isnan", &z))
1044 return NULL;
1045 return PyBool_FromLong(Py_IS_INFINITY(z.real) ||
1046 Py_IS_INFINITY(z.imag));
1047 }
1048
1049 PyDoc_STRVAR(cmath_isinf_doc,
1050 "isinf(z) -> bool\n\
1051 Checks if the real or imaginary part of z is infinite.");
1052
1053
1054 PyDoc_STRVAR(module_doc,
1055 "This module is always available. It provides access to mathematical\n"
1056 "functions for complex numbers.");
1057
1058 static PyMethodDef cmath_methods[] = {
1059 {"acos", cmath_acos, METH_VARARGS, c_acos_doc},
1060 {"acosh", cmath_acosh, METH_VARARGS, c_acosh_doc},
1061 {"asin", cmath_asin, METH_VARARGS, c_asin_doc},
1062 {"asinh", cmath_asinh, METH_VARARGS, c_asinh_doc},
1063 {"atan", cmath_atan, METH_VARARGS, c_atan_doc},
1064 {"atanh", cmath_atanh, METH_VARARGS, c_atanh_doc},
1065 {"cos", cmath_cos, METH_VARARGS, c_cos_doc},
1066 {"cosh", cmath_cosh, METH_VARARGS, c_cosh_doc},
1067 {"exp", cmath_exp, METH_VARARGS, c_exp_doc},
1068 {"isinf", cmath_isinf, METH_VARARGS, cmath_isinf_doc},
1069 {"isnan", cmath_isnan, METH_VARARGS, cmath_isnan_doc},
1070 {"log", cmath_log, METH_VARARGS, cmath_log_doc},
1071 {"log10", cmath_log10, METH_VARARGS, c_log10_doc},
1072 {"phase", cmath_phase, METH_VARARGS, cmath_phase_doc},
1073 {"polar", cmath_polar, METH_VARARGS, cmath_polar_doc},
1074 {"rect", cmath_rect, METH_VARARGS, cmath_rect_doc},
1075 {"sin", cmath_sin, METH_VARARGS, c_sin_doc},
1076 {"sinh", cmath_sinh, METH_VARARGS, c_sinh_doc},
1077 {"sqrt", cmath_sqrt, METH_VARARGS, c_sqrt_doc},
1078 {"tan", cmath_tan, METH_VARARGS, c_tan_doc},
1079 {"tanh", cmath_tanh, METH_VARARGS, c_tanh_doc},
1080 {NULL, NULL} /* sentinel */
1081 };
1082
1083 PyMODINIT_FUNC
1084 initcmath(void)
1085 {
1086 PyObject *m;
1087
1088 m = Py_InitModule3("cmath", cmath_methods, module_doc);
1089 if (m == NULL)
1090 return;
1091
1092 PyModule_AddObject(m, "pi",
1093 PyFloat_FromDouble(Py_MATH_PI));
1094 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
1095
1096 /* initialize special value tables */
1097
1098 #define INIT_SPECIAL_VALUES(NAME, BODY) { Py_complex* p = (Py_complex*)NAME; BODY }
1099 #define C(REAL, IMAG) p->real = REAL; p->imag = IMAG; ++p;
1100
1101 INIT_SPECIAL_VALUES(acos_special_values, {
1102 C(P34,INF) C(P,INF) C(P,INF) C(P,-INF) C(P,-INF) C(P34,-INF) C(N,INF)
1103 C(P12,INF) C(U,U) C(U,U) C(U,U) C(U,U) C(P12,-INF) C(N,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(P12,0.) C(P12,-0.) C(U,U) C(P12,-INF) C(P12,N)
1106 C(P12,INF) C(U,U) C(U,U) C(U,U) C(U,U) C(P12,-INF) C(N,N)
1107 C(P14,INF) C(0.,INF) C(0.,INF) C(0.,-INF) C(0.,-INF) C(P14,-INF) C(N,INF)
1108 C(N,INF) C(N,N) C(N,N) C(N,N) C(N,N) C(N,-INF) C(N,N)
1109 })
1110
1111 INIT_SPECIAL_VALUES(acosh_special_values, {
1112 C(INF,-P34) C(INF,-P) C(INF,-P) C(INF,P) C(INF,P) C(INF,P34) C(INF,N)
1113 C(INF,-P12) C(U,U) C(U,U) C(U,U) 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(0.,-P12) C(0.,P12) C(U,U) C(INF,P12) C(N,N)
1116 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1117 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
1118 C(INF,N) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,N) C(N,N)
1119 })
1120
1121 INIT_SPECIAL_VALUES(asinh_special_values, {
1122 C(-INF,-P14) C(-INF,-0.) C(-INF,-0.) C(-INF,0.) C(-INF,0.) C(-INF,P14) C(-INF,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,-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(0.,-0.) C(0.,0.) C(U,U) C(INF,P12) C(N,N)
1126 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1127 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
1128 C(INF,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(INF,N) C(N,N)
1129 })
1130
1131 INIT_SPECIAL_VALUES(atanh_special_values, {
1132 C(-0.,-P12) C(-0.,-P12) C(-0.,-P12) C(-0.,P12) C(-0.,P12) C(-0.,P12) C(-0.,N)
1133 C(-0.,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(-0.,P12) C(N,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(0.,-0.) C(0.,0.) C(U,U) C(0.,P12) C(0.,N)
1136 C(0.,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(0.,P12) C(N,N)
1137 C(0.,-P12) C(0.,-P12) C(0.,-P12) C(0.,P12) C(0.,P12) C(0.,P12) C(0.,N)
1138 C(0.,-P12) C(N,N) C(N,N) C(N,N) C(N,N) C(0.,P12) C(N,N)
1139 })
1140
1141 INIT_SPECIAL_VALUES(cosh_special_values, {
1142 C(INF,N) C(U,U) C(INF,0.) C(INF,-0.) C(U,U) C(INF,N) C(INF,N)
1143 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
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,0.) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(N,0.) C(N,0.)
1146 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1147 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
1148 C(N,N) C(N,N) C(N,0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1149 })
1150
1151 INIT_SPECIAL_VALUES(exp_special_values, {
1152 C(0.,0.) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,0.) C(0.,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(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(1.,-0.) C(1.,0.) C(U,U) C(N,N) C(N,N)
1156 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1157 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
1158 C(N,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1159 })
1160
1161 INIT_SPECIAL_VALUES(log_special_values, {
1162 C(INF,-P34) C(INF,-P) C(INF,-P) C(INF,P) C(INF,P) C(INF,P34) C(INF,N)
1163 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1164 C(INF,-P12) C(U,U) C(-INF,-P) C(-INF,P) C(U,U) C(INF,P12) C(N,N)
1165 C(INF,-P12) C(U,U) C(-INF,-0.) C(-INF,0.) C(U,U) C(INF,P12) C(N,N)
1166 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N)
1167 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
1168 C(INF,N) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,N) C(N,N)
1169 })
1170
1171 INIT_SPECIAL_VALUES(sinh_special_values, {
1172 C(INF,N) C(U,U) C(-INF,-0.) C(-INF,0.) C(U,U) C(INF,N) C(INF,N)
1173 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,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(0.,N) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,N) C(0.,N)
1176 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1177 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
1178 C(N,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1179 })
1180
1181 INIT_SPECIAL_VALUES(sqrt_special_values, {
1182 C(INF,-INF) C(0.,-INF) C(0.,-INF) C(0.,INF) C(0.,INF) C(INF,INF) C(N,INF)
1183 C(INF,-INF) C(U,U) C(U,U) C(U,U) 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(0.,-0.) C(0.,0.) C(U,U) C(INF,INF) C(N,N)
1186 C(INF,-INF) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,INF) C(N,N)
1187 C(INF,-INF) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,INF) C(INF,N)
1188 C(INF,-INF) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,INF) C(N,N)
1189 })
1190
1191 INIT_SPECIAL_VALUES(tanh_special_values, {
1192 C(-1.,0.) C(U,U) C(-1.,-0.) C(-1.,0.) C(U,U) C(-1.,0.) C(-1.,0.)
1193 C(N,N) C(U,U) C(U,U) C(U,U) 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(0.,-0.) C(0.,0.) C(U,U) C(N,N) C(N,N)
1196 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1197 C(1.,0.) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(1.,0.) C(1.,0.)
1198 C(N,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1199 })
1200
1201 INIT_SPECIAL_VALUES(rect_special_values, {
1202 C(INF,N) C(U,U) C(-INF,0.) C(-INF,-0.) C(U,U) C(INF,N) C(INF,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(0.,0.) C(U,U) C(-0.,0.) C(-0.,-0.) C(U,U) C(0.,0.) C(0.,0.)
1205 C(0.,0.) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,0.) C(0.,0.)
1206 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N)
1207 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
1208 C(N,N) C(N,N) C(N,0.) C(N,0.) C(N,N) C(N,N) C(N,N)
1209 })
1210 }