1
/*
2
 * vim: syntax=c
3
 *
4
 * Implement some C99-compatible complex math functions
5
 *
6
 * Most of the code is taken from the msun library in FreeBSD (HEAD @ 4th
7
 * October 2013), under the following license:
8
 *
9
 * Copyright (c) 2007, 2011 David Schultz <das@FreeBSD.ORG>
10
 * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
11
 * All rights reserved.
12
 *
13
 * Redistribution and use in source and binary forms, with or without
14
 * modification, are permitted provided that the following conditions
15
 * are met:
16
 * 1. Redistributions of source code must retain the above copyright
17
 *    notice, this list of conditions and the following disclaimer.
18
 * 2. Redistributions in binary form must reproduce the above copyright
19
 *    notice, this list of conditions and the following disclaimer in the
20
 *    documentation and/or other materials provided with the distribution.
21
 *
22
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
23
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
24
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
25
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
26
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
27
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
28
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
29
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
30
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
31
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
32
 * SUCH DAMAGE.
33
 */
34
#include "npy_math_common.h"
35
#include "npy_math_private.h"
36
#include <numpy/utils.h>
37

38
/*
39
 * Hack inherited from BSD, the intent is to set the FPU inexact
40
 * flag in an efficient way. The flag is IEEE specific. See
41
 * https://github.com/freebsd/freebsd/blob/4c6378299/lib/msun/src/catrig.c#L42
42
 */
43
#if !defined(HAVE_CACOSF) || !defined(HAVE_CACOSL) || !defined(HAVE_CASINHF) || !defined(HAVE_CASINHL)
44
#define raise_inexact() do {                        \
45
    volatile npy_float NPY_UNUSED(junk) = 1 + tiny; \
46
} while (0)
47

48

49
static const volatile npy_float tiny = 3.9443045e-31f;
50
#endif
51

52
/**begin repeat
53
 * #type = npy_float, npy_double, npy_longdouble#
54
 * #ctype = npy_cfloat,npy_cdouble,npy_clongdouble#
55
 * #c = f, , l#
56
 * #C = F, , L#
57
 * #TMAX = FLT_MAX, DBL_MAX, LDBL_MAX#
58
 * #TMIN = FLT_MIN, DBL_MIN, LDBL_MIN#
59
 * #TMANT_DIG = FLT_MANT_DIG, DBL_MANT_DIG, LDBL_MANT_DIG#
60
 * #TEPS = FLT_EPSILON, DBL_EPSILON, LDBL_EPSILON#
61
 * #precision = 1, 2, 3#
62
 */
63

64
/*==========================================================
65
 * Constants
66
 *=========================================================*/
67
static const @ctype@ c_1@c@ = {1.0@C@, 0.0};
68

69
/*==========================================================
70
 * Helper functions
71
 *
72
 * These are necessary because we do not count on using a
73
 * C99 compiler.
74
 *=========================================================*/
75
static NPY_INLINE
76
@ctype@
77
cmul@c@(@ctype@ a, @ctype@ b)
78
{
79
    @type@ ar, ai, br, bi;
80 1
    ar = npy_creal@c@(a);
81 1
    ai = npy_cimag@c@(a);
82 1
    br = npy_creal@c@(b);
83 1
    bi = npy_cimag@c@(b);
84 1
    return npy_cpack@c@(ar*br - ai*bi, ar*bi + ai*br);
85
}
86

87
static NPY_INLINE
88
@ctype@
89 1
cdiv@c@(@ctype@ a, @ctype@ b)
90
{
91
    @type@ ar, ai, br, bi, abs_br, abs_bi;
92 1
    ar = npy_creal@c@(a);
93 1
    ai = npy_cimag@c@(a);
94 1
    br = npy_creal@c@(b);
95 1
    bi = npy_cimag@c@(b);
96 1
    abs_br = npy_fabs@c@(br);
97 1
    abs_bi = npy_fabs@c@(bi);
98

99 1
    if (abs_br >= abs_bi) {
100 1
        if (abs_br == 0 && abs_bi == 0) {
101
            /* divide by zeros should yield a complex inf or nan */
102 0
            return npy_cpack@c@(ar/abs_br, ai/abs_bi);
103
        }
104
        else {
105 1
            @type@ rat = bi/br;
106 1
            @type@ scl = 1.0@C@/(br+bi*rat);
107 1
            return npy_cpack@c@((ar + ai*rat)*scl, (ai - ar*rat)*scl);
108
        }
109
    }
110
    else {
111 1
        @type@ rat = br/bi;
112 1
        @type@ scl = 1.0@C@/(bi + br*rat);
113 1
        return npy_cpack@c@((ar*rat + ai)*scl, (ai*rat - ar)*scl);
114
    }
115
}
116

117
/*==========================================================
118
 * Custom implementation of missing complex C99 functions
119
 *=========================================================*/
120

121
#ifndef HAVE_CABS@C@
122
@type@
123
npy_cabs@c@(@ctype@ z)
124
{
125
    return npy_hypot@c@(npy_creal@c@(z), npy_cimag@c@(z));
126
}
127
#endif
128

129
#ifndef HAVE_CARG@C@
130
@type@
131
npy_carg@c@(@ctype@ z)
132
{
133
    return npy_atan2@c@(npy_cimag@c@(z), npy_creal@c@(z));
134
}
135
#endif
136

137
/*
138
 * cexp and (ccos, csin)h functions need to calculate exp scaled by another
139
 * number.  This can be difficult if exp(x) overflows.  By doing this way, we
140
 * don't risk overflowing exp. This likely raises floating-point exceptions,
141
 * if we decide that we care.
142
 *
143
 * This is only useful over a limited range, (see below) an expects that the
144
 * input values are in this range.
145
 *
146
 * This is based on the technique used in FreeBSD's __frexp_exp and
147
 * __ldexp_(c)exp functions by David Schultz.
148
 *
149
 * SCALED_CEXP_LOWER = log(FLT_MAX)
150
 * SCALED_CEXP_UPPER = log(2) + log(FLT_MAX) - log(FLT_TRUE_MIN),
151
 * where FLT_TRUE_MIN is the smallest possible subnormal number.
152
 */
153

154
#define SCALED_CEXP_LOWERF 88.722839f
155
#define SCALED_CEXP_UPPERF 192.69492f
156
#define SCALED_CEXP_LOWER 710.47586007394386
157
#define SCALED_CEXP_UPPER 1454.9159319953251
158
#define SCALED_CEXP_LOWERL 11357.216553474703895L
159
#define SCALED_CEXP_UPPERL 22756.021937783004509L
160

161
#if !defined(HAVE_CSINH@C@) || \
162
    !defined(HAVE_CCOSH@C@) || \
163
    !defined(HAVE_CEXP@C@)
164

165
static
166
@ctype@
167
_npy_scaled_cexp@c@(@type@ x, @type@ y, npy_int expt)
168
{
169
#if @precision@ == 1
170
    const npy_int k = 235;
171
#endif
172
#if @precision@ == 2
173
    const npy_int k = 1799;
174
#endif
175
#if @precision@ == 3
176
    const npy_int k = 19547;
177
#endif
178
    const @type@ kln2 = k * NPY_LOGE2@c@;
179
    @type@ mant, mantcos, mantsin;
180
    npy_int ex, excos, exsin;
181

182
    mant = npy_frexp@c@(npy_exp@c@(x - kln2), &ex);
183
    mantcos = npy_frexp@c@(npy_cos@c@(y), &excos);
184
    mantsin = npy_frexp@c@(npy_sin@c@(y), &exsin);
185

186
    expt += ex + k;
187
    return npy_cpack@c@( npy_ldexp@c@(mant * mantcos, expt + excos),
188
                         npy_ldexp@c@(mant * mantsin, expt + exsin));
189
}
190

191
#endif
192

193
#ifndef HAVE_CEXP@C@
194

195
@ctype@
196
npy_cexp@c@(@ctype@ z)
197
{
198
    @type@ x, c, s;
199
    @type@ r, i;
200
    @ctype@ ret;
201

202
    r = npy_creal@c@(z);
203
    i = npy_cimag@c@(z);
204

205
    if (npy_isfinite(r)) {
206
        if (r >= SCALED_CEXP_LOWER@C@ && r <= SCALED_CEXP_UPPER@C@) {
207
            ret = _npy_scaled_cexp@c@(r, i, 0);
208
        }
209
        else {
210
            x = npy_exp@c@(r);
211

212
            c = npy_cos@c@(i);
213
            s = npy_sin@c@(i);
214

215
            if (npy_isfinite(i)) {
216
                ret = npy_cpack@c@(x * c, x * s);
217
            }
218
            else {
219
                ret = npy_cpack@c@(NPY_NAN@C@, npy_copysign@c@(NPY_NAN@C@, i));
220
            }
221
        }
222

223
    }
224
    else  if (npy_isnan(r)) {
225
        /* r is nan */
226
        if (i == 0) {
227
            ret = z;
228
        }
229
        else {
230
            ret = npy_cpack@c@(r, npy_copysign@c@(NPY_NAN@C@, i));
231
        }
232
    }
233
    else {
234
        /* r is +- inf */
235
        if (r > 0) {
236
            if (i == 0) {
237
                ret = npy_cpack@c@(r, i);
238
            }
239
            else if (npy_isfinite(i)) {
240
                c = npy_cos@c@(i);
241
                s = npy_sin@c@(i);
242

243
                ret = npy_cpack@c@(r * c, r * s);
244
            }
245
            else {
246
                /* x = +inf, y = +-inf | nan */
247
                npy_set_floatstatus_invalid();
248
                ret = npy_cpack@c@(r, NPY_NAN@C@);
249
            }
250
        }
251
        else {
252
            if (npy_isfinite(i)) {
253
                x = npy_exp@c@(r);
254
                c = npy_cos@c@(i);
255
                s = npy_sin@c@(i);
256

257
                ret = npy_cpack@c@(x * c, x * s);
258
            }
259
            else {
260
                /* x = -inf, y = nan | +i inf */
261
                ret = npy_cpack@c@(0, 0);
262
            }
263
        }
264
    }
265

266
    return ret;
267
}
268
#endif
269

270
#ifndef HAVE_CLOG@C@
271
/* algorithm from cpython, rev. d86f5686cef9
272
 *
273
 * The usual formula for the real part is log(hypot(z.real, z.imag)).
274
 * There are four situations where this formula is potentially
275
 * problematic:
276
 *
277
 * (1) the absolute value of z is subnormal.  Then hypot is subnormal,
278
 * so has fewer than the usual number of bits of accuracy, hence may
279
 * have large relative error.  This then gives a large absolute error
280
 * in the log.  This can be solved by rescaling z by a suitable power
281
 * of 2.
282
 *
283
 * (2) the absolute value of z is greater than DBL_MAX (e.g. when both
284
 * z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX)
285
 * Again, rescaling solves this.
286
 *
287
 * (3) the absolute value of z is close to 1.  In this case it's
288
 * difficult to achieve good accuracy, at least in part because a
289
 * change of 1ulp in the real or imaginary part of z can result in a
290
 * change of billions of ulps in the correctly rounded answer.
291
 *
292
 * (4) z = 0.  The simplest thing to do here is to call the
293
 * floating-point log with an argument of 0, and let its behaviour
294
 * (returning -infinity, signaling a floating-point exception, setting
295
 * errno, or whatever) determine that of c_log.  So the usual formula
296
 * is fine here.
297
*/
298
@ctype@
299
npy_clog@c@(@ctype@ z)
300
{
301
    @type@ ax = npy_fabs@c@(npy_creal@c@(z));
302
    @type@ ay = npy_fabs@c@(npy_cimag@c@(z));
303
    @type@ rr, ri;
304

305
    if (ax > @TMAX@/4 || ay > @TMAX@/4) {
306
        rr = npy_log@c@(npy_hypot@c@(ax/2, ay/2)) + NPY_LOGE2@c@;
307
    }
308
    else if (ax < @TMIN@ && ay < @TMIN@) {
309
        if (ax > 0  || ay > 0) {
310
            /* catch cases where hypot(ax, ay) is subnormal */
311
            rr = npy_log@c@(npy_hypot@c@(npy_ldexp@c@(ax, @TMANT_DIG@),
312
                 npy_ldexp@c@(ay, @TMANT_DIG@))) - @TMANT_DIG@*NPY_LOGE2@c@;
313
        }
314
        else {
315
            /* log(+/-0 +/- 0i) */
316
            /* raise divide-by-zero floating point exception */
317
            rr = -1.0@c@ / npy_creal@c@(z);
318
            rr = npy_copysign@c@(rr, -1);
319
            ri = npy_carg@c@(z);
320
            return npy_cpack@c@(rr, ri);
321
        }
322
    }
323
    else {
324
        @type@ h = npy_hypot@c@(ax, ay);
325
        if (0.71 <= h && h <= 1.73) {
326
            @type@ am = ax > ay ? ax : ay; /* max(ax, ay) */
327
            @type@ an = ax > ay ? ay : ax; /* min(ax, ay) */
328
            rr = npy_log1p@c@((am-1)*(am+1)+an*an)/2;
329
        }
330
        else {
331
            rr = npy_log@c@(h);
332
        }
333
    }
334
    ri = npy_carg@c@(z);
335

336
    return npy_cpack@c@(rr, ri);
337
}
338
#endif
339

340
#ifndef HAVE_CSQRT@C@
341

342
/* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */
343
#define THRESH  (@TMAX@ / (1 + NPY_SQRT2@c@))
344

345
@ctype@
346
npy_csqrt@c@(@ctype@ z)
347
{
348
    @ctype@ result;
349
    @type@ a, b;
350
    @type@ t;
351
    int scale;
352

353
    a = npy_creal@c@(z);
354
    b = npy_cimag@c@(z);
355

356
    /* Handle special cases. */
357
    if (a == 0 && b == 0) {
358
        return (npy_cpack@c@(0, b));
359
    }
360
    if (npy_isinf(b)) {
361
        return (npy_cpack@c@(NPY_INFINITY@C@, b));
362
    }
363
    if (npy_isnan(a)) {
364
        t = (b - b) / (b - b);  /* raise invalid if b is not a NaN */
365
        return (npy_cpack@c@(a, t));    /* return NaN + NaN i */
366
    }
367
    if (npy_isinf(a)) {
368
        /*
369
         * csqrt(inf + NaN i)  = inf +  NaN i
370
         * csqrt(inf + y i)    = inf +  0 i
371
         * csqrt(-inf + NaN i) = NaN +- inf i
372
         * csqrt(-inf + y i)   = 0   +  inf i
373
         */
374
        if (npy_signbit(a)) {
375
            return (npy_cpack@c@(npy_fabs@c@(b - b), npy_copysign@c@(a, b)));
376
        }
377
        else {
378
            return (npy_cpack@c@(a, npy_copysign@c@(b - b, b)));
379
        }
380
    }
381
    /*
382
     * The remaining special case (b is NaN) is handled just fine by
383
     * the normal code path below.
384
     */
385

386
    /* Scale to avoid overflow. */
387
    if (npy_fabs@c@(a) >= THRESH || npy_fabs@c@(b) >= THRESH) {
388
        a *= 0.25;
389
        b *= 0.25;
390
        scale = 1;
391
    }
392
    else {
393
        scale = 0;
394
    }
395

396
    /* Algorithm 312, CACM vol 10, Oct 1967. */
397
    if (a >= 0) {
398
        t = npy_sqrt@c@((a + npy_hypot@c@(a, b)) * 0.5@c@);
399
        result = npy_cpack@c@(t, b / (2 * t));
400
    }
401
    else {
402
        t = npy_sqrt@c@((-a + npy_hypot@c@(a, b)) * 0.5@c@);
403
        result = npy_cpack@c@(npy_fabs@c@(b) / (2 * t), npy_copysign@c@(t, b));
404
    }
405

406
    /* Rescale. */
407
    if (scale) {
408
        return (npy_cpack@c@(npy_creal@c@(result) * 2, npy_cimag@c@(result)));
409
    }
410
    else {
411
        return (result);
412
    }
413
}
414
#undef THRESH
415
#endif
416

417
/*
418
 * Always use this function because of the multiplication for small
419
 * integer powers, but in the body use cpow if it is available.
420
 */
421

422
/* private function for use in npy_pow{f, ,l} */
423
#ifdef HAVE_CPOW@C@
424
static @ctype@
425
sys_cpow@c@(@ctype@ x, @ctype@ y)
426
{
427
    __@ctype@_to_c99_cast xcast;
428
    __@ctype@_to_c99_cast ycast;
429
    __@ctype@_to_c99_cast ret;
430 1
    xcast.npy_z = x;
431 1
    ycast.npy_z = y;
432 1
    ret.c99_z = cpow@c@(xcast.c99_z, ycast.c99_z);
433 1
    return ret.npy_z;
434
}
435
#endif
436

437

438
@ctype@
439 1
npy_cpow@c@ (@ctype@ a, @ctype@ b)
440
{
441
    npy_intp n;
442 1
    @type@ ar = npy_creal@c@(a);
443 1
    @type@ br = npy_creal@c@(b);
444 1
    @type@ ai = npy_cimag@c@(a);
445 1
    @type@ bi = npy_cimag@c@(b);
446
    @ctype@ r;
447

448 1
    if (br == 0. && bi == 0.) {
449
        return npy_cpack@c@(1., 0.);
450
    }
451 1
    if (ar == 0. && ai == 0.) {
452 1
        if (br > 0 && bi == 0) {
453
            return npy_cpack@c@(0., 0.);
454
        }
455
        else {
456 1
            volatile @type@ tmp = NPY_INFINITY@C@;
457
            /*
458
             * NB: there are four complex zeros; c0 = (+-0, +-0), so that
459
             * unlike for reals, c0**p, with `p` negative is in general
460
             * ill-defined.
461
             *
462
             *     c0**z with z complex is also ill-defined.
463
             */
464 1
            r = npy_cpack@c@(NPY_NAN@C@, NPY_NAN@C@);
465

466
            /* Raise invalid */
467 1
            tmp -= NPY_INFINITY@C@;
468 1
            ar = tmp;
469 1
            return r;
470
        }
471
    }
472 1
    if (bi == 0 && (n=(npy_intp)br) == br) {
473 1
        if (n == 1) {
474
            /* unroll: handle inf better */
475
            return npy_cpack@c@(ar, ai);
476
        }
477 1
        else if (n == 2) {
478
            /* unroll: handle inf better */
479
            return cmul@c@(a, a);
480
        }
481 1
        else if (n == 3) {
482
            /* unroll: handle inf better */
483
            return cmul@c@(a, cmul@c@(a, a));
484
        }
485 1
        else if (n > -100 && n < 100) {
486
            @ctype@ p, aa;
487 1
            npy_intp mask = 1;
488 1
            if (n < 0) {
489 1
                n = -n;
490
            }
491 1
            aa = c_1@c@;
492
            p = npy_cpack@c@(ar, ai);
493
            while (1) {
494 1
                if (n & mask) {
495
                    aa = cmul@c@(aa,p);
496
                }
497 1
                mask <<= 1;
498 1
                if (n < mask || mask <= 0) {
499
                    break;
500
                }
501 1
                p = cmul@c@(p,p);
502
            }
503 1
            r = npy_cpack@c@(npy_creal@c@(aa), npy_cimag@c@(aa));
504 1
            if (br < 0) {
505 1
                r = cdiv@c@(c_1@c@, r);
506
            }
507 1
            return r;
508
        }
509
    }
510

511
#ifdef HAVE_CPOW@C@
512
    return sys_cpow@c@(a, b);
513

514
#else
515
    {
516
        @ctype@ loga = npy_clog@c@(a);
517

518
        ar = npy_creal@c@(loga);
519
        ai = npy_cimag@c@(loga);
520
        return npy_cexp@c@(npy_cpack@c@(ar*br - ai*bi, ar*bi + ai*br));
521
    }
522

523
#endif
524
}
525

526

527
#ifndef HAVE_CCOS@C@
528
@ctype@
529
npy_ccos@c@(@ctype@ z)
530
{
531
    /* ccos(z) = ccosh(I * z) */
532
    return npy_ccosh@c@(npy_cpack@c@(-npy_cimag@c@(z), npy_creal@c@(z)));
533
}
534
#endif
535

536
#ifndef HAVE_CSIN@C@
537
@ctype@
538
npy_csin@c@(@ctype@ z)
539
{
540
    /* csin(z) = -I * csinh(I * z) */
541
    z = npy_csinh@c@(npy_cpack@c@(-npy_cimag@c@(z), npy_creal@c@(z)));
542
    return npy_cpack@c@(npy_cimag@c@(z), -npy_creal@c@(z));
543
}
544
#endif
545

546
#ifndef HAVE_CTAN@C@
547
@ctype@
548
npy_ctan@c@(@ctype@ z)
549
{
550
    /* ctan(z) = -I * ctanh(I * z) */
551
    z = npy_ctanh@c@(npy_cpack@c@(-npy_cimag@c@(z), npy_creal@c@(z)));
552
    return (npy_cpack@c@(npy_cimag@c@(z), -npy_creal@c@(z)));
553
}
554
#endif
555

556
#ifndef HAVE_CCOSH@C@
557
/*
558
 * Taken from the msun library in FreeBSD, rev 226599.
559
 *
560
 * Hyperbolic cosine of a complex argument z = x + i y.
561
 *
562
 * cosh(z) = cosh(x+iy)
563
 *         = cosh(x) cos(y) + i sinh(x) sin(y).
564
 *
565
 * Exceptional values are noted in the comments within the source code.
566
 * These values and the return value were taken from n1124.pdf.
567
 *
568
 * CCOSH_BIG is chosen such that
569
 * spacing(0.5 * exp(CCOSH_BIG)) > 0.5*exp(-CCOSH_BIG)
570
 * although the exact value assigned to CCOSH_BIG is not so important
571
 */
572
@ctype@
573
npy_ccosh@c@(@ctype@ z)
574
{
575
#if @precision@ == 1
576
    const npy_float CCOSH_BIG = 9.0f;
577
    const npy_float CCOSH_HUGE = 1.70141183e+38f;
578
#endif
579
#if @precision@ == 2
580
    const npy_double CCOSH_BIG = 22.0;
581
    const npy_double CCOSH_HUGE = 8.9884656743115795e+307;
582
#endif
583
#if @precision@ >= 3
584
#if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE
585
    const npy_longdouble CCOSH_BIG = 22.0L;
586
    const npy_longdouble CCOSH_HUGE = 8.9884656743115795e+307L;
587
#else
588
    const npy_longdouble CCOSH_BIG = 24.0L;
589
    const npy_longdouble CCOSH_HUGE = 5.94865747678615882543e+4931L;
590
#endif
591
#endif
592

593
    @type@  x, y, h, absx;
594
    npy_int xfinite, yfinite;
595

596
    x = npy_creal@c@(z);
597
    y = npy_cimag@c@(z);
598

599
    xfinite = npy_isfinite(x);
600
    yfinite = npy_isfinite(y);
601

602
    /* Handle the nearly-non-exceptional cases where x and y are finite. */
603
    if (xfinite && yfinite) {
604
        if (y == 0) {
605
            return npy_cpack@c@(npy_cosh@c@(x), x * y);
606
        }
607
        absx = npy_fabs@c@(x);
608
        if (absx < CCOSH_BIG) {
609
            /* small x: normal case */
610
            return npy_cpack@c@(npy_cosh@c@(x) * npy_cos@c@(y),
611
                                npy_sinh@c@(x) * npy_sin@c@(y));
612
        }
613

614
        /* |x| >= 22, so cosh(x) ~= exp(|x|) */
615
        if (absx < SCALED_CEXP_LOWER@C@) {
616
            /* x < 710: exp(|x|) won't overflow */
617
            h = npy_exp@c@(absx) * 0.5@c@;
618
            return npy_cpack@c@(h * npy_cos@c@(y),
619
                                npy_copysign@c@(h, x) * npy_sin@c@(y));
620
        }
621
        else if (absx < SCALED_CEXP_UPPER@C@) {
622
            /* x < 1455: scale to avoid overflow */
623
            z = _npy_scaled_cexp@c@(absx, y, -1);
624
            return npy_cpack@c@(npy_creal@c@(z),
625
                                npy_cimag@c@(z) * npy_copysign@c@(1, x));
626
        }
627
        else {
628
            /* x >= 1455: the result always overflows */
629
            h = CCOSH_HUGE * x;
630
            return npy_cpack@c@(h * h * npy_cos@c@(y), h * npy_sin@c@(y));
631
        }
632
    }
633

634
    /*
635
     * cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0.
636
     * The sign of 0 in the result is unspecified.  Choice = normally
637
     * the same as dNaN.  Raise the invalid floating-point exception.
638
     *
639
     * cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0.
640
     * The sign of 0 in the result is unspecified.  Choice = normally
641
     * the same as d(NaN).
642
     */
643
    if (x == 0 && !yfinite) {
644
        return npy_cpack@c@(y - y, npy_copysign@c@(0, x * (y - y)));
645
    }
646

647
    /*
648
     * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
649
     *
650
     * cosh(NaN +- I 0)   = d(NaN) + I sign(d(NaN, +-0))0.
651
     * The sign of 0 in the result is unspecified.
652
     */
653
    if (y == 0 && !xfinite) {
654
        return npy_cpack@c@(x * x, npy_copysign@c@(0, x) * y);
655
    }
656

657
    /*
658
     * cosh(x +- I Inf) = dNaN + I dNaN.
659
     * Raise the invalid floating-point exception for finite nonzero x.
660
     *
661
     * cosh(x + I NaN) = d(NaN) + I d(NaN).
662
     * Optionally raises the invalid floating-point exception for finite
663
     * nonzero x.  Choice = don't raise (except for signaling NaNs).
664
     */
665
    if (xfinite && !yfinite) {
666
        return npy_cpack@c@(y - y, x * (y - y));
667
    }
668

669
    /*
670
     * cosh(+-Inf + I NaN)  = +Inf + I d(NaN).
671
     *
672
     * cosh(+-Inf +- I Inf) = +Inf + I dNaN.
673
     * The sign of Inf in the result is unspecified.  Choice = always +.
674
     * Raise the invalid floating-point exception.
675
     *
676
     * cosh(+-Inf + I y)   = +Inf cos(y) +- I Inf sin(y)
677
     */
678
    if (npy_isinf(x)) {
679
        if (!yfinite) {
680
            return npy_cpack@c@(x * x, x * (y - y));
681
        }
682
        return npy_cpack@c@((x * x) * npy_cos@c@(y), x * npy_sin@c@(y));
683
    }
684

685
    /*
686
     * cosh(NaN + I NaN)  = d(NaN) + I d(NaN).
687
     *
688
     * cosh(NaN +- I Inf) = d(NaN) + I d(NaN).
689
     * Optionally raises the invalid floating-point exception.
690
     * Choice = raise.
691
     *
692
     * cosh(NaN + I y)    = d(NaN) + I d(NaN).
693
     * Optionally raises the invalid floating-point exception for finite
694
     * nonzero y.  Choice = don't raise (except for signaling NaNs).
695
     */
696
    return npy_cpack@c@((x * x) * (y - y), (x + x) * (y - y));
697
}
698
#undef CCOSH_BIG
699
#undef CCOSH_HUGE
700
#endif
701

702
#ifndef HAVE_CSINH@C@
703
/*
704
 * Taken from the msun library in FreeBSD, rev 226599.
705
 *
706
 * Hyperbolic sine of a complex argument z = x + i y.
707
 *
708
 * sinh(z) = sinh(x+iy)
709
 *         = sinh(x) cos(y) + i cosh(x) sin(y).
710
 *
711
 * Exceptional values are noted in the comments within the source code.
712
 * These values and the return value were taken from n1124.pdf.
713
 */
714
@ctype@
715
npy_csinh@c@(@ctype@ z)
716
{
717
#if @precision@ == 1
718
    const npy_float CSINH_BIG = 9.0f;
719
    const npy_float CSINH_HUGE = 1.70141183e+38f;
720
#endif
721
#if @precision@ == 2
722
    const npy_double CSINH_BIG = 22.0;
723
    const npy_double CSINH_HUGE = 8.9884656743115795e+307;
724
#endif
725
#if @precision@ >= 3
726
#if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE
727
    const npy_longdouble CSINH_BIG = 22.0L;
728
    const npy_longdouble CSINH_HUGE = 8.9884656743115795e+307L;
729
#else
730
    const npy_longdouble CSINH_BIG = 24.0L;
731
    const npy_longdouble CSINH_HUGE = 5.94865747678615882543e+4931L;
732
#endif
733
#endif
734

735
    @type@ x, y, h, absx;
736
    npy_int xfinite, yfinite;
737

738
    x = npy_creal@c@(z);
739
    y = npy_cimag@c@(z);
740

741
    xfinite = npy_isfinite(x);
742
    yfinite = npy_isfinite(y);
743

744
    /* Handle the nearly-non-exceptional cases where x and y are finite. */
745
    if (xfinite && yfinite) {
746
        if (y == 0) {
747
            return npy_cpack@c@(npy_sinh@c@(x), y);
748
        }
749
        absx = npy_fabs@c@(x);
750
        if (absx < CSINH_BIG) {
751
            /* small x: normal case */
752
            return npy_cpack@c@(npy_sinh@c@(x) * npy_cos@c@(y),
753
                                npy_cosh@c@(x) * npy_sin@c@(y));
754
        }
755

756
        /* |x| >= 22, so cosh(x) ~= exp(|x|) */
757
        if (absx < SCALED_CEXP_LOWER@C@) {
758
            /* x < 710: exp(|x|) won't overflow */
759
            h = npy_exp@c@(npy_fabs@c@(x)) * 0.5@c@;
760
            return npy_cpack@c@(npy_copysign@c@(h, x) * npy_cos@c@(y),
761
                                h * npy_sin@c@(y));
762
        }
763
        else if (x < SCALED_CEXP_UPPER@C@) {
764
            /* x < 1455: scale to avoid overflow */
765
            z = _npy_scaled_cexp@c@(absx, y, -1);
766
            return npy_cpack@c@(npy_creal@c@(z) * npy_copysign@c@(1, x),
767
                                npy_cimag@c@(z));
768
        }
769
        else {
770
            /* x >= 1455: the result always overflows */
771
            h = CSINH_HUGE * x;
772
            return npy_cpack@c@(h * npy_cos@c@(y), h * h * npy_sin@c@(y));
773
        }
774
    }
775

776
    /*
777
     * sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN.
778
     * The sign of 0 in the result is unspecified.  Choice = normally
779
     * the same as dNaN.  Raise the invalid floating-point exception.
780
     *
781
     * sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN).
782
     * The sign of 0 in the result is unspecified.  Choice = normally
783
     * the same as d(NaN).
784
     */
785
    if (x == 0 && !yfinite) {
786
        return npy_cpack@c@(npy_copysign@c@(0, x * (y - y)), y - y);
787
    }
788

789
    /*
790
     * sinh(+-Inf +- I 0) = +-Inf + I +-0.
791
     *
792
     * sinh(NaN +- I 0)   = d(NaN) + I +-0.
793
     */
794
    if (y == 0 && !xfinite) {
795
        if (npy_isnan(x)) {
796
            return z;
797
        }
798
        return npy_cpack@c@(x, npy_copysign@c@(0, y));
799
    }
800

801
    /*
802
     * sinh(x +- I Inf) = dNaN + I dNaN.
803
     * Raise the invalid floating-point exception for finite nonzero x.
804
     *
805
     * sinh(x + I NaN) = d(NaN) + I d(NaN).
806
     * Optionally raises the invalid floating-point exception for finite
807
     * nonzero x.  Choice = don't raise (except for signaling NaNs).
808
     */
809
    if (xfinite && !yfinite) {
810
        return npy_cpack@c@(y - y, x * (y - y));
811
    }
812

813
    /*
814
     * sinh(+-Inf + I NaN)  = +-Inf + I d(NaN).
815
     * The sign of Inf in the result is unspecified.  Choice = normally
816
     * the same as d(NaN).
817
     *
818
     * sinh(+-Inf +- I Inf) = +Inf + I dNaN.
819
     * The sign of Inf in the result is unspecified.  Choice = always +.
820
     * Raise the invalid floating-point exception.
821
     *
822
     * sinh(+-Inf + I y)   = +-Inf cos(y) + I Inf sin(y)
823
     */
824
    if (!xfinite && !npy_isnan(x)) {
825
        if (!yfinite) {
826
            return npy_cpack@c@(x * x, x * (y - y));
827
        }
828
        return npy_cpack@c@(x * npy_cos@c@(y),
829
                            NPY_INFINITY@C@ * npy_sin@c@(y));
830
    }
831

832
    /*
833
     * sinh(NaN + I NaN)  = d(NaN) + I d(NaN).
834
     *
835
     * sinh(NaN +- I Inf) = d(NaN) + I d(NaN).
836
     * Optionally raises the invalid floating-point exception.
837
     * Choice = raise.
838
     *
839
     * sinh(NaN + I y)    = d(NaN) + I d(NaN).
840
     * Optionally raises the invalid floating-point exception for finite
841
     * nonzero y.  Choice = don't raise (except for signaling NaNs).
842
     */
843
    return npy_cpack@c@((x * x) * (y - y), (x + x) * (y - y));
844
}
845
#undef CSINH_BIG
846
#undef CSINH_HUGE
847
#endif
848

849
#ifndef HAVE_CTANH@C@
850
/*
851
 * Taken from the msun library in FreeBSD, rev 226600.
852
 *
853
 * Hyperbolic tangent of a complex argument z = x + i y.
854
 *
855
 * The algorithm is from:
856
 *
857
 *   W. Kahan.  Branch Cuts for Complex Elementary Functions or Much
858
 *   Ado About Nothing's Sign Bit.  In The State of the Art in
859
 *   Numerical Analysis, pp. 165 ff.  Iserles and Powell, eds., 1987.
860
 *
861
 * Method:
862
 *
863
 *   Let t    = tan(x)
864
 *       beta = 1/cos^2(y)
865
 *       s    = sinh(x)
866
 *       rho  = cosh(x)
867
 *
868
 *   We have:
869
 *
870
 *   tanh(z) = sinh(z) / cosh(z)
871
 *
872
 *             sinh(x) cos(y) + i cosh(x) sin(y)
873
 *           = ---------------------------------
874
 *             cosh(x) cos(y) + i sinh(x) sin(y)
875
 *
876
 *             cosh(x) sinh(x) / cos^2(y) + i tan(y)
877
 *           = -------------------------------------
878
 *                    1 + sinh^2(x) / cos^2(y)
879
 *
880
 *             beta rho s + i t
881
 *           = ----------------
882
 *               1 + beta s^2
883
 *
884
 * Modifications:
885
 *
886
 *   I omitted the original algorithm's handling of overflow in tan(x) after
887
 *   verifying with nearpi.c that this can't happen in IEEE single or double
888
 *   precision.  I also handle large x differently.
889
 */
890

891
#define TANH_HUGE 22.0
892
#define TANHF_HUGE 11.0F
893
#define TANHL_HUGE 42.0L
894

895
@ctype@
896
npy_ctanh@c@(@ctype@ z)
897
{
898
    @type@ x, y;
899
    @type@ t, beta, s, rho, denom;
900

901
    x = npy_creal@c@(z);
902
    y = npy_cimag@c@(z);
903

904
    /*
905
     * ctanh(NaN + i 0) = NaN + i 0
906
     *
907
     * ctanh(NaN + i y) = NaN + i NaN        for y != 0
908
     *
909
     * The imaginary part has the sign of x*sin(2*y), but there's no
910
     * special effort to get this right.
911
     *
912
     * ctanh(+-Inf +- i Inf) = +-1 +- 0
913
     *
914
     * ctanh(+-Inf + i y) = +-1 + 0 sin(2y)        for y finite
915
     *
916
     * The imaginary part of the sign is unspecified.  This special
917
     * case is only needed to avoid a spurious invalid exception when
918
     * y is infinite.
919
     */
920
        if (!npy_isfinite(x)) {
921
            if (npy_isnan(x)) {
922
                return npy_cpack@c@(x, (y == 0 ? y : x * y));
923
            }
924
            return npy_cpack@c@(npy_copysign@c@(1,x),
925
                                npy_copysign@c@(0,
926
                                npy_isinf(y) ?
927
                                    y : npy_sin@c@(y) * npy_cos@c@(y)));
928
        }
929

930
    /*
931
     * ctanh(x + i NAN) = NaN + i NaN
932
     * ctanh(x +- i Inf) = NaN + i NaN
933
     */
934
    if (!npy_isfinite(y)) {
935
        return (npy_cpack@c@(y - y, y - y));
936
    }
937

938
    /*
939
     * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the
940
     * approximation sinh^2(huge) ~= exp(2*huge) / 4.
941
     * We use a modified formula to avoid spurious overflow.
942
     */
943
    if (npy_fabs@c@(x) >= TANH@C@_HUGE) {
944
        @type@ exp_mx = npy_exp@c@(-npy_fabs@c@(x));
945
        return npy_cpack@c@(npy_copysign@c@(1, x),
946
                            4 * npy_sin@c@(y) * npy_cos@c@(y) *
947
                                exp_mx * exp_mx);
948
    }
949

950
    /* Kahan's algorithm */
951
    t = npy_tan@c@(y);
952
    beta = 1 + t * t;    /* = 1 / cos^2(y) */
953
    s = npy_sinh@c@(x);
954
    rho = npy_sqrt@c@(1 + s * s);    /* = cosh(x) */
955
    denom = 1 + beta * s * s;
956
    return (npy_cpack@c@((beta * rho * s) / denom, t / denom));
957
}
958
#undef TANH_HUGE
959
#undef TANHF_HUGE
960
#undef TANHL_HUGE
961
#endif
962

963
#if !defined (HAVE_CACOS@C@) || !defined (HAVE_CASINH@C@)
964
/*
965
 * Complex inverse trig functions taken from the msum library in FreeBSD
966
 * revision 251404
967
 *
968
 * The algorithm is very close to that in "Implementing the complex arcsine
969
 * and arccosine functions using exception handling" by T. E. Hull, Thomas F.
970
 * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on
971
 * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
972
 * http://dl.acm.org/citation.cfm?id=275324.
973
 *
974
 * Throughout we use the convention z = x + I*y.
975
 *
976
 * casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B)
977
 * where
978
 * A = (|z+I| + |z-I|) / 2
979
 * B = (|z+I| - |z-I|) / 2 = y/A
980
 *
981
 * These formulas become numerically unstable:
982
 *   (a) for Re(casinh(z)) when z is close to the line segment [-I, I] (that
983
 *       is, Re(casinh(z)) is close to 0);
984
 *   (b) for Im(casinh(z)) when z is close to either of the intervals
985
 *       [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is
986
 *       close to PI/2).
987
 *
988
 * These numerical problems are overcome by defining
989
 * f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2
990
 * Then if A < A_crossover, we use
991
 *   log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1)))
992
 *   A-1 = f(x, 1+y) + f(x, 1-y)
993
 * and if B > B_crossover, we use
994
 *   asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y)))
995
 *   A-y = f(x, y+1) + f(x, y-1)
996
 * where without loss of generality we have assumed that x and y are
997
 * non-negative.
998
 *
999
 * Much of the difficulty comes because the intermediate computations may
1000
 * produce overflows or underflows.  This is dealt with in the paper by Hull
1001
 * et al by using exception handling.  We do this by detecting when
1002
 * computations risk underflow or overflow.  The hardest part is handling the
1003
 * underflows when computing f(a, b).
1004
 *
1005
 * Note that the function f(a, b) does not appear explicitly in the paper by
1006
 * Hull et al, but the idea may be found on pages 308 and 309.  Introducing the
1007
 * function f(a, b) allows us to concentrate many of the clever tricks in this
1008
 * paper into one function.
1009
 */
1010

1011
/*
1012
 * Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2.
1013
 * Pass hypot(a, b) as the third argument.
1014
 */
1015
static NPY_INLINE @type@
1016
_f@c@(@type@ a, @type@ b, @type@ hypot_a_b)
1017
{
1018
    if (b < 0) {
1019
        return ((hypot_a_b - b) / 2);
1020
    }
1021
    if (b == 0) {
1022
        return (a / 2);
1023
    }
1024
    return (a * a / (hypot_a_b + b) / 2);
1025
}
1026

1027
/*
1028
 * All the hard work is contained in this function.
1029
 * x and y are assumed positive or zero, and less than RECIP_EPSILON.
1030
 * Upon return:
1031
 * rx = Re(casinh(z)) = -Im(cacos(y + I*x)).
1032
 * B_is_usable is set to 1 if the value of B is usable.
1033
 * If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y.
1034
 * If returning sqrt_A2my2 has potential to result in an underflow, it is
1035
 * rescaled, and new_y is similarly rescaled.
1036
 */
1037
static NPY_INLINE void
1038
_do_hard_work@c@(@type@ x, @type@ y, @type@ *rx,
1039
    npy_int *B_is_usable, @type@ *B, @type@ *sqrt_A2my2, @type@ *new_y)
1040
{
1041
#if @precision@ == 1
1042
    const npy_float A_crossover = 10.0f;
1043
    const npy_float B_crossover = 0.6417f;
1044
    const npy_float FOUR_SQRT_MIN = 4.3368086899420177e-19f;
1045
#endif
1046
#if @precision@ == 2
1047
    const npy_double A_crossover = 10.0;
1048
    const npy_double B_crossover = 0.6417;
1049
    const npy_double FOUR_SQRT_MIN = 5.9666725849601654e-154;
1050
#endif
1051
#if @precision@ == 3
1052
    const npy_longdouble A_crossover = 10.0l;
1053
    const npy_longdouble B_crossover = 0.6417l;
1054
#if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE
1055
    const npy_longdouble FOUR_SQRT_MIN = 5.9666725849601654e-154;
1056
#else
1057
    const npy_longdouble FOUR_SQRT_MIN = 7.3344154702193886625e-2466l;
1058
#endif
1059
#endif
1060
    @type@ R, S, A; /* A, B, R, and S are as in Hull et al. */
1061
    @type@ Am1, Amy; /* A-1, A-y. */
1062

1063
    R = npy_hypot@c@(x, y + 1);        /* |z+I| */
1064
    S = npy_hypot@c@(x, y - 1);        /* |z-I| */
1065

1066
    /* A = (|z+I| + |z-I|) / 2 */
1067
    A = (R + S) / 2;
1068
    /*
1069
     * Mathematically A >= 1.  There is a small chance that this will not
1070
     * be so because of rounding errors.  So we will make certain it is
1071
     * so.
1072
     */
1073
    if (A < 1) {
1074
        A = 1;
1075
    }
1076

1077
    if (A < A_crossover) {
1078
        /*
1079
         * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y).
1080
         * rx = log1p(Am1 + sqrt(Am1*(A+1)))
1081
         */
1082
        if (y == 1 && x < @TEPS@ * @TEPS@ / 128) {
1083
            /*
1084
             * fp is of order x^2, and fm = x/2.
1085
             * A = 1 (inexactly).
1086
             */
1087
            *rx = npy_sqrt@c@(x);
1088
        }
1089
        else if (x >= @TEPS@ * npy_fabs@c@(y - 1)) {
1090
            /*
1091
             * Underflow will not occur because
1092
             * x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN
1093
             */
1094
            Am1 = _f@c@(x, 1 + y, R) + _f@c@(x, 1 - y, S);
1095
            *rx = npy_log1p@c@(Am1 + npy_sqrt@c@(Am1 * (A + 1)));
1096
        }
1097
        else if (y < 1) {
1098
            /*
1099
             * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and
1100
             * A = 1 (inexactly).
1101
             */
1102
            *rx = x / npy_sqrt@c@((1 - y) * (1 + y));
1103
        }
1104
        else {        /* if (y > 1) */
1105
            /*
1106
             * A-1 = y-1 (inexactly).
1107
             */
1108
            *rx = npy_log1p@c@((y - 1) + npy_sqrt@c@((y - 1) * (y + 1)));
1109
        }
1110
    }
1111
    else {
1112
        *rx = npy_log@c@(A + npy_sqrt@c@(A * A - 1));
1113
    }
1114

1115
    *new_y = y;
1116

1117
    if (y < FOUR_SQRT_MIN) {
1118
        /*
1119
         * Avoid a possible underflow caused by y/A.  For casinh this
1120
         * would be legitimate, but will be picked up by invoking atan2
1121
         * later on.  For cacos this would not be legitimate.
1122
         */
1123
        *B_is_usable = 0;
1124
        *sqrt_A2my2 = A * (2 / @TEPS@);
1125
        *new_y = y * (2 / @TEPS@);
1126
        return;
1127
    }
1128

1129
    /* B = (|z+I| - |z-I|) / 2 = y/A */
1130
    *B = y / A;
1131
    *B_is_usable = 1;
1132

1133
    if (*B > B_crossover) {
1134
        *B_is_usable = 0;
1135
        /*
1136
         * Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1).
1137
         * sqrt_A2my2 = sqrt(Amy*(A+y))
1138
         */
1139
        if (y == 1 && x < @TEPS@ / 128) {
1140
            /*
1141
             * fp is of order x^2, and fm = x/2.
1142
             * A = 1 (inexactly).
1143
             */
1144
            *sqrt_A2my2 = npy_sqrt@c@(x) * npy_sqrt@c@((A + y) / 2);
1145
        }
1146
        else if (x >= @TEPS@ * npy_fabs@c@(y - 1)) {
1147
            /*
1148
             * Underflow will not occur because
1149
             * x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN
1150
             * and
1151
             * x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN
1152
             */
1153
            Amy = _f@c@(x, y + 1, R) + _f@c@(x, y - 1, S);
1154
            *sqrt_A2my2 = npy_sqrt@c@(Amy * (A + y));
1155
        }
1156
        else if (y > 1) {
1157
            /*
1158
             * fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and
1159
             * A = y (inexactly).
1160
             *
1161
             * y < RECIP_EPSILON.  So the following
1162
             * scaling should avoid any underflow problems.
1163
             */
1164
            *sqrt_A2my2 = x * (4 / @TEPS@ / @TEPS@) * y /
1165
                npy_sqrt@c@((y + 1) * (y - 1));
1166
            *new_y = y * (4 / @TEPS@ / @TEPS@);
1167
        }
1168
        else {        /* if (y < 1) */
1169
            /*
1170
             * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and
1171
             * A = 1 (inexactly).
1172
             */
1173
            *sqrt_A2my2 = npy_sqrt@c@((1 - y) * (1 + y));
1174
        }
1175
    }
1176
}
1177

1178
/*
1179
 * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON.
1180
 */
1181
static NPY_INLINE void
1182
_clog_for_large_values@c@(@type@ x, @type@ y,
1183
    @type@ *rr, @type@ *ri)
1184
{
1185
#if @precision@ == 1
1186
    const npy_float QUARTER_SQRT_MAX = 4.611685743549481e+18f;
1187
    const npy_float SQRT_MIN = 1.0842021724855044e-19f;
1188
 #endif
1189
#if @precision@ == 2
1190
    const npy_double QUARTER_SQRT_MAX = 3.3519519824856489e+153;
1191
    const npy_double SQRT_MIN = 1.4916681462400413e-154;
1192
 #endif
1193
#if @precision@ == 3
1194
#if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE
1195
    const npy_longdouble QUARTER_SQRT_MAX = 3.3519519824856489e+153;
1196
    const npy_longdouble SQRT_MIN = 1.4916681462400413e-154;
1197
#else
1198
    const npy_longdouble QUARTER_SQRT_MAX = 2.7268703390485398235e+2465l;
1199
    const npy_longdouble SQRT_MIN = 1.8336038675548471656e-2466l;
1200
#endif
1201
#endif
1202
    @type@ ax, ay, t;
1203

1204
    ax = npy_fabs@c@(x);
1205
    ay = npy_fabs@c@(y);
1206
    if (ax < ay) {
1207
        t = ax;
1208
        ax = ay;
1209
        ay = t;
1210
    }
1211

1212
    /*
1213
     * Avoid overflow in hypot() when x and y are both very large.
1214
     * Divide x and y by E, and then add 1 to the logarithm.  This depends
1215
     * on E being larger than sqrt(2).
1216
     * Dividing by E causes an insignificant loss of accuracy; however
1217
     * this method is still poor since it is unnecessarily slow.
1218
     */
1219
    if (ax > @TMAX@ / 2) {
1220
        *rr = npy_log@c@(npy_hypot@c@(x / NPY_E@c@, y / NPY_E@c@)) + 1;
1221
    }
1222
    /*
1223
     * Avoid overflow when x or y is large.  Avoid underflow when x or
1224
     * y is small.
1225
     */
1226
    else if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) {
1227
        *rr = npy_log@c@(npy_hypot@c@(x, y));
1228
    }
1229
    else {
1230
        *rr = npy_log@c@(ax * ax + ay * ay) / 2;
1231
    }
1232
    *ri = npy_atan2@c@(y, x);
1233
}
1234
#endif
1235

1236
#ifndef HAVE_CACOS@C@
1237
@ctype@
1238
npy_cacos@c@(@ctype@ z)
1239
{
1240
#if @precision@ == 1
1241
    /* this is sqrt(6*EPS) */
1242
    const npy_float SQRT_6_EPSILON = 8.4572793338e-4f;
1243
    /* chosen such that pio2_hi + pio2_lo == pio2_hi but causes FE_INEXACT. */
1244
    const volatile npy_float pio2_lo = 7.5497899549e-9f;
1245
#endif
1246
#if @precision@ == 2
1247
    const npy_double SQRT_6_EPSILON = 3.65002414998885671e-08;
1248
    const volatile npy_double pio2_lo = 6.1232339957367659e-17;
1249
#endif
1250
#if @precision@ == 3
1251
    const npy_longdouble SQRT_6_EPSILON = 8.0654900873493277169e-10l;
1252
    const volatile npy_longdouble pio2_lo = 2.710505431213761085e-20l;
1253
#endif
1254
    const @type@ RECIP_EPSILON = 1.0@c@ / @TEPS@;
1255
    const @type@ pio2_hi = NPY_PI_2@c@;
1256
    @type@ x, y, ax, ay, wx, wy, rx, ry, B, sqrt_A2mx2, new_x;
1257
    npy_int sx, sy;
1258
    npy_int B_is_usable;
1259

1260
    x = npy_creal@c@(z);
1261
    y = npy_cimag@c@(z);
1262
    sx = npy_signbit(x);
1263
    sy = npy_signbit(y);
1264
    ax = npy_fabs@c@(x);
1265
    ay = npy_fabs@c@(y);
1266

1267
    if (npy_isnan(x) || npy_isnan(y)) {
1268
        /* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
1269
        if (npy_isinf(x)) {
1270
            return npy_cpack@c@(y + y, -NPY_INFINITY@C@);
1271
        }
1272
        /* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */
1273
        if (npy_isinf(y)) {
1274
            return npy_cpack@c@(x + x, -y);
1275
        }
1276
        /* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */
1277
        if (x == 0) {
1278
            return npy_cpack@c@(pio2_hi + pio2_lo, y + y);
1279
        }
1280
        /*
1281
         * All other cases involving NaN return NaN + I*NaN.
1282
         * C99 leaves it optional whether to raise invalid if one of
1283
         * the arguments is not NaN, so we opt not to raise it.
1284
         */
1285
        return npy_cpack@c@(NPY_NAN@C@, NPY_NAN@C@);
1286
    }
1287

1288
    if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
1289
        /* clog...() will raise inexact unless x or y is infinite. */
1290
        _clog_for_large_values@c@(x, y, &wx, &wy);
1291
        rx = npy_fabs@c@(wy);
1292
        ry = wx + NPY_LOGE2@c@;
1293
        if (sy == 0) {
1294
            ry = -ry;
1295
        }
1296
        return npy_cpack@c@(rx, ry);
1297
    }
1298

1299
    /* Avoid spuriously raising inexact for z = 1. */
1300
    if (x == 1 && y == 0) {
1301
        return npy_cpack@c@(0, -y);
1302
    }
1303

1304
    /* All remaining cases are inexact. */
1305
    raise_inexact();
1306

1307
    if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) {
1308
        return npy_cpack@c@(pio2_hi - (x - pio2_lo), -y);
1309
    }
1310

1311
    _do_hard_work@c@(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
1312
    if (B_is_usable) {
1313
        if (sx == 0) {
1314
            rx = npy_acos@c@(B);
1315
        }
1316
        else {
1317
            rx = npy_acos@c@(-B);
1318
        }
1319
    }
1320
    else {
1321
        if (sx == 0) {
1322
            rx = npy_atan2@c@(sqrt_A2mx2, new_x);
1323
        }
1324
        else {
1325
            rx = npy_atan2@c@(sqrt_A2mx2, -new_x);
1326
        }
1327
    }
1328
    if (sy == 0) {
1329
        ry = -ry;
1330
    }
1331
    return npy_cpack@c@(rx, ry);
1332
}
1333
#endif
1334

1335
#ifndef HAVE_CASIN@C@
1336
@ctype@
1337
npy_casin@c@(@ctype@ z)
1338
{
1339
    /* casin(z) = I * conj( casinh(I * conj(z)) ) */
1340
    z = npy_casinh@c@(npy_cpack@c@(npy_cimag@c@(z), npy_creal@c@(z)));
1341
    return npy_cpack@c@(npy_cimag@c@(z), npy_creal@c@(z));
1342
}
1343
#endif
1344

1345
#ifndef HAVE_CATAN@C@
1346
@ctype@
1347
npy_catan@c@(@ctype@ z)
1348
{
1349
    /* catan(z) = I * conj( catanh(I * conj(z)) ) */
1350
    z = npy_catanh@c@(npy_cpack@c@(npy_cimag@c@(z), npy_creal@c@(z)));
1351
    return npy_cpack@c@(npy_cimag@c@(z), npy_creal@c@(z));
1352
}
1353
#endif
1354

1355
#ifndef HAVE_CACOSH@C@
1356
@ctype@
1357
npy_cacosh@c@(@ctype@ z)
1358
{
1359
    /*
1360
     * cacosh(z) = I*cacos(z) or -I*cacos(z)
1361
     * where the sign is chosen so Re(cacosh(z)) >= 0.
1362
     */
1363
    @ctype@  w;
1364
    @type@ rx, ry;
1365

1366
    w = npy_cacos@c@(z);
1367
    rx = npy_creal@c@(w);
1368
    ry = npy_cimag@c@(w);
1369
    /* cacosh(NaN + I*NaN) = NaN + I*NaN */
1370
    if (npy_isnan(rx) && npy_isnan(ry)) {
1371
        return npy_cpack@c@(ry, rx);
1372
    }
1373
    /* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */
1374
    /* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */
1375
    if (npy_isnan(rx)) {
1376
        return npy_cpack@c@(npy_fabs@c@(ry), rx);
1377
    }
1378
    /* cacosh(0 + I*NaN) = NaN + I*NaN */
1379
    if (npy_isnan(ry)) {
1380
        return npy_cpack@c@(ry, ry);
1381
    }
1382
    return npy_cpack@c@(npy_fabs@c@(ry), npy_copysign@c@(rx, npy_cimag@c@(z)));
1383
}
1384
#endif
1385

1386
#ifndef HAVE_CASINH@C@
1387
/*
1388
 * casinh(z) = z + O(z^3)   as z -> 0
1389
 *
1390
 * casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2)   as z -> infinity
1391
 * The above formula works for the imaginary part as well, because
1392
 * Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3)
1393
 *    as z -> infinity, uniformly in y
1394
 */
1395
@ctype@
1396
npy_casinh@c@(@ctype@ z)
1397
{
1398
#if @precision@ == 1
1399
    /* this is sqrt(6*EPS) */
1400
    const npy_float SQRT_6_EPSILON = 8.4572793338e-4f;
1401
#endif
1402
#if @precision@ == 2
1403
    const npy_double SQRT_6_EPSILON = 3.65002414998885671e-08;
1404
#endif
1405
#if @precision@ == 3
1406
    const npy_longdouble SQRT_6_EPSILON = 8.0654900873493277169e-10l;
1407
#endif
1408
    const @type@ RECIP_EPSILON = 1.0@c@ / @TEPS@;
1409
    @type@ x, y, ax, ay, wx, wy, rx, ry, B, sqrt_A2my2, new_y;
1410
    npy_int B_is_usable;
1411

1412
    x = npy_creal@c@(z);
1413
    y = npy_cimag@c@(z);
1414
    ax = npy_fabs@c@(x);
1415
    ay = npy_fabs@c@(y);
1416

1417
    if (npy_isnan(x) || npy_isnan(y)) {
1418
        /* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */
1419
        if (npy_isinf(x)) {
1420
            return npy_cpack@c@(x, y + y);
1421
        }
1422
        /* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
1423
        if (npy_isinf(y)) {
1424
            return npy_cpack@c@(y, x + x);
1425
        }
1426
        /* casinh(NaN + I*0) = NaN + I*0 */
1427
        if (y == 0) {
1428
            return npy_cpack@c@(x + x, y);
1429
        }
1430
        /*
1431
         * All other cases involving NaN return NaN + I*NaN.
1432
         * C99 leaves it optional whether to raise invalid if one of
1433
         * the arguments is not NaN, so we opt not to raise it.
1434
         */
1435
        return npy_cpack@c@(NPY_NAN@C@, NPY_NAN@C@);
1436
    }
1437

1438
    if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
1439
        /* clog...() will raise inexact unless x or y is infinite. */
1440
        if (npy_signbit(x) == 0) {
1441
            _clog_for_large_values@c@(x, y, &wx, &wy);
1442
            wx += NPY_LOGE2@c@;
1443
        }
1444
        else {
1445
            _clog_for_large_values@c@(-x, -y, &wx, &wy);
1446
            wx += NPY_LOGE2@c@;
1447
        }
1448
        return npy_cpack@c@(npy_copysign@c@(wx, x), npy_copysign@c@(wy, y));
1449
    }
1450

1451
    /* Avoid spuriously raising inexact for z = 0. */
1452
    if (x == 0 && y == 0) {
1453
        return (z);
1454
    }
1455

1456
    /* All remaining cases are inexact. */
1457
    raise_inexact();
1458

1459
    if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) {
1460
        return (z);
1461
    }
1462

1463
    _do_hard_work@c@(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
1464
    if (B_is_usable) {
1465
        ry = npy_asin@c@(B);
1466
    }
1467
    else {
1468
        ry = npy_atan2@c@(new_y, sqrt_A2my2);
1469
    }
1470
    return npy_cpack@c@(npy_copysign@c@(rx, x), npy_copysign@c@(ry, y));
1471
}
1472
#endif
1473

1474
#ifndef HAVE_CATANH@C@
1475
/*
1476
 * sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow).
1477
 * Assumes x*x and y*y will not overflow.
1478
 * Assumes x and y are finite.
1479
 * Assumes y is non-negative.
1480
 * Assumes fabs(x) >= DBL_EPSILON.
1481
 */
1482
static NPY_INLINE @type@
1483
_sum_squares@c@(@type@ x, @type@ y)
1484
{
1485
#if @precision@ == 1
1486
const npy_float SQRT_MIN = 1.0842022e-19f;
1487
#endif
1488
#if @precision@ == 2
1489
const npy_double SQRT_MIN = 1.4916681462400413e-154; /* sqrt(DBL_MIN) */
1490
#endif
1491
#if @precision@ == 3
1492
#if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE
1493
const npy_longdouble SQRT_MIN = 1.4916681462400413e-154; /* sqrt(DBL_MIN) */
1494
#else
1495
/* this is correct for 80 bit long doubles */
1496
const npy_longdouble SQRT_MIN = 1.8336038675548471656e-2466l;
1497
#endif
1498
#endif
1499
    /* Avoid underflow when y is small. */
1500
    if (y < SQRT_MIN) {
1501
        return (x * x);
1502
    }
1503

1504
    return (x * x + y * y);
1505
}
1506

1507
/*
1508
 * real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y).
1509
 * Assumes x and y are not NaN, and one of x and y is larger than
1510
 * RECIP_EPSILON.  We avoid unwarranted underflow.  It is important to not use
1511
 * the code creal(1/z), because the imaginary part may produce an unwanted
1512
 * underflow.
1513
 * This is only called in a context where inexact is always raised before
1514
 * the call, so no effort is made to avoid or force inexact.
1515
 */
1516
#if @precision@ == 1
1517
#define BIAS (FLT_MAX_EXP - 1)
1518
#define CUTOFF (FLT_MANT_DIG / 2 + 1)
1519
static NPY_INLINE npy_float
1520
_real_part_reciprocalf(npy_float x, npy_float y)
1521
{
1522
    npy_float scale;
1523
    npy_uint32 hx, hy;
1524
    npy_int32 ix, iy;
1525

1526
    GET_FLOAT_WORD(hx, x);
1527
    ix = hx & 0x7f800000;
1528
    GET_FLOAT_WORD(hy, y);
1529
    iy = hy & 0x7f800000;
1530
    if (ix - iy >= CUTOFF << 23 || npy_isinf(x)) {
1531
        return (1 / x);
1532
    }
1533
    if (iy - ix >= CUTOFF << 23) {
1534
        return (x / y / y);
1535
    }
1536
    if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23) {
1537
        return (x / (x * x + y * y));
1538
    }
1539
    SET_FLOAT_WORD(scale, 0x7f800000 - ix);
1540
    x *= scale;
1541
    y *= scale;
1542
    return (x / (x * x + y * y) * scale);
1543
}
1544
#undef BIAS
1545
#undef CUTOFF
1546
#endif
1547

1548
#if @precision@ == 2
1549
#define BIAS (DBL_MAX_EXP - 1)
1550
/*  more guard digits are useful iff there is extra precision. */
1551
#define CUTOFF (DBL_MANT_DIG / 2 + 1)  /* just half or 1 guard digit */
1552
static NPY_INLINE npy_double
1553
_real_part_reciprocal(npy_double x, npy_double y)
1554
{
1555
    npy_double scale;
1556
    npy_uint32 hx, hy;
1557
    npy_int32 ix, iy;
1558

1559
    /*
1560
     * This code is inspired by the C99 document n1124.pdf, Section G.5.1,
1561
     * example 2.
1562
     */
1563
    GET_HIGH_WORD(hx, x);
1564
    ix = hx & 0x7ff00000;
1565
    GET_HIGH_WORD(hy, y);
1566
    iy = hy & 0x7ff00000;
1567
    if (ix - iy >= CUTOFF << 20 || npy_isinf(x)) {
1568
        /* +-Inf -> +-0 is special */
1569
        return (1 / x);
1570
    }
1571
    if (iy - ix >= CUTOFF << 20) {
1572
        /* should avoid double div, but hard */
1573
        return (x / y / y);
1574
    }
1575
    if (ix <= (BIAS + DBL_MAX_EXP / 2 - CUTOFF) << 20) {
1576
        return (x / (x * x + y * y));
1577
    }
1578
    scale = 1;
1579
    SET_HIGH_WORD(scale, 0x7ff00000 - ix);  /* 2**(1-ilogb(x)) */
1580
    x *= scale;
1581
    y *= scale;
1582
    return (x / (x * x + y * y) * scale);
1583
}
1584
#undef BIAS
1585
#undef CUTOFF
1586
#endif
1587

1588
#if @precision@ == 3
1589
#if !defined(HAVE_LDOUBLE_DOUBLE_DOUBLE_BE) && \
1590
    !defined(HAVE_LDOUBLE_DOUBLE_DOUBLE_LE)
1591

1592
#define BIAS (LDBL_MAX_EXP - 1)
1593
#define CUTOFF (LDBL_MANT_DIG / 2 + 1)
1594
static NPY_INLINE npy_longdouble
1595
_real_part_reciprocall(npy_longdouble x,
1596
    npy_longdouble y)
1597
{
1598
    npy_longdouble scale;
1599
    union IEEEl2bitsrep ux, uy, us;
1600
    npy_int32 ix, iy;
1601

1602
    ux.e = x;
1603
    ix = GET_LDOUBLE_EXP(ux);
1604
    uy.e = y;
1605
    iy = GET_LDOUBLE_EXP(uy);
1606
    if (ix - iy >= CUTOFF || npy_isinf(x)) {
1607
        return (1/x);
1608
    }
1609
    if (iy - ix >= CUTOFF) {
1610
        return (x/y/y);
1611
    }
1612
    if (ix <= BIAS + LDBL_MAX_EXP / 2 - CUTOFF) {
1613
        return (x/(x*x + y*y));
1614
    }
1615
    us.e = 1;
1616
    SET_LDOUBLE_EXP(us, 0x7fff - ix);
1617
    scale = us.e;
1618
    x *= scale;
1619
    y *= scale;
1620
    return (x/(x*x + y*y) * scale);
1621
}
1622
#undef BIAS
1623
#undef CUTOFF
1624

1625
#else
1626

1627
static NPY_INLINE npy_longdouble
1628
_real_part_reciprocall(npy_longdouble x,
1629
    npy_longdouble y)
1630
{
1631
    return x/(x*x + y*y);
1632
}
1633

1634
#endif
1635
#endif
1636

1637
@ctype@
1638
npy_catanh@c@(@ctype@ z)
1639
{
1640
#if @precision@ == 1
1641
    /* this is sqrt(3*EPS) */
1642
    const npy_float SQRT_3_EPSILON = 5.9801995673e-4f;
1643
    /* chosen such that pio2_hi + pio2_lo == pio2_hi but causes FE_INEXACT. */
1644
    const volatile npy_float pio2_lo = 7.5497899549e-9f;
1645
#endif
1646
#if @precision@ == 2
1647
    const npy_double SQRT_3_EPSILON = 2.5809568279517849e-8;
1648
    const volatile npy_double pio2_lo = 6.1232339957367659e-17;
1649
#endif
1650
#if @precision@ == 3
1651
    const npy_longdouble SQRT_3_EPSILON = 5.70316273435758915310e-10l;
1652
    const volatile npy_longdouble pio2_lo = 2.710505431213761085e-20l;
1653
#endif
1654
    const @type@ RECIP_EPSILON = 1.0@c@ / @TEPS@;
1655
    const @type@ pio2_hi = NPY_PI_2@c@;
1656
    @type@ x, y, ax, ay, rx, ry;
1657

1658
    x = npy_creal@c@(z);
1659
    y = npy_cimag@c@(z);
1660
    ax = npy_fabs@c@(x);
1661
    ay = npy_fabs@c@(y);
1662

1663
    /* This helps handle many cases. */
1664
    if (y == 0 && ax <= 1) {
1665
        return npy_cpack@c@(npy_atanh@c@(x), y);
1666
    }
1667

1668
    /* To ensure the same accuracy as atan(), and to filter out z = 0. */
1669
    if (x == 0) {
1670
        return npy_cpack@c@(x, npy_atan@c@(y));
1671
    }
1672

1673
    if (npy_isnan(x) || npy_isnan(y)) {
1674
        /* catanh(+-Inf + I*NaN) = +-0 + I*NaN */
1675
        if (npy_isinf(x)) {
1676
            return npy_cpack@c@(npy_copysign@c@(0, x), y + y);
1677
        }
1678
        /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
1679
        if (npy_isinf(y)) {
1680
            return npy_cpack@c@(npy_copysign@c@(0, x),
1681
                npy_copysign@c@(pio2_hi + pio2_lo, y));
1682
        }
1683
        /*
1684
         * All other cases involving NaN return NaN + I*NaN.
1685
         * C99 leaves it optional whether to raise invalid if one of
1686
         * the arguments is not NaN, so we opt not to raise it.
1687
         */
1688
        return npy_cpack@c@(NPY_NAN@C@, NPY_NAN@C@);
1689
    }
1690

1691
    if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
1692
        return npy_cpack@c@(_real_part_reciprocal@c@(x, y),
1693
            npy_copysign@c@(pio2_hi + pio2_lo, y));
1694
    }
1695

1696
    if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
1697
        /*
1698
         * z = 0 was filtered out above.  All other cases must raise
1699
         * inexact, but this is the only only that needs to do it
1700
         * explicitly.
1701
         */
1702
        raise_inexact();
1703
        return (z);
1704
    }
1705

1706
    if (ax == 1 && ay < @TEPS@) {
1707
        rx = (NPY_LOGE2@c@ - npy_log@c@(ay)) / 2;
1708
    }
1709
    else {
1710
        rx = npy_log1p@c@(4 * ax / _sum_squares@c@(ax - 1, ay)) / 4;
1711
    }
1712

1713
    if (ax == 1) {
1714
        ry = npy_atan2@c@(2, -ay) / 2;
1715
    }
1716
    else if (ay < @TEPS@) {
1717
        ry = npy_atan2@c@(2 * ay, (1 - ax) * (1 + ax)) / 2;
1718
    }
1719
    else {
1720
        ry = npy_atan2@c@(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
1721
    }
1722

1723
    return npy_cpack@c@(npy_copysign@c@(rx, x), npy_copysign@c@(ry, y));
1724
}
1725
#endif
1726
/**end repeat**/
1727

1728
/*==========================================================
1729
 * Decorate all the functions which are available natively
1730
 *=========================================================*/
1731

1732
/**begin repeat
1733
 * #type = npy_float, npy_double, npy_longdouble#
1734
 * #ctype = npy_cfloat, npy_cdouble, npy_clongdouble#
1735
 * #c = f, , l#
1736
 * #C = F, , L#
1737
 */
1738

1739
/**begin repeat1
1740
 * #kind = cabs,carg#
1741
 * #KIND = CABS,CARG#
1742
 */
1743
#ifdef HAVE_@KIND@@C@
1744
@type@
1745 1
npy_@kind@@c@(@ctype@ z)
1746
{
1747
    __@ctype@_to_c99_cast z1;
1748 1
    z1.npy_z = z;
1749 1
    return @kind@@c@(z1.c99_z);
1750
}
1751
#endif
1752
/**end repeat1**/
1753

1754
/**begin repeat1
1755
 * #kind = cexp,clog,csqrt,ccos,csin,ctan,ccosh,csinh,ctanh,
1756
 *         cacos,casin,catan,cacosh,casinh,catanh#
1757
 * #KIND = CEXP,CLOG,CSQRT,CCOS,CSIN,CTAN,CCOSH,CSINH,CTANH,
1758
 *         CACOS,CASIN,CATAN,CACOSH,CASINH,CATANH#
1759
 */
1760
#ifdef HAVE_@KIND@@C@
1761
@ctype@
1762 1
npy_@kind@@c@(@ctype@ z)
1763
{
1764
    __@ctype@_to_c99_cast z1;
1765
    __@ctype@_to_c99_cast ret;
1766 1
    z1.npy_z = z;
1767 1
    ret.c99_z = @kind@@c@(z1.c99_z);
1768 1
    return ret.npy_z;
1769
}
1770
#endif
1771
/**end repeat1**/
1772

1773

1774
/**end repeat**/

Read our documentation on viewing source code .

Loading