diff --git a/.gitignore b/.gitignore
index 06ba5eb893..a262bd806d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,4 @@
 docs/build
 Manifest.toml
+benchmark/tune.json
+benchmark/results.md
diff --git a/Project.toml b/Project.toml
index aa5b6a6179..91148c8368 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,6 +14,7 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
 
 [extras]
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
@@ -22,4 +23,4 @@ ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["Test", "DoubleFloats", "ForwardDiff", "ReverseDiff"]
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
new file mode 100644
index 0000000000..52bb05095b
--- /dev/null
+++ b/benchmark/benchmarks.jl
@@ -0,0 +1,131 @@
+using Manifolds
+using BenchmarkTools
+using StaticArrays
+using LinearAlgebra
+using Random
+
+# Define a parent BenchmarkGroup to contain our suite
+SUITE = BenchmarkGroup()
+
+Random.seed!(12334)
+
+function add_manifold(M::Manifold, pts, name;
+    test_tangent_vector_broadcasting = true,
+    retraction_methods = [],
+    inverse_retraction_methods = [],
+    point_distributions = [],
+    tvector_distributions = [])
+
+    SUITE["manifolds"][name] = BenchmarkGroup()
+    tv = log(M, pts[1], pts[2])
+    tv1 = log(M, pts[1], pts[2])
+    tv2 = log(M, pts[1], pts[3])
+    p = similar(pts[1])
+    SUITE["manifolds"][name]["similar"] = @benchmarkable similar($(pts[1]))
+    SUITE["manifolds"][name]["log"] = @benchmarkable log($M, $(pts[1]), $(pts[2]))
+    SUITE["manifolds"][name]["log!"] = @benchmarkable log!($M, $tv, $(pts[1]), $(pts[2]))
+    for iretr ∈ inverse_retraction_methods
+        SUITE["manifolds"][name]["inverse_retract: "*string(iretr)] = @benchmarkable inverse_retract($M, $(pts[1]), $(pts[2]), $iretr)
+        SUITE["manifolds"][name]["inverse_retract!: "*string(iretr)] = @benchmarkable inverse_retract!($M, $tv, $(pts[1]), $(pts[2]), $iretr)
+    end
+    SUITE["manifolds"][name]["exp"] = @benchmarkable exp($M, $(pts[1]), $tv1)
+    SUITE["manifolds"][name]["exp!"] = @benchmarkable exp!($M, $p, $(pts[1]), $tv1)
+    for retr ∈ retraction_methods
+        SUITE["manifolds"][name]["retract: "*string(retr)] = @benchmarkable retract($M, $(pts[1]), $tv1, $retr)
+        SUITE["manifolds"][name]["retract!: "*string(retr)] = @benchmarkable retract!($M, $p, $(pts[1]), $tv1, $retr)
+    end
+    SUITE["manifolds"][name]["norm"] = @benchmarkable norm($M, $(pts[1]), $tv1)
+    SUITE["manifolds"][name]["inner"] = @benchmarkable inner($M, $(pts[1]), $tv1, $tv2)
+    SUITE["manifolds"][name]["distance"] = @benchmarkable distance($M, $(pts[1]), $(pts[2]))
+    SUITE["manifolds"][name]["isapprox (pt)"] = @benchmarkable isapprox($M, $(pts[1]), $(pts[2]))
+    SUITE["manifolds"][name]["isapprox (tv)"] = @benchmarkable isapprox($M, $(pts[1]), $tv1, $tv2)
+    SUITE["manifolds"][name]["2 * tv1 + 3 * tv2"] = @benchmarkable 2 * $tv1 + 3 * $tv2
+    SUITE["manifolds"][name]["tv = 2 * tv1 + 3 * tv2"] = @benchmarkable $tv = 2 * $tv1 + 3 * $tv2
+    if test_tangent_vector_broadcasting
+        SUITE["manifolds"][name]["tv = 2 .* tv1 .+ 3 .* tv2"] = @benchmarkable $tv = 2 .* $tv1 .+ 3 .* $tv2
+        SUITE["manifolds"][name]["tv .= 2 .* tv1 .+ 3 .* tv2"] = @benchmarkable $tv .= 2 .* $tv1 .+ 3 .* $tv2
+    end
+    for pd ∈ point_distributions
+        distr_name = string(pd)
+        distr_name = distr_name[1:min(length(distr_name), 50)]
+        SUITE["manifolds"][name]["point distribution "*distr_name] = @benchmarkable rand($pd)
+    end
+    for tvd ∈ point_distributions
+        distr_name = string(tvd)
+        distr_name = distr_name[1:min(length(distr_name), 50)]
+        SUITE["manifolds"][name]["tangent vector distribution "*distr_name] = @benchmarkable rand($tvd)
+    end
+end
+
+# General manifold benchmarks
+function add_manifold_benchmarks()
+
+    SUITE["manifolds"] = BenchmarkGroup()
+
+    s2 = Manifolds.Sphere(2)
+    array_s2 = ArrayManifold(s2)
+    r2 = Manifolds.Euclidean(2)
+
+    pts_r2 = [Size(2)([1.0, 1.0]),
+              Size(2)([-2.0, 3.0]),
+              Size(2)([3.0, -2.0])]
+
+    add_manifold(r2,
+                 pts_r2,
+                 "Euclidean{2} -- SizedArray")
+
+    add_manifold(r2,
+                 [MVector{2,Float64}([1.0, 1.0]),
+                  MVector{2,Float64}([-2.0, 3.0]),
+                  MVector{2,Float64}([3.0, -2.0])],
+                  "Euclidean{2} -- MVector")
+
+    ud_sphere = Manifolds.uniform_distribution(s2, Size(3)([1.0, 0.0, 0.0]))
+    gtd_sphere = Manifolds.normal_tvector_distribution(s2, Size(3)([1.0, 0.0, 0.0]), 1.0)
+
+    pts_s2 = [Size(3)([1.0, 0.0, 0.0]),
+              Size(3)([0.0, 1.0, 0.0]),
+              Size(3)([0.0, 0.0, 1.0])]
+
+    add_manifold(s2,
+                 pts_s2,
+                 "Sphere{2} -- SizedArray";
+                 point_distributions = [ud_sphere],
+                 tvector_distributions = [gtd_sphere])
+
+    add_manifold(array_s2,
+                 [Size(3)([1.0, 0.0, 0.0]),
+                  Size(3)([0.0, 1.0, 0.0]),
+                  Size(3)([0.0, 0.0, 1.0])],
+                 "ArrayManifold{Sphere{2}} -- SizedArray";
+                 test_tangent_vector_broadcasting = false)
+
+    retraction_methods_rot = [Manifolds.PolarRetraction(),
+                              Manifolds.QRRetraction()]
+
+    inverse_retraction_methods_rot = [Manifolds.PolarInverseRetraction(),
+                                      Manifolds.QRInverseRetraction()]
+
+    so2 = Manifolds.Rotations(2)
+    angles = (0.0, π/2, 2π/3)
+    add_manifold(so2,
+                 [Size(2, 2)([cos(ϕ) sin(ϕ); -sin(ϕ) cos(ϕ)]) for ϕ in angles],
+                 "Rotations(2) -- SizedArray",
+                 retraction_methods = retraction_methods_rot,
+                 inverse_retraction_methods = inverse_retraction_methods_rot)
+
+    m_prod = Manifolds.ProductManifold(s2, r2)
+    shape_s2r2 = Manifolds.ShapeSpecification(s2, r2)
+
+    pts_prd_base = [[1.0, 0.0, 0.0, 0.0, 0.0],
+                    [0.0, 1.0, 0.0, 1.0, 0.0],
+                    [0.0, 0.0, 1.0, 0.0, 0.1]]
+    pts_prod = map(p -> Manifolds.ProductArray(shape_s2r2, p), pts_prd_base)
+    add_manifold(m_prod, pts_prod, "ProductManifold with ProductArray")
+
+    pts_prod_mpoints = [Manifolds.ProductMPoint(p[1], p[2]) for p in zip(pts_s2, pts_r2)]
+    add_manifold(m_prod, pts_prod_mpoints, "ProductManifold with MPoint";
+        test_tangent_vector_broadcasting = false)
+end
+
+add_manifold_benchmarks()
diff --git a/benchmark/runbenchmarks.jl b/benchmark/runbenchmarks.jl
new file mode 100644
index 0000000000..9b76f7b335
--- /dev/null
+++ b/benchmark/runbenchmarks.jl
@@ -0,0 +1,9 @@
+using PkgBenchmark
+
+#--track-allocation=all
+config = BenchmarkConfig(id = nothing,
+    juliacmd = `julia -O3`,
+    env = Dict("JULIA_NUM_THREADS" => 4))
+
+results = benchmarkpkg("Manifolds", config)
+export_markdown("benchmark/results.md", results)
diff --git a/docs/Project.toml b/docs/Project.toml
new file mode 100644
index 0000000000..e741138502
--- /dev/null
+++ b/docs/Project.toml
@@ -0,0 +1,5 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+
+[compat]
+Documenter = "0.22"
diff --git a/docs/make.jl b/docs/make.jl
index 3d9ed663f7..e0f61e4f97 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -8,9 +8,22 @@ makedocs(
     pages = [
         "Home" => "index.md",
         "Manifolds" => [
-            "Sphere" => "manifolds/sphere.md",
-            "Rotations" => "manifolds/rotations.md",
+            "Basic manifolds" => [
+                "Euclidean" => "manifolds/euclidean.md",
+                "Rotations" => "manifolds/rotations.md",
+                "Sphere" => "manifolds/sphere.md"
+            ],
+            "Combined manifolds" => [
+                "Product manifold" => "manifolds/product.md"
+            ],
+            "Manifold decorators" => [
+                "Array manifold" => "manifolds/array.md",
+                "Metric manifold" => "manifolds/metric.md"
+            ]
         ],
-        "Distributions" => "distributions.md"
+        "Distributions" => "distributions.md",
+        "Library" => [
+            "Internals" => "lib/internals.md"
+        ]
     ]
 )
diff --git a/docs/src/distributions.md b/docs/src/distributions.md
index b9090374cf..2ba7994892 100644
--- a/docs/src/distributions.md
+++ b/docs/src/distributions.md
@@ -7,3 +7,7 @@ Modules = [Manifolds]
 Pages = ["DistributionsBase.jl"]
 Order = [:type, :function]
 ```
+```@docs
+Manifolds.ProjectedPointDistribution
+Manifolds.ProjectedTVectorDistribution
+```
diff --git a/docs/src/index.md b/docs/src/index.md
index 1a18678b6e..6e10a648fa 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -12,17 +12,8 @@ Order = [:type, :function]
 ```
 
 ## Decorators
-### ArrayManifold
-```@autodocs
-Modules = [Manifolds]
-Pages = ["ArrayManifold.jl"]
-Order = [:type, :function]
-```
-### Lie Group
 
-### Metric
-```@autodocs
-Modules = [Manifolds]
-Pages = ["Metric.jl"]
-Order = [:type, :function]
-```
+* [Array manifold](@ref)
+* [Metric manifold](@ref)
+
+### Lie Group
diff --git a/docs/src/lib/internals.md b/docs/src/lib/internals.md
new file mode 100644
index 0000000000..01a8a3cf6f
--- /dev/null
+++ b/docs/src/lib/internals.md
@@ -0,0 +1,9 @@
+# Internal Documentation
+
+Documentation for `Manifolds.jl`'s internal methods.
+
+```@docs
+Manifolds.SizedAbstractArray
+Manifolds.ziptuples
+Manifolds.size_to_tuple
+```
diff --git a/docs/src/manifolds/array.md b/docs/src/manifolds/array.md
new file mode 100644
index 0000000000..1ee496f343
--- /dev/null
+++ b/docs/src/manifolds/array.md
@@ -0,0 +1,7 @@
+# Array manifold
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["ArrayManifold.jl"]
+Order = [:type, :function]
+```
diff --git a/docs/src/manifolds/euclidean.md b/docs/src/manifolds/euclidean.md
new file mode 100644
index 0000000000..8f445f2c08
--- /dev/null
+++ b/docs/src/manifolds/euclidean.md
@@ -0,0 +1,7 @@
+# Sphere
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["Euclidean.jl"]
+Order = [:type, :function]
+```
diff --git a/docs/src/manifolds/metric.md b/docs/src/manifolds/metric.md
new file mode 100644
index 0000000000..c24c12709f
--- /dev/null
+++ b/docs/src/manifolds/metric.md
@@ -0,0 +1,7 @@
+# Metric manifold
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["Metric.jl"]
+Order = [:type, :function]
+```
diff --git a/docs/src/manifolds/product.md b/docs/src/manifolds/product.md
new file mode 100644
index 0000000000..b502571b56
--- /dev/null
+++ b/docs/src/manifolds/product.md
@@ -0,0 +1,9 @@
+# Product Manifold
+
+Product manifold $M = M_1 \times M_2 \times \dots M_n$ of manifolds $M_1, M_2, \dots, M_n$. Points on the product manifold can be constructed using [`ProductMPoint`](@ref) with canonical projections $\Pi_i \colon M \to M_i$ for $i \in 1, 2, \dots, n$ provided by [`submanifold_component`](@ref).
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["ProductManifold.jl"]
+Order = [:type, :function]
+```
diff --git a/src/Euclidean.jl b/src/Euclidean.jl
index c37f950868..07ff585cb6 100644
--- a/src/Euclidean.jl
+++ b/src/Euclidean.jl
@@ -19,6 +19,14 @@ struct Euclidean{T<:Tuple} <: Manifold where {T} end
 Euclidean(n::Int) = Euclidean{Tuple{n}}()
 Euclidean(m::Int, n::Int) = Euclidean{Tuple{m,n}}()
 
+function representation_size(::Euclidean{Tuple{n}}, ::Type{T}) where {n,T<:Union{MPoint, TVector, CoTVector}}
+    return (n,)
+end
+
+function representation_size(::Euclidean{Tuple{m,n}}, ::Type{T}) where {m,n,T<:Union{MPoint, TVector, CoTVector}}
+    return (m,n)
+end
+
 @generated manifold_dimension(::Euclidean{T}) where {T} = *(T.parameters...)
 
 struct EuclideanMetric <: RiemannianMetric end
@@ -37,12 +45,50 @@ det_local_metric(M::MetricManifold{<:Manifold,EuclideanMetric}, x) = one(eltype(
 
 log_local_metric_density(M::MetricManifold{<:Manifold,EuclideanMetric}, x) = zero(eltype(x))
 
-inner(::Euclidean, x, v, w) = dot(v, w)
-inner(::MetricManifold{<:Manifold,EuclideanMetric}, x, v, w) = dot(v, w)
+@inline inner(::Euclidean, x, v, w) = dot(v, w)
+@inline inner(::MetricManifold{<:Manifold,EuclideanMetric}, x, v, w) = dot(v, w)
 
+distance(::Euclidean, x, y) = norm(x-y)
 norm(::Euclidean, x, v) = norm(v)
 norm(::MetricManifold{<:Manifold,EuclideanMetric}, x, v) = norm(v)
 
 exp!(M::Euclidean, y, x, v) = (y .= x + v)
 
 log!(M::Euclidean, v, x, y) = (v .= y - x)
+
+function zero_tangent_vector!(M::Euclidean, v, x)
+    fill!(v, 0)
+    return v
+end
+
+project_point!(M::Euclidean, x) = x
+
+function project_tangent!(M::Euclidean, w, x, v)
+    w .= v
+    return w
+end
+
+"""
+    projected_distribution(M::Euclidean, d, [x])
+
+Wraps standard distribution `d` into a manifold-valued distribution. Generated
+points will be of similar type to `x`. By default, the type is not changed.
+"""
+function projected_distribution(M::Euclidean, d, x)
+    return ProjectedPointDistribution(M, d, project_point!, x)
+end
+
+function projected_distribution(M::Euclidean, d)
+    return ProjectedPointDistribution(M, d, project_point!, rand(d))
+end
+
+"""
+    normal_tvector_distribution(S::Euclidean, x, σ)
+
+Normal distribution in ambient space with standard deviation `σ`
+projected to tangent space at `x`.
+"""
+function normal_tvector_distribution(M::Euclidean{Tuple{N}}, x, σ) where N
+    d = Distributions.MvNormal(zero(x), σ)
+    return ProjectedTVectorDistribution(M, x, d, project_tangent!, x)
+end
diff --git a/src/Manifolds.jl b/src/Manifolds.jl
index e1c7d7ab25..e373aaa0e7 100644
--- a/src/Manifolds.jl
+++ b/src/Manifolds.jl
@@ -6,13 +6,19 @@ import Base: isapprox,
     angle,
     eltype,
     similar,
+    getindex,
+    setindex!,
+    size,
+    copy,
     convert,
+    dataids,
     +,
     -,
     *
 import LinearAlgebra: dot,
     norm,
     det,
+    cross,
     I,
     UniformScaling,
     Diagonal
@@ -24,6 +30,7 @@ using LinearAlgebra
 using Random: AbstractRNG
 using SimpleTraits
 using ForwardDiff
+using UnsafeArrays
 import Einsum: @einsum
 import OrdinaryDiffEq: ODEProblem,
     AutoVern9,
@@ -115,6 +122,25 @@ of each point of the manifold is homeomorphic.
 """
 function manifold_dimension end
 
+@doc doc"""
+    representation_size(M::Manifold, ::Type{T}) where {T}
+
+The size of array representing an object of type `T` on manifold `M`,
+for example point, tangent vector or cotangent vector.
+The second argument should in these cases be equal to, respectively,
+`MPoint`, `TVector` and `CoTVector`, regardless of the type used to represent
+said objects.
+"""
+function representation_size end
+
+@traitfn function representation_size(M::MT, ::Type{T}) where {MT<:Manifold,T;!IsDecoratorManifold{MT}}
+    error("representation_size not implemented for manifold $(typeof(M)) and type $(T).")
+end
+
+@traitfn function representation_size(M::MT, ::Type{T}) where {MT<:Manifold,T;IsDecoratorManifold{MT}}
+    return representation_size(base_manifold(M), T)
+end
+
 @traitfn function manifold_dimension(M::MT) where {MT<:Manifold;!IsDecoratorManifold{MT}}
     error("manifold_dimension not implemented for a $(typeof(M)).")
 end
@@ -243,6 +269,14 @@ function inverse_retract(M::Manifold, x, y)
     return vr
 end
 
+project_point!(M::Manifold, x) = error("project onto tangent space not implemented for a $(typeof(M)) and point $(typeof(x)).")
+
+function project_point(M::Manifold, x)
+    y = similar_result(M, project_point, x)
+    project_tangent!(M, y, x)
+    return y
+end
+
 project_tangent!(M::Manifold, w, x, v) = error("project onto tangent space not implemented for a $(typeof(M)) and point $(typeof(x)) with input $(typeof(v)).")
 
 function project_tangent(M::Manifold, x, v)
@@ -458,19 +492,27 @@ the assumption is to be optimistic.
 is_tangent_vector(M::Manifold, x, v; kwargs...) = true
 is_tangent_vector(M::Manifold, x::MPoint, v::TVector) = error("A validation for a $(typeof(v)) in the tangent space of a $(typeof(x)) on $(typeof(M)) not implemented.")
 
+include("utils.jl")
+
 include("ArrayManifold.jl")
 
 include("DistributionsBase.jl")
 include("Metric.jl")
 include("Euclidean.jl")
+include("ProductManifold.jl")
 include("Rotations.jl")
 include("Sphere.jl")
 include("ProjectedDistribution.jl")
 
 export Manifold,
     IsDecoratorManifold,
-    Euclidean
-export manifold_dimension,
+    Euclidean,
+    Sphere,
+    ProductManifold,
+    ProductMPoint,
+    ProductTVector,
+    ProductCoTVector
+export ×,
     base_manifold,
     distance,
     exp,
@@ -488,8 +530,14 @@ export manifold_dimension,
     log!,
     manifold_dimension,
     norm,
+    project_point,
+    project_point!,
+    project_tangent,
+    project_tangent!,
     retract,
     retract!,
+    submanifold,
+    submanifold_component,
     zero_tangent_vector,
     zero_tangent_vector!
 export Metric,
diff --git a/src/ProductManifold.jl b/src/ProductManifold.jl
new file mode 100644
index 0000000000..cc2afc4025
--- /dev/null
+++ b/src/ProductManifold.jl
@@ -0,0 +1,505 @@
+@doc doc"""
+    ProductManifold{TM<:Tuple, TRanges<:Tuple, TSizes<:Tuple} <: Manifold
+
+Product manifold $M_1 \times M_2 \times \dots \times M_n$ with product geometry.
+`TRanges` and `TSizes` statically define the relationship between representation
+of the product manifold and representations of point, tangent vectors
+and cotangent vectors of respective manifolds.
+
+# Constructor
+
+    ProductManifold(M_1, M_2, ..., M_n)
+
+generates the product manifold $M_1 \times M_2 \times \dots \times M_n$.
+Alternatively, the same manifold can be contructed using the `×` operator:
+`M_1 × M_2 × M_3`.
+"""
+struct ProductManifold{TM<:Tuple} <: Manifold
+    manifolds::TM
+end
+
+ProductManifold(manifolds::Manifold...) = ProductManifold{typeof(manifolds)}(manifolds)
+
+function cross(M1::Manifold, M2::Manifold)
+    return ProductManifold(M1, M2)
+end
+
+function cross(M1::ProductManifold, M2::Manifold)
+    return ProductManifold(M1.manifolds..., M2)
+end
+
+function cross(M1::Manifold, M2::ProductManifold)
+    return ProductManifold(M1, M2.manifolds...)
+end
+
+function cross(M1::ProductManifold, M2::ProductManifold)
+    return ProductManifold(M1.manifolds..., M2.manifolds...)
+end
+
+"""
+    submanifold(M::ProductManifold, i::Integer)
+
+Extract the `i`th factor of the product manifold `M`.
+"""
+function submanifold(M::ProductManifold, i::Integer)
+    return M.manifolds[i]
+end
+
+"""
+    submanifold(M::ProductManifold, i::Val)
+
+Extract the factor of the product manifold `M` indicated by indices in `i`.
+For example, for `i` equal to `Val((1, 3))` the product manifold constructed
+from the first and the third factor is returned.
+"""
+function submanifold(M::ProductManifold, i::Val)
+    return ProductManifold(select_from_tuple(M.manifolds, i))
+end
+
+"""
+    submanifold(M::ProductManifold, i::AbstractVector)
+
+Extract the factor of the product manifold `M` indicated by indices in `i`.
+For example, for `i` equal to `[1, 3]` the product manifold constructed
+from the first and the third factor is returned.
+
+This function is not type-stable, for better preformance use
+`submanifold(M::ProductManifold, i::Val)`.
+"""
+submanifold(M::ProductManifold, i::AbstractVector) = submanifold(M, Val(tuple(i...)))
+
+"""
+    ShapeSpecification(manifolds::Manifold...)
+
+A structure for specifying array size and offset information for linear
+storage of points and tangent vectors on the product manifold of `manifolds`.
+
+For example, consider the shape specification for the product of
+a sphere and group of rotations:
+
+```julia-repl
+julia> M1 = Sphere(2)
+Sphere{2}()
+
+julia> M3 = Rotations(2)
+Rotations{2}()
+
+julia> shape = Manifolds.ShapeSpecification(M1, M3)
+Manifolds.ShapeSpecification{(1:3, 4:7),Tuple{Tuple{3},Tuple{2,2}}}()
+```
+
+`TRanges` contains ranges in the linear storage that correspond to a specific
+manifold. `Sphere(2)` needs three numbers and is first, so it is allocated the
+first three elements of the linear storage (`1:3`). `Rotations(2)` needs four
+numbers and is second, so the next four numbers are allocated to it (`4:7`).
+`TSizes` describe how the linear storage must be reshaped to correctly
+represent points. In this case, `Sphere(2)` expects a three-element vector, so
+the corresponding size is `Tuple{3}`. On the other hand, `Rotations(2)`
+expects two-by-two matrices, so its size specification is `Tuple{2,2}`.
+"""
+struct ShapeSpecification{TRanges, TSizes} end
+
+function ShapeSpecification(manifolds::Manifold...)
+    sizes = map(m -> representation_size(m, MPoint), manifolds)
+    lengths = map(prod, sizes)
+    ranges = UnitRange{Int64}[]
+    k = 1
+    for len ∈ lengths
+        push!(ranges, k:(k+len-1))
+        k += len
+    end
+    TRanges = tuple(ranges...)
+    TSizes = Tuple{map(s -> Tuple{s...}, sizes)...}
+    return ShapeSpecification{TRanges, TSizes}()
+end
+
+"""
+    ProductArray(shape::ShapeSpecification, data)
+
+An array-based representation for points and tangent vectors on the
+product manifold. `data` contains underlying representation of points
+arranged according to `TRanges` and `TSizes` from `shape`.
+Internal views for each specific sub-point are created and stored in `parts`.
+"""
+struct ProductArray{TM<:ShapeSpecification,T,N,TData<:AbstractArray{T,N},TV<:Tuple} <: AbstractArray{T,N}
+    data::TData
+    parts::TV
+end
+
+# The two-argument version of this constructor is substantially faster than
+# the generic one.
+function ProductArray(M::Type{ShapeSpecification{TRanges, Tuple{Size1, Size2}}}, data::TData) where {TRanges, Size1, Size2, T, N, TData<:AbstractArray{T,N}}
+    views = (SizedAbstractArray{Size1}(view(data, TRanges[1])),
+             SizedAbstractArray{Size2}(view(data, TRanges[2])))
+    return ProductArray{M, T, N, TData, typeof(views)}(data, views)
+end
+
+function ProductArray(M::Type{ShapeSpecification{TRanges, Tuple{Size1, Size2, Size3}}}, data::TData) where {TRanges, Size1, Size2, Size3, T, N, TData<:AbstractArray{T,N}}
+    views = (SizedAbstractArray{Size1}(view(data, TRanges[1])),
+             SizedAbstractArray{Size2}(view(data, TRanges[2])),
+             SizedAbstractArray{Size3}(view(data, TRanges[3])))
+    return ProductArray{M, T, N, TData, typeof(views)}(data, views)
+end
+
+function ProductArray(M::Type{ShapeSpecification{TRanges, TSizes}}, data::TData) where {TRanges, TSizes, T, N, TData<:AbstractArray{T,N}}
+    views = map((size, range) -> SizedAbstractArray{size}(view(data, range)), size_to_tuple(TSizes), TRanges)
+    return ProductArray{M, T, N, TData, typeof(views)}(data, views)
+end
+
+function ProductArray(M::ShapeSpecification{TRanges, TSizes}, data::TData) where {TM, TRanges, TSizes, TData<:AbstractArray}
+    return ProductArray(typeof(M), data)
+end
+
+@doc doc"""
+    prod_point(M::ShapeSpecification, pts...)
+
+Construct a product point from product manifold `M` based on point `pts`
+represented by [`ProductArray`](@ref).
+
+# Example
+To construct a point on the product manifold $S^2 \times \mathbb{R}^2$
+from points on the sphere and in the euclidean space represented by,
+respectively, `[1.0, 0.0, 0.0]` and `[-3.0, 2.0]` you need to construct shape
+specification first. It describes how linear storage of `ProductArray`
+corresponds to array representations expected by `Sphere(2)` and `Euclidean(2)`.
+
+    M1 = Sphere(2)
+    M2 = Euclidean(2)
+    Mshape = Manifolds.ShapeSpecification(M1, M2)
+
+Next, the desired point on the product manifold can be obtained by calling
+`Manifolds.prod_point(Mshape, [1.0, 0.0, 0.0], [-3.0, 2.0])`.
+"""
+function prod_point(M::ShapeSpecification, pts...)
+    data = mapreduce(vcat, pts) do pt
+        reshape(pt, :)
+    end
+    # Array(data) is used to ensure that the data is mutable
+    # `mapreduce` can return `SArray` for some arguments
+    return ProductArray(M, Array(data))
+end
+
+"""
+    submanifold_component(x::ProductArray, i::Integer)
+
+Project the product array `x` to its `i`th component. A new array is returned.
+"""
+function submanifold_component(x::ProductArray, i::Integer)
+    return x.parts[i]
+end
+
+Base.BroadcastStyle(::Type{<:ProductArray{ShapeSpec}}) where ShapeSpec<:ShapeSpecification = Broadcast.ArrayStyle{ProductArray{ShapeSpec}}()
+
+function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{ProductArray{ShapeSpec}}}, ::Type{ElType}) where {ShapeSpec, ElType}
+    A = find_pv(bc)
+    return ProductArray(ShapeSpec, similar(A.data, ElType))
+end
+
+Base.dataids(x::ProductArray) = Base.dataids(x.data)
+
+"""
+    find_pv(x...)
+
+`A = find_pv(x...)` returns the first `ProductArray` among the arguments.
+"""
+@inline find_pv(bc::Base.Broadcast.Broadcasted) = find_pv(bc.args)
+@inline find_pv(args::Tuple) = find_pv(find_pv(args[1]), Base.tail(args))
+@inline find_pv(x) = x
+@inline find_pv(a::ProductArray, rest) = a
+@inline find_pv(::Any, rest) = find_pv(rest)
+
+size(x::ProductArray) = size(x.data)
+Base.@propagate_inbounds getindex(x::ProductArray, i) = getindex(x.data, i)
+Base.@propagate_inbounds setindex!(x::ProductArray, val, i) = setindex!(x.data, val, i)
+
+(+)(v1::ProductArray{ShapeSpec}, v2::ProductArray{ShapeSpec}) where ShapeSpec<:ShapeSpecification = ProductArray(ShapeSpec, v1.data + v2.data)
+(-)(v1::ProductArray{ShapeSpec}, v2::ProductArray{ShapeSpec}) where ShapeSpec<:ShapeSpecification = ProductArray(ShapeSpec, v1.data - v2.data)
+(-)(v::ProductArray{ShapeSpec}) where ShapeSpec<:ShapeSpecification = ProductArray(ShapeSpec, -v.data)
+(*)(a::Number, v::ProductArray{ShapeSpec}) where ShapeSpec<:ShapeSpecification = ProductArray(ShapeSpec, a*v.data)
+
+eltype(::Type{ProductArray{TM, TData, TV}}) where {TM, TData, TV} = eltype(TData)
+similar(x::ProductArray{ShapeSpec}) where ShapeSpec<:ShapeSpecification = ProductArray(ShapeSpec, similar(x.data))
+similar(x::ProductArray{ShapeSpec}, ::Type{T}) where {ShapeSpec<:ShapeSpecification, T} = ProductArray(ShapeSpec, similar(x.data, T))
+
+"""
+    ProductMPoint(parts)
+
+A more general but slower representation of points on a product manifold.
+
+# Example:
+
+A product point on a product manifold `Sphere(2) × Euclidean(2)` might be
+created as
+
+    ProductMPoint([1.0, 0.0, 0.0], [2.0, 3.0])
+
+where `[1.0, 0.0, 0.0]` is the part corresponding to the sphere factor
+and `[2.0, 3.0]` is the part corresponding to the euclidean manifold.
+"""
+struct ProductMPoint{TM<:Tuple} <: MPoint
+    parts::TM
+end
+
+ProductMPoint(points...) = ProductMPoint{typeof(points)}(points)
+eltype(x::ProductMPoint) = eltype(Tuple{map(eltype, x.parts)...})
+similar(x::ProductMPoint) = ProductMPoint(map(similar, x.parts)...)
+similar(x::ProductMPoint, ::Type{T}) where T = ProductMPoint(map(t -> similar(t, T), x.parts)...)
+
+function submanifold_component(x::ProductMPoint, i::Integer)
+    return x.parts[i]
+end
+
+"""
+    ProductTVector(parts)
+
+A more general but slower representation of tangent vectors on a product
+manifold.
+"""
+struct ProductTVector{TM<:Tuple} <: TVector
+    parts::TM
+end
+
+ProductTVector(vectors...) = ProductTVector{typeof(vectors)}(vectors)
+eltype(x::ProductTVector) = eltype(Tuple{map(eltype, x.parts)...})
+similar(x::ProductTVector) = ProductTVector(map(similar, x.parts)...)
+similar(x::ProductTVector, ::Type{T}) where T = ProductTVector(map(t -> similar(t, T), x.parts)...)
+
+(+)(v1::ProductTVector, v2::ProductTVector) = ProductTVector(map(+, v1.parts, v2.parts)...)
+(-)(v1::ProductTVector, v2::ProductTVector) = ProductTVector(map(-, v1.parts, v2.parts)...)
+(-)(v::ProductTVector) = ProductTVector(map(-, v.parts))
+(*)(a::Number, v::ProductTVector) = ProductTVector(map(t -> a*t, v.parts))
+
+function submanifold_component(x::ProductTVector, i::Integer)
+    return x.parts[i]
+end
+
+"""
+    ProductCoTVector(parts)
+
+A more general but slower representation of cotangent vectors on a product
+manifold.
+"""
+struct ProductCoTVector{TM<:Tuple} <: CoTVector
+    parts::TM
+end
+
+ProductCoTVector(covectors...) = ProductCoTVector{typeof(covectors)}(covectors)
+eltype(x::ProductCoTVector) = eltype(Tuple{map(eltype, x.parts)...})
+similar(x::ProductCoTVector) = ProductCoTVector(map(similar, x.parts)...)
+similar(x::ProductCoTVector, ::Type{T}) where T = ProductCoTVector(map(t -> similar(t, T), x.parts)...)
+
+function submanifold_component(x::ProductCoTVector, i::Integer)
+    return x.parts[i]
+end
+
+function isapprox(M::ProductManifold, x, y; kwargs...)
+    return all(t -> isapprox(t...; kwargs...), ziptuples(M.manifolds, x.parts, y.parts))
+end
+
+function isapprox(M::ProductManifold, x, v, w; kwargs...)
+    return all(t -> isapprox(t...; kwargs...), ziptuples(M.manifolds, x.parts, v.parts, w.parts))
+end
+
+function representation_size(M::ProductManifold, ::Type{T}) where {T}
+    return (mapreduce(m -> sum(representation_size(m, T)), +, M.manifolds),)
+end
+
+manifold_dimension(M::ProductManifold) = mapreduce(manifold_dimension, +, M.manifolds)
+
+struct ProductMetric <: Metric end
+
+function det_local_metric(M::MetricManifold{ProductManifold, ProductMetric}, x::ProductArray)
+    dets = map(det_local_metric, M.manifolds, x.parts)
+    return prod(dets)
+end
+
+function inner(M::ProductManifold, x, v, w)
+    subproducts = map(inner, M.manifolds, x.parts, v.parts, w.parts)
+    return sum(subproducts)
+end
+
+function norm(M::ProductManifold, x, v)
+    norms_squared = map(norm, M.manifolds, x.parts, v.parts).^2
+    return sqrt(sum(norms_squared))
+end
+
+function distance(M::ProductManifold, x, y)
+    return sqrt(sum(map(distance, M.manifolds, x.parts, y.parts).^2))
+end
+
+function exp!(M::ProductManifold, y, x, v)
+    map(exp!, M.manifolds, y.parts, x.parts, v.parts)
+    return y
+end
+
+function exp(M::ProductManifold, x::ProductMPoint, v::ProductTVector)
+    return ProductMPoint(map(exp, M.manifolds, x.parts, v.parts)...)
+end
+
+function log!(M::ProductManifold, v, x, y)
+    map(log!, M.manifolds, v.parts, x.parts, y.parts)
+    return v
+end
+
+function log(M::ProductManifold, x::ProductMPoint, y::ProductMPoint)
+    return ProductTVector(map(log, M.manifolds, x.parts, y.parts)...)
+end
+
+function injectivity_radius(M::ProductManifold, x)
+    return min(map(injectivity_radius, M.manifolds, x.parts)...)
+end
+
+function injectivity_radius(M::ProductManifold)
+    return min(map(injectivity_radius, M.manifolds)...)
+end
+
+"""
+    ProductRetraction(retractions::AbstractRetractionMethod...)
+
+Product retraction of `retractions`. Works on [`ProductManifold`](@ref).
+"""
+struct ProductRetraction{TR<:Tuple} <: AbstractRetractionMethod
+    retractions::TR
+end
+
+ProductRetraction(retractions::AbstractRetractionMethod...) = ProductRetraction{typeof(retractions)}(retractions)
+
+function retract!(M::ProductManifold, y, x, v, method::ProductRetraction)
+    map(retract!, M.manifolds, y.parts, x.parts, v.parts, method.retractions)
+    return y
+end
+
+struct InverseProductRetraction{TR<:Tuple} <: AbstractInverseRetractionMethod
+    inverse_retractions::TR
+end
+
+"""
+    InverseProductRetraction(inverse_retractions::AbstractInverseRetractionMethod...)
+
+Product inverse retraction of `inverse_retractions`.
+Works on [`ProductManifold`](@ref).
+"""
+InverseProductRetraction(inverse_retractions::AbstractInverseRetractionMethod...) = InverseProductRetraction{typeof(inverse_retractions)}(inverse_retractions)
+
+function inverse_retract!(M::ProductManifold, v, x, y, method::InverseProductRetraction)
+    map(inverse_retract!, M.manifolds, v.parts, x.parts, y.parts, method.inverse_retractions)
+    return v
+end
+
+"""
+    is_manifold_point(M::ProductManifold, x; kwargs...)
+
+Check whether `x` is a valid point on the [`ProductManifold`](@ref) `M`.
+
+The tolerance for the last test can be set using the ´kwargs...`.
+"""
+function is_manifold_point(M::ProductManifold, x::MPoint; kwargs...)
+    return all(t -> is_manifold_point(t...; kwargs...), ziptuples(M.manifolds, x.parts))
+end
+
+function is_manifold_point(M::ProductManifold, x::ProductArray; kwargs...)
+    return all(t -> is_manifold_point(t...; kwargs...), ziptuples(M.manifolds, x.parts))
+end
+
+"""
+    is_tangent_vector(M::ProductManifold, x, v; kwargs... )
+
+Check whether `v` is a tangent vector to `x` on the [`ProductManifold`](@ref)
+`M`, i.e. atfer [`is_manifold_point`](@ref)`(M, x)`, and all projections to
+base manifolds must be respective tangent vectors.
+
+The tolerance for the last test can be set using the ´kwargs...`.
+"""
+function is_tangent_vector(M::ProductManifold, x::MPoint, v::TVector; kwargs...)
+    is_manifold_point(M, x)
+    return all(t -> is_tangent_vector(t...; kwargs...), ziptuples(M.manifolds, x.parts, v.parts))
+end
+
+function is_tangent_vector(M::ProductManifold, x::ProductArray, v::ProductArray; kwargs...)
+    is_manifold_point(M, x)
+    return all(t -> is_tangent_vector(t...; kwargs...), ziptuples(M.manifolds, x.parts, v.parts))
+end
+
+"""
+    ProductPointDistribution(M::ProductManifold, distributions)
+
+Product distribution on manifold `M`, combined from `distributions`.
+"""
+struct ProductPointDistribution{TM<:ProductManifold, TD<:(NTuple{N,Distribution} where N)} <: MPointDistribution{TM}
+    manifold::TM
+    distributions::TD
+end
+
+function ProductPointDistribution(M::ProductManifold, distributions::MPointDistribution...)
+    return ProductPointDistribution{typeof(M), typeof(distributions)}(M, distributions)
+end
+
+function ProductPointDistribution(distributions::MPointDistribution...)
+    M = ProductManifold(map(d -> support(d).manifold, distributions)...)
+    return ProductPointDistribution(M, distributions...)
+end
+
+function support(d::ProductPointDistribution)
+    return MPointSupport(d.manifold)
+end
+
+function rand(rng::AbstractRNG, d::ProductPointDistribution)
+    return ProductMPoint(map(d -> rand(rng, d), d.distributions)...)
+end
+
+function _rand!(rng::AbstractRNG, d::ProductPointDistribution, x::AbstractArray{<:Number})
+    x .= rand(rng, d)
+    return x
+end
+
+function _rand!(rng::AbstractRNG, d::ProductPointDistribution, x::ProductMPoint)
+    map(t -> _rand!(rng, t[1], t[2]), d.distributions, x.parts)
+    return x
+end
+
+"""
+    ProductTVectorDistribution([m::ProductManifold], [x], distrs...)
+
+Generates a random tangent vector at point `x` from manifold `m` using the
+product distribution of given distributions.
+
+Manifold and `x` can be automatically inferred from distributions `distrs`.
+"""
+struct ProductTVectorDistribution{TM<:ProductManifold, TD<:(NTuple{N,Distribution} where N), TX} <: TVectorDistribution{TM, TX}
+    manifold::TM
+    x::TX
+    distributions::TD
+end
+
+function ProductTVectorDistribution(M::ProductManifold, x::Union{AbstractArray, MPoint}, distributions::TVectorDistribution...)
+    return ProductTVectorDistribution{typeof(M), typeof(distributions), typeof(x)}(M, x, distributions)
+end
+
+function ProductTVectorDistribution(M::ProductManifold, distributions::TVectorDistribution...)
+    x = ProductMPoint(map(d -> support(d).x, distributions))
+    return ProductTVectorDistribution(M, x, distributions...)
+end
+
+function ProductTVectorDistribution(distributions::TVectorDistribution...)
+    M = ProductManifold(map(d -> support(d).manifold, distributions)...)
+    x = ProductMPoint(map(d -> support(d).x, distributions)...)
+    return ProductTVectorDistribution(M, x, distributions...)
+end
+
+function support(tvd::ProductTVectorDistribution)
+    return TVectorSupport(tvd.manifold, ProductMPoint(map(d -> support(d).x, tvd.distributions)...))
+end
+
+function rand(rng::AbstractRNG, d::ProductTVectorDistribution)
+    return ProductTVector(map(d -> rand(rng, d), d.distributions)...)
+end
+
+function _rand!(rng::AbstractRNG, d::ProductTVectorDistribution, v::AbstractArray{<:Number})
+    v .= rand(rng, d)
+    return v
+end
+
+function _rand!(rng::AbstractRNG, d::ProductTVectorDistribution, v::ProductTVector)
+    map(t -> _rand!(rng, t[1], t[2]), d.distributions, v.parts)
+    return v
+end
diff --git a/src/ProjectedDistribution.jl b/src/ProjectedDistribution.jl
index fa1045d09f..5e7a4268c6 100644
--- a/src/ProjectedDistribution.jl
+++ b/src/ProjectedDistribution.jl
@@ -16,6 +16,10 @@ function ProjectedPointDistribution(M::Manifold, d::Distribution, proj!, x::TRes
     return ProjectedPointDistribution{TResult, typeof(M), typeof(d), typeof(proj!)}(M, d, proj!)
 end
 
+function support(d::ProjectedPointDistribution)
+    return MPointSupport(d.manifold)
+end
+
 function rand(rng::AbstractRNG, d::ProjectedPointDistribution{TResult}) where TResult
     x = convert(TResult, rand(rng, d.d))
     d.proj!(d.manifold, x)
@@ -36,7 +40,7 @@ Generates a random tangent vector in ambient space of `m` and projects it
 to `m` using function `project_tangent!`.
 Generated arrays are of type `TResult`.
 """
-struct ProjectedTVectorDistribution{TResult, TM<:Manifold, TX, TD<:Distribution, TProj} <: MPointDistribution{TM}
+struct ProjectedTVectorDistribution{TResult, TM<:Manifold, TX, TD<:Distribution, TProj} <: TVectorDistribution{TM, TX}
     manifold::TM
     x::TX
     d::TD
diff --git a/src/Rotations.jl b/src/Rotations.jl
index c13054415c..18e96f19db 100644
--- a/src/Rotations.jl
+++ b/src/Rotations.jl
@@ -13,6 +13,11 @@ generates the $\mathrm{SO}(n) \subset \mathbb R^{n\times n}$
 struct Rotations{N} <: Manifold end
 Rotations(n::Int) = Rotations{n}()
 
+function representation_size(::Rotations{N}, ::Type{T}) where {N,T<:Union{MPoint, TVector, CoTVector}}
+    return (N, N)
+end
+
+
 @doc doc"""
     manifold_dimension(S::Rotations)
 
@@ -33,13 +38,40 @@ Tangent vectors are represented by matrices.
 """
 inner(S::Rotations, x, w, v) = dot(w, v)
 
-project_tangent!(S::Rotations, w, x, v) = w .= (transpose(v).-v)./2
+project_tangent!(S::Rotations, w, x, v) = w .= (v .- transpose(v))./2
 
 function exp!(S::Rotations, y, x, v)
     y .= x * exp(v)
     return y
 end
 
+function exp!(S::Rotations{2}, y, x, v)
+    α = sign(v[1,2]) * norm(S, x, v) / sqrt(2)
+    @assert size(y) == (2, 2)
+    @assert size(x) == (2, 2)
+    @inbounds begin
+        sinα, cosα = sincos(α)
+        y[1] = x[1]*cosα - x[3]*sinα
+        y[2] = x[2]*cosα - x[4]*sinα
+        y[3] = x[1]*sinα + x[3]*cosα
+        y[4] = x[2]*sinα + x[4]*cosα
+    end
+    return y
+end
+
+function exp!(S::Rotations{3}, y, x, v)
+    α = norm(S, x, v) / sqrt(2)
+    if α ≈ 0
+        a = 1 - α^2 / 6
+        b = α / 2
+    else
+        a = sin(α) / α
+        b = (1-cos(α)) / α^2
+    end
+    y .= x .+ x * (a .* v .+ b .* (v^2))
+    return y
+end
+
 @doc doc"""
     log!(M, v, x, y)
 
@@ -56,11 +88,28 @@ and save the result to `v`.
 For antipodal rotations the function returns one of the tangent vectors that
 point at `y`.
 """
-function log!(S::Rotations, v::TV, x, y) where TV
+function log!(S::Rotations, v, x, y)
+    U = transpose(x) * y
+    v .= real(log(Matrix(U)))
+    project_tangent!(S, v, x, v)
+    return v
+end
+
+function log!(S::Rotations{2}, v, x, y)
     U = transpose(x) * y
-    # MMatrix doesn't have `log` defined
-    U1 = TV(real(log(Array(U))))
-    v .= (U1 .- transpose(U1))./2
+    α = asin(U[1,2])
+    v .= [0 α; -α 0]
+    return v
+end
+
+function log!(S::Rotations{3}, v, x, y)
+    U = transpose(x) * y
+    θ = acos(clamp((tr(U)-1)/2, -1, 1))
+    if θ ≈ 0
+        fill!(v, 0)
+    else
+        v .= θ/(2*sin(θ)) .* (U .- transpose(U))
+    end
     return v
 end
 
diff --git a/src/Sphere.jl b/src/Sphere.jl
index f86e9a7e85..d55f5aa12e 100644
--- a/src/Sphere.jl
+++ b/src/Sphere.jl
@@ -15,6 +15,10 @@ Sphere(n::Int) = Sphere{n}()
 
 @traitimpl HasMetric{Sphere,EuclideanMetric}
 
+function representation_size(::Sphere{N}, ::Type{T}) where {N,T<:Union{MPoint, TVector, CoTVector}}
+    return (N+1,)
+end
+
 @doc doc"""
     manifold_dimension(S::Sphere)
 
@@ -22,7 +26,7 @@ Return the dimension of the manifold $\mathbb S^n$, i.e. $n$.
 """
 manifold_dimension(S::Sphere{N}) where {N} = N
 
-proj!(S::Sphere, x) = (x ./= norm(x))
+project_point!(S::Sphere, x) = (x ./= norm(x))
 
 project_tangent!(S::Sphere, w, x, v) = (w .= v .- dot(x, v) .* x)
 
@@ -33,7 +37,7 @@ compute the inner product of the two tangent vectors `w,v` from the tangent
 plane at `x` on the sphere `S=`$\mathbb S^n$ using the restriction of the
 metric from the embedding, i.e. $ (v,w)_x = v^\mathrm{T}w $.
 """
-inner(S::Sphere, x, w, v) = dot(w, v)
+@inline inner(S::Sphere, x, w, v) = dot(w, v)
 
 norm(S::Sphere, x, v) = norm(v)
 
@@ -62,18 +66,20 @@ end
 
 injectivity_radius(S::Sphere, args...) = π
 
-zero_tangent_vector(S::Sphere, x) = zero(x)
-zero_tangent_vector!(S::Sphere, v, x) = (v .= zero(x))
+function zero_tangent_vector!(S::Sphere, v, x)
+    fill!(v, 0)
+    return v
+end
 
 """
-    uniform_sphere_distribution(S::Sphere, x)
+    uniform_distribution(S::Sphere, x)
 
 Uniform distribution on given sphere. Generated points will be of similar
 type to `x`.
 """
 function uniform_distribution(S::Sphere, x)
     d = Distributions.MvNormal(zero(x), 1.0)
-    return ProjectedPointDistribution(S, d, proj!, x)
+    return ProjectedPointDistribution(S, d, project_point!, x)
 end
 
 """
@@ -116,7 +122,7 @@ function is_tangent_vector(S::Sphere{N},x,v; kwargs...) where N
 end
 
 """
-    gaussian_sphere_tvector_distribution(S::Sphere, x, σ)
+    normal_tvector_distribution(S::Sphere, x, σ)
 
 Normal distribution in ambient space with standard deviation `σ`
 projected to tangent space at `x`.
diff --git a/src/utils.jl b/src/utils.jl
new file mode 100644
index 0000000000..d9e533fd47
--- /dev/null
+++ b/src/utils.jl
@@ -0,0 +1,174 @@
+
+"""
+    select_from_tuple(t::NTuple{N, Any}, positions::Val{P})
+
+Selects elements of tuple `t` at positions specified by the second argument.
+For example `select_from_tuple(("a", "b", "c"), Val((3, 1, 1)))` returns
+`("c", "a", "a")`.
+"""
+@generated function select_from_tuple(t::NTuple{N, Any}, positions::Val{P}) where {N, P}
+    for k in P
+        (k < 0 || k > N) && error("positions must be between 1 and $N")
+    end
+    return Expr(:tuple, [Expr(:ref, :t, k) for k in P]...)
+end
+
+"""
+    ziptuples(a, b)
+
+Zips tuples `a` and `b` in a fast, type-stable way. If they have different
+lengths, the result is trimmed to the length of the shorter tuple.
+"""
+@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}) where {N,M}
+    ex = Expr(:tuple)
+    for i = 1:min(N, M)
+        push!(ex.args, :((a[$i], b[$i])))
+    end
+    ex
+end
+
+"""
+    ziptuples(a, b, c)
+
+Zips tuples `a`, `b` and `c` in a fast, type-stable way. If they have different
+lengths, the result is trimmed to the length of the shorter tuple.
+"""
+@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}, c::NTuple{L,Any}) where {N,M,L}
+    ex = Expr(:tuple)
+    for i = 1:min(N, M, L)
+        push!(ex.args, :((a[$i], b[$i], c[$i])))
+    end
+    ex
+end
+
+"""
+    ziptuples(a, b, c, d)
+
+Zips tuples `a`, `b`, `c` and `d` in a fast, type-stable way. If they have
+different lengths, the result is trimmed to the length of the shorter tuple.
+"""
+@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}, c::NTuple{L,Any}, d::NTuple{K,Any}) where {N,M,L,K}
+    ex = Expr(:tuple)
+    for i = 1:min(N, M, L, K)
+        push!(ex.args, :((a[$i], b[$i], c[$i], d[$i])))
+    end
+    ex
+end
+
+"""
+    ziptuples(a, b, c, d, e)
+
+Zips tuples `a`, `b`, `c`, `d` and `e` in a fast, type-stable way. If they have
+different lengths, the result is trimmed to the length of the shorter tuple.
+"""
+@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}, c::NTuple{L,Any}, d::NTuple{K,Any}, e::NTuple{J,Any}) where {N,M,L,K,J}
+    ex = Expr(:tuple)
+    for i = 1:min(N, M, L, K, J)
+        push!(ex.args, :((a[$i], b[$i], c[$i], d[$i], e[$i])))
+    end
+    ex
+end
+
+# conversion of SizedArray from StaticArrays.jl
+"""
+    SizedAbstractArray{Tuple{dims...}}(array)
+
+Wraps an `Array` with a static size, so to take advantage of the (faster)
+methods defined by the static array package. The size is checked once upon
+construction to determine if the number of elements (`length`) match, but the
+array may be reshaped.
+(Also, `Size(dims...)(array)` acheives the same thing)
+"""
+struct SizedAbstractArray{S<:Tuple, T, N, M, TData<:AbstractArray{T,M}} <: StaticArray{S, T, N}
+    data::TData
+
+    function SizedAbstractArray{S, T, N, M, TData}(a::TData) where {S, T, N, M, TData<:AbstractArray}
+        if length(a) != StaticArrays.tuple_prod(S)
+            error("Dimensions $(size(a)) don't match static size $S")
+        end
+        new{S,T,N,M,TData}(a)
+    end
+
+    function SizedAbstractArray{S, T, N, 1}(::UndefInitializer) where {S, T, N, TData<:AbstractArray}
+        new{S, T, N, 1, Array{T, 1}}(Array{T, 1}(undef, StaticArrays.tuple_prod(S)))
+    end
+    function SizedAbstractArray{S, T, N, N}(::UndefInitializer) where {S, T, N, TData<:AbstractArray}
+        new{S, T, N, N, Array{T, N}}(Array{T, N}(undef, size_to_tuple(S)...))
+    end
+end
+
+@inline SizedAbstractArray{S,T,N}(a::TData) where {S,T,N,M,TData<:AbstractArray{T,M}} = SizedAbstractArray{S,T,N,M,TData}(a)
+@inline SizedAbstractArray{S,T}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}} = SizedAbstractArray{S,T,StaticArrays.tuple_length(S),M,TData}(a)
+@inline SizedAbstractArray{S}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}} = SizedAbstractArray{S,T,StaticArrays.tuple_length(S),M,TData}(a)
+
+@inline SizedAbstractArray{S,T,N}(::UndefInitializer) where {S,T,N} = SizedAbstractArray{S,T,N,N}(undef)
+@inline SizedAbstractArray{S,T}(::UndefInitializer) where {S,T} = SizedAbstractArray{S,T,StaticArrays.tuple_length(S),StaticArrays.tuple_length(S)}(undef)
+@generated function (::Type{SizedAbstractArray{S,T,N,M,TData}})(x::NTuple{L,Any}) where {S,T,N,M,TData,L}
+    if L != StaticArrays.tuple_prod(S)
+        error("Dimension mismatch")
+    end
+    exprs = [:(a[$i] = x[$i]) for i = 1:L]
+    return quote
+        $(Expr(:meta, :inline))
+        a = SizedAbstractArray{S,T,N,M}(undef)
+        @inbounds $(Expr(:block, exprs...))
+        return a
+    end
+end
+
+@inline SizedAbstractArray{S,T,N}(x::Tuple) where {S,T,N} = SizedAbstractArray{S,T,N,N}(x)
+@inline SizedAbstractArray{S,T}(x::Tuple) where {S,T} = SizedAbstractArray{S,T,StaticArrays.tuple_length(S),StaticArrays.tuple_length(S)}(x)
+@inline SizedAbstractArray{S}(x::NTuple{L,T}) where {S,T,L} = SizedAbstractArray{S,T,StaticArrays.tuple_length(S),StaticArrays.tuple_length(S)}(x)
+
+# Overide some problematic default behaviour
+@inline convert(::Type{SA}, sa::SizedAbstractArray) where {SA<:SizedAbstractArray} = SA(sa.data)
+@inline convert(::Type{SA}, sa::SA) where {SA<:SizedAbstractArray} = sa
+
+# Back to Array (unfortunately need both convert and construct to overide other methods)
+import Base.Array
+@inline Array(sa::SizedAbstractArray) = Array(sa.data)
+@inline Array{T}(sa::SizedAbstractArray{S,T}) where {T,S} = Array(sa.data)
+@inline Array{T,N}(sa::SizedAbstractArray{S,T,N}) where {T,S,N} = Array(sa.data)
+
+@inline convert(::Type{AbstractArray}, sa::SizedAbstractArray) = sa.data
+@inline convert(::Type{AbstractArray{T}}, sa::SizedAbstractArray{S,T}) where {T,S} = sa.data
+@inline convert(::Type{AbstractArray{T,N}}, sa::SizedAbstractArray{S,T,N}) where {T,S,N} = sa.data
+
+Base.@propagate_inbounds getindex(a::SizedAbstractArray, i::Int) = getindex(a.data, i)
+Base.@propagate_inbounds setindex!(a::SizedAbstractArray, v, i::Int) = setindex!(a.data, v, i)
+
+SizedAbstractVector{S,T,M} = SizedAbstractArray{Tuple{S},T,1,M}
+@inline SizedAbstractVector{S}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}} = SizedAbstractArray{Tuple{S},T,1,M,TData}(a)
+@inline SizedAbstractVector{S}(x::NTuple{L,T}) where {S,T,L} = SizedAbstractArray{Tuple{S},T,1,1,Vector{T}}(x)
+
+SizedAbstractMatrix{S1,S2,T,M} = SizedAbstractArray{Tuple{S1,S2},T,2,M}
+@inline SizedAbstractMatrix{S1,S2}(a::TData) where {S1,S2,T,M,TData<:AbstractArray{T,M}} = SizedAbstractArray{Tuple{S1,S2},T,2,M,TData}(a)
+@inline SizedAbstractMatrix{S1,S2}(x::NTuple{L,T}) where {S1,S2,T,L} = SizedAbstractArray{Tuple{S1,S2},T,2,2,Matrix{2}}(x)
+
+"""
+    Size(dims)(array)
+
+Creates a `SizedAbstractArray` wrapping `array` with the specified statically-known
+`dims`, so to take advantage of the (faster) methods defined by the static array
+package.
+"""
+(::Size{S})(a::AbstractArray) where {S} = SizedAbstractArray{Tuple{S...}}(a)
+
+function promote_rule(::Type{<:SizedAbstractArray{S,T,N,M,TDataA}}, ::Type{<:SizedAbstractArray{S,U,N,M,TDataB}}) where {S,T,U,N,M,TDataA,TDataB}
+    SizedAbstractArray{S,promote_type(T,U),N,M,T,promote_type(TDataA, TDataB)}
+end
+
+@inline copy(a::SizedAbstractArray) = typeof(a)(copy(a.data))
+
+similar(::Type{<:SizedAbstractArray{S,T,N,M}},::Type{T2}) where {S,T,N,M,T2} = SizedAbstractArray{S,T2,N,M}(undef)
+similar(::Type{SA},::Type{T},s::Size{S}) where {SA<:SizedAbstractArray,T,S} = sizedabstractarray_similar_type(T,s,StaticArrays.length_val(s))(undef)
+sizedabstractarray_similar_type(::Type{T},s::Size{S},::Type{Val{D}}) where {T,S,D} = SizedAbstractArray{Tuple{S...},T,D,length(s)}
+
+"""
+    size_to_tuple(::Type{S}) where S<:Tuple
+
+Converts a size given by `Tuple{N, M, ...}` into a tuple `(N, M, ...)`.
+"""
+@generated function size_to_tuple(::Type{S}) where S<:Tuple
+    return tuple(S.parameters...)
+end
diff --git a/test/euclidean.jl b/test/euclidean.jl
new file mode 100644
index 0000000000..1482b14cae
--- /dev/null
+++ b/test/euclidean.jl
@@ -0,0 +1,28 @@
+include("utils.jl")
+
+@testset "Euclidean" begin
+    M = Manifolds.Euclidean(3)
+    types = [Vector{Float64},
+             SizedVector{3, Float64},
+             MVector{3, Float64},
+             Vector{Float32},
+             SizedVector{3, Float32},
+             MVector{3, Float32},
+             Vector{Double64},
+             MVector{3, Double64},
+             SizedVector{3, Double64}]
+    for T in types
+        @testset "Type $T" begin
+            pts = [convert(T, [1.0, 0.0, 0.0]),
+                   convert(T, [0.0, 1.0, 0.0]),
+                   convert(T, [0.0, 0.0, 1.0])]
+            test_manifold(M,
+                          pts,
+                          test_reverse_diff = isa(T, Vector),
+                          test_project_tangent = true,
+                          point_distributions = [Manifolds.projected_distribution(M, Distributions.MvNormal(zero(pts[1]), 1.0))],
+                          tvector_distributions = [Manifolds.normal_tvector_distribution(M, pts[1], 1.0)])
+        end
+    end
+
+end
diff --git a/test/metric_test.jl b/test/metric_test.jl
index 736cfcec94..17b659046e 100644
--- a/test/metric_test.jl
+++ b/test/metric_test.jl
@@ -1,7 +1,8 @@
 
+struct TestEuclidean{N} <: Manifold end
+struct TestEuclideanMetric <: Metric end
+
 @testset "scaled Euclidean metric" begin
-    struct TestEuclidean{N} <: Manifold end
-    struct TestEuclideanMetric <: Metric end
     Manifolds.manifold_dimension(::TestEuclidean{N}) where {N} = N
     function Manifolds.local_metric(M::MetricManifold{<:TestEuclidean,<:TestEuclideanMetric}, x)
         return Diagonal(1.0:manifold_dimension(M))
@@ -41,11 +42,13 @@
     end
 end
 
+struct TestSphere{N,T} <: Manifold
+    r::T
+end
+
+struct TestSphericalMetric <: Metric end
+
 @testset "scaled Sphere metric" begin
-    struct TestSphere{N,T} <: Manifold
-        r::T
-    end
-    struct TestSphericalMetric <: Metric end
     Manifolds.manifold_dimension(::TestSphere{N}) where {N} = N
     function Manifolds.local_metric(M::MetricManifold{<:TestSphere,<:TestSphericalMetric}, x)
         r = base_manifold(M).r
@@ -136,9 +139,10 @@ end
     end
 end
 
+struct BaseManifold{N} <: Manifold end
+struct BaseManifoldMetric{M} <: Metric end
+
 @testset "HasMetric trait" begin
-    struct BaseManifold{N} <: Manifold end
-    struct BaseManifoldMetric{M} <: Metric end
     Manifolds.manifold_dimension(::BaseManifold{N}) where {N} = N
     @traitimpl HasMetric{BaseManifold,BaseManifoldMetric}
     Manifolds.inner(::BaseManifold, x, v, w) = 2 * dot(v,w)
@@ -162,12 +166,12 @@ end
     @test exp!(MM, y, x, v) === exp!(M, y, x, v)
     @test log(M, x, y) == (y - x) / 2
     @test log!(MM, v, x, y) === log!(M, v, x, y)
-    @test Manifolds.retract!(MM, y, x, v) === Manifolds.retract!(M, y, x, v)
-    @test Manifolds.retract!(MM, y, x, v, 1) === Manifolds.retract!(M, y, x, v, 1)
-    @test Manifolds.project_tangent!(MM, w, x, v) === Manifolds.project_tangent!(M, w, x, v)
+    @test retract!(MM, y, x, v) === retract!(M, y, x, v)
+    @test retract!(MM, y, x, v, 1) === retract!(M, y, x, v, 1)
+    @test project_tangent!(MM, w, x, v) === project_tangent!(M, w, x, v)
     @test zero_tangent_vector!(MM, v, x) === zero_tangent_vector!(M, v, x)
     @test injectivity_radius(MM, x) === injectivity_radius(M, x)
     @test injectivity_radius(MM) === injectivity_radius(M)
-    @test Manifolds.is_manifold_point(MM, x) === Manifolds.is_manifold_point(M, x)
-    @test Manifolds.is_tangent_vector(MM, x, v) === Manifolds.is_tangent_vector(M, x, v)
+    @test is_manifold_point(MM, x) === is_manifold_point(M, x)
+    @test is_tangent_vector(MM, x, v) === is_tangent_vector(M, x, v)
 end
diff --git a/test/product_manifold.jl b/test/product_manifold.jl
new file mode 100644
index 0000000000..2baede50a7
--- /dev/null
+++ b/test/product_manifold.jl
@@ -0,0 +1,110 @@
+include("utils.jl")
+
+@testset "Product manifold" begin
+    M1 = Sphere(2)
+    M2 = Euclidean(2)
+    Mse = ProductManifold(M1, M2)
+    @test Mse == M1 × M2
+    shape_se = Manifolds.ShapeSpecification(M1, M2)
+
+    types = [Vector{Float64},
+             SizedVector{5, Float64},
+             MVector{5, Float64},
+             Vector{Float32},
+             SizedVector{5, Float32},
+             MVector{5, Float32}]
+
+    retraction_methods = [Manifolds.ProductRetraction(Manifolds.ExponentialRetraction(), Manifolds.ExponentialRetraction())]
+    inverse_retraction_methods = [Manifolds.InverseProductRetraction(Manifolds.LogarithmicInverseRetraction(), Manifolds.LogarithmicInverseRetraction())]
+    for T in types
+        @testset "Type $T" begin
+            pts_base = [convert(T, [1.0, 0.0, 0.0, 0.0, 0.0]),
+                        convert(T, [0.0, 1.0, 0.0, 1.0, 0.0]),
+                        convert(T, [0.0, 0.0, 1.0, 0.0, 0.1])]
+            pts = map(p -> Manifolds.ProductArray(shape_se, p), pts_base)
+            distr_M1 = Manifolds.uniform_distribution(M1, pts_base[1][1:3])
+            distr_M2 = Manifolds.projected_distribution(M2, Distributions.MvNormal(zero(pts_base[1][4:5]), 1.0))
+            distr_tv_M1 = Manifolds.normal_tvector_distribution(M1, pts_base[1][1:3], 1.0)
+            distr_tv_M2 = Manifolds.normal_tvector_distribution(M2, pts_base[1][4:5], 1.0)
+            test_manifold(Mse,
+                          pts;
+                          test_reverse_diff = isa(T, Vector),
+                          retraction_methods = retraction_methods,
+                          inverse_retraction_methods = inverse_retraction_methods,
+                          point_distributions = [Manifolds.ProductPointDistribution(distr_M1, distr_M2)],
+                          tvector_distributions = [Manifolds.ProductTVectorDistribution(distr_tv_M1, distr_tv_M2)])
+        end
+    end
+
+    M3 = Manifolds.Rotations(2)
+    Mser = ProductManifold(M1, M2, M3)
+    shape_ser = Manifolds.ShapeSpecification(M1, M2, M3)
+
+    @test submanifold(Mser, 2) == M2
+    @test submanifold(Mser, Val((1, 3))) == M1 × M3
+    @test submanifold(Mser, 2:3) == M2 × M3
+    @test submanifold(Mser, [1, 3]) == M1 × M3
+    @inferred submanifold(Mser, Val((1, 3)))
+
+    # testing the slower generic constructor
+    Mprod4 = ProductManifold(M2, M2, M2, M2)
+    shape4 = Manifolds.ShapeSpecification(Mprod4.manifolds...)
+    data_x4 = rand(8)
+    @test Manifolds.ProductArray(shape4, data_x4).data === data_x4
+
+    pts_sphere = [[1.0, 0.0, 0.0],
+                  [0.0, 1.0, 0.0],
+                  [0.0, 0.0, 1.0]]
+    pts_r2 = [[0.0, 0.0],
+              [1.0, 0.0],
+              [0.0, 0.1]]
+    angles = (0.0, π/2, 2π/3)
+    pts_rot = [[cos(ϕ) sin(ϕ); -sin(ϕ) cos(ϕ)] for ϕ in angles]
+    pts = [Manifolds.prod_point(shape_ser, p[1], p[2], p[3]) for p in zip(pts_sphere, pts_r2, pts_rot)]
+    test_manifold(Mser,
+                  pts,
+                  test_forward_diff = false,
+                  test_reverse_diff = false)
+
+    @testset "prod_point" begin
+        Ts = SizedVector{3, Float64}
+        Tr2 = SizedVector{2, Float64}
+        T = SizedVector{5, Float64}
+        pts_base = [convert(T, [1.0, 0.0, 0.0, 0.0, 0.0]),
+                    convert(T, [0.0, 1.0, 0.0, 1.0, 0.0]),
+                    convert(T, [0.0, 0.0, 1.0, 0.0, 0.1])]
+        pts = map(p -> Manifolds.ProductArray(shape_se, p), pts_base)
+        pts_sphere = [convert(Ts, [1.0, 0.0, 0.0]),
+                      convert(Ts, [0.0, 1.0, 0.0]),
+                      convert(Ts, [0.0, 0.0, 1.0])]
+        pts_r2 = [convert(Tr2, [0.0, 0.0]),
+                  convert(Tr2, [1.0, 0.0]),
+                  convert(Tr2, [0.0, 0.1])]
+        pts_prod = [Manifolds.prod_point(shape_se, p[1], p[2]) for p in zip(pts_sphere, pts_r2)]
+        for p in zip(pts, pts_prod)
+            @test isapprox(Mse, p[1], p[2])
+        end
+        for p in zip(pts_sphere, pts_r2, pts_prod)
+            @test isapprox(M1, p[1], submanifold_component(p[3], 1))
+            @test isapprox(M2, p[2], submanifold_component(p[3], 2))
+        end
+    end
+
+    @testset "ProductMPoint" begin
+        Ts = SizedVector{3, Float64}
+        Tr2 = SizedVector{2, Float64}
+        pts_sphere = [convert(Ts, [1.0, 0.0, 0.0]),
+                      convert(Ts, [0.0, 1.0, 0.0]),
+                      convert(Ts, [0.0, 0.0, 1.0])]
+        pts_r2 = [convert(Tr2, [0.0, 0.0]),
+                  convert(Tr2, [1.0, 0.0]),
+                  convert(Tr2, [0.0, 0.1])]
+
+        pts = [ProductMPoint(p[1], p[2]) for p in zip(pts_sphere, pts_r2)]
+        test_manifold(Mse,
+                      pts,
+                      test_tangent_vector_broadcasting = false,
+                      test_forward_diff = false,
+                      test_reverse_diff = false)
+    end
+end
diff --git a/test/rotations.jl b/test/rotations.jl
new file mode 100644
index 0000000000..fdb77527b0
--- /dev/null
+++ b/test/rotations.jl
@@ -0,0 +1,68 @@
+include("utils.jl")
+
+@testset "Rotations" begin
+    M = Manifolds.Rotations(2)
+
+    types = [Matrix{Float64},
+             SizedMatrix{2, 2, Float64},
+             MMatrix{2, 2, Float64},
+             Matrix{Float32},
+             SizedMatrix{2, 2, Float32},
+             MMatrix{2, 2, Float32}]
+
+    retraction_methods = [Manifolds.PolarRetraction(),
+                          Manifolds.QRRetraction()]
+
+    inverse_retraction_methods = [Manifolds.PolarInverseRetraction(),
+                                  Manifolds.QRInverseRetraction()]
+
+    for T in types
+        angles = (0.0, π/2, 2π/3, π/4)
+        pts = [convert(T, [cos(ϕ) sin(ϕ); -sin(ϕ) cos(ϕ)]) for ϕ in angles]
+        test_manifold(M, pts;
+            test_reverse_diff = false,
+            test_project_tangent = true,
+            retraction_methods = retraction_methods,
+            inverse_retraction_methods = inverse_retraction_methods,
+            point_distributions = [Manifolds.normal_rotation_distribution(M, pts[1], 1.0)],
+            tvector_distributions = [Manifolds.normal_tvector_distribution(M, pts[1], 1.0)])
+
+        v = log(M, pts[1], pts[2])
+        @test norm(M, pts[1], v) ≈ (angles[2] - angles[1])*sqrt(2)
+
+        v14_polar = inverse_retract(M, pts[1], pts[4], Manifolds.PolarInverseRetraction())
+        p4_polar = retract(M, pts[1], v14_polar, Manifolds.PolarRetraction())
+        @test isapprox(M, pts[4], p4_polar)
+
+        v14_qr = inverse_retract(M, pts[1], pts[4], Manifolds.QRInverseRetraction())
+        p4_qr = retract(M, pts[1], v14_qr, Manifolds.QRRetraction())
+        @test isapprox(M, pts[4], p4_qr)
+    end
+
+    @testset "Distribution tests" begin
+        usd_mmatrix = Manifolds.normal_rotation_distribution(M, (@MMatrix [1.0 0.0; 0.0 1.0]), 1.0)
+        @test isa(rand(usd_mmatrix), MMatrix)
+
+        gtsd_mvector = Manifolds.normal_tvector_distribution(M, (@MMatrix [1.0 0.0; 0.0 1.0]), 1.0)
+        @test isa(rand(gtsd_mvector), MMatrix)
+    end
+
+    Random.seed!(42)
+    for n ∈ (3, 4, 5)
+        @testset "Rotations: SO($n)" begin
+            SOn = Manifolds.Rotations(n)
+            ptd = Manifolds.normal_rotation_distribution(SOn, Matrix(1.0I, n, n), 1.0)
+            tvd = Manifolds.normal_tvector_distribution(SOn, Matrix(1.0I, n, n), 1.0)
+            pts = [rand(ptd) for _ in 1:3]
+            test_manifold(SOn, pts;
+                test_forward_diff = n==3,
+                test_reverse_diff = false,
+                retraction_methods = retraction_methods,
+                inverse_retraction_methods = inverse_retraction_methods,
+                point_distributions = [ptd],
+                tvector_distributions = [tvd],
+                exp_log_atol_multiplier = 6)
+        end
+    end
+
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index b47cdf7667..f3ee61fc97 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,262 +1,10 @@
-using Manifolds
 
-using LinearAlgebra
-using DoubleFloats
-using ForwardDiff
-using ReverseDiff
-using StaticArrays
-using SimpleTraits
-using Test
+include("utils.jl")
 
-"""
-    test_manifold(m::Manifold, pts::AbstractVector;
-        test_forward_diff = true,
-        retraction_methods = [],
-        inverse_retraction_methods = [])
-
-Tests general properties of manifold `m`, given at least three different points
-that lie on it (contained in `pts`).
-
-# Arguments
-- `test_forward_diff = true`: if true, automatic differentiation using ForwardDiff is
-tested.
-- `retraction_methods = []`: retraction methods that will be tested.
-"""
-function test_manifold(M::Manifold, pts::AbstractVector;
-    test_forward_diff = true,
-    test_reverse_diff = true,
-    test_tangent_vector_broadcasting = true,
-    retraction_methods = [],
-    inverse_retraction_methods = [],
-    point_distributions = [],
-    tvector_distributions = [])
-    # log/exp
-    length(pts) ≥ 3 || error("Not enough points (at least three expected)")
-    isapprox(M, pts[1], pts[2]) && error("Points 1 and 2 are equal")
-    isapprox(M, pts[1], pts[3]) && error("Points 1 and 3 are equal")
-
-    @testset "dimension" begin
-        @test isa(manifold_dimension(M), Integer)
-        @test manifold_dimension(M) ≥ 0
-    end
-
-    @testset "injectivity radius" begin
-        @test injectivity_radius(M, pts[1]) > 0
-        for rm ∈ retraction_methods
-            @test injectivity_radius(M, pts[1], rm) > 0
-            @test injectivity_radius(M, pts[1], rm) ≤ injectivity_radius(M, pts[1])
-        end
-    end
-
-    tv1 = log(M, pts[1], pts[2])
-
-    @testset "is_manifold_point / is_tangent_vector" begin
-        for pt ∈ pts
-            @test is_manifold_point(M, pt)
-        end
-        @test is_tangent_vector(M, pts[1], tv1)
-    end
-
-    @testset "log/exp tests" begin
-        @test isapprox(M, pts[2], exp(M, pts[1], tv1))
-        @test isapprox(M, pts[1], exp(M, pts[1], tv1, 0))
-        @test isapprox(M, pts[2], exp(M, pts[1], tv1, 1))
-        @test is_manifold_point(M, retract(M, pts[1], tv1))
-        @test isapprox(M, pts[1], retract(M, pts[1], tv1, 0))
-        for retr_method ∈ retraction_methods
-            @test is_manifold_point(M, retract(M, pts[1], tv1, retr_method))
-            @test isapprox(M, pts[1], retract(M, pts[1], tv1, 0, retr_method))
-        end
-        new_pt = exp(M, pts[1], tv1)
-        retract!(M, new_pt, pts[1], tv1)
-        @test is_manifold_point(M, new_pt)
-        for x ∈ pts
-            @test isapprox(M, zero_tangent_vector(M, x), log(M, x, x); atol = eps(eltype(x)))
-            @test isapprox(M, zero_tangent_vector(M, x), inverse_retract(M, x, x); atol = eps(eltype(x)))
-            for inv_retr_method ∈ inverse_retraction_methods
-                @test isapprox(M, zero_tangent_vector(M, x), inverse_retract(M, x, x, inv_retr_method); atol = eps(eltype(x)))
-            end
-        end
-        zero_tangent_vector!(M, tv1, pts[1])
-        @test isapprox(M, pts[1], tv1, zero_tangent_vector(M, pts[1]))
-        log!(M, tv1, pts[1], pts[2])
-        @test norm(M, pts[1], tv1) ≈ sqrt(inner(M, pts[1], tv1, tv1))
-
-        @test isapprox(M, exp(M, pts[1], tv1, 1), pts[2])
-        @test isapprox(M, exp(M, pts[1], tv1, 0), pts[1])
-
-        @test distance(M, pts[1], pts[2]) ≈ norm(M, pts[1], tv1)
-    end
-
-    @testset "basic linear algebra in tangent space" begin
-        @test isapprox(M, pts[1], 0*tv1, zero_tangent_vector(M, pts[1]))
-        @test isapprox(M, pts[1], 2*tv1, tv1+tv1)
-        @test isapprox(M, pts[1], 0*tv1, tv1-tv1)
-        @test isapprox(M, pts[1], (-1)*tv1, -tv1)
-    end
-
-    test_tangent_vector_broadcasting && @testset "broadcasted linear algebra in tangent space" begin
-        @test isapprox(M, pts[1], 3*tv1, 2 .* tv1 .+ tv1)
-        @test isapprox(M, pts[1], -tv1, tv1 .- 2 .* tv1)
-        @test isapprox(M, pts[1], -tv1, .-tv1)
-        v = similar(tv1)
-        v .= 2 .* tv1 .+ tv1
-        @test v ≈ 3*tv1
-    end
-
-    test_forward_diff && @testset "ForwardDiff support" begin
-        exp_f(t) = distance(M, pts[1], exp(M, pts[1], t*tv1))
-        d12 = distance(M, pts[1], pts[2])
-        for t ∈ 0.1:0.1:1.0
-            @test d12 ≈ ForwardDiff.derivative(exp_f, t)
-        end
-
-        retract_f(t) = distance(M, pts[1], retract(M, pts[1], t*tv1))
-        for t ∈ 0.1:0.1:1.0
-            @test d12 ≈ ForwardDiff.derivative(retract_f, t)
-        end
-    end
-
-    test_reverse_diff && @testset "ReverseDiff support" begin
-        exp_f(t) = distance(M, pts[1], exp(M, pts[1], t[1]*tv1))
-        d12 = distance(M, pts[1], pts[2])
-        for t ∈ 0.1:0.1:1.0
-            @test d12 ≈ ReverseDiff.gradient(exp_f, [t])[1]
-        end
-
-        retract_f(t) = distance(M, pts[1], retract(M, pts[1], t[1]*tv1))
-        for t ∈ 0.1:0.1:1.0
-            @test d12 ≈ ReverseDiff.gradient(retract_f, [t])[1]
-        end
-    end
-
-    @testset "eltype" begin
-        tv1 = log(M, pts[1], pts[2])
-        @test eltype(tv1) == eltype(pts[1])
-        @test eltype(exp(M, pts[1], tv1)) == eltype(pts[1])
-    end
-
-    @testset "point distributions" begin
-        for pd ∈ point_distributions
-            for _ in 1:10
-                @test is_manifold_point(M, rand(pd))
-            end
-        end
-    end
-
-    @testset "tangent vector distributions" begin
-        for tvd ∈ tvector_distributions
-            supp = Manifolds.support(tvd)
-            for _ in 1:10
-                @test is_tangent_vector(M, supp.x, rand(tvd))
-            end
-        end
-    end
-end
-
-function test_arraymanifold()
-    M = Manifolds.Sphere(2)
-    A = ArrayManifold(M)
-    x = [1., 0., 0.]
-    y = 1/sqrt(2)*[1., 1., 0.]
-    z = [0., 1., 0.]
-    v = log(M,x,y)
-    v2 = log(A,x,y)
-    y2 = exp(A,x,v2)
-    w = log(M,x,z)
-    w2 = log(A,x,z; atol=10^(-15))
-    @test isapprox(y2.value,y)
-    @test distance(A,x,y) == distance(M,x,y)
-    @test norm(A,x,v) == norm(M,x,v)
-    @test inner(A,x,v2,w2; atol=10^(-15)) == inner(M,x,v,w)
-    @test_throws DomainError Manifolds.is_manifold_point(M,2*y)
-    @test_throws DomainError Manifolds.is_tangent_vector(M,y,v; atol=10^(-15))
-
-    test_manifold(A, [x, y, z],
-        test_tangent_vector_broadcasting = false)
-end
-
-@testset "Sphere" begin
-    M = Manifolds.Sphere(2)
-    types = [Vector{Float64},
-             SizedVector{3, Float64},
-             MVector{3, Float64},
-             Vector{Float32},
-             SizedVector{3, Float32},
-             MVector{3, Float32},
-             Vector{Double64},
-             MVector{3, Double64},
-             SizedVector{3, Double64}]
-    for T in types
-        @testset "Type $T" begin
-            pts = [convert(T, [1.0, 0.0, 0.0]),
-                   convert(T, [0.0, 1.0, 0.0]),
-                   convert(T, [0.0, 0.0, 1.0])]
-            test_manifold(M,
-                          pts,
-                          test_reverse_diff = isa(T, Vector),
-                          point_distributions = [Manifolds.uniform_distribution(M, pts[1])],
-                          tvector_distributions = [Manifolds.normal_tvector_distribution(M, pts[1], 1.0)])
-        end
-    end
-
-    @testset "Distribution tests" begin
-        usd_mvector = Manifolds.uniform_distribution(M, @MVector [1.0, 0.0, 0.0])
-        @test isa(rand(usd_mvector), MVector)
-
-        gtsd_mvector = Manifolds.normal_tvector_distribution(M, (@MVector [1.0, 0.0, 0.0]), 1.0)
-        @test isa(rand(gtsd_mvector), MVector)
-    end
-
-    test_arraymanifold()
-end
-
-@testset "Rotations" begin
-    M = Manifolds.Rotations(2)
-
-    types = [Matrix{Float64},
-             SizedMatrix{2, 2, Float64},
-             MMatrix{2, 2, Float64},
-             Matrix{Float32},
-             SizedMatrix{2, 2, Float32},
-             MMatrix{2, 2, Float32}]
-
-    retraction_methods = [Manifolds.PolarRetraction(),
-                          Manifolds.QRRetraction()]
-
-    inverse_retraction_methods = [Manifolds.PolarInverseRetraction(),
-                                  Manifolds.QRInverseRetraction()]
-
-    for T in types
-        angles = (0.0, π/2, 2π/3, π/4)
-        pts = [convert(T, [cos(ϕ) sin(ϕ); -sin(ϕ) cos(ϕ)]) for ϕ in angles]
-        test_manifold(M, pts;
-            test_forward_diff = false,
-            test_reverse_diff = false,
-            retraction_methods = retraction_methods,
-            inverse_retraction_methods = inverse_retraction_methods,
-            point_distributions = [Manifolds.normal_rotation_distribution(M, pts[1], 1.0)],
-            tvector_distributions = [Manifolds.normal_tvector_distribution(M, pts[1], 1.0)])
-
-        v = log(M, pts[1], pts[2])
-        @test norm(M, pts[1], v) ≈ (angles[2] - angles[1])*sqrt(2)
-
-        v14_polar = inverse_retract(M, pts[1], pts[4], Manifolds.PolarInverseRetraction())
-        p4_polar = retract(M, pts[1], v14_polar, Manifolds.PolarRetraction())
-        @test isapprox(M, pts[4], p4_polar)
-
-        v14_qr = inverse_retract(M, pts[1], pts[4], Manifolds.QRInverseRetraction())
-        p4_qr = retract(M, pts[1], v14_qr, Manifolds.QRRetraction())
-        @test isapprox(M, pts[4], p4_qr)
-    end
-
-    @testset "Distribution tests" begin
-        usd_mmatrix = Manifolds.normal_rotation_distribution(M, (@MMatrix [1.0 0.0; 0.0 1.0]), 1.0)
-        @test isa(rand(usd_mmatrix), MMatrix)
-
-        gtsd_mvector = Manifolds.normal_tvector_distribution(M, (@MMatrix [1.0 0.0; 0.0 1.0]), 1.0)
-        @test isa(rand(gtsd_mvector), MMatrix)
-    end
-end
+# starting with tests of simple manifolds
+include("euclidean.jl")
+include("sphere.jl")
+include("rotations.jl")
+include("product_manifold.jl")
 
 include("metric_test.jl")
diff --git a/test/sphere.jl b/test/sphere.jl
new file mode 100644
index 0000000000..ae38367dd2
--- /dev/null
+++ b/test/sphere.jl
@@ -0,0 +1,59 @@
+include("utils.jl")
+
+function test_arraymanifold()
+    M = Sphere(2)
+    A = ArrayManifold(M)
+    x = [1., 0., 0.]
+    y = 1/sqrt(2)*[1., 1., 0.]
+    z = [0., 1., 0.]
+    v = log(M,x,y)
+    v2 = log(A,x,y)
+    y2 = exp(A,x,v2)
+    w = log(M,x,z)
+    w2 = log(A,x,z; atol=10^(-15))
+    @test isapprox(y2.value,y)
+    @test distance(A,x,y) == distance(M,x,y)
+    @test norm(A,x,v) == norm(M,x,v)
+    @test inner(A,x,v2,w2; atol=10^(-15)) == inner(M,x,v,w)
+    @test_throws DomainError is_manifold_point(M,2*y)
+    @test_throws DomainError is_tangent_vector(M,y,v; atol=10^(-15))
+
+    test_manifold(A, [x, y, z],
+        test_tangent_vector_broadcasting = false)
+end
+
+@testset "Sphere" begin
+    M = Sphere(2)
+    types = [Vector{Float64},
+             SizedVector{3, Float64},
+             MVector{3, Float64},
+             Vector{Float32},
+             SizedVector{3, Float32},
+             MVector{3, Float32},
+             Vector{Double64},
+             MVector{3, Double64},
+             SizedVector{3, Double64}]
+    for T in types
+        @testset "Type $T" begin
+            pts = [convert(T, [1.0, 0.0, 0.0]),
+                   convert(T, [0.0, 1.0, 0.0]),
+                   convert(T, [0.0, 0.0, 1.0])]
+            test_manifold(M,
+                          pts,
+                          test_reverse_diff = isa(T, Vector),
+                          test_project_tangent = true,
+                          point_distributions = [Manifolds.uniform_distribution(M, pts[1])],
+                          tvector_distributions = [Manifolds.normal_tvector_distribution(M, pts[1], 1.0)])
+        end
+    end
+
+    @testset "Distribution tests" begin
+        usd_mvector = Manifolds.uniform_distribution(M, @MVector [1.0, 0.0, 0.0])
+        @test isa(rand(usd_mvector), MVector)
+
+        gtsd_mvector = Manifolds.normal_tvector_distribution(M, (@MVector [1.0, 0.0, 0.0]), 1.0)
+        @test isa(rand(gtsd_mvector), MVector)
+    end
+
+    test_arraymanifold()
+end
diff --git a/test/utils.jl b/test/utils.jl
new file mode 100644
index 0000000000..bfb940689c
--- /dev/null
+++ b/test/utils.jl
@@ -0,0 +1,178 @@
+using Manifolds
+
+using LinearAlgebra
+using Distributions
+using DoubleFloats
+using ForwardDiff
+using Random
+using ReverseDiff
+using StaticArrays
+using SimpleTraits
+using Test
+
+"""
+    test_manifold(m::Manifold, pts::AbstractVector;
+        test_forward_diff = true,
+        retraction_methods = [],
+        inverse_retraction_methods = [])
+
+Tests general properties of manifold `m`, given at least three different points
+that lie on it (contained in `pts`).
+
+# Arguments
+- `test_forward_diff = true`: if true, automatic differentiation using ForwardDiff is
+tested.
+- `retraction_methods = []`: retraction methods that will be tested.
+"""
+function test_manifold(M::Manifold, pts::AbstractVector;
+    test_forward_diff = true,
+    test_reverse_diff = true,
+    test_tangent_vector_broadcasting = true,
+    test_project_tangent = false,
+    retraction_methods = [],
+    inverse_retraction_methods = [],
+    point_distributions = [],
+    tvector_distributions = [],
+    exp_log_atol_multiplier = 1)
+    # log/exp
+    length(pts) ≥ 3 || error("Not enough points (at least three expected)")
+    isapprox(M, pts[1], pts[2]) && error("Points 1 and 2 are equal")
+    isapprox(M, pts[1], pts[3]) && error("Points 1 and 3 are equal")
+
+    @testset "dimension" begin
+        @test isa(manifold_dimension(M), Integer)
+        @test manifold_dimension(M) ≥ 0
+    end
+
+    @testset "representation" begin
+        for T ∈ (Manifolds.MPoint, Manifolds.TVector, Manifolds.CoTVector)
+            repr = Manifolds.representation_size(M, T)
+            @test isa(repr, Tuple)
+            for rs ∈ repr
+                @test rs > 0
+            end
+        end
+    end
+
+    @testset "injectivity radius" begin
+        @test injectivity_radius(M, pts[1]) > 0
+        for rm ∈ retraction_methods
+            @test injectivity_radius(M, pts[1], rm) > 0
+            @test injectivity_radius(M, pts[1], rm) ≤ injectivity_radius(M, pts[1])
+        end
+    end
+
+    tv1 = log(M, pts[1], pts[2])
+
+    @testset "is_manifold_point / is_tangent_vector" begin
+        for pt ∈ pts
+            @test is_manifold_point(M, pt)
+        end
+        @test is_tangent_vector(M, pts[1], tv1)
+    end
+
+    @testset "log/exp tests" begin
+        tv2 = log(M, pts[2], pts[1])
+        @test isapprox(M, pts[2], exp(M, pts[1], tv1))
+        @test isapprox(M, pts[1], exp(M, pts[1], tv1, 0))
+        @test isapprox(M, pts[2], exp(M, pts[1], tv1, 1))
+        @test isapprox(M, pts[1], exp(M, pts[2], tv2))
+        @test is_manifold_point(M, retract(M, pts[1], tv1))
+        @test isapprox(M, pts[1], retract(M, pts[1], tv1, 0))
+        for retr_method ∈ retraction_methods
+            @test is_manifold_point(M, retract(M, pts[1], tv1, retr_method))
+            @test isapprox(M, pts[1], retract(M, pts[1], tv1, 0, retr_method))
+        end
+        new_pt = exp(M, pts[1], tv1)
+        retract!(M, new_pt, pts[1], tv1)
+        @test is_manifold_point(M, new_pt)
+        for x ∈ pts
+            @test isapprox(M, zero_tangent_vector(M, x), log(M, x, x); atol = eps(eltype(x)) * exp_log_atol_multiplier)
+            @test isapprox(M, zero_tangent_vector(M, x), inverse_retract(M, x, x); atol = eps(eltype(x)) * exp_log_atol_multiplier)
+            for inv_retr_method ∈ inverse_retraction_methods
+                @test isapprox(M, zero_tangent_vector(M, x), inverse_retract(M, x, x, inv_retr_method); atol = eps(eltype(x)) * exp_log_atol_multiplier)
+            end
+        end
+        zero_tangent_vector!(M, tv1, pts[1])
+        @test isapprox(M, pts[1], tv1, zero_tangent_vector(M, pts[1]))
+        log!(M, tv1, pts[1], pts[2])
+        @test norm(M, pts[1], tv1) ≈ sqrt(inner(M, pts[1], tv1, tv1))
+
+        @test isapprox(M, exp(M, pts[1], tv1, 1), pts[2])
+        @test isapprox(M, exp(M, pts[1], tv1, 0), pts[1])
+
+        @test distance(M, pts[1], pts[2]) ≈ norm(M, pts[1], tv1)
+    end
+
+    @testset "basic linear algebra in tangent space" begin
+        @test isapprox(M, pts[1], 0*tv1, zero_tangent_vector(M, pts[1]))
+        @test isapprox(M, pts[1], 2*tv1, tv1+tv1)
+        @test isapprox(M, pts[1], 0*tv1, tv1-tv1)
+        @test isapprox(M, pts[1], (-1)*tv1, -tv1)
+    end
+
+    test_tangent_vector_broadcasting && @testset "broadcasted linear algebra in tangent space" begin
+        @test isapprox(M, pts[1], 3*tv1, 2 .* tv1 .+ tv1)
+        @test isapprox(M, pts[1], -tv1, tv1 .- 2 .* tv1)
+        @test isapprox(M, pts[1], -tv1, .-tv1)
+        v = similar(tv1)
+        v .= 2 .* tv1 .+ tv1
+        @test v ≈ 3*tv1
+    end
+
+    test_project_tangent && @testset "project_tangent test" begin
+        @test isapprox(M, pts[1], tv1, project_tangent(M, pts[1], tv1))
+        tv = similar(tv1)
+        project_tangent!(M, tv, pts[1], tv1)
+        @test isapprox(M, pts[1], tv, tv1)
+    end
+
+    test_forward_diff && @testset "ForwardDiff support" begin
+        exp_f(t) = distance(M, pts[1], exp(M, pts[1], t*tv1))
+        d12 = distance(M, pts[1], pts[2])
+        for t ∈ 0.1:0.1:0.9
+            @test d12 ≈ ForwardDiff.derivative(exp_f, t)
+        end
+
+        retract_f(t) = distance(M, pts[1], retract(M, pts[1], t*tv1))
+        for t ∈ 0.1:0.1:0.9
+            @test ForwardDiff.derivative(retract_f, t) ≥ 0
+        end
+    end
+
+    test_reverse_diff && @testset "ReverseDiff support" begin
+        exp_f(t) = distance(M, pts[1], exp(M, pts[1], t[1]*tv1))
+        d12 = distance(M, pts[1], pts[2])
+        for t ∈ 0.1:0.1:0.9
+            @test d12 ≈ ReverseDiff.gradient(exp_f, [t])[1]
+        end
+
+        retract_f(t) = distance(M, pts[1], retract(M, pts[1], t[1]*tv1))
+        for t ∈ 0.1:0.1:0.9
+            @test ReverseDiff.gradient(retract_f, [t])[1] ≥ 0
+        end
+    end
+
+    @testset "eltype" begin
+        tv1 = log(M, pts[1], pts[2])
+        @test eltype(tv1) == eltype(pts[1])
+        @test eltype(exp(M, pts[1], tv1)) == eltype(pts[1])
+    end
+
+    @testset "point distributions" begin
+        for pd ∈ point_distributions
+            for _ in 1:10
+                @test is_manifold_point(M, rand(pd))
+            end
+        end
+    end
+
+    @testset "tangent vector distributions" begin
+        for tvd ∈ tvector_distributions
+            supp = Manifolds.support(tvd)
+            for _ in 1:10
+                @test is_tangent_vector(M, supp.x, rand(tvd))
+            end
+        end
+    end
+end