1# Implement sin(pi*z) and cos(pi*z) for complex z. Since the periods
2# of these functions are integral (and thus better representable in
3# floating point), it's possible to compute them with greater accuracy
4# than sin(z), cos(z).
5#
6from libc.math cimport sin, cos, sinh, cosh, exp, fabs, fmod, M_PI
7
8from ._cephes cimport sinpi as dsinpi, cospi as dcospi
9from ._complexstuff cimport number_t, double_complex, zpack
10
11cdef extern from "numpy/npy_math.h":
12    double npy_copysign(double x, double y) nogil
13    double NPY_INFINITY
14
15
16cdef inline double complex csinpi(double complex z) nogil:
17    """Compute sin(pi*z) for complex arguments."""
18    cdef:
19        double x = z.real
20        double piy = M_PI*z.imag
21        double abspiy = fabs(piy)
22        double sinpix = sinpi(x)
23        double cospix = cospi(x)
24        double exphpiy, coshfac, sinhfac
25
26    if abspiy < 700:
27        return zpack(sinpix*cosh(piy), cospix*sinh(piy))
28
29    # Have to be careful--sinh/cosh could overflow while cos/sin are
30    # small. At this large of values
31    #
32    # cosh(y) ~ exp(y)/2
33    # sinh(y) ~ sgn(y)*exp(y)/2
34    #
35    # so we can compute exp(y/2), scale by the right factor of sin/cos
36    # and then multiply by exp(y/2) to avoid overflow.
37    exphpiy = exp(abspiy/2)
38    if exphpiy == NPY_INFINITY:
39        if sinpix == 0:
40            # Preserve the sign of zero
41            coshfac = npy_copysign(0.0, sinpix)
42        else:
43            coshfac = npy_copysign(NPY_INFINITY, sinpix)
44        if cospix == 0:
45            sinhfac = npy_copysign(0.0, cospix)
46        else:
47            sinhfac = npy_copysign(NPY_INFINITY, cospix)
48        return zpack(coshfac, sinhfac)
49
50    coshfac = 0.5*sinpix*exphpiy
51    sinhfac = 0.5*cospix*exphpiy
52    return zpack(coshfac*exphpiy, sinhfac*exphpiy)
53
54
55cdef inline double complex ccospi(double complex z) nogil:
56    """Compute cos(pi*z) for complex arguments."""
57    cdef:
58        double x = z.real
59        double piy = M_PI*z.imag
60        double abspiy = fabs(piy)
61        double sinpix = sinpi(x)
62        double cospix = cospi(x)
63        double exphpiy, coshfac, sinhfac
64
65    if abspiy < 700:
66        return zpack(cospix*cosh(piy), -sinpix*sinh(piy))
67
68    # See csinpi(z) for an idea of what's going on here
69    exphpiy = exp(abspiy/2)
70    if exphpiy == NPY_INFINITY:
71        if sinpix == 0:
72            coshfac = npy_copysign(0.0, cospix)
73        else:
74            coshfac = npy_copysign(NPY_INFINITY, cospix)
75        if cospix == 0:
76            sinhfac = npy_copysign(0.0, sinpix)
77        else:
78            sinhfac = npy_copysign(NPY_INFINITY, sinpix)
79        return zpack(coshfac, sinhfac)
80
81    coshfac = 0.5*cospix*exphpiy
82    sinhfac = 0.5*sinpix*exphpiy
83    return zpack(coshfac*exphpiy, sinhfac*exphpiy)
84
85
86cdef inline number_t sinpi(number_t z) nogil:
87    if number_t is double:
88        return dsinpi(z)
89    else:
90        return csinpi(z)
91
92
93cdef inline number_t cospi(number_t z) nogil:
94    if number_t is double:
95        return dcospi(z)
96    else:
97        return ccospi(z)
98