Fixes to norm and friends (closes #915) (#929)

* various fixes to `norm` & more tests

* fix test failure; widen signatures

* restrict type-stability tests to v1.2+

* fixes/improvements from code-review

* address remaining issue for `_inner_eltype` of empty arrays

* use `_inner_eltype(x)` for `_norm_p0` and adjust signature to `AbstractArray`

- also fix copying errors from last commit

* bump patch version

Author: thchr

Project.toml

@@ -1,6 +1,6 @@
 name = "StaticArrays"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.2.7"
+version = "1.2.8"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/linalg.jl

@@ -208,17 +208,23 @@ end
 
 #--------------------------------------------------
 # Norms
-@inline LinearAlgebra.norm_sqr(v::StaticVector) = mapreduce(abs2, +, v; init=zero(real(eltype(v))))
+_inner_eltype(v::AbstractArray) = isempty(v) ? eltype(v) : _inner_eltype(first(v))
+_inner_eltype(x::Number) = typeof(x)
+@inline _init_zero(v::StaticArray) = float(norm(zero(_inner_eltype(v))))
+
+@inline function LinearAlgebra.norm_sqr(v::StaticArray)
+    return mapreduce(LinearAlgebra.norm_sqr, +, v; init=_init_zero(v))
+end
 
 @inline norm(a::StaticArray) = _norm(Size(a), a)
 @generated function _norm(::Size{S}, a::StaticArray) where {S}
     if prod(S) == 0
-        return :(zero(real(eltype(a))))
+        return :(_init_zero(a))
     end
 
-    expr = :(abs2(a[1]))
+    expr = :(LinearAlgebra.norm_sqr(a[1]))
     for j = 2:prod(S)
-        expr = :($expr + abs2(a[$j]))
+        expr = :($expr + LinearAlgebra.norm_sqr(a[$j]))
     end
 
     return quote
@@ -227,28 +233,31 @@ end
     end
 end
 
-_norm_p0(x) = x == 0 ? zero(x) : one(x)
+function _norm_p0(x)
+    T = _inner_eltype(x)
+    return float(norm(iszero(x) ? zero(T) : one(T)))
+end
 
 @inline norm(a::StaticArray, p::Real) = _norm(Size(a), a, p)
 @generated function _norm(::Size{S}, a::StaticArray, p::Real) where {S}
     if prod(S) == 0
-        return :(zero(real(eltype(a))))
+        return :(_init_zero(a))
     end
 
-    expr = :(abs(a[1])^p)
+    expr = :(norm(a[1])^p)
     for j = 2:prod(S)
-        expr = :($expr + abs(a[$j])^p)
+        expr = :($expr + norm(a[$j])^p)
     end
 
-    expr_p1 = :(abs(a[1]))
+    expr_p1 = :(norm(a[1]))
     for j = 2:prod(S)
-        expr_p1 = :($expr_p1 + abs(a[$j]))
+        expr_p1 = :($expr_p1 + norm(a[$j]))
     end
 
     return quote
         $(Expr(:meta, :inline))
         if p == Inf
-            return mapreduce(abs, max, a)
+            return mapreduce(norm, max, a)
         elseif p == 1
             @inbounds return $expr_p1
         elseif p == 2

test/linalg.jl

@@ -202,6 +202,57 @@ StaticArrays.similar_type(::Union{RotMat2,Type{RotMat2}}) = SMatrix{2,2,Float64,
         @test normalize(SVector(1,2,3), 1) ≈ normalize([1,2,3], 1)
         @test normalize!(MVector(1.,2.,3.)) ≈ normalize([1.,2.,3.])
         @test normalize!(MVector(1.,2.,3.), 1) ≈ normalize([1.,2.,3.], 1)
+
+        # nested vectors
+        a  = SA[SA[1, 2], SA[3, 4]]
+        av = convert(Vector{Vector{Int}}, a)
+        aa = SA[a,a]
+        @test norm(a) ≈ norm(av)
+        @test norm(aa) ≈ norm([a,a])
+        @test norm(aa) ≈ norm([av,av])
+        @test norm(SVector{0,Int}()) === norm(Vector{Float64}()) === 0.0
+
+        # do not overflow for Int
+        c = SA[typemax(Int)÷2, typemax(Int)÷3]
+        @test norm(c) ≈ norm(Vector(c))
+        @test norm(SA[c,c]) ≈ norm([Vector(c), Vector(c)])
+
+        # 0-norm of vectors w/ zero-vectors
+        @test norm(SA[0,0], 0) == norm([0,0], 0)
+        @test norm(SA[[0,0],[1,1]], 0) == norm([[0,0],[1,1]], 0) == 1.0
+        @test norm(SA[[0,1],[1,1]], 0) == norm([[0,1],[1,1]], 0) == 2.0
+
+        # complex numbers
+        @test norm(SA[1+im, 2+3im]) ≈ norm([1+im, 2+3im])
+        @test norm(SA[22.0+0.1im, 2.0-23.0im]) ≈ norm([22.0+0.1im, 2.0-23.0im])
+        a_c, av_c = a+1im*a, av+1im*av
+        @test norm(a_c) ≈ norm(av_c)
+
+        # p-norms for nested vectors
+        for (x,xv) in ((a,av), (a_c, av_c))
+            @test norm(x, 2) ≈ norm(xv,2)
+            @test norm(x, Inf) ≈ norm(xv, Inf)
+            @test norm(x, 1) ≈ norm(xv, 1)
+            @test norm(x, 0) ≈ norm(xv, 0)
+            @test norm(SA[Int[], [1,2]], 0) ≈ norm([Int[], [1,2]], 0)
+        end
+
+        # type-stability
+        if VERSION ≥ v"1.2"
+            # only test strict type-stability on v1.2+, since there were Base-related type
+            # instabilities in `norm` prior to https://github.com/JuliaLang/julia/pull/30481
+            @test (@inferred norm(a[1])) == (@inferred norm(a[1], 2))
+            @test (@inferred norm(a)) == (@inferred norm(a, 2))
+            @test (@inferred norm(aa)) == (@inferred norm(aa, 2))
+            @test (@inferred norm(float.(aa))) == (@inferred norm(float.(aa), 2))
+            @test (@inferred norm(SVector{0,Int}())) == (@inferred norm(SVector{0,Int}(), 2))
+        end
+
+        # norm of empty SVector
+        @test norm(SVector{0,Int}()) isa float(Int)
+        @test norm(SVector{0,Float64}()) isa Float64
+        @test norm(SA[SVector{0,Int}(),SVector{0,Int}()]) isa float(Int)
+        @test norm(SA[SVector{0,Int}(),SVector{0,Int}()]) == norm([Int[], Int[]])
     end
 
     @testset "trace" begin