• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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     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 *
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=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 *
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 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
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