~ubuntu-branches/ubuntu/maverick/python3.1/maverick

« back to all changes in this revision

Viewing changes to Modules/cmathmodule.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2009-03-23 00:01:27 UTC
  • Revision ID: james.westby@ubuntu.com-20090323000127-5fstfxju4ufrhthq
Tags: upstream-3.1~a1+20090322
ImportĀ upstreamĀ versionĀ 3.1~a1+20090322

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* Complex math module */
 
2
 
 
3
/* much code borrowed from mathmodule.c */
 
4
 
 
5
#include "Python.h"
 
6
/* we need DBL_MAX, DBL_MIN, DBL_EPSILON, DBL_MANT_DIG and FLT_RADIX from
 
7
   float.h.  We assume that FLT_RADIX is either 2 or 16. */
 
8
#include <float.h>
 
9
 
 
10
#if (FLT_RADIX != 2 && FLT_RADIX != 16)
 
11
#error "Modules/cmathmodule.c expects FLT_RADIX to be 2 or 16"
 
12
#endif
 
13
 
 
14
#ifndef M_LN2
 
15
#define M_LN2 (0.6931471805599453094) /* natural log of 2 */
 
16
#endif
 
17
 
 
18
#ifndef M_LN10
 
19
#define M_LN10 (2.302585092994045684) /* natural log of 10 */
 
20
#endif
 
21
 
 
22
/*
 
23
   CM_LARGE_DOUBLE is used to avoid spurious overflow in the sqrt, log,
 
24
   inverse trig and inverse hyperbolic trig functions.  Its log is used in the
 
25
   evaluation of exp, cos, cosh, sin, sinh, tan, and tanh to avoid unecessary
 
26
   overflow.
 
27
 */
 
28
 
 
29
#define CM_LARGE_DOUBLE (DBL_MAX/4.)
 
30
#define CM_SQRT_LARGE_DOUBLE (sqrt(CM_LARGE_DOUBLE))
 
31
#define CM_LOG_LARGE_DOUBLE (log(CM_LARGE_DOUBLE))
 
32
#define CM_SQRT_DBL_MIN (sqrt(DBL_MIN))
 
33
 
 
34
/* 
 
35
   CM_SCALE_UP is an odd integer chosen such that multiplication by
 
36
   2**CM_SCALE_UP is sufficient to turn a subnormal into a normal.
 
37
   CM_SCALE_DOWN is (-(CM_SCALE_UP+1)/2).  These scalings are used to compute
 
38
   square roots accurately when the real and imaginary parts of the argument
 
39
   are subnormal.
 
40
*/
 
41
 
 
42
#if FLT_RADIX==2
 
43
#define CM_SCALE_UP (2*(DBL_MANT_DIG/2) + 1)
 
44
#elif FLT_RADIX==16
 
45
#define CM_SCALE_UP (4*DBL_MANT_DIG+1)
 
46
#endif
 
47
#define CM_SCALE_DOWN (-(CM_SCALE_UP+1)/2)
 
48
 
 
49
/* forward declarations */
 
50
static Py_complex c_asinh(Py_complex);
 
51
static Py_complex c_atanh(Py_complex);
 
52
static Py_complex c_cosh(Py_complex);
 
53
static Py_complex c_sinh(Py_complex);
 
54
static Py_complex c_sqrt(Py_complex);
 
55
static Py_complex c_tanh(Py_complex);
 
56
static PyObject * math_error(void);
 
57
 
 
58
/* Code to deal with special values (infinities, NaNs, etc.). */
 
59
 
 
60
/* special_type takes a double and returns an integer code indicating
 
61
   the type of the double as follows:
 
62
*/
 
63
 
 
64
enum special_types {
 
65
        ST_NINF,        /* 0, negative infinity */
 
66
        ST_NEG,         /* 1, negative finite number (nonzero) */
 
67
        ST_NZERO,       /* 2, -0. */
 
68
        ST_PZERO,       /* 3, +0. */
 
69
        ST_POS,         /* 4, positive finite number (nonzero) */
 
70
        ST_PINF,        /* 5, positive infinity */
 
71
        ST_NAN,         /* 6, Not a Number */
 
72
};
 
73
 
 
74
static enum special_types
 
75
special_type(double d)
 
76
{
 
77
        if (Py_IS_FINITE(d)) {
 
78
                if (d != 0) {
 
79
                        if (copysign(1., d) == 1.)
 
80
                                return ST_POS;
 
81
                        else
 
82
                                return ST_NEG;
 
83
                }
 
84
                else {
 
85
                        if (copysign(1., d) == 1.)
 
86
                                return ST_PZERO;
 
87
                        else
 
88
                                return ST_NZERO;
 
89
                }
 
90
        }
 
91
        if (Py_IS_NAN(d))
 
92
                return ST_NAN;
 
93
        if (copysign(1., d) == 1.)
 
94
                return ST_PINF;
 
95
        else
 
96
                return ST_NINF;
 
97
}
 
98
 
 
99
#define SPECIAL_VALUE(z, table)                                         \
 
100
        if (!Py_IS_FINITE((z).real) || !Py_IS_FINITE((z).imag)) {       \
 
101
                errno = 0;                                              \
 
102
                return table[special_type((z).real)]                    \
 
103
                            [special_type((z).imag)];                   \
 
104
        }
 
105
 
 
106
#define P Py_MATH_PI
 
107
#define P14 0.25*Py_MATH_PI
 
108
#define P12 0.5*Py_MATH_PI
 
109
#define P34 0.75*Py_MATH_PI
 
110
#define INF Py_HUGE_VAL
 
111
#define N Py_NAN
 
112
#define U -9.5426319407711027e33 /* unlikely value, used as placeholder */
 
113
 
 
114
/* First, the C functions that do the real work.  Each of the c_*
 
115
   functions computes and returns the C99 Annex G recommended result
 
116
   and also sets errno as follows: errno = 0 if no floating-point
 
117
   exception is associated with the result; errno = EDOM if C99 Annex
 
118
   G recommends raising divide-by-zero or invalid for this result; and
 
119
   errno = ERANGE where the overflow floating-point signal should be
 
120
   raised.
 
121
*/
 
122
 
 
123
static Py_complex acos_special_values[7][7];
 
124
 
 
125
static Py_complex
 
126
c_acos(Py_complex z)
 
127
{
 
128
        Py_complex s1, s2, r;
 
129
 
 
130
        SPECIAL_VALUE(z, acos_special_values);
 
131
 
 
132
        if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
 
133
                /* avoid unnecessary overflow for large arguments */
 
134
                r.real = atan2(fabs(z.imag), z.real);
 
135
                /* split into cases to make sure that the branch cut has the
 
136
                   correct continuity on systems with unsigned zeros */
 
137
                if (z.real < 0.) {
 
138
                        r.imag = -copysign(log(hypot(z.real/2., z.imag/2.)) +
 
139
                                           M_LN2*2., z.imag);
 
140
                } else {
 
141
                        r.imag = copysign(log(hypot(z.real/2., z.imag/2.)) +
 
142
                                          M_LN2*2., -z.imag);
 
143
                }
 
144
        } else {
 
145
                s1.real = 1.-z.real;
 
146
                s1.imag = -z.imag;
 
147
                s1 = c_sqrt(s1);
 
148
                s2.real = 1.+z.real;
 
149
                s2.imag = z.imag;
 
150
                s2 = c_sqrt(s2);
 
151
                r.real = 2.*atan2(s1.real, s2.real);
 
152
                r.imag = asinh(s2.real*s1.imag - s2.imag*s1.real);
 
153
        }
 
154
        errno = 0;
 
155
        return r;
 
156
}
 
157
 
 
158
PyDoc_STRVAR(c_acos_doc,
 
159
"acos(x)\n"
 
160
"\n"
 
161
"Return the arc cosine of x.");
 
162
 
 
163
 
 
164
static Py_complex acosh_special_values[7][7];
 
165
 
 
166
static Py_complex
 
167
c_acosh(Py_complex z)
 
168
{
 
169
        Py_complex s1, s2, r;
 
170
 
 
171
        SPECIAL_VALUE(z, acosh_special_values);
 
172
 
 
173
        if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
 
174
                /* avoid unnecessary overflow for large arguments */
 
175
                r.real = log(hypot(z.real/2., z.imag/2.)) + M_LN2*2.;
 
176
                r.imag = atan2(z.imag, z.real);
 
177
        } else {
 
178
                s1.real = z.real - 1.;
 
179
                s1.imag = z.imag;
 
180
                s1 = c_sqrt(s1);
 
181
                s2.real = z.real + 1.;
 
182
                s2.imag = z.imag;
 
183
                s2 = c_sqrt(s2);
 
184
                r.real = asinh(s1.real*s2.real + s1.imag*s2.imag);
 
185
                r.imag = 2.*atan2(s1.imag, s2.real);
 
186
        }
 
187
        errno = 0;
 
188
        return r;
 
189
}
 
190
 
 
191
PyDoc_STRVAR(c_acosh_doc,
 
192
"acosh(x)\n"
 
193
"\n"
 
194
"Return the hyperbolic arccosine of x.");
 
195
 
 
196
 
 
197
static Py_complex
 
198
c_asin(Py_complex z)
 
199
{
 
200
        /* asin(z) = -i asinh(iz) */
 
201
        Py_complex s, r;
 
202
        s.real = -z.imag;
 
203
        s.imag = z.real;
 
204
        s = c_asinh(s);
 
205
        r.real = s.imag;
 
206
        r.imag = -s.real;
 
207
        return r;
 
208
}
 
209
 
 
210
PyDoc_STRVAR(c_asin_doc,
 
211
"asin(x)\n"
 
212
"\n"
 
213
"Return the arc sine of x.");
 
214
 
 
215
 
 
216
static Py_complex asinh_special_values[7][7];
 
217
 
 
218
static Py_complex
 
219
c_asinh(Py_complex z)
 
220
{
 
221
        Py_complex s1, s2, r;
 
222
 
 
223
        SPECIAL_VALUE(z, asinh_special_values);
 
224
 
 
225
        if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
 
226
                if (z.imag >= 0.) {
 
227
                        r.real = copysign(log(hypot(z.real/2., z.imag/2.)) +
 
228
                                          M_LN2*2., z.real);
 
229
                } else {
 
230
                        r.real = -copysign(log(hypot(z.real/2., z.imag/2.)) +
 
231
                                           M_LN2*2., -z.real);
 
232
                }
 
233
                r.imag = atan2(z.imag, fabs(z.real));
 
234
        } else {
 
235
                s1.real = 1.+z.imag;
 
236
                s1.imag = -z.real;
 
237
                s1 = c_sqrt(s1);
 
238
                s2.real = 1.-z.imag;
 
239
                s2.imag = z.real;
 
240
                s2 = c_sqrt(s2);
 
241
                r.real = asinh(s1.real*s2.imag-s2.real*s1.imag);
 
242
                r.imag = atan2(z.imag, s1.real*s2.real-s1.imag*s2.imag);
 
243
        }
 
244
        errno = 0;
 
245
        return r;
 
246
}
 
247
 
 
248
PyDoc_STRVAR(c_asinh_doc,
 
249
"asinh(x)\n"
 
250
"\n"
 
251
"Return the hyperbolic arc sine of x.");
 
252
 
 
253
 
 
254
static Py_complex
 
255
c_atan(Py_complex z)
 
256
{
 
257
        /* atan(z) = -i atanh(iz) */
 
258
        Py_complex s, r;
 
259
        s.real = -z.imag;
 
260
        s.imag = z.real;
 
261
        s = c_atanh(s);
 
262
        r.real = s.imag;
 
263
        r.imag = -s.real;
 
264
        return r;
 
265
}
 
266
 
 
267
/* Windows screws up atan2 for inf and nan, and alpha Tru64 5.1 doesn't follow
 
268
   C99 for atan2(0., 0.). */
 
269
static double
 
270
c_atan2(Py_complex z)
 
271
{
 
272
        if (Py_IS_NAN(z.real) || Py_IS_NAN(z.imag))
 
273
                return Py_NAN;
 
274
        if (Py_IS_INFINITY(z.imag)) {
 
275
                if (Py_IS_INFINITY(z.real)) {
 
276
                        if (copysign(1., z.real) == 1.)
 
277
                                /* atan2(+-inf, +inf) == +-pi/4 */
 
278
                                return copysign(0.25*Py_MATH_PI, z.imag);
 
279
                        else
 
280
                                /* atan2(+-inf, -inf) == +-pi*3/4 */
 
281
                                return copysign(0.75*Py_MATH_PI, z.imag);
 
282
                }
 
283
                /* atan2(+-inf, x) == +-pi/2 for finite x */
 
284
                return copysign(0.5*Py_MATH_PI, z.imag);
 
285
        }
 
286
        if (Py_IS_INFINITY(z.real) || z.imag == 0.) {
 
287
                if (copysign(1., z.real) == 1.)
 
288
                        /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
 
289
                        return copysign(0., z.imag);
 
290
                else
 
291
                        /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
 
292
                        return copysign(Py_MATH_PI, z.imag);
 
293
        }
 
294
        return atan2(z.imag, z.real);
 
295
}
 
296
 
 
297
PyDoc_STRVAR(c_atan_doc,
 
298
"atan(x)\n"
 
299
"\n"
 
300
"Return the arc tangent of x.");
 
301
 
 
302
 
 
303
static Py_complex atanh_special_values[7][7];
 
304
 
 
305
static Py_complex
 
306
c_atanh(Py_complex z)
 
307
{
 
308
        Py_complex r;
 
309
        double ay, h;
 
310
 
 
311
        SPECIAL_VALUE(z, atanh_special_values);
 
312
 
 
313
        /* Reduce to case where z.real >= 0., using atanh(z) = -atanh(-z). */
 
314
        if (z.real < 0.) {
 
315
                return c_neg(c_atanh(c_neg(z)));
 
316
        }
 
317
 
 
318
        ay = fabs(z.imag);
 
319
        if (z.real > CM_SQRT_LARGE_DOUBLE || ay > CM_SQRT_LARGE_DOUBLE) {
 
320
                /*
 
321
                   if abs(z) is large then we use the approximation
 
322
                   atanh(z) ~ 1/z +/- i*pi/2 (+/- depending on the sign
 
323
                   of z.imag)
 
324
                */
 
325
                h = hypot(z.real/2., z.imag/2.);  /* safe from overflow */
 
326
                r.real = z.real/4./h/h;
 
327
                /* the two negations in the next line cancel each other out
 
328
                   except when working with unsigned zeros: they're there to
 
329
                   ensure that the branch cut has the correct continuity on
 
330
                   systems that don't support signed zeros */
 
331
                r.imag = -copysign(Py_MATH_PI/2., -z.imag);
 
332
                errno = 0;
 
333
        } else if (z.real == 1. && ay < CM_SQRT_DBL_MIN) {
 
334
                /* C99 standard says:  atanh(1+/-0.) should be inf +/- 0i */
 
335
                if (ay == 0.) {
 
336
                        r.real = INF;
 
337
                        r.imag = z.imag;
 
338
                        errno = EDOM;
 
339
                } else {
 
340
                        r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.)));
 
341
                        r.imag = copysign(atan2(2., -ay)/2, z.imag);
 
342
                        errno = 0;
 
343
                }
 
344
        } else {
 
345
                r.real = log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.;
 
346
                r.imag = -atan2(-2.*z.imag, (1-z.real)*(1+z.real) - ay*ay)/2.;
 
347
                errno = 0;
 
348
        }
 
349
        return r;
 
350
}
 
351
 
 
352
PyDoc_STRVAR(c_atanh_doc,
 
353
"atanh(x)\n"
 
354
"\n"
 
355
"Return the hyperbolic arc tangent of x.");
 
356
 
 
357
 
 
358
static Py_complex
 
359
c_cos(Py_complex z)
 
360
{
 
361
        /* cos(z) = cosh(iz) */
 
362
        Py_complex r;
 
363
        r.real = -z.imag;
 
364
        r.imag = z.real;
 
365
        r = c_cosh(r);
 
366
        return r;
 
367
}
 
368
 
 
369
PyDoc_STRVAR(c_cos_doc,
 
370
"cos(x)\n"
 
371
"\n"
 
372
"Return the cosine of x.");
 
373
 
 
374
 
 
375
/* cosh(infinity + i*y) needs to be dealt with specially */
 
376
static Py_complex cosh_special_values[7][7];
 
377
 
 
378
static Py_complex
 
379
c_cosh(Py_complex z)
 
380
{
 
381
        Py_complex r;
 
382
        double x_minus_one;
 
383
 
 
384
        /* special treatment for cosh(+/-inf + iy) if y is not a NaN */
 
385
        if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
 
386
                if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) &&
 
387
                    (z.imag != 0.)) {
 
388
                        if (z.real > 0) {
 
389
                                r.real = copysign(INF, cos(z.imag));
 
390
                                r.imag = copysign(INF, sin(z.imag));
 
391
                        }
 
392
                        else {
 
393
                                r.real = copysign(INF, cos(z.imag));
 
394
                                r.imag = -copysign(INF, sin(z.imag));
 
395
                        }
 
396
                }
 
397
                else {
 
398
                        r = cosh_special_values[special_type(z.real)]
 
399
                                               [special_type(z.imag)];
 
400
                }
 
401
                /* need to set errno = EDOM if y is +/- infinity and x is not
 
402
                   a NaN */
 
403
                if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
 
404
                        errno = EDOM;
 
405
                else
 
406
                        errno = 0;
 
407
                return r;
 
408
        }
 
409
 
 
410
        if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
 
411
                /* deal correctly with cases where cosh(z.real) overflows but
 
412
                   cosh(z) does not. */
 
413
                x_minus_one = z.real - copysign(1., z.real);
 
414
                r.real = cos(z.imag) * cosh(x_minus_one) * Py_MATH_E;
 
415
                r.imag = sin(z.imag) * sinh(x_minus_one) * Py_MATH_E;
 
416
        } else {
 
417
                r.real = cos(z.imag) * cosh(z.real);
 
418
                r.imag = sin(z.imag) * sinh(z.real);
 
419
        }
 
420
        /* detect overflow, and set errno accordingly */
 
421
        if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
 
422
                errno = ERANGE;
 
423
        else
 
424
                errno = 0;
 
425
        return r;
 
426
}
 
427
 
 
428
PyDoc_STRVAR(c_cosh_doc,
 
429
"cosh(x)\n"
 
430
"\n"
 
431
"Return the hyperbolic cosine of x.");
 
432
 
 
433
 
 
434
/* exp(infinity + i*y) and exp(-infinity + i*y) need special treatment for
 
435
   finite y */
 
436
static Py_complex exp_special_values[7][7];
 
437
 
 
438
static Py_complex
 
439
c_exp(Py_complex z)
 
440
{
 
441
        Py_complex r;
 
442
        double l;
 
443
 
 
444
        if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
 
445
                if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
 
446
                    && (z.imag != 0.)) {
 
447
                        if (z.real > 0) {
 
448
                                r.real = copysign(INF, cos(z.imag));
 
449
                                r.imag = copysign(INF, sin(z.imag));
 
450
                        }
 
451
                        else {
 
452
                                r.real = copysign(0., cos(z.imag));
 
453
                                r.imag = copysign(0., sin(z.imag));
 
454
                        }
 
455
                }
 
456
                else {
 
457
                        r = exp_special_values[special_type(z.real)]
 
458
                                              [special_type(z.imag)];
 
459
                }
 
460
                /* need to set errno = EDOM if y is +/- infinity and x is not
 
461
                   a NaN and not -infinity */
 
462
                if (Py_IS_INFINITY(z.imag) &&
 
463
                    (Py_IS_FINITE(z.real) ||
 
464
                     (Py_IS_INFINITY(z.real) && z.real > 0)))
 
465
                        errno = EDOM;
 
466
                else
 
467
                        errno = 0;
 
468
                return r;
 
469
        }
 
470
 
 
471
        if (z.real > CM_LOG_LARGE_DOUBLE) {
 
472
                l = exp(z.real-1.);
 
473
                r.real = l*cos(z.imag)*Py_MATH_E;
 
474
                r.imag = l*sin(z.imag)*Py_MATH_E;
 
475
        } else {
 
476
                l = exp(z.real);
 
477
                r.real = l*cos(z.imag);
 
478
                r.imag = l*sin(z.imag);
 
479
        }
 
480
        /* detect overflow, and set errno accordingly */
 
481
        if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
 
482
                errno = ERANGE;
 
483
        else
 
484
                errno = 0;
 
485
        return r;
 
486
}
 
487
 
 
488
PyDoc_STRVAR(c_exp_doc,
 
489
"exp(x)\n"
 
490
"\n"
 
491
"Return the exponential value e**x.");
 
492
 
 
493
 
 
494
static Py_complex log_special_values[7][7];
 
495
 
 
496
static Py_complex
 
497
c_log(Py_complex z)
 
498
{
 
499
        /*
 
500
           The usual formula for the real part is log(hypot(z.real, z.imag)).
 
501
           There are four situations where this formula is potentially
 
502
           problematic:
 
503
 
 
504
           (1) the absolute value of z is subnormal.  Then hypot is subnormal,
 
505
           so has fewer than the usual number of bits of accuracy, hence may
 
506
           have large relative error.  This then gives a large absolute error
 
507
           in the log.  This can be solved by rescaling z by a suitable power
 
508
           of 2.
 
509
 
 
510
           (2) the absolute value of z is greater than DBL_MAX (e.g. when both
 
511
           z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX)
 
512
           Again, rescaling solves this.
 
513
 
 
514
           (3) the absolute value of z is close to 1.  In this case it's
 
515
           difficult to achieve good accuracy, at least in part because a
 
516
           change of 1ulp in the real or imaginary part of z can result in a
 
517
           change of billions of ulps in the correctly rounded answer.
 
518
 
 
519
           (4) z = 0.  The simplest thing to do here is to call the
 
520
           floating-point log with an argument of 0, and let its behaviour
 
521
           (returning -infinity, signaling a floating-point exception, setting
 
522
           errno, or whatever) determine that of c_log.  So the usual formula
 
523
           is fine here.
 
524
 
 
525
         */
 
526
 
 
527
        Py_complex r;
 
528
        double ax, ay, am, an, h;
 
529
 
 
530
        SPECIAL_VALUE(z, log_special_values);
 
531
 
 
532
        ax = fabs(z.real);
 
533
        ay = fabs(z.imag);
 
534
 
 
535
        if (ax > CM_LARGE_DOUBLE || ay > CM_LARGE_DOUBLE) {
 
536
                r.real = log(hypot(ax/2., ay/2.)) + M_LN2;
 
537
        } else if (ax < DBL_MIN && ay < DBL_MIN) {
 
538
                if (ax > 0. || ay > 0.) {
 
539
                        /* catch cases where hypot(ax, ay) is subnormal */
 
540
                        r.real = log(hypot(ldexp(ax, DBL_MANT_DIG),
 
541
                                 ldexp(ay, DBL_MANT_DIG))) - DBL_MANT_DIG*M_LN2;
 
542
                }
 
543
                else {
 
544
                        /* log(+/-0. +/- 0i) */
 
545
                        r.real = -INF;
 
546
                        r.imag = atan2(z.imag, z.real);
 
547
                        errno = EDOM;
 
548
                        return r;
 
549
                }
 
550
        } else {
 
551
                h = hypot(ax, ay);
 
552
                if (0.71 <= h && h <= 1.73) {
 
553
                        am = ax > ay ? ax : ay;  /* max(ax, ay) */
 
554
                        an = ax > ay ? ay : ax;  /* min(ax, ay) */
 
555
                        r.real = log1p((am-1)*(am+1)+an*an)/2.;
 
556
                } else {
 
557
                        r.real = log(h);
 
558
                }
 
559
        }
 
560
        r.imag = atan2(z.imag, z.real);
 
561
        errno = 0;
 
562
        return r;
 
563
}
 
564
 
 
565
 
 
566
static Py_complex
 
567
c_log10(Py_complex z)
 
568
{
 
569
        Py_complex r;
 
570
        int errno_save;
 
571
 
 
572
        r = c_log(z);
 
573
        errno_save = errno; /* just in case the divisions affect errno */
 
574
        r.real = r.real / M_LN10;
 
575
        r.imag = r.imag / M_LN10;
 
576
        errno = errno_save;
 
577
        return r;
 
578
}
 
579
 
 
580
PyDoc_STRVAR(c_log10_doc,
 
581
"log10(x)\n"
 
582
"\n"
 
583
"Return the base-10 logarithm of x.");
 
584
 
 
585
 
 
586
static Py_complex
 
587
c_sin(Py_complex z)
 
588
{
 
589
        /* sin(z) = -i sin(iz) */
 
590
        Py_complex s, r;
 
591
        s.real = -z.imag;
 
592
        s.imag = z.real;
 
593
        s = c_sinh(s);
 
594
        r.real = s.imag;
 
595
        r.imag = -s.real;
 
596
        return r;
 
597
}
 
598
 
 
599
PyDoc_STRVAR(c_sin_doc,
 
600
"sin(x)\n"
 
601
"\n"
 
602
"Return the sine of x.");
 
603
 
 
604
 
 
605
/* sinh(infinity + i*y) needs to be dealt with specially */
 
606
static Py_complex sinh_special_values[7][7];
 
607
 
 
608
static Py_complex
 
609
c_sinh(Py_complex z)
 
610
{
 
611
        Py_complex r;
 
612
        double x_minus_one;
 
613
 
 
614
        /* special treatment for sinh(+/-inf + iy) if y is finite and
 
615
           nonzero */
 
616
        if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
 
617
                if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
 
618
                    && (z.imag != 0.)) {
 
619
                        if (z.real > 0) {
 
620
                                r.real = copysign(INF, cos(z.imag));
 
621
                                r.imag = copysign(INF, sin(z.imag));
 
622
                        }
 
623
                        else {
 
624
                                r.real = -copysign(INF, cos(z.imag));
 
625
                                r.imag = copysign(INF, sin(z.imag));
 
626
                        }
 
627
                }
 
628
                else {
 
629
                        r = sinh_special_values[special_type(z.real)]
 
630
                                               [special_type(z.imag)];
 
631
                }
 
632
                /* need to set errno = EDOM if y is +/- infinity and x is not
 
633
                   a NaN */
 
634
                if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
 
635
                        errno = EDOM;
 
636
                else
 
637
                        errno = 0;
 
638
                return r;
 
639
        }
 
640
 
 
641
        if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
 
642
                x_minus_one = z.real - copysign(1., z.real);
 
643
                r.real = cos(z.imag) * sinh(x_minus_one) * Py_MATH_E;
 
644
                r.imag = sin(z.imag) * cosh(x_minus_one) * Py_MATH_E;
 
645
        } else {
 
646
                r.real = cos(z.imag) * sinh(z.real);
 
647
                r.imag = sin(z.imag) * cosh(z.real);
 
648
        }
 
649
        /* detect overflow, and set errno accordingly */
 
650
        if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
 
651
                errno = ERANGE;
 
652
        else
 
653
                errno = 0;
 
654
        return r;
 
655
}
 
656
 
 
657
PyDoc_STRVAR(c_sinh_doc,
 
658
"sinh(x)\n"
 
659
"\n"
 
660
"Return the hyperbolic sine of x.");
 
661
 
 
662
 
 
663
static Py_complex sqrt_special_values[7][7];
 
664
 
 
665
static Py_complex
 
666
c_sqrt(Py_complex z)
 
667
{
 
668
        /*
 
669
           Method: use symmetries to reduce to the case when x = z.real and y
 
670
           = z.imag are nonnegative.  Then the real part of the result is
 
671
           given by
 
672
 
 
673
             s = sqrt((x + hypot(x, y))/2)
 
674
 
 
675
           and the imaginary part is
 
676
 
 
677
             d = (y/2)/s
 
678
 
 
679
           If either x or y is very large then there's a risk of overflow in
 
680
           computation of the expression x + hypot(x, y).  We can avoid this
 
681
           by rewriting the formula for s as:
 
682
 
 
683
             s = 2*sqrt(x/8 + hypot(x/8, y/8))
 
684
 
 
685
           This costs us two extra multiplications/divisions, but avoids the
 
686
           overhead of checking for x and y large.
 
687
 
 
688
           If both x and y are subnormal then hypot(x, y) may also be
 
689
           subnormal, so will lack full precision.  We solve this by rescaling
 
690
           x and y by a sufficiently large power of 2 to ensure that x and y
 
691
           are normal.
 
692
        */
 
693
 
 
694
 
 
695
        Py_complex r;
 
696
        double s,d;
 
697
        double ax, ay;
 
698
 
 
699
        SPECIAL_VALUE(z, sqrt_special_values);
 
700
 
 
701
        if (z.real == 0. && z.imag == 0.) {
 
702
                r.real = 0.;
 
703
                r.imag = z.imag;
 
704
                return r;
 
705
        }
 
706
 
 
707
        ax = fabs(z.real);
 
708
        ay = fabs(z.imag);
 
709
 
 
710
        if (ax < DBL_MIN && ay < DBL_MIN && (ax > 0. || ay > 0.)) {
 
711
                /* here we catch cases where hypot(ax, ay) is subnormal */
 
712
                ax = ldexp(ax, CM_SCALE_UP);
 
713
                s = ldexp(sqrt(ax + hypot(ax, ldexp(ay, CM_SCALE_UP))),
 
714
                          CM_SCALE_DOWN);
 
715
        } else {
 
716
                ax /= 8.;
 
717
                s = 2.*sqrt(ax + hypot(ax, ay/8.));
 
718
        }
 
719
        d = ay/(2.*s);
 
720
 
 
721
        if (z.real >= 0.) {
 
722
                r.real = s;
 
723
                r.imag = copysign(d, z.imag);
 
724
        } else {
 
725
                r.real = d;
 
726
                r.imag = copysign(s, z.imag);
 
727
        }
 
728
        errno = 0;
 
729
        return r;
 
730
}
 
731
 
 
732
PyDoc_STRVAR(c_sqrt_doc,
 
733
"sqrt(x)\n"
 
734
"\n"
 
735
"Return the square root of x.");
 
736
 
 
737
 
 
738
static Py_complex
 
739
c_tan(Py_complex z)
 
740
{
 
741
        /* tan(z) = -i tanh(iz) */
 
742
        Py_complex s, r;
 
743
        s.real = -z.imag;
 
744
        s.imag = z.real;
 
745
        s = c_tanh(s);
 
746
        r.real = s.imag;
 
747
        r.imag = -s.real;
 
748
        return r;
 
749
}
 
750
 
 
751
PyDoc_STRVAR(c_tan_doc,
 
752
"tan(x)\n"
 
753
"\n"
 
754
"Return the tangent of x.");
 
755
 
 
756
 
 
757
/* tanh(infinity + i*y) needs to be dealt with specially */
 
758
static Py_complex tanh_special_values[7][7];
 
759
 
 
760
static Py_complex
 
761
c_tanh(Py_complex z)
 
762
{
 
763
        /* Formula:
 
764
 
 
765
           tanh(x+iy) = (tanh(x)(1+tan(y)^2) + i tan(y)(1-tanh(x))^2) /
 
766
           (1+tan(y)^2 tanh(x)^2)
 
767
 
 
768
           To avoid excessive roundoff error, 1-tanh(x)^2 is better computed
 
769
           as 1/cosh(x)^2.  When abs(x) is large, we approximate 1-tanh(x)^2
 
770
           by 4 exp(-2*x) instead, to avoid possible overflow in the
 
771
           computation of cosh(x).
 
772
 
 
773
        */
 
774
 
 
775
        Py_complex r;
 
776
        double tx, ty, cx, txty, denom;
 
777
 
 
778
        /* special treatment for tanh(+/-inf + iy) if y is finite and
 
779
           nonzero */
 
780
        if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
 
781
                if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
 
782
                    && (z.imag != 0.)) {
 
783
                        if (z.real > 0) {
 
784
                                r.real = 1.0;
 
785
                                r.imag = copysign(0.,
 
786
                                                  2.*sin(z.imag)*cos(z.imag));
 
787
                        }
 
788
                        else {
 
789
                                r.real = -1.0;
 
790
                                r.imag = copysign(0.,
 
791
                                                  2.*sin(z.imag)*cos(z.imag));
 
792
                        }
 
793
                }
 
794
                else {
 
795
                        r = tanh_special_values[special_type(z.real)]
 
796
                                               [special_type(z.imag)];
 
797
                }
 
798
                /* need to set errno = EDOM if z.imag is +/-infinity and
 
799
                   z.real is finite */
 
800
                if (Py_IS_INFINITY(z.imag) && Py_IS_FINITE(z.real))
 
801
                        errno = EDOM;
 
802
                else
 
803
                        errno = 0;
 
804
                return r;
 
805
        }
 
806
 
 
807
        /* danger of overflow in 2.*z.imag !*/
 
808
        if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
 
809
                r.real = copysign(1., z.real);
 
810
                r.imag = 4.*sin(z.imag)*cos(z.imag)*exp(-2.*fabs(z.real));
 
811
        } else {
 
812
                tx = tanh(z.real);
 
813
                ty = tan(z.imag);
 
814
                cx = 1./cosh(z.real);
 
815
                txty = tx*ty;
 
816
                denom = 1. + txty*txty;
 
817
                r.real = tx*(1.+ty*ty)/denom;
 
818
                r.imag = ((ty/denom)*cx)*cx;
 
819
        }
 
820
        errno = 0;
 
821
        return r;
 
822
}
 
823
 
 
824
PyDoc_STRVAR(c_tanh_doc,
 
825
"tanh(x)\n"
 
826
"\n"
 
827
"Return the hyperbolic tangent of x.");
 
828
 
 
829
 
 
830
static PyObject *
 
831
cmath_log(PyObject *self, PyObject *args)
 
832
{
 
833
        Py_complex x;
 
834
        Py_complex y;
 
835
 
 
836
        if (!PyArg_ParseTuple(args, "D|D", &x, &y))
 
837
                return NULL;
 
838
 
 
839
        errno = 0;
 
840
        PyFPE_START_PROTECT("complex function", return 0)
 
841
        x = c_log(x);
 
842
        if (PyTuple_GET_SIZE(args) == 2) {
 
843
                y = c_log(y);
 
844
                x = c_quot(x, y);
 
845
        }
 
846
        PyFPE_END_PROTECT(x)
 
847
        if (errno != 0)
 
848
                return math_error();
 
849
        return PyComplex_FromCComplex(x);
 
850
}
 
851
 
 
852
PyDoc_STRVAR(cmath_log_doc,
 
853
"log(x[, base]) -> the logarithm of x to the given base.\n\
 
854
If the base not specified, returns the natural logarithm (base e) of x.");
 
855
 
 
856
 
 
857
/* And now the glue to make them available from Python: */
 
858
 
 
859
static PyObject *
 
860
math_error(void)
 
861
{
 
862
        if (errno == EDOM)
 
863
                PyErr_SetString(PyExc_ValueError, "math domain error");
 
864
        else if (errno == ERANGE)
 
865
                PyErr_SetString(PyExc_OverflowError, "math range error");
 
866
        else    /* Unexpected math error */
 
867
                PyErr_SetFromErrno(PyExc_ValueError);
 
868
        return NULL;
 
869
}
 
870
 
 
871
static PyObject *
 
872
math_1(PyObject *args, Py_complex (*func)(Py_complex))
 
873
{
 
874
        Py_complex x,r ;
 
875
        if (!PyArg_ParseTuple(args, "D", &x))
 
876
                return NULL;
 
877
        errno = 0;
 
878
        PyFPE_START_PROTECT("complex function", return 0);
 
879
        r = (*func)(x);
 
880
        PyFPE_END_PROTECT(r);
 
881
        if (errno == EDOM) {
 
882
                PyErr_SetString(PyExc_ValueError, "math domain error");
 
883
                return NULL;
 
884
        }
 
885
        else if (errno == ERANGE) {
 
886
                PyErr_SetString(PyExc_OverflowError, "math range error");
 
887
                return NULL;
 
888
        }
 
889
        else {
 
890
                return PyComplex_FromCComplex(r);
 
891
        }
 
892
}
 
893
 
 
894
#define FUNC1(stubname, func) \
 
895
        static PyObject * stubname(PyObject *self, PyObject *args) { \
 
896
                return math_1(args, func); \
 
897
        }
 
898
 
 
899
FUNC1(cmath_acos, c_acos)
 
900
FUNC1(cmath_acosh, c_acosh)
 
901
FUNC1(cmath_asin, c_asin)
 
902
FUNC1(cmath_asinh, c_asinh)
 
903
FUNC1(cmath_atan, c_atan)
 
904
FUNC1(cmath_atanh, c_atanh)
 
905
FUNC1(cmath_cos, c_cos)
 
906
FUNC1(cmath_cosh, c_cosh)
 
907
FUNC1(cmath_exp, c_exp)
 
908
FUNC1(cmath_log10, c_log10)
 
909
FUNC1(cmath_sin, c_sin)
 
910
FUNC1(cmath_sinh, c_sinh)
 
911
FUNC1(cmath_sqrt, c_sqrt)
 
912
FUNC1(cmath_tan, c_tan)
 
913
FUNC1(cmath_tanh, c_tanh)
 
914
 
 
915
static PyObject *
 
916
cmath_phase(PyObject *self, PyObject *args)
 
917
{
 
918
        Py_complex z;
 
919
        double phi;
 
920
        if (!PyArg_ParseTuple(args, "D:phase", &z))
 
921
                return NULL;
 
922
        errno = 0;
 
923
        PyFPE_START_PROTECT("arg function", return 0)
 
924
        phi = c_atan2(z);
 
925
        PyFPE_END_PROTECT(phi)
 
926
        if (errno != 0)
 
927
                return math_error();
 
928
        else
 
929
                return PyFloat_FromDouble(phi);
 
930
}
 
931
 
 
932
PyDoc_STRVAR(cmath_phase_doc,
 
933
"phase(z) -> float\n\n\
 
934
Return argument, also known as the phase angle, of a complex.");
 
935
 
 
936
static PyObject *
 
937
cmath_polar(PyObject *self, PyObject *args)
 
938
{
 
939
        Py_complex z;
 
940
        double r, phi;
 
941
        if (!PyArg_ParseTuple(args, "D:polar", &z))
 
942
                return NULL;
 
943
        PyFPE_START_PROTECT("polar function", return 0)
 
944
        phi = c_atan2(z); /* should not cause any exception */
 
945
        r = c_abs(z); /* sets errno to ERANGE on overflow;  otherwise 0 */
 
946
        PyFPE_END_PROTECT(r)
 
947
        if (errno != 0)
 
948
                return math_error();
 
949
        else
 
950
                return Py_BuildValue("dd", r, phi);
 
951
}
 
952
 
 
953
PyDoc_STRVAR(cmath_polar_doc,
 
954
"polar(z) -> r: float, phi: float\n\n\
 
955
Convert a complex from rectangular coordinates to polar coordinates. r is\n\
 
956
the distance from 0 and phi the phase angle.");
 
957
 
 
958
/*
 
959
  rect() isn't covered by the C99 standard, but it's not too hard to
 
960
  figure out 'spirit of C99' rules for special value handing:
 
961
 
 
962
    rect(x, t) should behave like exp(log(x) + it) for positive-signed x
 
963
    rect(x, t) should behave like -exp(log(-x) + it) for negative-signed x
 
964
    rect(nan, t) should behave like exp(nan + it), except that rect(nan, 0)
 
965
      gives nan +- i0 with the sign of the imaginary part unspecified.
 
966
 
 
967
*/
 
968
 
 
969
static Py_complex rect_special_values[7][7];
 
970
 
 
971
static PyObject *
 
972
cmath_rect(PyObject *self, PyObject *args)
 
973
{
 
974
        Py_complex z;
 
975
        double r, phi;
 
976
        if (!PyArg_ParseTuple(args, "dd:rect", &r, &phi))
 
977
                return NULL;
 
978
        errno = 0;
 
979
        PyFPE_START_PROTECT("rect function", return 0)
 
980
 
 
981
        /* deal with special values */
 
982
        if (!Py_IS_FINITE(r) || !Py_IS_FINITE(phi)) {
 
983
                /* if r is +/-infinity and phi is finite but nonzero then
 
984
                   result is (+-INF +-INF i), but we need to compute cos(phi)
 
985
                   and sin(phi) to figure out the signs. */
 
986
                if (Py_IS_INFINITY(r) && (Py_IS_FINITE(phi)
 
987
                                          && (phi != 0.))) {
 
988
                        if (r > 0) {
 
989
                                z.real = copysign(INF, cos(phi));
 
990
                                z.imag = copysign(INF, sin(phi));
 
991
                        }
 
992
                        else {
 
993
                                z.real = -copysign(INF, cos(phi));
 
994
                                z.imag = -copysign(INF, sin(phi));
 
995
                        }
 
996
                }
 
997
                else {
 
998
                        z = rect_special_values[special_type(r)]
 
999
                                               [special_type(phi)];
 
1000
                }
 
1001
                /* need to set errno = EDOM if r is a nonzero number and phi
 
1002
                   is infinite */
 
1003
                if (r != 0. && !Py_IS_NAN(r) && Py_IS_INFINITY(phi))
 
1004
                        errno = EDOM;
 
1005
                else
 
1006
                        errno = 0;
 
1007
        }
 
1008
        else {
 
1009
                z.real = r * cos(phi);
 
1010
                z.imag = r * sin(phi);
 
1011
                errno = 0;
 
1012
        }
 
1013
 
 
1014
        PyFPE_END_PROTECT(z)
 
1015
        if (errno != 0)
 
1016
                return math_error();
 
1017
        else
 
1018
                return PyComplex_FromCComplex(z);
 
1019
}
 
1020
 
 
1021
PyDoc_STRVAR(cmath_rect_doc,
 
1022
"rect(r, phi) -> z: complex\n\n\
 
1023
Convert from polar coordinates to rectangular coordinates.");
 
1024
 
 
1025
static PyObject *
 
1026
cmath_isnan(PyObject *self, PyObject *args)
 
1027
{
 
1028
        Py_complex z;
 
1029
        if (!PyArg_ParseTuple(args, "D:isnan", &z))
 
1030
                return NULL;
 
1031
        return PyBool_FromLong(Py_IS_NAN(z.real) || Py_IS_NAN(z.imag));
 
1032
}
 
1033
 
 
1034
PyDoc_STRVAR(cmath_isnan_doc,
 
1035
"isnan(z) -> bool\n\
 
1036
Checks if the real or imaginary part of z not a number (NaN)");
 
1037
 
 
1038
static PyObject *
 
1039
cmath_isinf(PyObject *self, PyObject *args)
 
1040
{
 
1041
        Py_complex z;
 
1042
        if (!PyArg_ParseTuple(args, "D:isnan", &z))
 
1043
                return NULL;
 
1044
        return PyBool_FromLong(Py_IS_INFINITY(z.real) ||
 
1045
                               Py_IS_INFINITY(z.imag));
 
1046
}
 
1047
 
 
1048
PyDoc_STRVAR(cmath_isinf_doc,
 
1049
"isinf(z) -> bool\n\
 
1050
Checks if the real or imaginary part of z is infinite.");
 
1051
 
 
1052
 
 
1053
PyDoc_STRVAR(module_doc,
 
1054
"This module is always available. It provides access to mathematical\n"
 
1055
"functions for complex numbers.");
 
1056
 
 
1057
static PyMethodDef cmath_methods[] = {
 
1058
        {"acos",   cmath_acos,  METH_VARARGS, c_acos_doc},
 
1059
        {"acosh",  cmath_acosh, METH_VARARGS, c_acosh_doc},
 
1060
        {"asin",   cmath_asin,  METH_VARARGS, c_asin_doc},
 
1061
        {"asinh",  cmath_asinh, METH_VARARGS, c_asinh_doc},
 
1062
        {"atan",   cmath_atan,  METH_VARARGS, c_atan_doc},
 
1063
        {"atanh",  cmath_atanh, METH_VARARGS, c_atanh_doc},
 
1064
        {"cos",    cmath_cos,   METH_VARARGS, c_cos_doc},
 
1065
        {"cosh",   cmath_cosh,  METH_VARARGS, c_cosh_doc},
 
1066
        {"exp",    cmath_exp,   METH_VARARGS, c_exp_doc},
 
1067
        {"isinf",  cmath_isinf, METH_VARARGS, cmath_isinf_doc},
 
1068
        {"isnan",  cmath_isnan, METH_VARARGS, cmath_isnan_doc},
 
1069
        {"log",    cmath_log,   METH_VARARGS, cmath_log_doc},
 
1070
        {"log10",  cmath_log10, METH_VARARGS, c_log10_doc},
 
1071
        {"phase",  cmath_phase, METH_VARARGS, cmath_phase_doc},
 
1072
        {"polar",  cmath_polar, METH_VARARGS, cmath_polar_doc},
 
1073
        {"rect",   cmath_rect,  METH_VARARGS, cmath_rect_doc},
 
1074
        {"sin",    cmath_sin,   METH_VARARGS, c_sin_doc},
 
1075
        {"sinh",   cmath_sinh,  METH_VARARGS, c_sinh_doc},
 
1076
        {"sqrt",   cmath_sqrt,  METH_VARARGS, c_sqrt_doc},
 
1077
        {"tan",    cmath_tan,   METH_VARARGS, c_tan_doc},
 
1078
        {"tanh",   cmath_tanh,  METH_VARARGS, c_tanh_doc},
 
1079
        {NULL,          NULL}           /* sentinel */
 
1080
};
 
1081
 
 
1082
 
 
1083
static struct PyModuleDef cmathmodule = {
 
1084
        PyModuleDef_HEAD_INIT,
 
1085
        "cmath",
 
1086
        module_doc,
 
1087
        -1,
 
1088
        cmath_methods,
 
1089
        NULL,
 
1090
        NULL,
 
1091
        NULL,
 
1092
        NULL
 
1093
};
 
1094
 
 
1095
PyMODINIT_FUNC
 
1096
PyInit_cmath(void)
 
1097
{
 
1098
        PyObject *m;
 
1099
 
 
1100
        m = PyModule_Create(&cmathmodule);
 
1101
        if (m == NULL)
 
1102
                return NULL;
 
1103
 
 
1104
        PyModule_AddObject(m, "pi",
 
1105
                           PyFloat_FromDouble(Py_MATH_PI));
 
1106
        PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
 
1107
 
 
1108
        /* initialize special value tables */
 
1109
 
 
1110
#define INIT_SPECIAL_VALUES(NAME, BODY) { Py_complex* p = (Py_complex*)NAME; BODY }
 
1111
#define C(REAL, IMAG) p->real = REAL; p->imag = IMAG; ++p;
 
1112
 
 
1113
        INIT_SPECIAL_VALUES(acos_special_values, {
 
1114
          C(P34,INF) C(P,INF)  C(P,INF)  C(P,-INF)  C(P,-INF)  C(P34,-INF) C(N,INF)
 
1115
          C(P12,INF) C(U,U)    C(U,U)    C(U,U)     C(U,U)     C(P12,-INF) C(N,N)
 
1116
          C(P12,INF) C(U,U)    C(P12,0.) C(P12,-0.) C(U,U)     C(P12,-INF) C(P12,N)
 
1117
          C(P12,INF) C(U,U)    C(P12,0.) C(P12,-0.) C(U,U)     C(P12,-INF) C(P12,N)
 
1118
          C(P12,INF) C(U,U)    C(U,U)    C(U,U)     C(U,U)     C(P12,-INF) C(N,N)
 
1119
          C(P14,INF) C(0.,INF) C(0.,INF) C(0.,-INF) C(0.,-INF) C(P14,-INF) C(N,INF)
 
1120
          C(N,INF)   C(N,N)    C(N,N)    C(N,N)     C(N,N)     C(N,-INF)   C(N,N)
 
1121
        })
 
1122
 
 
1123
        INIT_SPECIAL_VALUES(acosh_special_values, {
 
1124
          C(INF,-P34) C(INF,-P)  C(INF,-P)  C(INF,P)  C(INF,P)  C(INF,P34) C(INF,N)
 
1125
          C(INF,-P12) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,P12) C(N,N)
 
1126
          C(INF,-P12) C(U,U)     C(0.,-P12) C(0.,P12) C(U,U)    C(INF,P12) C(N,N)
 
1127
          C(INF,-P12) C(U,U)     C(0.,-P12) C(0.,P12) C(U,U)    C(INF,P12) C(N,N)
 
1128
          C(INF,-P12) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,P12) C(N,N)
 
1129
          C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)
 
1130
          C(INF,N)    C(N,N)     C(N,N)     C(N,N)    C(N,N)    C(INF,N)   C(N,N)
 
1131
        })
 
1132
 
 
1133
        INIT_SPECIAL_VALUES(asinh_special_values, {
 
1134
          C(-INF,-P14) C(-INF,-0.) C(-INF,-0.) C(-INF,0.) C(-INF,0.) C(-INF,P14) C(-INF,N)
 
1135
          C(-INF,-P12) C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(-INF,P12) C(N,N)
 
1136
          C(-INF,-P12) C(U,U)      C(-0.,-0.)  C(-0.,0.)  C(U,U)     C(-INF,P12) C(N,N)
 
1137
          C(INF,-P12)  C(U,U)      C(0.,-0.)   C(0.,0.)   C(U,U)     C(INF,P12)  C(N,N)
 
1138
          C(INF,-P12)  C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(INF,P12)  C(N,N)
 
1139
          C(INF,-P14)  C(INF,-0.)  C(INF,-0.)  C(INF,0.)  C(INF,0.)  C(INF,P14)  C(INF,N)
 
1140
          C(INF,N)     C(N,N)      C(N,-0.)    C(N,0.)    C(N,N)     C(INF,N)    C(N,N)
 
1141
        })
 
1142
 
 
1143
        INIT_SPECIAL_VALUES(atanh_special_values, {
 
1144
          C(-0.,-P12) C(-0.,-P12) C(-0.,-P12) C(-0.,P12) C(-0.,P12) C(-0.,P12) C(-0.,N)
 
1145
          C(-0.,-P12) C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(-0.,P12) C(N,N)
 
1146
          C(-0.,-P12) C(U,U)      C(-0.,-0.)  C(-0.,0.)  C(U,U)     C(-0.,P12) C(-0.,N)
 
1147
          C(0.,-P12)  C(U,U)      C(0.,-0.)   C(0.,0.)   C(U,U)     C(0.,P12)  C(0.,N)
 
1148
          C(0.,-P12)  C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(0.,P12)  C(N,N)
 
1149
          C(0.,-P12)  C(0.,-P12)  C(0.,-P12)  C(0.,P12)  C(0.,P12)  C(0.,P12)  C(0.,N)
 
1150
          C(0.,-P12)  C(N,N)      C(N,N)      C(N,N)     C(N,N)     C(0.,P12)  C(N,N)
 
1151
        })
 
1152
 
 
1153
        INIT_SPECIAL_VALUES(cosh_special_values, {
 
1154
          C(INF,N) C(U,U) C(INF,0.)  C(INF,-0.) C(U,U) C(INF,N) C(INF,N)
 
1155
          C(N,N)   C(U,U) C(U,U)     C(U,U)     C(U,U) C(N,N)   C(N,N)
 
1156
          C(N,0.)  C(U,U) C(1.,0.)   C(1.,-0.)  C(U,U) C(N,0.)  C(N,0.)
 
1157
          C(N,0.)  C(U,U) C(1.,-0.)  C(1.,0.)   C(U,U) C(N,0.)  C(N,0.)
 
1158
          C(N,N)   C(U,U) C(U,U)     C(U,U)     C(U,U) C(N,N)   C(N,N)
 
1159
          C(INF,N) C(U,U) C(INF,-0.) C(INF,0.)  C(U,U) C(INF,N) C(INF,N)
 
1160
          C(N,N)   C(N,N) C(N,0.)    C(N,0.)    C(N,N) C(N,N)   C(N,N)
 
1161
        })
 
1162
 
 
1163
        INIT_SPECIAL_VALUES(exp_special_values, {
 
1164
          C(0.,0.) C(U,U) C(0.,-0.)  C(0.,0.)  C(U,U) C(0.,0.) C(0.,0.)
 
1165
          C(N,N)   C(U,U) C(U,U)     C(U,U)    C(U,U) C(N,N)   C(N,N)
 
1166
          C(N,N)   C(U,U) C(1.,-0.)  C(1.,0.)  C(U,U) C(N,N)   C(N,N)
 
1167
          C(N,N)   C(U,U) C(1.,-0.)  C(1.,0.)  C(U,U) C(N,N)   C(N,N)
 
1168
          C(N,N)   C(U,U) C(U,U)     C(U,U)    C(U,U) C(N,N)   C(N,N)
 
1169
          C(INF,N) C(U,U) C(INF,-0.) C(INF,0.) C(U,U) C(INF,N) C(INF,N)
 
1170
          C(N,N)   C(N,N) C(N,-0.)   C(N,0.)   C(N,N) C(N,N)   C(N,N)
 
1171
        })
 
1172
 
 
1173
        INIT_SPECIAL_VALUES(log_special_values, {
 
1174
          C(INF,-P34) C(INF,-P)  C(INF,-P)   C(INF,P)   C(INF,P)  C(INF,P34)  C(INF,N)
 
1175
          C(INF,-P12) C(U,U)     C(U,U)      C(U,U)     C(U,U)    C(INF,P12)  C(N,N)
 
1176
          C(INF,-P12) C(U,U)     C(-INF,-P)  C(-INF,P)  C(U,U)    C(INF,P12)  C(N,N)
 
1177
          C(INF,-P12) C(U,U)     C(-INF,-0.) C(-INF,0.) C(U,U)    C(INF,P12)  C(N,N)
 
1178
          C(INF,-P12) C(U,U)     C(U,U)      C(U,U)     C(U,U)    C(INF,P12)  C(N,N)
 
1179
          C(INF,-P14) C(INF,-0.) C(INF,-0.)  C(INF,0.)  C(INF,0.) C(INF,P14)  C(INF,N)
 
1180
          C(INF,N)    C(N,N)     C(N,N)      C(N,N)     C(N,N)    C(INF,N)    C(N,N)
 
1181
        })
 
1182
 
 
1183
        INIT_SPECIAL_VALUES(sinh_special_values, {
 
1184
          C(INF,N) C(U,U) C(-INF,-0.) C(-INF,0.) C(U,U) C(INF,N) C(INF,N)
 
1185
          C(N,N)   C(U,U) C(U,U)      C(U,U)     C(U,U) C(N,N)   C(N,N)
 
1186
          C(0.,N)  C(U,U) C(-0.,-0.)  C(-0.,0.)  C(U,U) C(0.,N)  C(0.,N)
 
1187
          C(0.,N)  C(U,U) C(0.,-0.)   C(0.,0.)   C(U,U) C(0.,N)  C(0.,N)
 
1188
          C(N,N)   C(U,U) C(U,U)      C(U,U)     C(U,U) C(N,N)   C(N,N)
 
1189
          C(INF,N) C(U,U) C(INF,-0.)  C(INF,0.)  C(U,U) C(INF,N) C(INF,N)
 
1190
          C(N,N)   C(N,N) C(N,-0.)    C(N,0.)    C(N,N) C(N,N)   C(N,N)
 
1191
        })
 
1192
 
 
1193
        INIT_SPECIAL_VALUES(sqrt_special_values, {
 
1194
          C(INF,-INF) C(0.,-INF) C(0.,-INF) C(0.,INF) C(0.,INF) C(INF,INF) C(N,INF)
 
1195
          C(INF,-INF) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,INF) C(N,N)
 
1196
          C(INF,-INF) C(U,U)     C(0.,-0.)  C(0.,0.)  C(U,U)    C(INF,INF) C(N,N)
 
1197
          C(INF,-INF) C(U,U)     C(0.,-0.)  C(0.,0.)  C(U,U)    C(INF,INF) C(N,N)
 
1198
          C(INF,-INF) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,INF) C(N,N)
 
1199
          C(INF,-INF) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,INF) C(INF,N)
 
1200
          C(INF,-INF) C(N,N)     C(N,N)     C(N,N)    C(N,N)    C(INF,INF) C(N,N)
 
1201
        })
 
1202
 
 
1203
        INIT_SPECIAL_VALUES(tanh_special_values, {
 
1204
          C(-1.,0.) C(U,U) C(-1.,-0.) C(-1.,0.) C(U,U) C(-1.,0.) C(-1.,0.)
 
1205
          C(N,N)    C(U,U) C(U,U)     C(U,U)    C(U,U) C(N,N)    C(N,N)
 
1206
          C(N,N)    C(U,U) C(-0.,-0.) C(-0.,0.) C(U,U) C(N,N)    C(N,N)
 
1207
          C(N,N)    C(U,U) C(0.,-0.)  C(0.,0.)  C(U,U) C(N,N)    C(N,N)
 
1208
          C(N,N)    C(U,U) C(U,U)     C(U,U)    C(U,U) C(N,N)    C(N,N)
 
1209
          C(1.,0.)  C(U,U) C(1.,-0.)  C(1.,0.)  C(U,U) C(1.,0.)  C(1.,0.)
 
1210
          C(N,N)    C(N,N) C(N,-0.)   C(N,0.)   C(N,N) C(N,N)    C(N,N)
 
1211
        })
 
1212
 
 
1213
        INIT_SPECIAL_VALUES(rect_special_values, {
 
1214
          C(INF,N) C(U,U) C(-INF,0.) C(-INF,-0.) C(U,U) C(INF,N) C(INF,N)
 
1215
          C(N,N)   C(U,U) C(U,U)     C(U,U)      C(U,U) C(N,N)   C(N,N)
 
1216
          C(0.,0.) C(U,U) C(-0.,0.)  C(-0.,-0.)  C(U,U) C(0.,0.) C(0.,0.)
 
1217
          C(0.,0.) C(U,U) C(0.,-0.)  C(0.,0.)    C(U,U) C(0.,0.) C(0.,0.)
 
1218
          C(N,N)   C(U,U) C(U,U)     C(U,U)      C(U,U) C(N,N)   C(N,N)
 
1219
          C(INF,N) C(U,U) C(INF,-0.) C(INF,0.)   C(U,U) C(INF,N) C(INF,N)
 
1220
          C(N,N)   C(N,N) C(N,0.)    C(N,0.)     C(N,N) C(N,N)   C(N,N)
 
1221
        })
 
1222
        return m;
 
1223
}