diff --git a/NEWS.md b/NEWS.md
index 41eb7873..578556d0 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,15 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.15.18] 18/10/2024
+
+### Added
+
+* `distance(M, p, q, r)` to compute `r`-norms on manifolds that have components.
+* `distance(M, p, q, m, r)` to compute (approximate) `r`-norms on manifolds that have components
+  using an `AbstractInverseRetractionMethod m` within every (inner) distance call.
+* `norm(M, p, X, r)` to compute `r`-norms on manifolds that have components.
+
 ## [0.15.17] 04/10/2024
 
 ### Changed
diff --git a/Project.toml b/Project.toml
index ad144c1a..3e73cc26 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ManifoldsBase"
 uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
 authors = ["Seth Axen <seth.axen@gmail.com>", "Mateusz Baran <mateuszbaran89@gmail.com>", "Ronny Bergmann <manopt@ronnybergmann.net>", "Antoine Levitt <antoine.levitt@gmail.com>"]
-version = "0.15.17"
+version = "0.15.18"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/DefaultManifold.jl b/src/DefaultManifold.jl
index 6bf2f41d..7817730f 100644
--- a/src/DefaultManifold.jl
+++ b/src/DefaultManifold.jl
@@ -57,7 +57,7 @@ function check_approx(M::DefaultManifold, p, X, Y; kwargs...)
     return ApproximatelyError(v, s)
 end
 
-distance(::DefaultManifold, p, q) = norm(p - q)
+distance(::DefaultManifold, p, q, r::Real = 2) = norm(p - q, r)
 
 embed!(::DefaultManifold, q, p) = copyto!(q, p)
 
@@ -123,6 +123,8 @@ function get_vector_orthonormal!(M::DefaultManifold{ℂ}, Y, p, c, ::RealNumbers
     return copyto!(Y, reshape(c[1:n] + c[(n + 1):(2n)] * 1im, representation_size(M)))
 end
 
+has_components(::DefaultManifold) = true
+
 injectivity_radius(::DefaultManifold) = Inf
 
 @inline inner(::DefaultManifold, p, X, Y) = dot(X, Y)
@@ -138,7 +140,7 @@ end
 
 number_system(::DefaultManifold{𝔽}) where {𝔽} = 𝔽
 
-norm(::DefaultManifold, p, X) = norm(X)
+norm(::DefaultManifold, p, X, r::Real = 2) = norm(X, r)
 
 project!(::DefaultManifold, q, p) = copyto!(q, p)
 project!(::DefaultManifold, Y, p, X) = copyto!(Y, X)
diff --git a/src/ManifoldsBase.jl b/src/ManifoldsBase.jl
index 80d5af26..cdecdd92 100644
--- a/src/ManifoldsBase.jl
+++ b/src/ManifoldsBase.jl
@@ -594,6 +594,15 @@ function embed_project!(M::AbstractManifold, Y, p, X)
     return project!(M, Y, p, embed(M, p, X))
 end
 
+"""
+    has_components(M::AbstractManifold)
+
+Return whether the [`AbstractManifold`](@ref)`(M)` consists of components,
+like the [`PowerManifold`](@ref) or the [`ProductManifold`](@ref), that one can iterate over.
+By default, this function returns `false`.
+"""
+has_components(M::AbstractManifold) = false
+
 @doc raw"""
     injectivity_radius(M::AbstractManifold)
 
@@ -1354,6 +1363,7 @@ export ×,
     get_vector,
     get_vector!,
     get_vectors,
+    has_components,
     hat,
     hat!,
     injectivity_radius,
diff --git a/src/PowerManifold.jl b/src/PowerManifold.jl
index c6938256..e0cec526 100644
--- a/src/PowerManifold.jl
+++ b/src/PowerManifold.jl
@@ -530,22 +530,132 @@ function default_vector_transport_method(M::PowerManifold, t::Type)
     return default_vector_transport_method(M.manifold, eltype(t))
 end
 
-@doc raw"""
-    distance(M::AbstractPowerManifold, p, q)
+_doc_distance_pow = """
+    distance(M::AbstractPowerManifold, p, q, r::Real=2)
+    distance(M::AbstractPowerManifold, p, q, m::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(), r::Real=2)
+
+Compute the distance between `q` and `p` on an [`AbstractPowerManifold`](@ref).
+
+First, the componentwise distances are computed using the Riemannian distance function
+on `M.manifold`. These can be approximated using the
+`norm` of an [`AbstractInverseRetractionMethod`](@ref) `m`.
+This yields an array of distance values.
 
-Compute the distance between `q` and `p` on an [`AbstractPowerManifold`](@ref),
-i.e. from the element wise distances the Forbenius norm is computed.
+Second, we compute the `r`-norm on this array of distances.
+This is also the only place, there the `r` is used.
 """
+
 function distance(M::AbstractPowerManifold, p, q)
-    sum_squares = zero(number_eltype(p))
+    return _distance_r(M, p, q, 2)
+end
+
+@doc "$(_doc_distance_pow)"
+function distance(M::AbstractPowerManifold, p, q, r::Real)
+    (isinf(r) && r > 0) && return _distance_max(M, p, q)
+    (isinf(r) && r < 0) && return _distance_min(M, p, q)
+    (r == 1) && return _distance_1(M, p, q)
+    return _distance_r(M, p, q, r)
+end
+function distance(
+    M::AbstractPowerManifold,
+    p,
+    q,
+    ::LogarithmicInverseRetraction,
+    r::Real = 2,
+)
+    return distance(M, p, q, r)
+end
+
+@doc "$(_doc_distance_pow)"
+function distance(
+    M::AbstractPowerManifold,
+    p,
+    q,
+    m::AbstractInverseRetractionMethod,
+    r::Real = 2,
+)
+    (isinf(r) && r > 0) && return _distance_max(M, p, q, m)
+    (isinf(r) && r < 0) && return _distance_min(M, p, q, m)
+    (r == 1) && return _distance_1(M, p, q, m)
+    return _distance_r(M, p, q, m, r)
+end
+#
+#
+# The three special cases
+function _distance_r(
+    M::AbstractPowerManifold,
+    p,
+    q,
+    m::AbstractInverseRetractionMethod,
+    r::Real,
+)
+    rep_size = representation_size(M.manifold)
+    values = [
+        distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i), m) for
+        i in get_iterator(M)
+    ]
+    return norm(values, r)
+end
+function _distance_r(M::AbstractPowerManifold, p, q, r::Real)
+    rep_size = representation_size(M.manifold)
+    values = [
+        distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i)) for
+        i in get_iterator(M)
+    ]
+    return norm(values, r)
+end
+function _distance_1(M::AbstractPowerManifold, p, q, m::AbstractInverseRetractionMethod)
+    s = zero(number_eltype(p))
     rep_size = representation_size(M.manifold)
     for i in get_iterator(M)
-        sum_squares +=
-            distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i))^2
+        s += distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i), m)
     end
-    return sqrt(sum_squares)
+    return s
+end
+function _distance_1(M::AbstractPowerManifold, p, q)
+    s = zero(number_eltype(p))
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        s += distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i))
+    end
+    return s
+end
+function _distance_max(M::AbstractPowerManifold, p, q, m::AbstractInverseRetractionMethod)
+    d = float(zero(number_eltype(p)))
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        v = distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i), m)
+        d = (v > d) ? v : d
+    end
+    return d
+end
+function _distance_max(M::AbstractPowerManifold, p, q)
+    d = float(zero(number_eltype(p)))
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        v = distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i))
+        d = (v > d) ? v : d
+    end
+    return d
+end
+function _distance_min(M::AbstractPowerManifold, p, q, m::AbstractInverseRetractionMethod)
+    d = Inf
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        v = distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i), m)
+        d = (v < d) ? v : d
+    end
+    return d
+end
+function _distance_min(M::AbstractPowerManifold, p, q)
+    d = Inf
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        v = distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i))
+        d = (v < d) ? v : d
+    end
+    return d
 end
-
 @doc raw"""
     exp(M::AbstractPowerManifold, p, X)
 
@@ -881,6 +991,13 @@ function Base.getindex(
     return TangentSpace(M.manifold, p[M, I...])
 end
 
+"""
+    has_components(::AbstractPowerManifold)
+
+Return `true`, since points on an [`AbstractPowerManifold`](@ref) consist of components.
+"""
+has_components(::AbstractPowerManifold) = true
+
 @doc raw"""
     injectivity_radius(M::AbstractPowerManifold[, p])
 
@@ -1092,20 +1209,51 @@ function mid_point!(M::PowerManifoldNestedReplacing, q, p1, p2)
 end
 
 @doc raw"""
-    norm(M::AbstractPowerManifold, p, X)
+    norm(M::AbstractPowerManifold, p, X, r::Real=2)
 
 Compute the norm of `X` from the tangent space of `p` on an
-[`AbstractPowerManifold`](@ref) `M`, i.e. from the element wise norms the
-Frobenius norm is computed.
+[`AbstractPowerManifold`](@ref) `M`, i.e. from the element wise norms `r`-norm is computed,
+where the default `r=2` yields the Frobenius norm is computed.
 """
-function LinearAlgebra.norm(M::AbstractPowerManifold, p, X)
-    sum_squares = zero(number_eltype(X))
+function LinearAlgebra.norm(M::AbstractPowerManifold, p, X, r::Real = 2)
+    (isinf(r) && r > 0) && return _norm_max(M, p, X)
+    (isinf(r) && r < 0) && return _norm_min(M, p, X)
+    (r == 1) && return _norm_1(M, p, X)
+    return _norm_r(M, p, X, r)
+end
+function _norm_r(M::AbstractPowerManifold, p, X, r::Real)
+    rep_size = representation_size(M.manifold)
+    values = [
+        norm(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, X, i)) for
+        i in get_iterator(M)
+    ]
+    return norm(values, r)
+end
+function _norm_1(M::AbstractPowerManifold, p, X)
+    s = zero(number_eltype(p))
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        s += norm(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, X, i))
+    end
+    return s
+end
+function _norm_max(M::AbstractPowerManifold, p, X)
+    d = 0.0
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        v = norm(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, X, i))
+        d = (v > d) ? v : d
+    end
+    return d
+end
+function _norm_min(M::AbstractPowerManifold, p, X)
+    d = Inf
     rep_size = representation_size(M.manifold)
     for i in get_iterator(M)
-        sum_squares +=
-            norm(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, X, i))^2
+        v = norm(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, X, i))
+        d = (v < d) ? v : d
     end
-    return sqrt(sum_squares)
+    return d
 end
 
 
diff --git a/src/ProductManifold.jl b/src/ProductManifold.jl
index cb8eae69..fdd86521 100644
--- a/src/ProductManifold.jl
+++ b/src/ProductManifold.jl
@@ -372,25 +372,57 @@ function default_vector_transport_method(M::ProductManifold)
     return ProductVectorTransport(map(default_vector_transport_method, M.manifolds)...)
 end
 
-@doc raw"""
-    distance(M::ProductManifold, p, q)
-
-Compute the distance between two points `p` and `q` on the [`ProductManifold`](@ref) `M`, which is
-the 2-norm of the elementwise distances on the internal manifolds that build `M`.
-"""
-function distance(M::ProductManifold, p, q)
-    return sqrt(
-        sum(
-            map(
-                distance,
-                M.manifolds,
-                submanifold_components(M, p),
-                submanifold_components(M, q),
-            ) .^ 2,
+_doc_distance_prod = """
+    distance(M::ProductManifold, p, q, r::Real=2)
+    distance(M::ProductManifold, p, q, m::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(), r::Real=2)
+
+Compute the distance between `q` and `p` on an [`ProductManifold`](@ref).
+
+First, the componentwise distances are computed. These can be approximated using the
+`norm` of an [`AbstractInverseRetractionMethod`](@ref) `m`.
+Then, the `r`-norm of the tuple of these elements is computed.
+"""
+
+@doc "$(_doc_distance_prod)"
+function distance(M::ProductManifold, p, q, r::Real = 2)
+    return norm(
+        map(
+            distance,
+            M.manifolds,
+            submanifold_components(M, p),
+            submanifold_components(M, q),
         ),
+        r,
     )
 end
+function distance(M::ProductManifold, p, q, ::LogarithmicInverseRetraction, r::Real = 2)
+    return distance(M, p, q, r)
+end
 
+@doc "$(_doc_distance_prod)"
+function distance(M::ProductManifold, p, q, m::AbstractInverseRetractionMethod, r::Real = 2)
+    return norm(
+        map(
+            (M, p, q) -> distance(M, p, q, m),
+            M.manifolds,
+            submanifold_components(M, p),
+            submanifold_components(M, q),
+        ),
+        r,
+    )
+end
+function distance(M::ProductManifold, p, q, method::InverseProductRetraction, r::Real = 2)
+    return norm(
+        map(
+            (M, p, q, m) -> distance(M, p, q, m),
+            M.manifolds,
+            submanifold_components(M, p),
+            submanifold_components(M, q),
+            method.inverse_retractions,
+        ),
+        r,
+    )
+end
 @doc raw"""
     exp(M::ProductManifold, p, X)
 
@@ -557,6 +589,13 @@ function get_vector!(
     return X
 end
 
+"""
+    has_components(::ProductManifold)
+
+Return `true` since points on an [`ProductManifold`](@ref) consist of components.
+"""
+has_components(::ProductManifold) = true
+
 @doc raw"""
     injectivity_radius(M::ProductManifold)
     injectivity_radius(M::ProductManifold, x)
@@ -725,21 +764,15 @@ function mid_point!(M::ProductManifold, q, p1, p2)
 end
 
 @doc raw"""
-    norm(M::ProductManifold, p, X)
+    norm(M::ProductManifold, p, X, r::Real=2)
 
-Compute the norm of `X` from the tangent space of `p` on the [`ProductManifold`](@ref),
+Compute the (`r`-)norm of `X` from the tangent space of `p` on the [`ProductManifold`](@ref),
 i.e. from the element wise norms the 2-norm is computed.
 """
-function LinearAlgebra.norm(M::ProductManifold, p, X)
-    norms_squared = (
-        map(
-            norm,
-            M.manifolds,
-            submanifold_components(M, p),
-            submanifold_components(M, X),
-        ) .^ 2
-    )
-    return sqrt(sum(norms_squared))
+function LinearAlgebra.norm(M::ProductManifold, p, X, r::Real = 2)
+    norms =
+        (map(norm, M.manifolds, submanifold_components(M, p), submanifold_components(M, X)))
+    return norm(norms, r)
 end
 
 """
diff --git a/test/ManifoldsBaseTestUtils.jl b/test/ManifoldsBaseTestUtils.jl
index ce9fe2aa..9bad4f2c 100644
--- a/test/ManifoldsBaseTestUtils.jl
+++ b/test/ManifoldsBaseTestUtils.jl
@@ -575,6 +575,11 @@ end
 Base.getindex(x::MatrixVectorTransport, i) = x.m[:, i]
 Base.size(x::MatrixVectorTransport) = (size(x.m, 2),)
 
+struct TestArrayRepresentation <: AbstractPowerRepresentation end
+
+const TestPowerManifoldMultidimensional =
+    AbstractPowerManifold{𝔽,<:AbstractManifold{𝔽},TestArrayRepresentation} where {𝔽}
+
 export CustomDefinedInverseRetraction, CustomDefinedKeywordInverseRetraction
 export CustomDefinedKeywordRetraction, CustomDefinedRetraction, CustomUndefinedRetraction
 export DefaultPoint, DefaultTVector
diff --git a/test/default_manifold.jl b/test/default_manifold.jl
index 44bc97d7..b802fe97 100644
--- a/test/default_manifold.jl
+++ b/test/default_manifold.jl
@@ -31,6 +31,7 @@ using ManifoldsBaseTestUtils
     @test isa(manifold_dimension(M), Integer)
     @test manifold_dimension(M) ≥ 0
     @test base_manifold(M) == M
+    @test has_components(M)
     @test number_system(M) == ManifoldsBase.ℝ
     @test ManifoldsBase.representation_size(M) == (3,)
 
diff --git a/test/empty_manifold.jl b/test/empty_manifold.jl
index 6bc437b9..4f93744a 100644
--- a/test/empty_manifold.jl
+++ b/test/empty_manifold.jl
@@ -11,6 +11,7 @@ using Test
     @test base_manifold(M) === M
     @test number_system(M) === ℝ
     @test representation_size(M) === nothing
+    @test !has_components(M)
 
     @test_throws MethodError manifold_dimension(M)
 
diff --git a/test/power.jl b/test/power.jl
index ab69e05c..6d5a2d21 100644
--- a/test/power.jl
+++ b/test/power.jl
@@ -16,8 +16,6 @@ s = @__DIR__
 !(s in LOAD_PATH) && (push!(LOAD_PATH, s))
 using ManifoldsBaseTestUtils
 
-struct TestVectorSpaceType <: VectorSpaceType end
-
 power_array_wrapper(::Type{NestedPowerRepresentation}, ::Int) = identity
 power_array_wrapper(::Type{NestedReplacingPowerRepresentation}, i::Int) = SVector{i}
 
@@ -28,17 +26,14 @@ function ManifoldsBase.allocate(
     return similar(x)
 end
 
-struct TestArrayRepresentation <: AbstractPowerRepresentation end
-
-const TestPowerManifoldMultidimensional =
-    AbstractPowerManifold{𝔽,<:AbstractManifold{𝔽},TestArrayRepresentation} where {𝔽}
-
-function ManifoldsBase.representation_size(M::TestPowerManifoldMultidimensional)
+function ManifoldsBase.representation_size(
+    M::ManifoldsBaseTestUtils.TestPowerManifoldMultidimensional,
+)
     return (representation_size(M.manifold)..., ManifoldsBase.get_parameter(M.size)...)
 end
 
 @inline function ManifoldsBase._write(
-    ::TestPowerManifoldMultidimensional,
+    ::ManifoldsBaseTestUtils.TestPowerManifoldMultidimensional,
     rep_size::Tuple,
     x::AbstractArray,
     i::Tuple,
@@ -51,12 +46,14 @@ end
 
     @testset "Power Manifold with a test representation" begin
         M = ManifoldsBase.DefaultManifold(3)
-        N = PowerManifold(M, TestArrayRepresentation(), 2)
-        O = PowerManifold(N, TestArrayRepresentation(), 3) # joins instead of nesting.
+        N = PowerManifold(M, ManifoldsBaseTestUtils.TestArrayRepresentation(), 2)
+        O = PowerManifold(N, ManifoldsBaseTestUtils.TestArrayRepresentation(), 3) # joins instead of nesting.
         @test repr(O) ==
-              "PowerManifold(DefaultManifold(3; field = ℝ), TestArrayRepresentation(), 2, 3)"
+              "PowerManifold(DefaultManifold(3; field = ℝ), ManifoldsBaseTestUtils.TestArrayRepresentation(), 2, 3)"
         p = zeros(6)
         X = zeros(6)
+        @test has_components(N)
+        @test has_components(O)
         @test ManifoldsBase.check_power_size(N, p) === nothing
         @test ManifoldsBase.check_power_size(O, p) isa DomainError
         @test ManifoldsBase.check_power_size(N, p, X) === nothing
@@ -178,7 +175,16 @@ end
                     end
                 end
                 @testset "specific functions" begin
-                    @test distance(N, p, q) == sqrt(sum(distance.(Ref(M), p, q) .^ 2))
+                    nsq = sqrt(sum(distance.(Ref(M), p, q) .^ 2))
+                    @test distance(N, p, q) == nsq
+                    @test distance(N, p, q, LogarithmicInverseRetraction(), 2) == nsq
+                    @test distance(N, p, q, 2) == sqrt(sum(distance.(Ref(M), p, q) .^ 2))
+                    absn = sum(distance.(Ref(M), p, q))
+                    @test distance(N, p, q, 1) == absn
+                    @test distance(N, p, q, LogarithmicInverseRetraction(), 1) == absn
+                    Infn = maximum(distance.(Ref(M), p, q))
+                    @test distance(N, p, q, Inf) == Infn
+                    @test distance(N, p, q, LogarithmicInverseRetraction(), Inf) == Infn
                     @test exp(N, p, q) == p .+ q
                     @test exp(N, p, q, 2) == p .+ 2 .* q
 
@@ -226,6 +232,10 @@ end
                     @test manifold_dimension(N) == prod(pow_size) * manifold_dimension(M)
                     @test mid_point(N, p, q) == mid_point.(Ref(M), p, q)
                     @test sqrt(inner(N, p, q, q)) ≈ norm(N, p, q)
+                    norms = norm.(Ref(M), p, q)
+                    @test norm(N, p, q, 1) == sum(norms)
+                    @test norm(N, p, q, 2) == sqrt(sum(norms .^ 2))
+                    @test norm(N, p, q, Inf) == maximum(norms)
                     @test project(N, p) == p
                     @test project(N, p, q) == q
                     @test power_dimensions(N) == pow_size
@@ -464,7 +474,7 @@ end
         @test P2[NR, 1] === p
         @test P2[NR, 2] === p
 
-        NAR = PowerManifold(M, TestArrayRepresentation(), 2)
+        NAR = PowerManifold(M, ManifoldsBaseTestUtils.TestArrayRepresentation(), 2)
         P1 = fill(p, NAR)
         @test P1 isa Matrix{Float64}
         @test P1 == [1.0 1.0; 2.0 2.0; 3.0 3.0]
@@ -485,4 +495,28 @@ end
             @test norm(M, p, v) == 1
         end
     end
+
+    @testset "r-norm with an inverse retraction" begin
+        M = ManifoldsBaseTestUtils.TestSphere(2)
+        m = ProjectionInverseRetraction()
+        N = PowerManifold(M, NestedPowerRepresentation(), 2)
+        p1 = [1 / sqrt(2) .* [1.0, 1.0, 0.0], [1.0, 0.0, 0.0]]
+        p2 = [[0.0, 1.0, 0.0], 1 / sqrt(2) .* [1.0, 0.0, 1.0]]
+        X = log(N, p1, p2)
+        Xns = norm.(Ref(M), p1, X)
+        dRns = norm.(Ref(M), p1, inverse_retract.(Ref(M), p1, p2, Ref(m)))
+        @test distance(N, p1, p2) == norm(Xns)
+        @test distance(N, p1, p2, 1) == norm(Xns, 1)
+        @test distance(N, p1, p2, Inf) == norm(Xns, Inf)
+        @test distance(N, p1, p2, -Inf) == norm(Xns, -Inf)
+        @test distance(N, p1, p2, m) == norm(dRns)
+        @test distance(N, p1, p2, m, 1) == norm(dRns, 1)
+        @test distance(N, p1, p2, m, Inf) == norm(dRns, Inf)
+        @test distance(N, p1, p2, m, -Inf) == norm(dRns, -Inf)
+
+        @test norm(N, p1, X) == norm(Xns)
+        @test norm(N, p1, X, 1) == norm(Xns, 1)
+        @test norm(N, p1, X, Inf) == norm(Xns, Inf)
+        @test norm(N, p1, X, -Inf) == norm(Xns, -Inf)
+    end
 end
diff --git a/test/product_manifold.jl b/test/product_manifold.jl
index 6478b5a1..e566a48c 100644
--- a/test/product_manifold.jl
+++ b/test/product_manifold.jl
@@ -6,7 +6,7 @@ using ManifoldsBase:
 using LinearAlgebra
 using Random
 
-s = (@__DIR__) * "/test/"
+s = (@__DIR__)
 !(s in LOAD_PATH) && (push!(LOAD_PATH, s))
 using ManifoldsBaseTestUtils
 
@@ -34,6 +34,7 @@ using RecursiveArrayTools
     X1 = ArrayPartition([0.0, 1.0, 0.2], [4.0 0.0; 2.0 7.0])
 
     @test !is_flat(M)
+    @test has_components(M)
     @test M[1] == M1
     @test M[2] == M2
     @test injectivity_radius(M) ≈ π
@@ -160,20 +161,34 @@ using RecursiveArrayTools
     end
 
     p1 = ArrayPartition([1.0, 0.0, 0.0], [4.0 5.0; 6.0 7.0])
-    p2 = ArrayPartition([0.0, 1.0, 0.0], [4.0 8.0; 3.0 7.5])
+    p2 = ArrayPartition(1 / sqrt(2) .* [1.0, 1.0, 0.0], [4.0 8.0; 3.0 7.5])
     X1 = ArrayPartition([0.0, 1.0, 0.2], [4.0 0.0; 2.0 7.0])
     X2 = ArrayPartition([0.0, -1.0, 0.4], [3.0 1.0; -2.0 2.0])
-
+    d1 = distance(M1, p1.x[1], p2.x[1])
+    d1r = distance(M1, p1.x[1], p2.x[1], ProjectionInverseRetraction())
+    d2 = distance(M2, p1.x[2], p2.x[2])
     @testset "Basic operations" begin
         @test manifold_dimension(M) == 6
         @test representation_size(M) === nothing
-        @test distance(M, p1, p2) ≈ 4.551637188998299
+        @test distance(M, p1, p2) ≈ norm([d1, d2])
+        @test distance(M, p1, p2, 1) ≈ norm([d1, d2], 1)
+        @test distance(M, p1, p2, Inf) ≈ norm([d1, d2], Inf)
+        # Fallback log
+        m1 = LogarithmicInverseRetraction()
+        @test distance(M, p1, p2, m1) ≈ norm([d1, d2])
+        @test distance(M, p1, p2, m1, 1) ≈ norm([d1, d2], 1)
+        @test distance(M, p1, p2, m1, Inf) ≈ norm([d1, d2], Inf)
+        m2 = InverseProductRetraction(
+            ProjectionInverseRetraction(),
+            LogarithmicInverseRetraction(),
+        )
+        @test distance(M, p1, p2, m2) ≈ norm([d1r, d2])
+        @test distance(M, p1, p2, m2, 1) ≈ norm([d1r, d2], 1)
+        @test distance(M, p1, p2, m2, Inf) ≈ norm([d1r, d2], Inf)
         qr = similar(p1)
         exp!(M, qr, p1, X1)
-        @test exp(M, p1, X1) ≈ ArrayPartition(
-            [0.5235330372543839, 0.8354600062374664, 0.1670920012474933],
-            [8.0 5.0; 8.0 14.0],
-        )
+        @test exp(M, p1, X1) ≈
+              ArrayPartition(exp(M1, p1[M, 1], X1[M, 1]), exp(M2, p1[M, 2], X1[M, 2]))
         @test exp(M, p1, X1) ≈ qr
         @test exp(M, p1, X1, 2.0) ≈ exp(M, p1, 2 * X1)
         exp!(M, qr, p1, X1, 2.0)
@@ -200,8 +215,7 @@ using RecursiveArrayTools
         @test qr3 ≈ qr
         Xr = similar(X1)
         log!(M, Xr, p1, p2)
-        @test log(M, p1, p2) ≈
-              ArrayPartition([0.0, 1.5707963267948966, 0.0], [0.0 3.0; -3.0 0.5])
+        @test log(M, p1, p2) ≈ ArrayPartition([0.0, π / 4, 0.0], [0.0 3.0; -3.0 0.5])
         @test log(M, p1, p2) ≈ Xr
         @test inverse_retract(
             M,
@@ -270,11 +284,11 @@ using RecursiveArrayTools
     @testset "Broadcasting" begin
         br_result = p1 .+ 2.0 .* p2
         @test br_result isa ArrayPartition
-        @test br_result.x[1] ≈ [1.0, 2.0, 0.0]
+        @test br_result.x[1] ≈ [1 + sqrt(2), sqrt(2), 0.0]
         @test br_result.x[2] ≈ [12.0 21.0; 12.0 22.0]
 
         br_result .= 2.0 .* p1 .+ p2
-        @test br_result.x[1] ≈ [2.0, 1.0, 0.0]
+        @test br_result.x[1] ≈ [2.0 + 1 / sqrt(2), 1 / sqrt(2), 0.0]
         @test br_result.x[2] ≈ [12.0 18.0; 15.0 21.5]
 
         br_result .= p1
@@ -327,8 +341,9 @@ using RecursiveArrayTools
     end
 
     @testset "arithmetic" begin
-        @test isapprox(M, p1 + p2, ArrayPartition([1.0, 1.0, 0.0], [8.0 13.0; 9.0 14.5]))
-        @test isapprox(M, p1 - p2, ArrayPartition([1.0, -1.0, 0.0], [0.0 -3.0; 3.0 -0.5]))
+        a = 1 / sqrt(2)
+        @test isapprox(M, p1 + p2, ArrayPartition([1.0 + a, a, 0.0], [8.0 13.0; 9.0 14.5]))
+        @test isapprox(M, p1 - p2, ArrayPartition([1.0 - a, -a, 0.0], [0.0 -3.0; 3.0 -0.5]))
         @test isapprox(M, -p1, ArrayPartition([-1.0, -0.0, -0.0], [-4.0 -5.0; -6.0 -7.0]))
         @test isapprox(M, p1 * 2, ArrayPartition([2.0, 0.0, 0.0], [8.0 10.0; 12.0 14.0]))
         @test isapprox(M, 2 * p1, ArrayPartition([2.0, 0.0, 0.0], [8.0 10.0; 12.0 14.0]))
@@ -639,4 +654,18 @@ using RecursiveArrayTools
         @test base_point(Tp1M1) === p1[M, 1]
         @test base_manifold(Tp1M1) === M[1]
     end
+
+    @testset "Special Case for r-distance" begin
+        M1 = TestSphere(2)
+        M = ProductManifold(M1, M1)
+        p1 = ArrayPartition([1.0, 0.0, 0.0], 1 / sqrt(2) .* [1.0, 1.0, 0.0])
+        p2 = ArrayPartition(1 / sqrt(2) .* [1.0, 1.0, 0.0], [0.0, 1.0, 0.0])
+        m = ProjectionInverseRetraction()
+        ds = [distance(M1, p1[M, i], p2[M, i], m) for i in [1, 2]]
+        @test distance(M, p1, p2, m) == norm(ds)
+        @test distance(M, p1, p2, m, 1) == norm(ds, 1)
+        @test distance(M, p1, p2, m, 2) == norm(ds, 2)
+        @test distance(M, p1, p2, m, Inf) == norm(ds, Inf)
+        @test distance(M, p1, p2, m, -Inf) == norm(ds, -Inf)
+    end
 end