Skip to content

Commit 521bda1

Browse files
committed
additional branch nudge tests
1 parent 3ebc128 commit 521bda1

1 file changed

Lines changed: 44 additions & 29 deletions

File tree

test/gamma.jl

Lines changed: 44 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -175,39 +175,54 @@
175175
@test loggamma(Complex{BigFloat}(big"3.099", big"0.0")) loggamma(big"3.099")
176176
@test loggamma(Complex{BigFloat}(big"1.15", big"0.0")) loggamma(big"1.15")
177177
@test gamma(Complex{BigFloat}(big"3.099", big"0.0")) gamma(big"3.099")
178+
# branch mapping
179+
ε = BigFloat(2)^(-1400)
180+
xs = BigFloat.(["-0.1", "-0.7", "-1.3", "-2.6", "-3.4", "-5.2"])
181+
for x in xs
182+
z_minus = Complex{BigFloat}(x, -ε)
183+
z_plus = Complex{BigFloat}(x, ε)
184+
185+
Lm = loggamma(z_minus)
186+
Lp = loggamma(z_plus)
187+
188+
zf_minus = _loggamma_oracle64_point(z_minus)
189+
zf_plus = _loggamma_oracle64_point(z_plus)
190+
Lf_minus = loggamma(zf_minus)
191+
Lf_plus = loggamma(zf_plus)
192+
193+
@test isapprox(Float64(real(Lm)), real(Lf_minus); rtol=0, atol=1e-14)
194+
@test isapprox(Float64(imag(Lm)), imag(Lf_minus); rtol=0, atol=1e-14)
195+
@test isapprox(Float64(real(Lp)), real(Lf_plus); rtol=0, atol=1e-14)
196+
@test isapprox(Float64(imag(Lp)), imag(Lf_plus); rtol=0, atol=1e-14)
197+
198+
kdiff = round(Int, (imag(Lm) - imag(Lp)) / (2*big(pi)))
199+
@test kdiff != 0
200+
@test isapprox(Lm - Lp, (2*big(pi))*BigFloat(kdiff)*im; rtol=0, atol=big"1e-70")
201+
end
202+
# Additional _loggamma_oracle64_point nudge tests at Float64 pole
203+
n = -3
204+
x64 = Float64(n)
205+
δp64 = nextfloat(x64) - x64
206+
δm64 = x64 - prevfloat(x64)
207+
ε64 = eps(Float64)
208+
209+
# right of pole: rez > n but Float64(rez) == n → nextfloat
210+
rezR = BigFloat(x64) + BigFloat(δp64)/4
211+
zR = Complex{BigFloat}(rezR, BigFloat(ε64)/4)
212+
zfR = _loggamma_oracle64_point(zR)
213+
@test real(zfR) == nextfloat(x64)
214+
@test abs(imag(zfR)) 2ε64
215+
216+
# left of pole: rez < n but Float64(rez) == n → prevfloat
217+
rezL = BigFloat(x64) - BigFloat(δm64)/4
218+
zL = Complex{BigFloat}(rezL, -BigFloat(ε64)/4)
219+
zfL = _loggamma_oracle64_point(zL)
220+
@test real(zfL) == prevfloat(x64)
221+
@test abs(imag(zfL)) 2ε64
178222
end
179223
end
180224
end
181225

182-
@testset "BigFloat branch mapping near the negative real axis" begin
183-
setprecision(256) do
184-
ε = BigFloat(2)^(-1400)
185-
xs = BigFloat.(["-0.1", "-0.7", "-1.3", "-2.6", "-3.4", "-5.2"])
186-
187-
for x in xs
188-
z_minus = Complex{BigFloat}(x, -ε)
189-
z_plus = Complex{BigFloat}(x, ε)
190-
191-
Lm = loggamma(z_minus)
192-
Lp = loggamma(z_plus)
193-
194-
zf_minus = _loggamma_oracle64_point(z_minus)
195-
zf_plus = _loggamma_oracle64_point(z_plus)
196-
Lf_minus = loggamma(zf_minus)
197-
Lf_plus = loggamma(zf_plus)
198-
199-
@test isapprox(Float64(real(Lm)), real(Lf_minus); rtol=0, atol=1e-14)
200-
@test isapprox(Float64(imag(Lm)), imag(Lf_minus); rtol=0, atol=1e-14)
201-
@test isapprox(Float64(real(Lp)), real(Lf_plus); rtol=0, atol=1e-14)
202-
@test isapprox(Float64(imag(Lp)), imag(Lf_plus); rtol=0, atol=1e-14)
203-
204-
kdiff = round(Int, (imag(Lm) - imag(Lp)) / (2*big(pi)))
205-
@test kdiff != 0
206-
@test isapprox(Lm - Lp, (2*big(pi))*BigFloat(kdiff)*im; rtol=0, atol=big"1e-70")
207-
end
208-
end
209-
end
210-
211226
@testset "Other float types" begin
212227
struct NotAFloat <: AbstractFloat
213228
end

0 commit comments

Comments
 (0)