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
m_inf(void)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
c_infj(void)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
m_nan(void)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
c_nanj(void)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
special_type(double d)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
cmath_acos_impl(PyObject * module,Py_complex z)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
cmath_acosh_impl(PyObject * module,Py_complex z)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
cmath_asin_impl(PyObject * module,Py_complex z)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
cmath_asinh_impl(PyObject * module,Py_complex z)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
cmath_atan_impl(PyObject * module,Py_complex z)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
c_atan2(Py_complex z)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
cmath_atanh_impl(PyObject * module,Py_complex z)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
cmath_cos_impl(PyObject * module,Py_complex z)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
cmath_cosh_impl(PyObject * module,Py_complex z)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
cmath_exp_impl(PyObject * module,Py_complex z)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
c_log(Py_complex z)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
cmath_log10_impl(PyObject * module,Py_complex z)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
cmath_sin_impl(PyObject * module,Py_complex z)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
cmath_sinh_impl(PyObject * module,Py_complex z)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
cmath_sqrt_impl(PyObject * module,Py_complex z)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
cmath_tan_impl(PyObject * module,Py_complex z)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
cmath_tanh_impl(PyObject * module,Py_complex z)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 z as x: Py_complex
948 base as y_obj: object = NULL
949 /
950
951 log(z[, base]) -> 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 *
cmath_log_impl(PyObject * module,Py_complex x,PyObject * y_obj)957 cmath_log_impl(PyObject *module, Py_complex x, PyObject *y_obj)
958 /*[clinic end generated code: output=4effdb7d258e0d94 input=230ed3a71ecd000a]*/
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 *
math_error(void)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 *
cmath_phase_impl(PyObject * module,Py_complex z)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 *
cmath_polar_impl(PyObject * module,Py_complex z)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 *
cmath_rect_impl(PyObject * module,double r,double phi)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 *
cmath_isfinite_impl(PyObject * module,Py_complex z)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 *
cmath_isnan_impl(PyObject * module,Py_complex z)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 *
cmath_isinf_impl(PyObject * module,Py_complex z)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
cmath_isclose_impl(PyObject * module,Py_complex a,Py_complex b,double rel_tol,double abs_tol)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 provides access to mathematical functions for complex\n"
1236 "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
PyInit_cmath(void)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