1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
|
# -*-cython-*-
#
# Common functions required when doing complex arithmetic with Cython.
#
import cython
cimport numpy as np
cimport libc.math
cdef extern from "_complexstuff.h":
double npy_cabs(np.npy_cdouble z) nogil
double npy_carg(np.npy_cdouble z) nogil
np.npy_cdouble npy_clog(np.npy_cdouble z) nogil
np.npy_cdouble npy_cexp(np.npy_cdouble z) nogil
np.npy_cdouble npy_csin(np.npy_cdouble z) nogil
np.npy_cdouble npy_ccos(np.npy_cdouble z) nogil
np.npy_cdouble npy_csqrt(np.npy_cdouble z) nogil
np.npy_cdouble npy_cpow(np.npy_cdouble x, np.npy_cdouble y) nogil
double npy_log1p(double x) nogil
int npy_isnan(double x) nogil
int npy_isinf(double x) nogil
int npy_isfinite(double x) nogil
double inf "NPY_INFINITY"
double pi "NPY_PI"
double nan "NPY_NAN"
DEF tol = 2.220446092504131e-16
ctypedef double complex double_complex
ctypedef fused number_t:
double
double_complex
ctypedef union _complex_pun:
np.npy_cdouble npy
double_complex c99
cdef inline np.npy_cdouble npy_cdouble_from_double_complex(
double_complex x) nogil:
cdef _complex_pun z
z.c99 = x
return z.npy
cdef inline double_complex double_complex_from_npy_cdouble(
np.npy_cdouble x) nogil:
cdef _complex_pun z
z.npy = x
return z.c99
cdef inline bint zisnan(number_t x) nogil:
if number_t is double_complex:
return npy_isnan(x.real) or npy_isnan(x.imag)
else:
return npy_isnan(x)
cdef inline bint zisfinite(number_t x) nogil:
if number_t is double_complex:
return npy_isfinite(x.real) and npy_isfinite(x.imag)
else:
return npy_isfinite(x)
cdef inline bint zisinf(number_t x) nogil:
if number_t is double_complex:
return not zisnan(x) and not zisfinite(x)
else:
return npy_isinf(x)
cdef inline double zreal(number_t x) nogil:
if number_t is double_complex:
return x.real
else:
return x
cdef inline double zabs(number_t x) nogil:
if number_t is double_complex:
return npy_cabs(npy_cdouble_from_double_complex(x))
else:
return libc.math.fabs(x)
cdef inline double zarg(double complex x) nogil:
return npy_carg(npy_cdouble_from_double_complex(x))
cdef inline number_t zlog(number_t x) nogil:
cdef np.npy_cdouble r
if number_t is double_complex:
r = npy_clog(npy_cdouble_from_double_complex(x))
return double_complex_from_npy_cdouble(r)
else:
return libc.math.log(x)
cdef inline number_t zexp(number_t x) nogil:
cdef np.npy_cdouble r
if number_t is double_complex:
r = npy_cexp(npy_cdouble_from_double_complex(x))
return double_complex_from_npy_cdouble(r)
else:
return libc.math.exp(x)
cdef inline number_t zsin(number_t x) nogil:
cdef np.npy_cdouble r
if number_t is double_complex:
r = npy_csin(npy_cdouble_from_double_complex(x))
return double_complex_from_npy_cdouble(r)
else:
return libc.math.sin(x)
cdef inline number_t zcos(number_t x) nogil:
cdef np.npy_cdouble r
if number_t is double_complex:
r = npy_ccos(npy_cdouble_from_double_complex(x))
return double_complex_from_npy_cdouble(r)
else:
return libc.math.cos(x)
cdef inline number_t zsqrt(number_t x) nogil:
cdef np.npy_cdouble r
if number_t is double_complex:
r = npy_csqrt(npy_cdouble_from_double_complex(x))
return double_complex_from_npy_cdouble(r)
else:
return libc.math.sqrt(x)
cdef inline number_t zpow(number_t x, double y) nogil:
cdef np.npy_cdouble r, z
# FIXME
if number_t is double_complex:
z.real = y
z.imag = 0.0
r = npy_cpow(npy_cdouble_from_double_complex(x), z)
return double_complex_from_npy_cdouble(r)
else:
return libc.math.pow(x, y)
cdef inline double_complex zpack(double zr, double zi) nogil:
cdef np.npy_cdouble z
z.real = zr
z.imag = zi
return double_complex_from_npy_cdouble(z)
@cython.cdivision(True)
cdef inline double complex zdiv(double complex x, double complex y) nogil:
"""
Cython 0.24 uses the naive complex division algorithm which
overflows far before it should. See
https://groups.google.com/forum/#!topic/cython-users/1oSGbfiX7qw
This function implements Smith's method to get around the
problem. Smith's method can be further improved, see
http://arxiv.org/pdf/1210.4539v2.pdf
This is an UGLY HACK and should be removed as soon as the problem
is fixed upstream.
"""
cdef:
double a = x.real
double b = x.imag
double c = y.real
double d = y.imag
double ratio, denom
double complex out
if libc.math.fabs(d) < libc.math.fabs(c):
ratio = d/c
denom = c + d*ratio
out.real = (a + b*ratio)/denom
out.imag = (b - a*ratio)/denom
else:
ratio = c/d
denom = c*ratio + d
out.real = (a*ratio + b)/denom
out.imag = (b*ratio - a)/denom
return out
@cython.cdivision(True)
cdef inline double complex zlog1(double complex z) nogil:
"""
Compute log, paying special attention to accuracy around 1. We
implement this ourselves because some systems (most notably the
Travis CI machines) are weak in this regime.
"""
cdef:
int n
double complex coeff = -1
double complex res = 0
if zabs(z - 1) > 0.1:
return zlog(z)
z = z - 1
if z == 0:
return 0
for n in range(1, 17):
coeff *= -z
res += coeff/n
if zabs(res/coeff) < tol:
break
return res
|
