Skip to content

Commit f0b0e08

Browse files
committed
Fix JET analysis (#502)
1 parent fddf5e0 commit f0b0e08

File tree

2 files changed

+36
-33
lines changed

2 files changed

+36
-33
lines changed

src/logabsgamma/e_lgamma_r.jl

Lines changed: 18 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -92,32 +92,37 @@ function _logabsgamma(x::Float64)
9292
lx = ux % UInt32
9393

9494
#= purge off +-inf, NaN, +-0, tiny and negative arguments =#
95-
signgam = 1
96-
isneg = hx < Int32(0)
9795
ix = hx & 0x7fffffff
98-
ix 0x7ff00000 && return x * x, signgam
99-
ix | lx == 0x00000000 && return Inf, signgam
96+
ix 0x7ff00000 && return x * x, 1
97+
ix | lx == 0x00000000 && return Inf, 1
98+
isneg = hx < Int32(0)
10099
if ix < 0x3b900000 #= |x|<2**-70, return -log(|x|) =#
101100
if isneg
102-
signgam = -1
103-
return -log(-x), signgam
101+
return -log(-x), -1
104102
else
105-
return -log(x), signgam
103+
return -log(x), 1
106104
end
107105
end
106+
107+
# remaining cases
108108
if isneg
109109
# ix ≥ 0x43300000 && return Inf, signgam #= |x|>=2**52, must be -integer =#
110110
t = sinpi(x)
111-
iszero(t) && return Inf, signgam #= -integer =#
112-
nadj = logπ - log(abs(t * x))
113-
if t < 0.0; signgam = -1; end
114-
x = -x
111+
iszero(t) && return Inf, 1 #= -integer =#
112+
r = logπ - log(abs(t * x)) - _logabsgamma_pos(-x, ix)
113+
signgam = t < 0.0 ? -1 : 1
114+
else
115+
r = _logabsgamma_pos(x, ix)
116+
signgam = 1
115117
end
118+
return r, signgam
119+
end
120+
function _logabsgamma_pos(x::Float64, ix::UInt32)
116121
if ix < 0x40000000 #= x < 2.0 =#
117122
i = round(x, RoundToZero)
118123
f = x - i
119124
if f == 0.0 #= purge off 1; 2 handled by x < 8.0 branch =#
120-
return 0.0, signgam
125+
return 0.0
121126
elseif i == 1.0
122127
r = 0.0
123128
c = 1.0
@@ -169,8 +174,5 @@ function _logabsgamma(x::Float64)
169174
else #= 2^58 ≤ x ≤ Inf =#
170175
r = muladd(x, log(x), -x)
171176
end
172-
if isneg
173-
r = nadj - r
174-
end
175-
return r, signgam
177+
return r
176178
end

src/logabsgamma/e_lgammaf_r.jl

Lines changed: 18 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -20,33 +20,37 @@ function _logabsgamma(x::Float32)
2020
hx = reinterpret(Int32, x)
2121

2222
#= purge off +-inf, NaN, +-0, tiny and negative arguments =#
23-
signgam = 1
24-
isneg = hx < Int32(0)
2523
ix = hx & 0x7fffffff
26-
ix 0x7f800000 && return x * x, signgam
27-
ix == 0x00000000 && return Inf32, signgam
24+
ix 0x7f800000 && return x * x, 1
25+
ix == 0x00000000 && return Inf32, 1
26+
isneg = hx < Int32(0)
2827
if ix < 0x35000000 #= |x|<2**-21, return -log(|x|) =#
2928
if isneg
30-
signgam = -1
31-
return -log(-x), signgam
29+
return -log(-x), -1
3230
else
33-
return -log(x), signgam
31+
return -log(x), 1
3432
end
3533
end
34+
35+
# remaining cases
3636
if isneg
3737
# ix ≥ 0x4b000000 && return Inf32, signgam #= |x|>=2**23, must be -integer =#
3838
t = sinpi(x)
39-
t == 0.0f0 && return Inf32, signgam #= -integer =#
40-
nadj = logπ - log(abs(t * x))
41-
if t < 0.0f0; signgam = -1; end
42-
x = -x
39+
t == 0.0f0 && return Inf32, 1 #= -integer =#
40+
r = logπ - log(abs(t * x)) - _logabsgamma_pos(-x, ix)
41+
signgam = copysign(1, t)
42+
else
43+
r = _logabsgamma_pos(x, ix)
44+
signgam = 1
4345
end
44-
46+
return r, signgam
47+
end
48+
function _logabsgamma_pos(x::Float32, ix::UInt32)
4549
if ix < 0x40000000 #= x < 2.0 =#
4650
i = round(x, RoundToZero)
4751
f = x - i
4852
if f == 0.0f0 #= purge off 1; 2 handled by x < 8.0 branch =#
49-
return 0.0f0, signgam
53+
return 0.0f0
5054
elseif i == 1.0f0
5155
r = 0.0f0
5256
c = 1.0f0
@@ -99,8 +103,5 @@ function _logabsgamma(x::Float32)
99103
#= 2^58 ≤ x ≤ Inf =#
100104
r = muladd(x, log(x), -x)
101105
end
102-
if isneg
103-
r = nadj - r
104-
end
105-
return r, signgam
106+
return r
106107
end

0 commit comments

Comments
 (0)