1 /* Complex math module */ 2 3 /* much code borrowed from mathmodule.c */ 4 5 #include "Python.h" 6 #include "_math.h" 7 /* we need DBL_MAX, DBL_MIN, DBL_EPSILON, DBL_MANT_DIG and FLT_RADIX from 8 float.h. We assume that FLT_RADIX is either 2 or 16. */ 9 #include <float.h> 10 11 #if (FLT_RADIX != 2 && FLT_RADIX != 16) 12 #error "Modules/cmathmodule.c expects FLT_RADIX to be 2 or 16" 13 #endif 14 15 #ifndef M_LN2 16 #define M_LN2 (0.6931471805599453094) /* natural log of 2 */ 17 #endif 18 19 #ifndef M_LN10 20 #define M_LN10 (2.302585092994045684) /* natural log of 10 */ 21 #endif 22 23 /* 24 CM_LARGE_DOUBLE is used to avoid spurious overflow in the sqrt, log, 25 inverse trig and inverse hyperbolic trig functions. Its log is used in the 26 evaluation of exp, cos, cosh, sin, sinh, tan, and tanh to avoid unnecessary 27 overflow. 28 */ 29 30 #define CM_LARGE_DOUBLE (DBL_MAX/4.) 31 #define CM_SQRT_LARGE_DOUBLE (sqrt(CM_LARGE_DOUBLE)) 32 #define CM_LOG_LARGE_DOUBLE (log(CM_LARGE_DOUBLE)) 33 #define CM_SQRT_DBL_MIN (sqrt(DBL_MIN)) 34 35 /* 36 CM_SCALE_UP is an odd integer chosen such that multiplication by 37 2**CM_SCALE_UP is sufficient to turn a subnormal into a normal. 38 CM_SCALE_DOWN is (-(CM_SCALE_UP+1)/2). These scalings are used to compute 39 square roots accurately when the real and imaginary parts of the argument 40 are subnormal. 41 */ 42 43 #if FLT_RADIX==2 44 #define CM_SCALE_UP (2*(DBL_MANT_DIG/2) + 1) 45 #elif FLT_RADIX==16 46 #define CM_SCALE_UP (4*DBL_MANT_DIG+1) 47 #endif 48 #define CM_SCALE_DOWN (-(CM_SCALE_UP+1)/2) 49 50 /* forward declarations */ 51 static Py_complex c_asinh(Py_complex); 52 static Py_complex c_atanh(Py_complex); 53 static Py_complex c_cosh(Py_complex); 54 static Py_complex c_sinh(Py_complex); 55 static Py_complex c_sqrt(Py_complex); 56 static Py_complex c_tanh(Py_complex); 57 static PyObject * math_error(void); 58 59 /* Code to deal with special values (infinities, NaNs, etc.). */ 60 61 /* special_type takes a double and returns an integer code indicating 62 the type of the double as follows: 63 */ 64 65 enum special_types { 66 ST_NINF, /* 0, negative infinity */ 67 ST_NEG, /* 1, negative finite number (nonzero) */ 68 ST_NZERO, /* 2, -0. */ 69 ST_PZERO, /* 3, +0. */ 70 ST_POS, /* 4, positive finite number (nonzero) */ 71 ST_PINF, /* 5, positive infinity */ 72 ST_NAN /* 6, Not a Number */ 73 }; 74 75 static enum special_types 76 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 } 1211
