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 #include "clinic/cmathmodule.c.h" 12 /*[clinic input] 13 module cmath 14 [clinic start generated code]*/ 15 /*[clinic end generated code: output=da39a3ee5e6b4b0d input=308d6839f4a46333]*/ 16 17 /*[python input] 18 class Py_complex_protected_converter(Py_complex_converter): 19 def modify(self): 20 return 'errno = 0; PyFPE_START_PROTECT("complex function", goto exit);' 21 22 23 class Py_complex_protected_return_converter(CReturnConverter): 24 type = "Py_complex" 25 26 def render(self, function, data): 27 self.declare(data) 28 data.return_conversion.append(""" 29 PyFPE_END_PROTECT(_return_value); 30 if (errno == EDOM) { 31 PyErr_SetString(PyExc_ValueError, "math domain error"); 32 goto exit; 33 } 34 else if (errno == ERANGE) { 35 PyErr_SetString(PyExc_OverflowError, "math range error"); 36 goto exit; 37 } 38 else { 39 return_value = PyComplex_FromCComplex(_return_value); 40 } 41 """.strip()) 42 [python start generated code]*/ 43 /*[python end generated code: output=da39a3ee5e6b4b0d input=345daa075b1028e7]*/ 44 45 #if (FLT_RADIX != 2 && FLT_RADIX != 16) 46 #error "Modules/cmathmodule.c expects FLT_RADIX to be 2 or 16" 47 #endif 48 49 #ifndef M_LN2 50 #define M_LN2 (0.6931471805599453094) /* natural log of 2 */ 51 #endif 52 53 #ifndef M_LN10 54 #define M_LN10 (2.302585092994045684) /* natural log of 10 */ 55 #endif 56 57 /* 58 CM_LARGE_DOUBLE is used to avoid spurious overflow in the sqrt, log, 59 inverse trig and inverse hyperbolic trig functions. Its log is used in the 60 evaluation of exp, cos, cosh, sin, sinh, tan, and tanh to avoid unnecessary 61 overflow. 62 */ 63 64 #define CM_LARGE_DOUBLE (DBL_MAX/4.) 65 #define CM_SQRT_LARGE_DOUBLE (sqrt(CM_LARGE_DOUBLE)) 66 #define CM_LOG_LARGE_DOUBLE (log(CM_LARGE_DOUBLE)) 67 #define CM_SQRT_DBL_MIN (sqrt(DBL_MIN)) 68 69 /* 70 CM_SCALE_UP is an odd integer chosen such that multiplication by 71 2**CM_SCALE_UP is sufficient to turn a subnormal into a normal. 72 CM_SCALE_DOWN is (-(CM_SCALE_UP+1)/2). These scalings are used to compute 73 square roots accurately when the real and imaginary parts of the argument 74 are subnormal. 75 */ 76 77 #if FLT_RADIX==2 78 #define CM_SCALE_UP (2*(DBL_MANT_DIG/2) + 1) 79 #elif FLT_RADIX==16 80 #define CM_SCALE_UP (4*DBL_MANT_DIG+1) 81 #endif 82 #define CM_SCALE_DOWN (-(CM_SCALE_UP+1)/2) 83 84 /* Constants cmath.inf, cmath.infj, cmath.nan, cmath.nanj. 85 cmath.nan and cmath.nanj are defined only when either 86 PY_NO_SHORT_FLOAT_REPR is *not* defined (which should be 87 the most common situation on machines using an IEEE 754 88 representation), or Py_NAN is defined. */ 89 90 static double 91 m_inf(void) 92 { 93 #ifndef PY_NO_SHORT_FLOAT_REPR 94 return _Py_dg_infinity(0); 95 #else 96 return Py_HUGE_VAL; 97 #endif 98 } 99 100 static Py_complex 101 c_infj(void) 102 { 103 Py_complex r; 104 r.real = 0.0; 105 r.imag = m_inf(); 106 return r; 107 } 108 109 #if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN) 110 111 static double 112 m_nan(void) 113 { 114 #ifndef PY_NO_SHORT_FLOAT_REPR 115 return _Py_dg_stdnan(0); 116 #else 117 return Py_NAN; 118 #endif 119 } 120 121 static Py_complex 122 c_nanj(void) 123 { 124 Py_complex r; 125 r.real = 0.0; 126 r.imag = m_nan(); 127 return r; 128 } 129 130 #endif 131 132 /* forward declarations */ 133 static Py_complex cmath_asinh_impl(PyObject *, Py_complex); 134 static Py_complex cmath_atanh_impl(PyObject *, Py_complex); 135 static Py_complex cmath_cosh_impl(PyObject *, Py_complex); 136 static Py_complex cmath_sinh_impl(PyObject *, Py_complex); 137 static Py_complex cmath_sqrt_impl(PyObject *, Py_complex); 138 static Py_complex cmath_tanh_impl(PyObject *, Py_complex); 139 static PyObject * math_error(void); 140 141 /* Code to deal with special values (infinities, NaNs, etc.). */ 142 143 /* special_type takes a double and returns an integer code indicating 144 the type of the double as follows: 145 */ 146 147 enum special_types { 148 ST_NINF, /* 0, negative infinity */ 149 ST_NEG, /* 1, negative finite number (nonzero) */ 150 ST_NZERO, /* 2, -0. */ 151 ST_PZERO, /* 3, +0. */ 152 ST_POS, /* 4, positive finite number (nonzero) */ 153 ST_PINF, /* 5, positive infinity */ 154 ST_NAN /* 6, Not a Number */ 155 }; 156 157 static enum special_types 158 special_type(double d) 159 { 160 if (Py_IS_FINITE(d)) { 161 if (d != 0) { 162 if (copysign(1., d) == 1.) 163 return ST_POS; 164 else 165 return ST_NEG; 166 } 167 else { 168 if (copysign(1., d) == 1.) 169 return ST_PZERO; 170 else 171 return ST_NZERO; 172 } 173 } 174 if (Py_IS_NAN(d)) 175 return ST_NAN; 176 if (copysign(1., d) == 1.) 177 return ST_PINF; 178 else 179 return ST_NINF; 180 } 181 182 #define SPECIAL_VALUE(z, table) \ 183 if (!Py_IS_FINITE((z).real) || !Py_IS_FINITE((z).imag)) { \ 184 errno = 0; \ 185 return table[special_type((z).real)] \ 186 [special_type((z).imag)]; \ 187 } 188 189 #define P Py_MATH_PI 190 #define P14 0.25*Py_MATH_PI 191 #define P12 0.5*Py_MATH_PI 192 #define P34 0.75*Py_MATH_PI 193 #define INF Py_HUGE_VAL 194 #define N Py_NAN 195 #define U -9.5426319407711027e33 /* unlikely value, used as placeholder */ 196 197 /* First, the C functions that do the real work. Each of the c_* 198 functions computes and returns the C99 Annex G recommended result 199 and also sets errno as follows: errno = 0 if no floating-point 200 exception is associated with the result; errno = EDOM if C99 Annex 201 G recommends raising divide-by-zero or invalid for this result; and 202 errno = ERANGE where the overflow floating-point signal should be 203 raised. 204 */ 205 206 static Py_complex acos_special_values[7][7]; 207 208 /*[clinic input] 209 cmath.acos -> Py_complex_protected 210 211 z: Py_complex_protected 212 / 213 214 Return the arc cosine of z. 215 [clinic start generated code]*/ 216 217 static Py_complex 218 cmath_acos_impl(PyObject *module, Py_complex z) 219 /*[clinic end generated code: output=40bd42853fd460ae input=bd6cbd78ae851927]*/ 220 { 221 Py_complex s1, s2, r; 222 223 SPECIAL_VALUE(z, acos_special_values); 224 225 if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) { 226 /* avoid unnecessary overflow for large arguments */ 227 r.real = atan2(fabs(z.imag), z.real); 228 /* split into cases to make sure that the branch cut has the 229 correct continuity on systems with unsigned zeros */ 230 if (z.real < 0.) { 231 r.imag = -copysign(log(hypot(z.real/2., z.imag/2.)) + 232 M_LN2*2., z.imag); 233 } else { 234 r.imag = copysign(log(hypot(z.real/2., z.imag/2.)) + 235 M_LN2*2., -z.imag); 236 } 237 } else { 238 s1.real = 1.-z.real; 239 s1.imag = -z.imag; 240 s1 = cmath_sqrt_impl(module, s1); 241 s2.real = 1.+z.real; 242 s2.imag = z.imag; 243 s2 = cmath_sqrt_impl(module, s2); 244 r.real = 2.*atan2(s1.real, s2.real); 245 r.imag = m_asinh(s2.real*s1.imag - s2.imag*s1.real); 246 } 247 errno = 0; 248 return r; 249 } 250 251 252 static Py_complex acosh_special_values[7][7]; 253 254 /*[clinic input] 255 cmath.acosh = cmath.acos 256 257 Return the inverse hyperbolic cosine of z. 258 [clinic start generated code]*/ 259 260 static Py_complex 261 cmath_acosh_impl(PyObject *module, Py_complex z) 262 /*[clinic end generated code: output=3e2454d4fcf404ca input=3f61bee7d703e53c]*/ 263 { 264 Py_complex s1, s2, r; 265 266 SPECIAL_VALUE(z, acosh_special_values); 267 268 if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) { 269 /* avoid unnecessary overflow for large arguments */ 270 r.real = log(hypot(z.real/2., z.imag/2.)) + M_LN2*2.; 271 r.imag = atan2(z.imag, z.real); 272 } else { 273 s1.real = z.real - 1.; 274 s1.imag = z.imag; 275 s1 = cmath_sqrt_impl(module, s1); 276 s2.real = z.real + 1.; 277 s2.imag = z.imag; 278 s2 = cmath_sqrt_impl(module, s2); 279 r.real = m_asinh(s1.real*s2.real + s1.imag*s2.imag); 280 r.imag = 2.*atan2(s1.imag, s2.real); 281 } 282 errno = 0; 283 return r; 284 } 285 286 /*[clinic input] 287 cmath.asin = cmath.acos 288 289 Return the arc sine of z. 290 [clinic start generated code]*/ 291 292 static Py_complex 293 cmath_asin_impl(PyObject *module, Py_complex z) 294 /*[clinic end generated code: output=3b264cd1b16bf4e1 input=be0bf0cfdd5239c5]*/ 295 { 296 /* asin(z) = -i asinh(iz) */ 297 Py_complex s, r; 298 s.real = -z.imag; 299 s.imag = z.real; 300 s = cmath_asinh_impl(module, s); 301 r.real = s.imag; 302 r.imag = -s.real; 303 return r; 304 } 305 306 307 static Py_complex asinh_special_values[7][7]; 308 309 /*[clinic input] 310 cmath.asinh = cmath.acos 311 312 Return the inverse hyperbolic sine of z. 313 [clinic start generated code]*/ 314 315 static Py_complex 316 cmath_asinh_impl(PyObject *module, Py_complex z) 317 /*[clinic end generated code: output=733d8107841a7599 input=5c09448fcfc89a79]*/ 318 { 319 Py_complex s1, s2, r; 320 321 SPECIAL_VALUE(z, asinh_special_values); 322 323 if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) { 324 if (z.imag >= 0.) { 325 r.real = copysign(log(hypot(z.real/2., z.imag/2.)) + 326 M_LN2*2., z.real); 327 } else { 328 r.real = -copysign(log(hypot(z.real/2., z.imag/2.)) + 329 M_LN2*2., -z.real); 330 } 331 r.imag = atan2(z.imag, fabs(z.real)); 332 } else { 333 s1.real = 1.+z.imag; 334 s1.imag = -z.real; 335 s1 = cmath_sqrt_impl(module, s1); 336 s2.real = 1.-z.imag; 337 s2.imag = z.real; 338 s2 = cmath_sqrt_impl(module, s2); 339 r.real = m_asinh(s1.real*s2.imag-s2.real*s1.imag); 340 r.imag = atan2(z.imag, s1.real*s2.real-s1.imag*s2.imag); 341 } 342 errno = 0; 343 return r; 344 } 345 346 347 /*[clinic input] 348 cmath.atan = cmath.acos 349 350 Return the arc tangent of z. 351 [clinic start generated code]*/ 352 353 static Py_complex 354 cmath_atan_impl(PyObject *module, Py_complex z) 355 /*[clinic end generated code: output=b6bfc497058acba4 input=3b21ff7d5eac632a]*/ 356 { 357 /* atan(z) = -i atanh(iz) */ 358 Py_complex s, r; 359 s.real = -z.imag; 360 s.imag = z.real; 361 s = cmath_atanh_impl(module, s); 362 r.real = s.imag; 363 r.imag = -s.real; 364 return r; 365 } 366 367 /* Windows screws up atan2 for inf and nan, and alpha Tru64 5.1 doesn't follow 368 C99 for atan2(0., 0.). */ 369 static double 370 c_atan2(Py_complex z) 371 { 372 if (Py_IS_NAN(z.real) || Py_IS_NAN(z.imag)) 373 return Py_NAN; 374 if (Py_IS_INFINITY(z.imag)) { 375 if (Py_IS_INFINITY(z.real)) { 376 if (copysign(1., z.real) == 1.) 377 /* atan2(+-inf, +inf) == +-pi/4 */ 378 return copysign(0.25*Py_MATH_PI, z.imag); 379 else 380 /* atan2(+-inf, -inf) == +-pi*3/4 */ 381 return copysign(0.75*Py_MATH_PI, z.imag); 382 } 383 /* atan2(+-inf, x) == +-pi/2 for finite x */ 384 return copysign(0.5*Py_MATH_PI, z.imag); 385 } 386 if (Py_IS_INFINITY(z.real) || z.imag == 0.) { 387 if (copysign(1., z.real) == 1.) 388 /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */ 389 return copysign(0., z.imag); 390 else 391 /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */ 392 return copysign(Py_MATH_PI, z.imag); 393 } 394 return atan2(z.imag, z.real); 395 } 396 397 398 static Py_complex atanh_special_values[7][7]; 399 400 /*[clinic input] 401 cmath.atanh = cmath.acos 402 403 Return the inverse hyperbolic tangent of z. 404 [clinic start generated code]*/ 405 406 static Py_complex 407 cmath_atanh_impl(PyObject *module, Py_complex z) 408 /*[clinic end generated code: output=e83355f93a989c9e input=2b3fdb82fb34487b]*/ 409 { 410 Py_complex r; 411 double ay, h; 412 413 SPECIAL_VALUE(z, atanh_special_values); 414 415 /* Reduce to case where z.real >= 0., using atanh(z) = -atanh(-z). */ 416 if (z.real < 0.) { 417 return _Py_c_neg(cmath_atanh_impl(module, _Py_c_neg(z))); 418 } 419 420 ay = fabs(z.imag); 421 if (z.real > CM_SQRT_LARGE_DOUBLE || ay > CM_SQRT_LARGE_DOUBLE) { 422 /* 423 if abs(z) is large then we use the approximation 424 atanh(z) ~ 1/z +/- i*pi/2 (+/- depending on the sign 425 of z.imag) 426 */ 427 h = hypot(z.real/2., z.imag/2.); /* safe from overflow */ 428 r.real = z.real/4./h/h; 429 /* the two negations in the next line cancel each other out 430 except when working with unsigned zeros: they're there to 431 ensure that the branch cut has the correct continuity on 432 systems that don't support signed zeros */ 433 r.imag = -copysign(Py_MATH_PI/2., -z.imag); 434 errno = 0; 435 } else if (z.real == 1. && ay < CM_SQRT_DBL_MIN) { 436 /* C99 standard says: atanh(1+/-0.) should be inf +/- 0i */ 437 if (ay == 0.) { 438 r.real = INF; 439 r.imag = z.imag; 440 errno = EDOM; 441 } else { 442 r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.))); 443 r.imag = copysign(atan2(2., -ay)/2, z.imag); 444 errno = 0; 445 } 446 } else { 447 r.real = m_log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.; 448 r.imag = -atan2(-2.*z.imag, (1-z.real)*(1+z.real) - ay*ay)/2.; 449 errno = 0; 450 } 451 return r; 452 } 453 454 455 /*[clinic input] 456 cmath.cos = cmath.acos 457 458 Return the cosine of z. 459 [clinic start generated code]*/ 460 461 static Py_complex 462 cmath_cos_impl(PyObject *module, Py_complex z) 463 /*[clinic end generated code: output=fd64918d5b3186db input=6022e39b77127ac7]*/ 464 { 465 /* cos(z) = cosh(iz) */ 466 Py_complex r; 467 r.real = -z.imag; 468 r.imag = z.real; 469 r = cmath_cosh_impl(module, r); 470 return r; 471 } 472 473 474 /* cosh(infinity + i*y) needs to be dealt with specially */ 475 static Py_complex cosh_special_values[7][7]; 476 477 /*[clinic input] 478 cmath.cosh = cmath.acos 479 480 Return the hyperbolic cosine of z. 481 [clinic start generated code]*/ 482 483 static Py_complex 484 cmath_cosh_impl(PyObject *module, Py_complex z) 485 /*[clinic end generated code: output=2e969047da601bdb input=d6b66339e9cc332b]*/ 486 { 487 Py_complex r; 488 double x_minus_one; 489 490 /* special treatment for cosh(+/-inf + iy) if y is not a NaN */ 491 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) { 492 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) && 493 (z.imag != 0.)) { 494 if (z.real > 0) { 495 r.real = copysign(INF, cos(z.imag)); 496 r.imag = copysign(INF, sin(z.imag)); 497 } 498 else { 499 r.real = copysign(INF, cos(z.imag)); 500 r.imag = -copysign(INF, sin(z.imag)); 501 } 502 } 503 else { 504 r = cosh_special_values[special_type(z.real)] 505 [special_type(z.imag)]; 506 } 507 /* need to set errno = EDOM if y is +/- infinity and x is not 508 a NaN */ 509 if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real)) 510 errno = EDOM; 511 else 512 errno = 0; 513 return r; 514 } 515 516 if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) { 517 /* deal correctly with cases where cosh(z.real) overflows but 518 cosh(z) does not. */ 519 x_minus_one = z.real - copysign(1., z.real); 520 r.real = cos(z.imag) * cosh(x_minus_one) * Py_MATH_E; 521 r.imag = sin(z.imag) * sinh(x_minus_one) * Py_MATH_E; 522 } else { 523 r.real = cos(z.imag) * cosh(z.real); 524 r.imag = sin(z.imag) * sinh(z.real); 525 } 526 /* detect overflow, and set errno accordingly */ 527 if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag)) 528 errno = ERANGE; 529 else 530 errno = 0; 531 return r; 532 } 533 534 535 /* exp(infinity + i*y) and exp(-infinity + i*y) need special treatment for 536 finite y */ 537 static Py_complex exp_special_values[7][7]; 538 539 /*[clinic input] 540 cmath.exp = cmath.acos 541 542 Return the exponential value e**z. 543 [clinic start generated code]*/ 544 545 static Py_complex 546 cmath_exp_impl(PyObject *module, Py_complex z) 547 /*[clinic end generated code: output=edcec61fb9dfda6c input=8b9e6cf8a92174c3]*/ 548 { 549 Py_complex r; 550 double l; 551 552 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) { 553 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) 554 && (z.imag != 0.)) { 555 if (z.real > 0) { 556 r.real = copysign(INF, cos(z.imag)); 557 r.imag = copysign(INF, sin(z.imag)); 558 } 559 else { 560 r.real = copysign(0., cos(z.imag)); 561 r.imag = copysign(0., sin(z.imag)); 562 } 563 } 564 else { 565 r = exp_special_values[special_type(z.real)] 566 [special_type(z.imag)]; 567 } 568 /* need to set errno = EDOM if y is +/- infinity and x is not 569 a NaN and not -infinity */ 570 if (Py_IS_INFINITY(z.imag) && 571 (Py_IS_FINITE(z.real) || 572 (Py_IS_INFINITY(z.real) && z.real > 0))) 573 errno = EDOM; 574 else 575 errno = 0; 576 return r; 577 } 578 579 if (z.real > CM_LOG_LARGE_DOUBLE) { 580 l = exp(z.real-1.); 581 r.real = l*cos(z.imag)*Py_MATH_E; 582 r.imag = l*sin(z.imag)*Py_MATH_E; 583 } else { 584 l = exp(z.real); 585 r.real = l*cos(z.imag); 586 r.imag = l*sin(z.imag); 587 } 588 /* detect overflow, and set errno accordingly */ 589 if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag)) 590 errno = ERANGE; 591 else 592 errno = 0; 593 return r; 594 } 595 596 static Py_complex log_special_values[7][7]; 597 598 static Py_complex 599 c_log(Py_complex z) 600 { 601 /* 602 The usual formula for the real part is log(hypot(z.real, z.imag)). 603 There are four situations where this formula is potentially 604 problematic: 605 606 (1) the absolute value of z is subnormal. Then hypot is subnormal, 607 so has fewer than the usual number of bits of accuracy, hence may 608 have large relative error. This then gives a large absolute error 609 in the log. This can be solved by rescaling z by a suitable power 610 of 2. 611 612 (2) the absolute value of z is greater than DBL_MAX (e.g. when both 613 z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX) 614 Again, rescaling solves this. 615 616 (3) the absolute value of z is close to 1. In this case it's 617 difficult to achieve good accuracy, at least in part because a 618 change of 1ulp in the real or imaginary part of z can result in a 619 change of billions of ulps in the correctly rounded answer. 620 621 (4) z = 0. The simplest thing to do here is to call the 622 floating-point log with an argument of 0, and let its behaviour 623 (returning -infinity, signaling a floating-point exception, setting 624 errno, or whatever) determine that of c_log. So the usual formula 625 is fine here. 626 627 */ 628 629 Py_complex r; 630 double ax, ay, am, an, h; 631 632 SPECIAL_VALUE(z, log_special_values); 633 634 ax = fabs(z.real); 635 ay = fabs(z.imag); 636 637 if (ax > CM_LARGE_DOUBLE || ay > CM_LARGE_DOUBLE) { 638 r.real = log(hypot(ax/2., ay/2.)) + M_LN2; 639 } else if (ax < DBL_MIN && ay < DBL_MIN) { 640 if (ax > 0. || ay > 0.) { 641 /* catch cases where hypot(ax, ay) is subnormal */ 642 r.real = log(hypot(ldexp(ax, DBL_MANT_DIG), 643 ldexp(ay, DBL_MANT_DIG))) - DBL_MANT_DIG*M_LN2; 644 } 645 else { 646 /* log(+/-0. +/- 0i) */ 647 r.real = -INF; 648 r.imag = atan2(z.imag, z.real); 649 errno = EDOM; 650 return r; 651 } 652 } else { 653 h = hypot(ax, ay); 654 if (0.71 <= h && h <= 1.73) { 655 am = ax > ay ? ax : ay; /* max(ax, ay) */ 656 an = ax > ay ? ay : ax; /* min(ax, ay) */ 657 r.real = m_log1p((am-1)*(am+1)+an*an)/2.; 658 } else { 659 r.real = log(h); 660 } 661 } 662 r.imag = atan2(z.imag, z.real); 663 errno = 0; 664 return r; 665 } 666 667 668 /*[clinic input] 669 cmath.log10 = cmath.acos 670 671 Return the base-10 logarithm of z. 672 [clinic start generated code]*/ 673 674 static Py_complex 675 cmath_log10_impl(PyObject *module, Py_complex z) 676 /*[clinic end generated code: output=2922779a7c38cbe1 input=cff5644f73c1519c]*/ 677 { 678 Py_complex r; 679 int errno_save; 680 681 r = c_log(z); 682 errno_save = errno; /* just in case the divisions affect errno */ 683 r.real = r.real / M_LN10; 684 r.imag = r.imag / M_LN10; 685 errno = errno_save; 686 return r; 687 } 688 689 690 /*[clinic input] 691 cmath.sin = cmath.acos 692 693 Return the sine of z. 694 [clinic start generated code]*/ 695 696 static Py_complex 697 cmath_sin_impl(PyObject *module, Py_complex z) 698 /*[clinic end generated code: output=980370d2ff0bb5aa input=2d3519842a8b4b85]*/ 699 { 700 /* sin(z) = -i sin(iz) */ 701 Py_complex s, r; 702 s.real = -z.imag; 703 s.imag = z.real; 704 s = cmath_sinh_impl(module, s); 705 r.real = s.imag; 706 r.imag = -s.real; 707 return r; 708 } 709 710 711 /* sinh(infinity + i*y) needs to be dealt with specially */ 712 static Py_complex sinh_special_values[7][7]; 713 714 /*[clinic input] 715 cmath.sinh = cmath.acos 716 717 Return the hyperbolic sine of z. 718 [clinic start generated code]*/ 719 720 static Py_complex 721 cmath_sinh_impl(PyObject *module, Py_complex z) 722 /*[clinic end generated code: output=38b0a6cce26f3536 input=d2d3fc8c1ddfd2dd]*/ 723 { 724 Py_complex r; 725 double x_minus_one; 726 727 /* special treatment for sinh(+/-inf + iy) if y is finite and 728 nonzero */ 729 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) { 730 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) 731 && (z.imag != 0.)) { 732 if (z.real > 0) { 733 r.real = copysign(INF, cos(z.imag)); 734 r.imag = copysign(INF, sin(z.imag)); 735 } 736 else { 737 r.real = -copysign(INF, cos(z.imag)); 738 r.imag = copysign(INF, sin(z.imag)); 739 } 740 } 741 else { 742 r = sinh_special_values[special_type(z.real)] 743 [special_type(z.imag)]; 744 } 745 /* need to set errno = EDOM if y is +/- infinity and x is not 746 a NaN */ 747 if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real)) 748 errno = EDOM; 749 else 750 errno = 0; 751 return r; 752 } 753 754 if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) { 755 x_minus_one = z.real - copysign(1., z.real); 756 r.real = cos(z.imag) * sinh(x_minus_one) * Py_MATH_E; 757 r.imag = sin(z.imag) * cosh(x_minus_one) * Py_MATH_E; 758 } else { 759 r.real = cos(z.imag) * sinh(z.real); 760 r.imag = sin(z.imag) * cosh(z.real); 761 } 762 /* detect overflow, and set errno accordingly */ 763 if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag)) 764 errno = ERANGE; 765 else 766 errno = 0; 767 return r; 768 } 769 770 771 static Py_complex sqrt_special_values[7][7]; 772 773 /*[clinic input] 774 cmath.sqrt = cmath.acos 775 776 Return the square root of z. 777 [clinic start generated code]*/ 778 779 static Py_complex 780 cmath_sqrt_impl(PyObject *module, Py_complex z) 781 /*[clinic end generated code: output=b6507b3029c339fc input=7088b166fc9a58c7]*/ 782 { 783 /* 784 Method: use symmetries to reduce to the case when x = z.real and y 785 = z.imag are nonnegative. Then the real part of the result is 786 given by 787 788 s = sqrt((x + hypot(x, y))/2) 789 790 and the imaginary part is 791 792 d = (y/2)/s 793 794 If either x or y is very large then there's a risk of overflow in 795 computation of the expression x + hypot(x, y). We can avoid this 796 by rewriting the formula for s as: 797 798 s = 2*sqrt(x/8 + hypot(x/8, y/8)) 799 800 This costs us two extra multiplications/divisions, but avoids the 801 overhead of checking for x and y large. 802 803 If both x and y are subnormal then hypot(x, y) may also be 804 subnormal, so will lack full precision. We solve this by rescaling 805 x and y by a sufficiently large power of 2 to ensure that x and y 806 are normal. 807 */ 808 809 810 Py_complex r; 811 double s,d; 812 double ax, ay; 813 814 SPECIAL_VALUE(z, sqrt_special_values); 815 816 if (z.real == 0. && z.imag == 0.) { 817 r.real = 0.; 818 r.imag = z.imag; 819 return r; 820 } 821 822 ax = fabs(z.real); 823 ay = fabs(z.imag); 824 825 if (ax < DBL_MIN && ay < DBL_MIN && (ax > 0. || ay > 0.)) { 826 /* here we catch cases where hypot(ax, ay) is subnormal */ 827 ax = ldexp(ax, CM_SCALE_UP); 828 s = ldexp(sqrt(ax + hypot(ax, ldexp(ay, CM_SCALE_UP))), 829 CM_SCALE_DOWN); 830 } else { 831 ax /= 8.; 832 s = 2.*sqrt(ax + hypot(ax, ay/8.)); 833 } 834 d = ay/(2.*s); 835 836 if (z.real >= 0.) { 837 r.real = s; 838 r.imag = copysign(d, z.imag); 839 } else { 840 r.real = d; 841 r.imag = copysign(s, z.imag); 842 } 843 errno = 0; 844 return r; 845 } 846 847 848 /*[clinic input] 849 cmath.tan = cmath.acos 850 851 Return the tangent of z. 852 [clinic start generated code]*/ 853 854 static Py_complex 855 cmath_tan_impl(PyObject *module, Py_complex z) 856 /*[clinic end generated code: output=7c5f13158a72eb13 input=fc167e528767888e]*/ 857 { 858 /* tan(z) = -i tanh(iz) */ 859 Py_complex s, r; 860 s.real = -z.imag; 861 s.imag = z.real; 862 s = cmath_tanh_impl(module, s); 863 r.real = s.imag; 864 r.imag = -s.real; 865 return r; 866 } 867 868 869 /* tanh(infinity + i*y) needs to be dealt with specially */ 870 static Py_complex tanh_special_values[7][7]; 871 872 /*[clinic input] 873 cmath.tanh = cmath.acos 874 875 Return the hyperbolic tangent of z. 876 [clinic start generated code]*/ 877 878 static Py_complex 879 cmath_tanh_impl(PyObject *module, Py_complex z) 880 /*[clinic end generated code: output=36d547ef7aca116c input=22f67f9dc6d29685]*/ 881 { 882 /* Formula: 883 884 tanh(x+iy) = (tanh(x)(1+tan(y)^2) + i tan(y)(1-tanh(x))^2) / 885 (1+tan(y)^2 tanh(x)^2) 886 887 To avoid excessive roundoff error, 1-tanh(x)^2 is better computed 888 as 1/cosh(x)^2. When abs(x) is large, we approximate 1-tanh(x)^2 889 by 4 exp(-2*x) instead, to avoid possible overflow in the 890 computation of cosh(x). 891 892 */ 893 894 Py_complex r; 895 double tx, ty, cx, txty, denom; 896 897 /* special treatment for tanh(+/-inf + iy) if y is finite and 898 nonzero */ 899 if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) { 900 if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) 901 && (z.imag != 0.)) { 902 if (z.real > 0) { 903 r.real = 1.0; 904 r.imag = copysign(0., 905 2.*sin(z.imag)*cos(z.imag)); 906 } 907 else { 908 r.real = -1.0; 909 r.imag = copysign(0., 910 2.*sin(z.imag)*cos(z.imag)); 911 } 912 } 913 else { 914 r = tanh_special_values[special_type(z.real)] 915 [special_type(z.imag)]; 916 } 917 /* need to set errno = EDOM if z.imag is +/-infinity and 918 z.real is finite */ 919 if (Py_IS_INFINITY(z.imag) && Py_IS_FINITE(z.real)) 920 errno = EDOM; 921 else 922 errno = 0; 923 return r; 924 } 925 926 /* danger of overflow in 2.*z.imag !*/ 927 if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) { 928 r.real = copysign(1., z.real); 929 r.imag = 4.*sin(z.imag)*cos(z.imag)*exp(-2.*fabs(z.real)); 930 } else { 931 tx = tanh(z.real); 932 ty = tan(z.imag); 933 cx = 1./cosh(z.real); 934 txty = tx*ty; 935 denom = 1. + txty*txty; 936 r.real = tx*(1.+ty*ty)/denom; 937 r.imag = ((ty/denom)*cx)*cx; 938 } 939 errno = 0; 940 return r; 941 } 942 943 944 /*[clinic input] 945 cmath.log 946 947 x: Py_complex 948 y_obj: object = NULL 949 / 950 951 The logarithm of z to the given base. 952 953 If the base not specified, returns the natural logarithm (base e) of z. 954 [clinic start generated code]*/ 955 956 static PyObject * 957 cmath_log_impl(PyObject *module, Py_complex x, PyObject *y_obj) 958 /*[clinic end generated code: output=4effdb7d258e0d94 input=ee0e823a7c6e68ea]*/ 959 { 960 Py_complex y; 961 962 errno = 0; 963 PyFPE_START_PROTECT("complex function", return 0) 964 x = c_log(x); 965 if (y_obj != NULL) { 966 y = PyComplex_AsCComplex(y_obj); 967 if (PyErr_Occurred()) { 968 return NULL; 969 } 970 y = c_log(y); 971 x = _Py_c_quot(x, y); 972 } 973 PyFPE_END_PROTECT(x) 974 if (errno != 0) 975 return math_error(); 976 return PyComplex_FromCComplex(x); 977 } 978 979 980 /* And now the glue to make them available from Python: */ 981 982 static PyObject * 983 math_error(void) 984 { 985 if (errno == EDOM) 986 PyErr_SetString(PyExc_ValueError, "math domain error"); 987 else if (errno == ERANGE) 988 PyErr_SetString(PyExc_OverflowError, "math range error"); 989 else /* Unexpected math error */ 990 PyErr_SetFromErrno(PyExc_ValueError); 991 return NULL; 992 } 993 994 995 /*[clinic input] 996 cmath.phase 997 998 z: Py_complex 999 / 1000 1001 Return argument, also known as the phase angle, of a complex. 1002 [clinic start generated code]*/ 1003 1004 static PyObject * 1005 cmath_phase_impl(PyObject *module, Py_complex z) 1006 /*[clinic end generated code: output=50725086a7bfd253 input=5cf75228ba94b69d]*/ 1007 { 1008 double phi; 1009 1010 errno = 0; 1011 PyFPE_START_PROTECT("arg function", return 0) 1012 phi = c_atan2(z); 1013 PyFPE_END_PROTECT(phi) 1014 if (errno != 0) 1015 return math_error(); 1016 else 1017 return PyFloat_FromDouble(phi); 1018 } 1019 1020 /*[clinic input] 1021 cmath.polar 1022 1023 z: Py_complex 1024 / 1025 1026 Convert a complex from rectangular coordinates to polar coordinates. 1027 1028 r is the distance from 0 and phi the phase angle. 1029 [clinic start generated code]*/ 1030 1031 static PyObject * 1032 cmath_polar_impl(PyObject *module, Py_complex z) 1033 /*[clinic end generated code: output=d0a8147c41dbb654 input=26c353574fd1a861]*/ 1034 { 1035 double r, phi; 1036 1037 errno = 0; 1038 PyFPE_START_PROTECT("polar function", return 0) 1039 phi = c_atan2(z); /* should not cause any exception */ 1040 r = _Py_c_abs(z); /* sets errno to ERANGE on overflow */ 1041 PyFPE_END_PROTECT(r) 1042 if (errno != 0) 1043 return math_error(); 1044 else 1045 return Py_BuildValue("dd", r, phi); 1046 } 1047 1048 /* 1049 rect() isn't covered by the C99 standard, but it's not too hard to 1050 figure out 'spirit of C99' rules for special value handing: 1051 1052 rect(x, t) should behave like exp(log(x) + it) for positive-signed x 1053 rect(x, t) should behave like -exp(log(-x) + it) for negative-signed x 1054 rect(nan, t) should behave like exp(nan + it), except that rect(nan, 0) 1055 gives nan +- i0 with the sign of the imaginary part unspecified. 1056 1057 */ 1058 1059 static Py_complex rect_special_values[7][7]; 1060 1061 /*[clinic input] 1062 cmath.rect 1063 1064 r: double 1065 phi: double 1066 / 1067 1068 Convert from polar coordinates to rectangular coordinates. 1069 [clinic start generated code]*/ 1070 1071 static PyObject * 1072 cmath_rect_impl(PyObject *module, double r, double phi) 1073 /*[clinic end generated code: output=385a0690925df2d5 input=24c5646d147efd69]*/ 1074 { 1075 Py_complex z; 1076 errno = 0; 1077 PyFPE_START_PROTECT("rect function", return 0) 1078 1079 /* deal with special values */ 1080 if (!Py_IS_FINITE(r) || !Py_IS_FINITE(phi)) { 1081 /* if r is +/-infinity and phi is finite but nonzero then 1082 result is (+-INF +-INF i), but we need to compute cos(phi) 1083 and sin(phi) to figure out the signs. */ 1084 if (Py_IS_INFINITY(r) && (Py_IS_FINITE(phi) 1085 && (phi != 0.))) { 1086 if (r > 0) { 1087 z.real = copysign(INF, cos(phi)); 1088 z.imag = copysign(INF, sin(phi)); 1089 } 1090 else { 1091 z.real = -copysign(INF, cos(phi)); 1092 z.imag = -copysign(INF, sin(phi)); 1093 } 1094 } 1095 else { 1096 z = rect_special_values[special_type(r)] 1097 [special_type(phi)]; 1098 } 1099 /* need to set errno = EDOM if r is a nonzero number and phi 1100 is infinite */ 1101 if (r != 0. && !Py_IS_NAN(r) && Py_IS_INFINITY(phi)) 1102 errno = EDOM; 1103 else 1104 errno = 0; 1105 } 1106 else if (phi == 0.0) { 1107 /* Workaround for buggy results with phi=-0.0 on OS X 10.8. See 1108 bugs.python.org/issue18513. */ 1109 z.real = r; 1110 z.imag = r * phi; 1111 errno = 0; 1112 } 1113 else { 1114 z.real = r * cos(phi); 1115 z.imag = r * sin(phi); 1116 errno = 0; 1117 } 1118 1119 PyFPE_END_PROTECT(z) 1120 if (errno != 0) 1121 return math_error(); 1122 else 1123 return PyComplex_FromCComplex(z); 1124 } 1125 1126 /*[clinic input] 1127 cmath.isfinite = cmath.polar 1128 1129 Return True if both the real and imaginary parts of z are finite, else False. 1130 [clinic start generated code]*/ 1131 1132 static PyObject * 1133 cmath_isfinite_impl(PyObject *module, Py_complex z) 1134 /*[clinic end generated code: output=ac76611e2c774a36 input=848e7ee701895815]*/ 1135 { 1136 return PyBool_FromLong(Py_IS_FINITE(z.real) && Py_IS_FINITE(z.imag)); 1137 } 1138 1139 /*[clinic input] 1140 cmath.isnan = cmath.polar 1141 1142 Checks if the real or imaginary part of z not a number (NaN). 1143 [clinic start generated code]*/ 1144 1145 static PyObject * 1146 cmath_isnan_impl(PyObject *module, Py_complex z) 1147 /*[clinic end generated code: output=e7abf6e0b28beab7 input=71799f5d284c9baf]*/ 1148 { 1149 return PyBool_FromLong(Py_IS_NAN(z.real) || Py_IS_NAN(z.imag)); 1150 } 1151 1152 /*[clinic input] 1153 cmath.isinf = cmath.polar 1154 1155 Checks if the real or imaginary part of z is infinite. 1156 [clinic start generated code]*/ 1157 1158 static PyObject * 1159 cmath_isinf_impl(PyObject *module, Py_complex z) 1160 /*[clinic end generated code: output=502a75a79c773469 input=363df155c7181329]*/ 1161 { 1162 return PyBool_FromLong(Py_IS_INFINITY(z.real) || 1163 Py_IS_INFINITY(z.imag)); 1164 } 1165 1166 /*[clinic input] 1167 cmath.isclose -> bool 1168 1169 a: Py_complex 1170 b: Py_complex 1171 * 1172 rel_tol: double = 1e-09 1173 maximum difference for being considered "close", relative to the 1174 magnitude of the input values 1175 abs_tol: double = 0.0 1176 maximum difference for being considered "close", regardless of the 1177 magnitude of the input values 1178 1179 Determine whether two complex numbers are close in value. 1180 1181 Return True if a is close in value to b, and False otherwise. 1182 1183 For the values to be considered close, the difference between them must be 1184 smaller than at least one of the tolerances. 1185 1186 -inf, inf and NaN behave similarly to the IEEE 754 Standard. That is, NaN is 1187 not close to anything, even itself. inf and -inf are only close to themselves. 1188 [clinic start generated code]*/ 1189 1190 static int 1191 cmath_isclose_impl(PyObject *module, Py_complex a, Py_complex b, 1192 double rel_tol, double abs_tol) 1193 /*[clinic end generated code: output=8a2486cc6e0014d1 input=df9636d7de1d4ac3]*/ 1194 { 1195 double diff; 1196 1197 /* sanity check on the inputs */ 1198 if (rel_tol < 0.0 || abs_tol < 0.0 ) { 1199 PyErr_SetString(PyExc_ValueError, 1200 "tolerances must be non-negative"); 1201 return -1; 1202 } 1203 1204 if ( (a.real == b.real) && (a.imag == b.imag) ) { 1205 /* short circuit exact equality -- needed to catch two infinities of 1206 the same sign. And perhaps speeds things up a bit sometimes. 1207 */ 1208 return 1; 1209 } 1210 1211 /* This catches the case of two infinities of opposite sign, or 1212 one infinity and one finite number. Two infinities of opposite 1213 sign would otherwise have an infinite relative tolerance. 1214 Two infinities of the same sign are caught by the equality check 1215 above. 1216 */ 1217 1218 if (Py_IS_INFINITY(a.real) || Py_IS_INFINITY(a.imag) || 1219 Py_IS_INFINITY(b.real) || Py_IS_INFINITY(b.imag)) { 1220 return 0; 1221 } 1222 1223 /* now do the regular computation 1224 this is essentially the "weak" test from the Boost library 1225 */ 1226 1227 diff = _Py_c_abs(_Py_c_diff(a, b)); 1228 1229 return (((diff <= rel_tol * _Py_c_abs(b)) || 1230 (diff <= rel_tol * _Py_c_abs(a))) || 1231 (diff <= abs_tol)); 1232 } 1233 1234 PyDoc_STRVAR(module_doc, 1235 "This module is always available. It provides access to mathematical\n" 1236 "functions for complex numbers."); 1237 1238 static PyMethodDef cmath_methods[] = { 1239 CMATH_ACOS_METHODDEF 1240 CMATH_ACOSH_METHODDEF 1241 CMATH_ASIN_METHODDEF 1242 CMATH_ASINH_METHODDEF 1243 CMATH_ATAN_METHODDEF 1244 CMATH_ATANH_METHODDEF 1245 CMATH_COS_METHODDEF 1246 CMATH_COSH_METHODDEF 1247 CMATH_EXP_METHODDEF 1248 CMATH_ISCLOSE_METHODDEF 1249 CMATH_ISFINITE_METHODDEF 1250 CMATH_ISINF_METHODDEF 1251 CMATH_ISNAN_METHODDEF 1252 CMATH_LOG_METHODDEF 1253 CMATH_LOG10_METHODDEF 1254 CMATH_PHASE_METHODDEF 1255 CMATH_POLAR_METHODDEF 1256 CMATH_RECT_METHODDEF 1257 CMATH_SIN_METHODDEF 1258 CMATH_SINH_METHODDEF 1259 CMATH_SQRT_METHODDEF 1260 CMATH_TAN_METHODDEF 1261 CMATH_TANH_METHODDEF 1262 {NULL, NULL} /* sentinel */ 1263 }; 1264 1265 1266 static struct PyModuleDef cmathmodule = { 1267 PyModuleDef_HEAD_INIT, 1268 "cmath", 1269 module_doc, 1270 -1, 1271 cmath_methods, 1272 NULL, 1273 NULL, 1274 NULL, 1275 NULL 1276 }; 1277 1278 PyMODINIT_FUNC 1279 PyInit_cmath(void) 1280 { 1281 PyObject *m; 1282 1283 m = PyModule_Create(&cmathmodule); 1284 if (m == NULL) 1285 return NULL; 1286 1287 PyModule_AddObject(m, "pi", 1288 PyFloat_FromDouble(Py_MATH_PI)); 1289 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E)); 1290 PyModule_AddObject(m, "tau", PyFloat_FromDouble(Py_MATH_TAU)); /* 2pi */ 1291 PyModule_AddObject(m, "inf", PyFloat_FromDouble(m_inf())); 1292 PyModule_AddObject(m, "infj", PyComplex_FromCComplex(c_infj())); 1293 #if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN) 1294 PyModule_AddObject(m, "nan", PyFloat_FromDouble(m_nan())); 1295 PyModule_AddObject(m, "nanj", PyComplex_FromCComplex(c_nanj())); 1296 #endif 1297 1298 /* initialize special value tables */ 1299 1300 #define INIT_SPECIAL_VALUES(NAME, BODY) { Py_complex* p = (Py_complex*)NAME; BODY } 1301 #define C(REAL, IMAG) p->real = REAL; p->imag = IMAG; ++p; 1302 1303 INIT_SPECIAL_VALUES(acos_special_values, { 1304 C(P34,INF) C(P,INF) C(P,INF) C(P,-INF) C(P,-INF) C(P34,-INF) C(N,INF) 1305 C(P12,INF) C(U,U) C(U,U) C(U,U) C(U,U) C(P12,-INF) C(N,N) 1306 C(P12,INF) C(U,U) C(P12,0.) C(P12,-0.) C(U,U) C(P12,-INF) C(P12,N) 1307 C(P12,INF) C(U,U) C(P12,0.) C(P12,-0.) C(U,U) C(P12,-INF) C(P12,N) 1308 C(P12,INF) C(U,U) C(U,U) C(U,U) C(U,U) C(P12,-INF) C(N,N) 1309 C(P14,INF) C(0.,INF) C(0.,INF) C(0.,-INF) C(0.,-INF) C(P14,-INF) C(N,INF) 1310 C(N,INF) C(N,N) C(N,N) C(N,N) C(N,N) C(N,-INF) C(N,N) 1311 }) 1312 1313 INIT_SPECIAL_VALUES(acosh_special_values, { 1314 C(INF,-P34) C(INF,-P) C(INF,-P) C(INF,P) C(INF,P) C(INF,P34) C(INF,N) 1315 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N) 1316 C(INF,-P12) C(U,U) C(0.,-P12) C(0.,P12) C(U,U) C(INF,P12) C(N,N) 1317 C(INF,-P12) C(U,U) C(0.,-P12) C(0.,P12) C(U,U) C(INF,P12) C(N,N) 1318 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N) 1319 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N) 1320 C(INF,N) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,N) C(N,N) 1321 }) 1322 1323 INIT_SPECIAL_VALUES(asinh_special_values, { 1324 C(-INF,-P14) C(-INF,-0.) C(-INF,-0.) C(-INF,0.) C(-INF,0.) C(-INF,P14) C(-INF,N) 1325 C(-INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(-INF,P12) C(N,N) 1326 C(-INF,-P12) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(-INF,P12) C(N,N) 1327 C(INF,-P12) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(INF,P12) C(N,N) 1328 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N) 1329 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N) 1330 C(INF,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(INF,N) C(N,N) 1331 }) 1332 1333 INIT_SPECIAL_VALUES(atanh_special_values, { 1334 C(-0.,-P12) C(-0.,-P12) C(-0.,-P12) C(-0.,P12) C(-0.,P12) C(-0.,P12) C(-0.,N) 1335 C(-0.,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(-0.,P12) C(N,N) 1336 C(-0.,-P12) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(-0.,P12) C(-0.,N) 1337 C(0.,-P12) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,P12) C(0.,N) 1338 C(0.,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(0.,P12) C(N,N) 1339 C(0.,-P12) C(0.,-P12) C(0.,-P12) C(0.,P12) C(0.,P12) C(0.,P12) C(0.,N) 1340 C(0.,-P12) C(N,N) C(N,N) C(N,N) C(N,N) C(0.,P12) C(N,N) 1341 }) 1342 1343 INIT_SPECIAL_VALUES(cosh_special_values, { 1344 C(INF,N) C(U,U) C(INF,0.) C(INF,-0.) C(U,U) C(INF,N) C(INF,N) 1345 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N) 1346 C(N,0.) C(U,U) C(1.,0.) C(1.,-0.) C(U,U) C(N,0.) C(N,0.) 1347 C(N,0.) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(N,0.) C(N,0.) 1348 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N) 1349 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N) 1350 C(N,N) C(N,N) C(N,0.) C(N,0.) C(N,N) C(N,N) C(N,N) 1351 }) 1352 1353 INIT_SPECIAL_VALUES(exp_special_values, { 1354 C(0.,0.) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,0.) C(0.,0.) 1355 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N) 1356 C(N,N) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(N,N) C(N,N) 1357 C(N,N) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(N,N) C(N,N) 1358 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N) 1359 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N) 1360 C(N,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(N,N) C(N,N) 1361 }) 1362 1363 INIT_SPECIAL_VALUES(log_special_values, { 1364 C(INF,-P34) C(INF,-P) C(INF,-P) C(INF,P) C(INF,P) C(INF,P34) C(INF,N) 1365 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N) 1366 C(INF,-P12) C(U,U) C(-INF,-P) C(-INF,P) C(U,U) C(INF,P12) C(N,N) 1367 C(INF,-P12) C(U,U) C(-INF,-0.) C(-INF,0.) C(U,U) C(INF,P12) C(N,N) 1368 C(INF,-P12) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,P12) C(N,N) 1369 C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N) 1370 C(INF,N) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,N) C(N,N) 1371 }) 1372 1373 INIT_SPECIAL_VALUES(sinh_special_values, { 1374 C(INF,N) C(U,U) C(-INF,-0.) C(-INF,0.) C(U,U) C(INF,N) C(INF,N) 1375 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N) 1376 C(0.,N) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(0.,N) C(0.,N) 1377 C(0.,N) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,N) C(0.,N) 1378 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N) 1379 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N) 1380 C(N,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(N,N) C(N,N) 1381 }) 1382 1383 INIT_SPECIAL_VALUES(sqrt_special_values, { 1384 C(INF,-INF) C(0.,-INF) C(0.,-INF) C(0.,INF) C(0.,INF) C(INF,INF) C(N,INF) 1385 C(INF,-INF) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,INF) C(N,N) 1386 C(INF,-INF) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(INF,INF) C(N,N) 1387 C(INF,-INF) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(INF,INF) C(N,N) 1388 C(INF,-INF) C(U,U) C(U,U) C(U,U) C(U,U) C(INF,INF) C(N,N) 1389 C(INF,-INF) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,INF) C(INF,N) 1390 C(INF,-INF) C(N,N) C(N,N) C(N,N) C(N,N) C(INF,INF) C(N,N) 1391 }) 1392 1393 INIT_SPECIAL_VALUES(tanh_special_values, { 1394 C(-1.,0.) C(U,U) C(-1.,-0.) C(-1.,0.) C(U,U) C(-1.,0.) C(-1.,0.) 1395 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N) 1396 C(N,N) C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(N,N) C(N,N) 1397 C(N,N) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(N,N) C(N,N) 1398 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N) 1399 C(1.,0.) C(U,U) C(1.,-0.) C(1.,0.) C(U,U) C(1.,0.) C(1.,0.) 1400 C(N,N) C(N,N) C(N,-0.) C(N,0.) C(N,N) C(N,N) C(N,N) 1401 }) 1402 1403 INIT_SPECIAL_VALUES(rect_special_values, { 1404 C(INF,N) C(U,U) C(-INF,0.) C(-INF,-0.) C(U,U) C(INF,N) C(INF,N) 1405 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N) 1406 C(0.,0.) C(U,U) C(-0.,0.) C(-0.,-0.) C(U,U) C(0.,0.) C(0.,0.) 1407 C(0.,0.) C(U,U) C(0.,-0.) C(0.,0.) C(U,U) C(0.,0.) C(0.,0.) 1408 C(N,N) C(U,U) C(U,U) C(U,U) C(U,U) C(N,N) C(N,N) 1409 C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N) 1410 C(N,N) C(N,N) C(N,0.) C(N,0.) C(N,N) C(N,N) C(N,N) 1411 }) 1412 return m; 1413 } 1414