← Back to branch summary

~ubuntu-branches/ubuntu/natty/python3.1/natty-security

« back to all changes in this revision

Viewing changes to Modules/cmathmodule.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2010-07-06 16:52:42 UTC
  • mfrom: (1.2.1 upstream) (2.1.11 sid)
  • Revision ID: james.westby@ubuntu.com-20100706165242-2xv4i019r3et6c0j
Tags: 3.1.2+20100706-1ubuntu1
* Merge with Debian; remaining changes:
  - Regenerate the control file.
  - Add debian/patches/overwrite-semaphore-check for Lucid buildds.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Modules/cmathmodule.c
31
31
#define CM_LOG_LARGE_DOUBLE (log(CM_LARGE_DOUBLE))
32
32
#define CM_SQRT_DBL_MIN (sqrt(DBL_MIN))
33
33
 
34
 
/* 
 
34
/*
35
35
   CM_SCALE_UP is an odd integer chosen such that multiplication by
36
36
   2**CM_SCALE_UP is sufficient to turn a subnormal into a normal.
37
37
   CM_SCALE_DOWN is (-(CM_SCALE_UP+1)/2).  These scalings are used to compute
62
62
*/
63
63
 
64
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 */
 
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
72
};
73
73
 
74
74
static enum special_types
75
75
special_type(double d)
76
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;
 
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
97
}
98
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
 
        }
 
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
105
 
106
106
#define P Py_MATH_PI
107
107
#define P14 0.25*Py_MATH_PI
125
125
static Py_complex
126
126
c_acos(Py_complex z)
127
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;
 
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
156
}
157
157
 
158
158
PyDoc_STRVAR(c_acos_doc,
166
166
static Py_complex
167
167
c_acosh(Py_complex z)
168
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;
 
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
189
}
190
190
 
191
191
PyDoc_STRVAR(c_acosh_doc,
197
197
static Py_complex
198
198
c_asin(Py_complex z)
199
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;
 
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
208
}
209
209
 
210
210
PyDoc_STRVAR(c_asin_doc,
218
218
static Py_complex
219
219
c_asinh(Py_complex z)
220
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;
 
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
246
}
247
247
 
248
248
PyDoc_STRVAR(c_asinh_doc,
254
254
static Py_complex
255
255
c_atan(Py_complex z)
256
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;
 
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
265
}
266
266
 
267
267
/* Windows screws up atan2 for inf and nan, and alpha Tru64 5.1 doesn't follow
269
269
static double
270
270
c_atan2(Py_complex z)
271
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);
 
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
295
}
296
296
 
297
297
PyDoc_STRVAR(c_atan_doc,
305
305
static Py_complex
306
306
c_atanh(Py_complex z)
307
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;
 
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
350
}
351
351
 
352
352
PyDoc_STRVAR(c_atanh_doc,
358
358
static Py_complex
359
359
c_cos(Py_complex z)
360
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;
 
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
367
}
368
368
 
369
369
PyDoc_STRVAR(c_cos_doc,
378
378
static Py_complex
379
379
c_cosh(Py_complex z)
380
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;
 
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
426
}
427
427
 
428
428
PyDoc_STRVAR(c_cosh_doc,
438
438
static Py_complex
439
439
c_exp(Py_complex z)
440
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;
 
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
486
}
487
487
 
488
488
PyDoc_STRVAR(c_exp_doc,
496
496
static Py_complex
497
497
c_log(Py_complex z)
498
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;
 
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
563
}
564
564
 
565
565
 
566
566
static Py_complex
567
567
c_log10(Py_complex z)
568
568
{
569
 
        Py_complex r;
570
 
        int errno_save;
 
569
    Py_complex r;
 
570
    int errno_save;
571
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;
 
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
578
}
579
579
 
580
580
PyDoc_STRVAR(c_log10_doc,
586
586
static Py_complex
587
587
c_sin(Py_complex z)
588
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;
 
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
597
}
598
598
 
599
599
PyDoc_STRVAR(c_sin_doc,
608
608
static Py_complex
609
609
c_sinh(Py_complex z)
610
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;
 
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
655
}
656
656
 
657
657
PyDoc_STRVAR(c_sinh_doc,
665
665
static Py_complex
666
666
c_sqrt(Py_complex z)
667
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;
 
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
730
}
731
731
 
732
732
PyDoc_STRVAR(c_sqrt_doc,
738
738
static Py_complex
739
739
c_tan(Py_complex z)
740
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;
 
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
749
}
750
750
 
751
751
PyDoc_STRVAR(c_tan_doc,
760
760
static Py_complex
761
761
c_tanh(Py_complex z)
762
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;
 
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
822
}
823
823
 
824
824
PyDoc_STRVAR(c_tanh_doc,
830
830
static PyObject *
831
831
cmath_log(PyObject *self, PyObject *args)
832
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);
 
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
850
}
851
851
 
852
852
PyDoc_STRVAR(cmath_log_doc,
859
859
static PyObject *
860
860
math_error(void)
861
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;
 
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
869
}
870
870
 
871
871
static PyObject *
872
872
math_1(PyObject *args, Py_complex (*func)(Py_complex))
873
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
 
        }
 
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
892
}
893
893
 
894
894
#define FUNC1(stubname, func) \
895
 
        static PyObject * stubname(PyObject *self, PyObject *args) { \
896
 
                return math_1(args, func); \
897
 
        }
 
895
    static PyObject * stubname(PyObject *self, PyObject *args) { \
 
896
        return math_1(args, func); \
 
897
    }
898
898
 
899
899
FUNC1(cmath_acos, c_acos)
900
900
FUNC1(cmath_acosh, c_acosh)
915
915
static PyObject *
916
916
cmath_phase(PyObject *self, PyObject *args)
917
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);
 
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
930
}
931
931
 
932
932
PyDoc_STRVAR(cmath_phase_doc,
936
936
static PyObject *
937
937
cmath_polar(PyObject *self, PyObject *args)
938
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);
 
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
951
}
952
952
 
953
953
PyDoc_STRVAR(cmath_polar_doc,
971
971
static PyObject *
972
972
cmath_rect(PyObject *self, PyObject *args)
973
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);
 
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
1019
}
1020
1020
 
1021
1021
PyDoc_STRVAR(cmath_rect_doc,
1025
1025
static PyObject *
1026
1026
cmath_isnan(PyObject *self, PyObject *args)
1027
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));
 
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
1032
}
1033
1033
 
1034
1034
PyDoc_STRVAR(cmath_isnan_doc,
1038
1038
static PyObject *
1039
1039
cmath_isinf(PyObject *self, PyObject *args)
1040
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));
 
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
1046
}
1047
1047
 
1048
1048
PyDoc_STRVAR(cmath_isinf_doc,
1055
1055
"functions for complex numbers.");
1056
1056
 
1057
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 */
 
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
1080
};
1081
1081
 
1082
1082
 
1083
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
 
1084
    PyModuleDef_HEAD_INIT,
 
1085
    "cmath",
 
1086
    module_doc,
 
1087
    -1,
 
1088
    cmath_methods,
 
1089
    NULL,
 
1090
    NULL,
 
1091
    NULL,
 
1092
    NULL
1093
1093
};
1094
1094
 
1095
1095
PyMODINIT_FUNC
1096
1096
PyInit_cmath(void)
1097
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 */
 
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
1109
 
1110
1110
#define INIT_SPECIAL_VALUES(NAME, BODY) { Py_complex* p = (Py_complex*)NAME; BODY }
1111
1111
#define C(REAL, IMAG) p->real = REAL; p->imag = IMAG; ++p;
1112
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;
 
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
1223
}