Skip to content

Commit aba8c8a

Browse files
musmmusm
authored
Use Refs for passing parameters by reference (#44)
1 parent 10368b0 commit aba8c8a

File tree

2 files changed

+94
-84
lines changed

2 files changed

+94
-84
lines changed

‎src/bessel.jl

Lines changed: 93 additions & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -28,37 +28,37 @@ end
2828
Base.showerror(io::IO, err::AmosException) = print(io, "AmosException with id $(err.id): $(err.msg)")
2929

3030
## Airy functions
31-
let
32-
const ai::Array{Float64,1} = Array{Float64}(2)
33-
const ae::Array{Int32,1} = Array{Int32}(2)
34-
global _airy, _biry
35-
function _airy(z::Complex128, id::Int32, kode::Int32)
36-
ccall((:zairy_,openspecfun), Void,
37-
(Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
38-
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
39-
&real(z), &imag(z),
40-
&id, &kode,
41-
pointer(ai,1), pointer(ai,2),
42-
pointer(ae,1), pointer(ae,2))
43-
if ae[2] == 0 || ae[2] == 3
44-
return complex(ai[1],ai[2])
45-
else
46-
throw(AmosException(ae[2]))
47-
end
31+
function _airy(z::Complex128, id::Int32, kode::Int32)
32+
ai1, ai2 = Ref{Float64}(), Ref{Float64}()
33+
ae1, ae2 = Ref{Int32}(), Ref{Int32}()
34+
35+
ccall((:zairy_,openspecfun), Void,
36+
(Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
37+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
38+
real(z), imag(z), id, kode,
39+
ai1, ai2, ae1, ae2)
40+
41+
if ae2[] == 0 || ae2[] == 3 # ignore underflow and less than half machine accuracy loss
42+
return complex(ai1[], ai2[])
43+
else
44+
throw(AmosException(ae2[]))
4845
end
49-
function _biry(z::Complex128, id::Int32, kode::Int32)
50-
ccall((:zbiry_,openspecfun), Void,
51-
(Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
52-
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}),
53-
&real(z), &imag(z),
54-
&id, &kode,
55-
pointer(ai,1), pointer(ai,2),
56-
pointer(ae,1))
57-
if ae[1] == 0 || ae[1] == 3 # ignore underflow
58-
return complex(ai[1],ai[2])
59-
else
60-
throw(AmosException(ae[2]))
61-
end
46+
end
47+
48+
function _biry(z::Complex128, id::Int32, kode::Int32)
49+
ai1, ai2 = Ref{Float64}(), Ref{Float64}()
50+
ae1 = Ref{Int32}()
51+
52+
ccall((:zbiry_,openspecfun), Void,
53+
(Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
54+
Ref{Float64}, Ref{Float64}, Ref{Int32}),
55+
real(z), imag(z), id, kode,
56+
ai1, ai2, ae1)
57+
58+
if ae1[] == 0 || ae1[] == 3 # ignore less than half machine accuracy loss
59+
return complex(ai1[], ai2[])
60+
else
61+
throw(AmosException(ae1[]))
6262
end
6363
end
6464

@@ -145,7 +145,7 @@ end
145145

146146
function airyai(x::BigFloat)
147147
z = BigFloat()
148-
ccall((:mpfr_ai, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
148+
ccall((:mpfr_ai, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
149149
return z
150150
end
151151

@@ -174,79 +174,89 @@ for jy in ("j","y"), nu in (0,1)
174174
end
175175

176176

177-
const cy = Array{Float64}(2)
178-
const ae = Array{Int32}(2)
179-
const wrk = Array{Float64}(2)
180-
181177
function _besselh(nu::Float64, k::Int32, z::Complex128, kode::Int32)
178+
ai1, ai2 = Ref{Float64}(), Ref{Float64}()
179+
ae1, ae2 = Ref{Int32}(), Ref{Int32}()
180+
182181
ccall((:zbesh_,openspecfun), Void,
183-
(Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
184-
Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
185-
&real(z), &imag(z), &nu, &kode, &k, &1,
186-
pointer(cy,1), pointer(cy,2),
187-
pointer(ae,1), pointer(ae,2))
188-
if ae[2] == 0 || ae[2] == 3
189-
return complex(cy[1],cy[2])
182+
(Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}, Ref{Int},
183+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
184+
real(z), imag(z), nu, kode, k, 1,
185+
ai1, ai2, ae1, ae2)
186+
187+
if ae2[] == 0 || ae2[] == 3
188+
return complex(ai1[], ai2[])
190189
else
191-
throw(AmosException(ae[2]))
190+
throw(AmosException(ae2[]))
192191
end
193192
end
194193

195194
function _besseli(nu::Float64, z::Complex128, kode::Int32)
195+
ai1, ai2 = Ref{Float64}(), Ref{Float64}()
196+
ae1, ae2 = Ref{Int32}(), Ref{Int32}()
197+
196198
ccall((:zbesi_,openspecfun), Void,
197-
(Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
198-
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
199-
&real(z), &imag(z), &nu, &kode, &1,
200-
pointer(cy,1), pointer(cy,2),
201-
pointer(ae,1), pointer(ae,2))
202-
if ae[2] == 0 || ae[2] == 3
203-
return complex(cy[1],cy[2])
199+
(Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
200+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
201+
real(z), imag(z), nu, kode, 1,
202+
ai1, ai2, ae1, ae2)
203+
204+
if ae2[] == 0 || ae2[] == 3
205+
return complex(ai1[], ai2[])
204206
else
205-
throw(AmosException(ae[2]))
207+
throw(AmosException(ae2[]))
206208
end
207209
end
208210

209211
function _besselj(nu::Float64, z::Complex128, kode::Int32)
212+
ai1, ai2 = Ref{Float64}(), Ref{Float64}()
213+
ae1, ae2 = Ref{Int32}(), Ref{Int32}()
214+
210215
ccall((:zbesj_,openspecfun), Void,
211-
(Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
212-
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
213-
&real(z), &imag(z), &nu, &kode, &1,
214-
pointer(cy,1), pointer(cy,2),
215-
pointer(ae,1), pointer(ae,2))
216-
if ae[2] == 0 || ae[2] == 3
217-
return complex(cy[1],cy[2])
216+
(Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
217+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
218+
real(z), imag(z), nu, kode, 1,
219+
ai1, ai2, ae1, ae2)
220+
221+
if ae2[] == 0 || ae2[] == 3
222+
return complex(ai1[], ai2[])
218223
else
219-
throw(AmosException(ae[2]))
224+
throw(AmosException(ae2[]))
220225
end
221226
end
222227

223228
function _besselk(nu::Float64, z::Complex128, kode::Int32)
229+
ai1, ai2 = Ref{Float64}(), Ref{Float64}()
230+
ae1, ae2 = Ref{Int32}(), Ref{Int32}()
231+
224232
ccall((:zbesk_,openspecfun), Void,
225-
(Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
226-
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
227-
&real(z), &imag(z), &nu, &kode, &1,
228-
pointer(cy,1), pointer(cy,2),
229-
pointer(ae,1), pointer(ae,2))
230-
if ae[2] == 0 || ae[2] == 3
231-
return complex(cy[1],cy[2])
233+
(Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
234+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
235+
real(z), imag(z), nu, kode, 1,
236+
ai1, ai2, ae1, ae2)
237+
238+
if ae2[] == 0 || ae2[] == 3
239+
return complex(ai1[], ai2[])
232240
else
233-
throw(AmosException(ae[2]))
241+
throw(AmosException(ae2[]))
234242
end
235243
end
236244

237245
function _bessely(nu::Float64, z::Complex128, kode::Int32)
246+
ai1, ai2 = Ref{Float64}(), Ref{Float64}()
247+
ae1, ae2 = Ref{Int32}(), Ref{Int32}()
248+
wrk1, wrk2 = Ref{Float64}(), Ref{Float64}()
249+
238250
ccall((:zbesy_,openspecfun), Void,
239-
(Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32},
240-
Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32},
241-
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}),
242-
&real(z), &imag(z), &nu, &kode, &1,
243-
pointer(cy,1), pointer(cy,2),
244-
pointer(ae,1), pointer(wrk,1),
245-
pointer(wrk,2), pointer(ae,2))
246-
if ae[2] == 0 || ae[2] == 3
247-
return complex(cy[1],cy[2])
251+
(Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
252+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Int32}),
253+
real(z), imag(z), nu, kode, 1,
254+
ai1, ai2, ae1, wrk1, wrk2, ae2)
255+
256+
if ae2[] == 0 || ae2[] == 3
257+
return complex(ai1[], ai2[])
248258
else
249-
throw(AmosException(ae[2]))
259+
throw(AmosException(ae2[]))
250260
end
251261
end
252262

@@ -520,7 +530,7 @@ Bessel function of the first kind of order 0, ``J_0(x)``.
520530
"""
521531
function besselj0(x::BigFloat)
522532
z = BigFloat()
523-
ccall((:mpfr_j0, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
533+
ccall((:mpfr_j0, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
524534
return z
525535
end
526536

@@ -531,13 +541,13 @@ Bessel function of the first kind of order 1, ``J_1(x)``.
531541
"""
532542
function besselj1(x::BigFloat)
533543
z = BigFloat()
534-
ccall((:mpfr_j1, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
544+
ccall((:mpfr_j1, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
535545
return z
536546
end
537547

538548
function besselj(n::Integer, x::BigFloat)
539549
z = BigFloat()
540-
ccall((:mpfr_jn, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, n, &x, ROUNDING_MODE[])
550+
ccall((:mpfr_jn, :libmpfr), Int32, (Ref{BigFloat}, Clong, Ref{BigFloat}, Int32), z, n, x, ROUNDING_MODE[])
541551
return z
542552
end
543553

@@ -551,7 +561,7 @@ function bessely0(x::BigFloat)
551561
throw(DomainError(x, "`x` must be nonnegative."))
552562
end
553563
z = BigFloat()
554-
ccall((:mpfr_y0, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
564+
ccall((:mpfr_y0, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
555565
return z
556566
end
557567

@@ -565,7 +575,7 @@ function bessely1(x::BigFloat)
565575
throw(DomainError(x, "`x` must be nonnegative."))
566576
end
567577
z = BigFloat()
568-
ccall((:mpfr_y1, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
578+
ccall((:mpfr_y1, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
569579
return z
570580
end
571581

@@ -574,7 +584,7 @@ function bessely(n::Integer, x::BigFloat)
574584
throw(DomainError(x, "`x` must be nonnegative."))
575585
end
576586
z = BigFloat()
577-
ccall((:mpfr_yn, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, n, &x, ROUNDING_MODE[])
587+
ccall((:mpfr_yn, :libmpfr), Int32, (Ref{BigFloat}, Clong, Ref{BigFloat}, Int32), z, n, x, ROUNDING_MODE[])
578588
return z
579589
end
580590

‎src/erf.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ for f in (:erf, :erfc)
1212
($f)(a::Complex32) = Complex32($f(Complex64(a)))
1313
function ($f)(x::BigFloat)
1414
z = BigFloat()
15-
ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
15+
ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
1616
return z
1717
end
1818
($f)(x::AbstractFloat) = error("not implemented for ", typeof(x))

0 commit comments

Comments
 (0)