Avoid overflow in generic_vecnorm()

No need to take the inverse too early. Fixes #11788.
JuliaLang · Jun 20, 2015 · f110415 · f110415
1 parent dc3f7d6
commit f110415
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 6 deletions.
diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl
@@ -107,12 +107,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))
@@ -128,11 +128,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

diff --git a/test/linalg2.jl b/test/linalg2.jl
@@ -333,6 +333,11 @@ end
 # overflow/underflow in norms:
 @test_approx_eq norm(Float64[1e-300, 1], -3)*1e300 1
 @test_approx_eq norm(Float64[1e300, 1], 3)*1e-300 1
+@test_approx_eq norm(Float64[1e300, 1], 3)*1e-300 1
+
+# 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
 
 # Uniform scaling
 @test I[1,1] == 1 # getindex