Skip to content

Commit

Permalink
Avoid overflow in generic_vecnorm()
Browse files Browse the repository at this point in the history
No need to take the inverse too early. Fixes #11788.
  • Loading branch information
nalimilan authored and andreasnoack committed Aug 11, 2015
1 parent fdf5295 commit 2996f74
Show file tree
Hide file tree
Showing 2 changed files with 10 additions and 6 deletions.
12 changes: 6 additions & 6 deletions base/linalg/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -109,12 +109,12 @@ function generic_vecnorm2(x)
s = start(x)
(v, s) = next(x, s)
T = typeof(maxabs)
scale::promote_type(Float64, T) = 1/maxabs
y = norm(v)*scale
scale::promote_type(Float64, T) = maxabs
y = norm(v)/scale
sum::promote_type(Float64, T) = y*y
while !done(x, s)
(v, s) = next(x, s)
y = norm(v)*scale
y = norm(v)/scale
sum += y*y
end
return convert(T, maxabs * sqrt(sum))
Expand All @@ -130,11 +130,11 @@ function generic_vecnormp(x, p)
(v, s) = next(x, s)
T = typeof(maxabs)
spp::promote_type(Float64, T) = p
scale::promote_type(Float64, T) = 1/maxabs
ssum::promote_type(Float64, T) = (norm(v)*scale)^spp
scale::promote_type(Float64, T) = maxabs
ssum::promote_type(Float64, T) = (norm(v)/scale)^spp
while !done(x, s)
(v, s) = next(x, s)
ssum += (norm(v)*scale)^spp
ssum += (norm(v)/scale)^spp
end
return convert(T, maxabs * ssum^inv(spp))
else # -1 ≤ p ≤ 1, no need for rescaling
Expand Down
4 changes: 4 additions & 0 deletions test/linalg/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -106,3 +106,7 @@ b = randn(Base.LinAlg.SCAL_CUTOFF) # make sure we try BLAS path
@test isequal(scale(Float64[1.0], big(2.0)im), Complex{BigFloat}[2.0im])
@test isequal(scale(BigFloat[1.0], 2.0im), Complex{BigFloat}[2.0im])
@test isequal(scale(BigFloat[1.0], 2.0f0im), Complex{BigFloat}[2.0im])

# Issue #11788
@test_approx_eq norm([2.4e-322, 4.4e-323]) 2.47e-322
@test_approx_eq norm([2.4e-322, 4.4e-323], 3) 2.4e-322

0 comments on commit 2996f74

Please sign in to comment.