diff --git a/NEWS.md b/NEWS.md
index 8b0714e5..3d7aed6c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,11 +5,18 @@ 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.9] unreleased
+## [0.15.9] 2024-05-02
 
 ### Added
 
 * Tests now also use `Aqua.jl` to spot problems in the code such as ambiguities.
+* introduce a `check_inverse_retraction` function to numerically check whether an inverse retraction method is a (correct) inverse retraction.
+* introduce a `check_retraction` function to numerically check whether a retraction method is a (correct) retraction.
+* introduce a `check_vector_transport` function to numerically check whether a vector transport is a (correct) vector transport.
+
+### Changed
+
+* introduced a `ManifoldsBaseTestUtils` module to encapsulate common types and function definitions in different parts of the tests.
 
 ## [0.15.8] 13/03/2024
 
diff --git a/Project.toml b/Project.toml
index 4dc08254..1c1fe3ca 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,44 +1,54 @@
 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.8"
+version = "0.15.9"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 
 [weakdeps]
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [extensions]
+ManifoldsBasePlotsExt = "Plots"
 ManifoldsBaseRecursiveArrayToolsExt = "RecursiveArrayTools"
+ManifoldsBaseStatisticsExt = "Statistics"
 
 [compat]
 Aqua = "0.8"
 DoubleFloats = ">= 0.9.2"
 ForwardDiff = "0.10"
-julia = "1.6"
 LinearAlgebra = "1.6"
 Markdown = "1.6"
 OrdinaryDiffEq = "6"
+Plots = "1"
+Printf = "1.6"
 Random = "1.6"
 RecursiveArrayTools = "2, 3"
 Requires = "1"
 ReverseDiff = "1"
 StaticArrays = "1"
+Statistics = "1.6"
 Test = "1.6"
+julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Aqua", "DoubleFloats", "ForwardDiff", "OrdinaryDiffEq", "ReverseDiff", "StaticArrays", "RecursiveArrayTools"]
+test = ["Test", "Aqua", "DoubleFloats", "ForwardDiff", "OrdinaryDiffEq", "Plots", "ReverseDiff", "StaticArrays", "Statistics", "RecursiveArrayTools"]
diff --git a/Readme.md b/Readme.md
index 8b3ccef5..3109d133 100644
--- a/Readme.md
+++ b/Readme.md
@@ -11,7 +11,7 @@
 [![codecov.io](http://codecov.io/github/JuliaManifolds/ManifoldsBase.jl/coverage.svg?branch=master)](https://codecov.io/gh/JuliaManifolds/ManifoldsBase.jl/)
 [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
-[![arXiv](https://img.shields.io/badge/arXiv%20CS.MS-2106.08777-blue.svg)](https://arxiv.org/abs/2106.08777)
+[![ACM TOMS](https://img.shields.io/badge/ACM%20TOMS-10.1145%2F3618296-blue.svg)](http://doi.org/10.1145/3618296)
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.5964340.svg)](https://doi.org/10.5281/zenodo.5964340)
 ## Installation
 
diff --git a/docs/make.jl b/docs/make.jl
index 736a6beb..fc2fecc3 100755
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -82,6 +82,7 @@ makedocs(;
         "Meta-Manifolds" => "metamanifolds.md",
         "Decorating/Extending a Manifold" => "decorator.md",
         "Bases for tangent spaces" => "bases.md",
+        "Numerical Verification" => "numerical_verification.md",
         "References" => "references.md",
     ],
     plugins = [bib],
diff --git a/docs/src/numerical_verification.md b/docs/src/numerical_verification.md
new file mode 100644
index 00000000..56817c4a
--- /dev/null
+++ b/docs/src/numerical_verification.md
@@ -0,0 +1,23 @@
+# Numerical Verification
+
+```@autodocs
+Modules = [ManifoldsBase]
+Pages = ["src/numerical_checks.jl"]
+Order = [:macro, :type, :function]
+Private = false
+Public = true
+```
+
+## Internal functions
+
+The following functions split the check into several parts, for
+example looking for the best fitting window and finding out the best slope,
+or plotting the slope.
+
+```@autodocs
+Modules = [ManifoldsBase]
+Pages = ["src/numerical_checks.jl"]
+Order = [:macro, :type, :function]
+Private = true
+Public = false
+```
\ No newline at end of file
diff --git a/docs/src/references.bib b/docs/src/references.bib
index c0e97eaf..169cef9f 100644
--- a/docs/src/references.bib
+++ b/docs/src/references.bib
@@ -25,6 +25,18 @@ @article{AxenBaranBergmannRzecki:2023
     VOLUME    = {49},
     YEAR      = {2023}
 }
+@book{Boumal:2023,
+    TITLE     = {An Introduction to Optimization on Smooth Manifolds},
+    AUTHOR    = {Boumal, Nicolas},
+    YEAR      = {2023},
+    MONTH     = mar,
+    EDITION   = {First},
+    PUBLISHER = {Cambridge University Press},
+    DOI       = {10.1017/9781009166164},
+    URL       = {https://www.nicolasboumal.net/book/index.html},
+    ABSTRACT  = {Optimization on Riemannian manifolds-the result of smooth geometry and optimization merging into one elegant modern framework-spans many areas of science and engineering, including machine learning, computer vision, signal processing, dynamical systems and scientific computing. This text introduces the differential geometry and Riemannian geometry concepts that will help students and researchers in applied mathematics, computer science and engineering gain a firm mathematical grounding to use these tools confidently in their research. Its charts-last approach will prove more intuitive from an optimizer's viewpoint, and all definitions and theorems are motivated to build time-tested optimization algorithms. Starting from first principles, the text goes on to cover current research on topics including worst-case complexity and geodesic convexity. Readers will appreciate the tricks of the trade for conducting research and for numerical implementations sprinkled throughout the book.},
+    ISBN      = {978-1-00-916616-4}
+}
 @article{EhlersPiraniSchild:1972,
     DOI       = {10.1007/s10714-012-1353-4},
     YEAR      = {1972},
diff --git a/ext/ManifoldsBasePlotsExt.jl b/ext/ManifoldsBasePlotsExt.jl
new file mode 100644
index 00000000..d9ca7406
--- /dev/null
+++ b/ext/ManifoldsBasePlotsExt.jl
@@ -0,0 +1,62 @@
+module ManifoldsBasePlotsExt
+
+if isdefined(Base, :get_extension)
+    using ManifoldsBase
+    using Plots
+    using Printf: @sprintf
+    import ManifoldsBase: plot_slope
+else
+    # imports need to be relative for Requires.jl-based workflows:
+    # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387
+    using ..ManifoldsBase
+    using ..Plots
+    using ..ManifoldsBase: @sprintf # is from Printf, but loaded in ManifoldsBase, and since Printf loading works here only on full moon days between 12 and noon, this trick might do it?
+    import ..ManifoldsBase: plot_slope
+end
+
+function ManifoldsBase.plot_slope(
+    x,
+    y;
+    slope = 2,
+    line_base = 0,
+    a = 0,
+    b = 2.0,
+    i = 1,
+    j = length(x),
+)
+    fig = plot(
+        x,
+        y;
+        xaxis = :log,
+        yaxis = :log,
+        label = "\$E(t)\$",
+        linewidth = 3,
+        legend = :topleft,
+        color = :lightblue,
+    )
+    s_line = [exp10(line_base + t * slope) for t in log10.(x)]
+    plot!(
+        fig,
+        x,
+        s_line;
+        label = "slope s=$slope",
+        linestyle = :dash,
+        color = :black,
+        linewidth = 2,
+    )
+    if (i != 0) && (j != 0)
+        best_line = [exp10(a + t * b) for t in log10.(x[i:j])]
+        plot!(
+            fig,
+            x[i:j],
+            best_line;
+            label = "best slope $(@sprintf("%.4f", b))",
+            color = :blue,
+            linestyle = :dot,
+            linewidth = 2,
+        )
+    end
+    return fig
+end
+
+end
diff --git a/ext/ManifoldsBaseStatisticsExt.jl b/ext/ManifoldsBaseStatisticsExt.jl
new file mode 100644
index 00000000..b5ba7172
--- /dev/null
+++ b/ext/ManifoldsBaseStatisticsExt.jl
@@ -0,0 +1,62 @@
+module ManifoldsBaseStatisticsExt
+
+if isdefined(Base, :get_extension)
+    using Statistics
+    using ManifoldsBase
+    import ManifoldsBase: find_best_slope_window
+else
+    # imports need to be relative for Requires.jl-based workflows:
+    # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387
+    using ..Statistics
+    using ..ManifoldsBase
+    import ..ManifoldsBase: find_best_slope_window
+end
+
+function ManifoldsBase.find_best_slope_window(
+    X,
+    Y,
+    window = nothing;
+    slope::Real = 2.0,
+    slope_tol::Real = 0.1,
+)
+    n = length(X)
+    if window !== nothing && (any(window .> n))
+        error(
+            "One of the window sizes ($(window)) is larger than the length of the signal (n=$n).",
+        )
+    end
+    a_best = 0
+    b_best = -Inf
+    i_best = 0
+    j_best = 0
+    r_best = 0 # longest interval
+    for w in (window === nothing ? (2:n) : [window...])
+        for j in 1:(n - w + 1)
+            x = X[j:(j + w - 1)]
+            y = Y[j:(j + w - 1)]
+            # fit a line a + bx
+            c = cor(x, y)
+            b = std(y) / std(x) * c
+            a = mean(y) - b * mean(x)
+            # look for the largest interval where b is within slope tolerance
+            r = (maximum(x) - minimum(x))
+            if (r > r_best) && abs(b - slope) < slope_tol #longer interval found.
+                r_best = r
+                a_best = a
+                b_best = b
+                i_best = j
+                j_best = j + w - 1 #last index (see x and y from before)
+            end
+            # not best interval - maybe it is still the (first) best slope?
+            if r_best == 0 && abs(b - slope) < abs(b_best - slope)
+                # but do not update `r` since this indicates only a best r
+                a_best = a
+                b_best = b
+                i_best = j
+                j_best = j + w - 1 #last index (see x and y from before)
+            end
+        end
+    end
+    return (a_best, b_best, i_best, j_best)
+end
+end
diff --git a/src/ManifoldsBase.jl b/src/ManifoldsBase.jl
index 10086ecc..09054298 100644
--- a/src/ManifoldsBase.jl
+++ b/src/ManifoldsBase.jl
@@ -22,6 +22,7 @@ import Random: rand, rand!
 
 using LinearAlgebra
 using Markdown: @doc_str
+using Printf: @sprintf
 using Random
 using Requires
 
@@ -1056,6 +1057,7 @@ include("vector_spaces.jl")
 include("point_vector_fallbacks.jl")
 include("nested_trait.jl")
 include("decorator_trait.jl")
+include("numerical_checks.jl")
 
 include("VectorFiber.jl")
 include("TangentSpace.jl")
@@ -1067,15 +1069,53 @@ include("PowerManifold.jl")
 
 #
 #
-# Requires
+# Init
 # -----
 function __init__()
+    #
+    # Error Hints
+    #
+    @static if isdefined(Base.Experimental, :register_error_hint) # COV_EXCL_LINE
+        Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs
+            if exc.f === plot_slope
+                print(
+                    io,
+                    """
+
+                    `plot_slope` has to be implemented using your favourite plotting package.
+                    A default is available when Plots.jl is added to the current environment.
+                    To then get the plotting functionality activated, do
+                    """,
+                )
+                printstyled(io, "`using Plots`"; color = :cyan)
+            end
+            if exc.f === find_best_slope_window
+                print(
+                    io,
+                    """
+
+                    `find_best_slope_window` has to be implemented using some statistics package
+                    A default is available when Statistics.jl is added to the current environment.
+                    To then get the functionality activated, do
+                    """,
+                )
+                printstyled(io, "`using Statistics`"; color = :cyan)
+            end
+        end
+    end
+    # Extensions in the pre 1.9 fallback using Requires.jl
     @static if !isdefined(Base, :get_extension)
         @require RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" begin
             include(
                 "../ext/ManifoldsBaseRecursiveArrayToolsExt/ManifoldsBaseRecursiveArrayToolsExt.jl",
             )
         end
+        @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
+            include("../ext/ManifoldsBasePlotsExt.jl")
+        end
+        @require Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" begin
+            include("../ext/ManifoldsBaseStatisticsExt.jl")
+        end
     end
 end
 #
@@ -1185,6 +1225,9 @@ export ×,
     change_metric!,
     change_representer,
     change_representer!,
+    check_inverse_retraction,
+    check_retraction,
+    check_vector_transport,
     copy,
     copyto!,
     default_approximation_method,
diff --git a/src/numerical_checks.jl b/src/numerical_checks.jl
new file mode 100644
index 00000000..c97dedad
--- /dev/null
+++ b/src/numerical_checks.jl
@@ -0,0 +1,404 @@
+
+@doc raw"""
+    check_inverse_retraction(
+        M::AbstractManifold,
+        inverse_rectraction_method::AbstractInverseRetractionMethod,
+        p=rand(M),
+        X=rand(M; vector_at=p);
+        #
+        exactness_tol::Real = 1e-12,
+        io::Union{IO,Nothing} = nothing,
+        limits::Tuple = (-8.0, 0.0),
+        log_range::AbstractVector = range(limits[1], limits[2]; length=N),
+        N::Int = 101,
+        name::String = "inverse retraction",
+        plot::Bool = false,
+        second_order::Bool = true
+        slope_tol::Real = 0.1,
+        error::Symbol = :none,
+        window = nothing,
+    )
+
+Check numerically wether the inverse retraction `inverse_retraction_method` is correct.
+This requires the [`exp`](@ref) and [`norm`](@ref) functions to be implemented for the [`AbstractManifold`](@ref) `M`.
+
+This implements a method similar to [Boumal:2023; Section 4.8 or Section 6.8](@cite).
+
+Note that if the errors are below the given tolerance and the method is exact,
+no plot is generated,
+
+# Keyword arguments
+
+* `exactness_tol`:     if all errors are below this tolerance, the inverse retraction is considered to be exact
+* `io`:                provide an `IO` to print the result to
+* `limits`:            specify the limits in the `log_range`, that is the exponent for the range
+* `log_range`:         specify the range of points (in log scale) to sample the length of the tangent vector `X`
+* `N`:                 number of points to verify within the `log_range` default range ``[10^{-8},10^{0}]``
+* `name`:              name to display in the plot
+* `plot`:              whether to plot the result (see [`plot_slope`](@ref))
+  The plot is in log-log-scale. This is returned and can then also be saved.
+* `second_order`:      check whether the retraction is of second order. if set to `false`, first order is checked.
+* `slope_tol`:         tolerance for the slope (global) of the approximation
+* `error`:             specify how to report errors: `:none`, `:info`, `:warn`, or `:error` are available
+* `window`:            specify window sizes within the `log_range` that are used for the slope estimation.
+  the default is, to use all window sizes `2:N`.
+"""
+function check_inverse_retraction(
+    M::AbstractManifold,
+    inverse_retraction_method::AbstractInverseRetractionMethod,
+    p = rand(M),
+    X = rand(M; vector_at = p);
+    exactness_tol::Real = 1e-12,
+    io::Union{IO,Nothing} = nothing,
+    limits = (-8.0, 0.0),
+    N::Int = 101,
+    second_order::Bool = true,
+    name::String = second_order ? "second order inverse retraction" : "inverse retraction",
+    log_range::AbstractVector = range(limits[1], limits[2]; length = N),
+    plot::Bool = false,
+    slope_tol::Real = 0.1,
+    error::Symbol = :none,
+    window = nothing,
+)
+    Xn = X ./ norm(M, p, X) # normalize tangent direction
+    # function for the directional derivative
+    #
+    T = exp10.(log_range)
+    # points `p_i` to evaluate the error function at
+    points = [exp(M, p, Xn, t) for t in T]
+    Xs = [t * Xn for t in T]
+    approx_Xs = [inverse_retract(M, p, q, inverse_retraction_method) for q in points]
+    errors = [norm(M, p, X - Y) for (X, Y) in zip(Xs, approx_Xs)]
+    return prepare_check_result(
+        log_range,
+        errors,
+        second_order ? 3.0 : 2.0;
+        exactness_tol = exactness_tol,
+        io = io,
+        name = name,
+        plot = plot,
+        slope_tol = slope_tol,
+        error = error,
+        window = window,
+    )
+end
+
+@doc raw"""
+    check_retraction(
+        M::AbstractManifold,
+        rectraction_method::AbstractRetractionMethod,
+        p=rand(M),
+        X=rand(M; vector_at=p);
+        #
+        exactness_tol::Real = 1e-12,
+        io::Union{IO,Nothing} = nothing,
+        limits::Tuple = (-8.0, 0.0),
+        log_range::AbstractVector = range(limits[1], limits[2]; length=N),
+        N::Int = 101,
+        name::String = "retraction",
+        plot::Bool = false,
+        second_order::Bool = true
+        slope_tol::Real = 0.1,
+        error::Symbol = :none,
+        window = nothing,
+    )
+
+Check numerically wether the retraction `vector_transport_to` is correct, by selecting
+a set of points ``q_i = \exp_p (t_i X)`` where ``t`` takes all values from `log_range`,
+to then compare [`parallel_transport_to`](@ref) to the `vector_transport_method`
+applied to the vector `Y`.
+
+This requires the [`exp`](@ref), [`parallel_transport_to`](@ref) and [`norm`](@ref) function
+to be implemented for the [`AbstractManifold`](@ref) `M`.
+
+This implements a method similar to [Boumal:2023; Section 4.8 or Section 6.8](@cite).
+
+Note that if the errors are below the given tolerance and the method is exact,
+no plot is generated,
+
+# Keyword arguments
+
+* `exactness_tol`:     if all errors are below this tolerance, the retraction is considered to be exact
+* `io`:                provide an `IO` to print the result to
+* `limits`:            specify the limits in the `log_range`, that is the exponent for the range
+* `log_range`:         specify the range of points (in log scale) to sample the length of the tangent vector `X`
+* `N`:                 number of points to verify within the `log_range` default range ``[10^{-8},10^{0}]``
+* `name`:              name to display in the plot
+* `plot`:              whether to plot the result (if `Plots.jl` is loaded).
+  The plot is in log-log-scale. This is returned and can then also be saved.
+* `second_order`:      check whether the retraction is of second order. if set to `false`, first order is checked.
+* `slope_tol`:         tolerance for the slope (global) of the approximation
+* `error`:             specify how to report errors: `:none`, `:info`, `:warn`, or `:error` are available
+* `window`:            specify window sizes within the `log_range` that are used for the slope estimation.
+  the default is, to use all window sizes `2:N`.
+"""
+function check_retraction(
+    M::AbstractManifold,
+    retraction_method::AbstractRetractionMethod,
+    p = rand(M),
+    X = rand(M; vector_at = p);
+    exactness_tol::Real = 1e-12,
+    io::Union{IO,Nothing} = nothing,
+    limits::Tuple = (-8.0, 0.0),
+    N::Int = 101,
+    second_order::Bool = true,
+    name::String = second_order ? "second order retraction" : "retraction",
+    log_range = range(limits[1], limits[2]; length = N),
+    plot::Bool = false,
+    slope_tol::Real = 0.1,
+    error::Symbol = :none,
+    window = nothing,
+)
+    Xn = X ./ norm(M, p, X) # normalize tangent direction
+    # function for the directional derivative
+    #
+    T = exp10.(log_range)
+    # points `p_i` to evaluate the error function at
+    points = [exp(M, p, Xn, t) for t in T]
+    approx_points = [retract(M, p, Xn, t, retraction_method) for t in T]
+    errors = [distance(M, p, q) for (p, q) in zip(points, approx_points)]
+    return prepare_check_result(
+        log_range,
+        errors,
+        second_order ? 3.0 : 2.0;
+        exactness_tol = exactness_tol,
+        io = io,
+        name = name,
+        plot = plot,
+        slope_tol = slope_tol,
+        error = error,
+        window = window,
+    )
+end
+
+@doc raw"""
+    check_vector_transport(
+        M::AbstractManifold,
+        vector_transport_method::AbstractVectorTransportMethod,
+        p=rand(M),
+        X=rand(M; vector_at=p),
+        Y=rand(M; vector_at=p);
+        #
+        exactness_tol::Real = 1e-12,
+        io::Union{IO,Nothing} = nothing,
+        limits::Tuple = (-8.0, 0.0),
+        log_range::AbstractVector = range(limits[1], limits[2]; length=N),
+        N::Int = 101,
+        name::String = "inverse retraction",
+        plot::Bool = false,
+        second_order::Bool = true
+        slope_tol::Real = 0.1,
+        error::Symbol = :none,
+        window = nothing,
+    )
+
+Check numerically wether the retraction `vector_transport_to` is correct, by selecting
+a set of points ``q_i = \exp_p (t_i X)`` where ``t`` takes all values from `log_range`,
+to then compare [`parallel_transport_to`](@ref) to the `vector_transport_method`
+applied to the vector `Y`.
+
+This requires the [`exp`](@ref), [`parallel_transport_to`](@ref) and [`norm`](@ref) function
+to be implemented for the [`AbstractManifold`](@ref) `M`.
+
+This implements a method similar to [Boumal:2023; Section 4.8 or Section 6.8](@cite).
+
+Note that if the errors are below the given tolerance and the method is exact,
+no plot is generated,
+
+# Keyword arguments
+
+* `exactness_tol`:     if all errors are below this tolerance, the differential is considered to be exact
+* `io`:                provide an `IO` to print the result to
+* `limits`:            specify the limits in the `log_range`, that is the exponent for the range
+* `log_range`:         specify the range of points (in log scale) to sample the differential line
+* `N`:                 number of points to verify within the `log_range` default range ``[10^{-8},10^{0}]``
+* `name`:              name to display in the plot
+* `plot`:              whether to plot the result (if `Plots.jl` is loaded).
+  The plot is in log-log-scale. This is returned and can then also be saved.
+* `second_order`:      check whether the retraction is of second order. if set to `false`, first order is checked.
+* `slope_tol`:         tolerance for the slope (global) of the approximation
+* `error`:             specify how to report errors: `:none`, `:info`, `:warn`, or `:error` are available
+* `window`:            specify window sizes within the `log_range` that are used for the slope estimation.
+  the default is, to use all window sizes `2:N`.
+"""
+function check_vector_transport(
+    M::AbstractManifold,
+    vector_transport_method::AbstractVectorTransportMethod,
+    p = rand(M),
+    X = rand(M; vector_at = p),
+    Y = rand(M; vector_at = p);
+    exactness_tol::Real = 1e-12,
+    io::Union{IO,Nothing} = nothing,
+    limits::Tuple = (-8.0, 0.0),
+    N::Int = 101,
+    second_order::Bool = true,
+    name::String = second_order ? "second order vector transport" : "vector transport",
+    log_range::AbstractVector = range(limits[1], limits[2]; length = N),
+    plot::Bool = false,
+    slope_tol::Real = 0.1,
+    error::Symbol = :none,
+    window = nothing,
+)
+    Xn = X ./ norm(M, p, X) # normalize tangent direction
+    # function for the directional derivative
+    #
+    T = exp10.(log_range)
+    # points `p_i` to evaluate the error function at
+    points = [exp(M, p, Xn, t) for t in T]
+    Yv = [vector_transport_to(M, p, Y, q, vector_transport_method) for q in points]
+    Yp = [parallel_transport_to(M, p, Y, q) for q in points]
+    errors = [norm(M, q, X - Y) for (q, X, Y) in zip(points, Yv, Yp)]
+    return prepare_check_result(
+        log_range,
+        errors,
+        second_order ? 3.0 : 2.0;
+        exactness_tol = exactness_tol,
+        io = io,
+        name = name,
+        plot = plot,
+        slope_tol = slope_tol,
+        error = error,
+        window = window,
+    )
+end
+
+function plot_slope end
+
+"""
+    plot_slope(x, y;
+        slope=2,
+        line_base=0,
+        a=0,
+        b=2.0,
+        i=1,
+        j=length(x)
+    )
+
+Plot the result from the verification functions on data `x,y` with two comparison lines
+
+1) `line_base` + t`slope`  as the global slope(s) the plot could have
+2) `a` + `b*t` on the interval [`x[i]`, `x[j]`] for some (best fitting) comparison slope
+
+!!! note
+    This function has to be implemented for a certain plotting package.
+    loading [Plots.jl](https://docs.juliaplots.org/stable/) provides a default implementation.
+"""
+plot_slope(x, y)
+
+"""
+    prepare_check_result(
+        log_range::AbstractVector,
+        errors::AbstractVector,
+        slope::Real;
+        exactness_to::Real = 1e3*eps(eltype(errors)),
+        io::Union{IO,Nothing} = nothing
+        name::String = "estimated slope",
+        plot::Bool = false,
+        slope_tol::Real = 0.1,
+        error::Symbol = :none,
+    )
+
+Given a range of values `log_range`, with computed `errors`,
+verify whether this yields a slope of `slope` in log-scale
+
+Note that if the errors are below the given tolerance and the method is exact,
+no plot is be generated,
+
+# Keyword arguments
+
+* `exactness_tol`: is all errors are below this tolerance, the verification is considered to be exact
+* `io`:            provide an `IO` to print the result to
+* `name`:          name to display in the plot title
+* `plot`:          whether to plot the result, see [`plot_slope`](@ref)
+  The plot is in log-log-scale. This is returned and can then also be saved.
+* `slope_tol`:     tolerance for the slope (global) of the approximation
+* `error`:         specify how to handle errors, `:none`, `:info`, `:warn`, `:error`
+"""
+
+function prepare_check_result(
+    log_range::AbstractVector,
+    errors::AbstractVector,
+    slope::Real;
+    io::Union{IO,Nothing} = nothing,
+    name::String = "estimated slope",
+    slope_tol::Real = 1e-1,
+    plot::Bool = false,
+    error::Symbol = :none,
+    window = nothing,
+    exactness_tol::Real = 1e3 * eps(eltype(errors)),
+)
+    if max(errors...) < exactness_tol
+        (io !== nothing) && print(
+            io,
+            "All errors are below the exactness tolerance $(exactness_tol). Your check can be considered exact, hence there is no use to check for a slope.\n",
+        )
+        return true
+    end
+    x = log_range[errors .> 0]
+    T = exp10.(x)
+    y = log10.(errors[errors .> 0])
+    (a, b) = find_best_slope_window(x, y, length(x))[1:2]
+    if isapprox(b, slope; atol = slope_tol)
+        plot && return plot_slope(
+            T,
+            errors[errors .> 0];
+            slope = slope,
+            line_base = errors[1],
+            a = a,
+            b = b,
+            i = 1,
+            j = length(y),
+        )
+        (io !== nothing) && print(
+            io,
+            "Your $name's slope is globally $(@sprintf("%.4f", b)), so within $slope ± $(slope_tol).\n",
+        )
+        return true
+    end
+    # otherwise
+    # find best contiguous window of length w
+    (ab, bb, ib, jb) = find_best_slope_window(x, y, window; slope_tol = slope_tol)
+    msg = "The $(name) fits best on [$(T[ib]),$(T[jb])] with slope  $(@sprintf("%.4f", bb)), but globally your slope $(@sprintf("%.4f", b)) is outside of the tolerance $slope ± $(slope_tol).\n"
+    (io !== nothing) && print(io, msg)
+    plot && return plot_slope(
+        T,
+        errors[errors .> 0];
+        slope = slope,
+        line_base = errors[1],
+        a = ab,
+        b = bb,
+        i = ib,
+        j = jb,
+    )
+    (error === :info) && @info msg
+    (error === :warn) && @warn msg
+    (error === :error) && throw(ErrorException(msg))
+    return false
+end
+
+function find_best_slope_window end
+
+"""
+    (a, b, i, j) = find_best_slope_window(X, Y, window=nothing; slope::Real=2.0, slope_tol::Real=0.1)
+
+Check data X,Y for the largest contiguous interval (window) with a regression line fitting “best”.
+Among all intervals with a slope within `slope_tol` to `slope` the longest one is taken.
+If no such interval exists, the one with the slope closest to `slope` is taken.
+
+If the window is set to `nothing` (default), all window sizes `2,...,length(X)` are checked.
+You can also specify a window size or an array of window sizes.
+
+For each window size, all its translates in the data is checked.
+For all these (shifted) windows the regression line is computed (with `a,b` in `a + t*b`)
+and the best line is computed.
+
+From the best line the following data is returned
+
+* `a`, `b` specifying the regression line `a + t*b`
+* `i`, `j` determining the window, i.e the regression line stems from data `X[i], ..., X[j]`
+
+!!! note
+    This function has to be implemented using some statistics package.
+    loading [Statistics.jl](https://github.com/JuliaStats/Statistics.jl) provides a default implementation.
+"""
+find_best_slope_window(X, Y, window = nothing)
diff --git a/src/point_vector_fallbacks.jl b/src/point_vector_fallbacks.jl
index 3e6f8382..ccb1b75b 100644
--- a/src/point_vector_fallbacks.jl
+++ b/src/point_vector_fallbacks.jl
@@ -113,11 +113,11 @@ macro default_manifold_fallbacks(TM, TP, TV, pfield::Symbol, vfield::Symbol)
         end
 
         function ManifoldsBase.angle(M::$TM, p::$TP, X::$TV, Y::$TV)
-            return angle(M, p.$pfield, X.$vfield, Y.$vfield)
+            return ManifoldsBase.angle(M, p.$pfield, X.$vfield, Y.$vfield)
         end
 
         function ManifoldsBase.check_point(M::$TM, p::$TP; kwargs...)
-            return check_point(M, p.$pfield; kwargs...)
+            return ManifoldsBase.check_point(M, p.$pfield; kwargs...)
         end
 
         function ManifoldsBase.check_vector(M::$TM, p::$TP, X::$TV; kwargs...)
@@ -132,26 +132,26 @@ macro default_manifold_fallbacks(TM, TP, TV, pfield::Symbol, vfield::Symbol)
         end
 
         function ManifoldsBase.distance(M::$TM, p::$TP, q::$TP)
-            return distance(M, p.$pfield, q.$pfield)
+            return ManifoldsBase.distance(M, p.$pfield, q.$pfield)
         end
 
         function ManifoldsBase.embed!(M::$TM, q::$TP, p::$TP)
-            embed!(M, q.$pfield, p.$pfield)
+            ManifoldsBase.embed!(M, q.$pfield, p.$pfield)
             return q
         end
 
         function ManifoldsBase.embed!(M::$TM, Y::$TV, p::$TP, X::$TV)
-            embed!(M, Y.$vfield, p.$pfield, X.$vfield)
+            ManifoldsBase.embed!(M, Y.$vfield, p.$pfield, X.$vfield)
             return Y
         end
 
         function ManifoldsBase.exp!(M::$TM, q::$TP, p::$TP, X::$TV)
-            exp!(M, q.$pfield, p.$pfield, X.$vfield)
+            ManifoldsBase.exp!(M, q.$pfield, p.$pfield, X.$vfield)
             return q
         end
 
         function ManifoldsBase.inner(M::$TM, p::$TP, X::$TV, Y::$TV)
-            return inner(M, p.$pfield, X.$vfield, Y.$vfield)
+            return ManifoldsBase.inner(M, p.$pfield, X.$vfield, Y.$vfield)
         end
 
         function ManifoldsBase.inverse_retract!(
@@ -161,25 +161,25 @@ macro default_manifold_fallbacks(TM, TP, TV, pfield::Symbol, vfield::Symbol)
             q::$TP,
             m::LogarithmicInverseRetraction,
         )
-            inverse_retract!(M, X.$vfield, p.$pfield, q.$pfield, m)
+            ManifoldsBase.inverse_retract!(M, X.$vfield, p.$pfield, q.$pfield, m)
             return X
         end
 
         function ManifoldsBase.isapprox(M::$TM, p::$TP, q::$TP; kwargs...)
-            return isapprox(M, p.$pfield, q.$pfield; kwargs...)
+            return ManifoldsBase.isapprox(M, p.$pfield, q.$pfield; kwargs...)
         end
 
         function ManifoldsBase.isapprox(M::$TM, p::$TP, X::$TV, Y::$TV; kwargs...)
-            return isapprox(M, p.$pfield, X.$vfield, Y.$vfield; kwargs...)
+            return ManifoldsBase.isapprox(M, p.$pfield, X.$vfield, Y.$vfield; kwargs...)
         end
 
         function ManifoldsBase.log!(M::$TM, X::$TV, p::$TP, q::$TP)
-            log!(M, X.$vfield, p.$pfield, q.$pfield)
+            ManifoldsBase.log!(M, X.$vfield, p.$pfield, q.$pfield)
             return X
         end
 
         function ManifoldsBase.norm(M::$TM, p::$TP, X::$TV)
-            return norm(M, p.$pfield, X.$vfield)
+            return ManifoldsBase.norm(M, p.$pfield, X.$vfield)
         end
 
         function ManifoldsBase.retract!(
@@ -189,7 +189,7 @@ macro default_manifold_fallbacks(TM, TP, TV, pfield::Symbol, vfield::Symbol)
             X::$TV,
             m::ExponentialRetraction,
         )
-            retract!(M, q.$pfield, p.$pfield, X.$vfield, m)
+            ManifoldsBase.retract!(M, q.$pfield, p.$pfield, X.$vfield, m)
             return q
         end
 
@@ -200,7 +200,7 @@ macro default_manifold_fallbacks(TM, TP, TV, pfield::Symbol, vfield::Symbol)
             X::$TV,
             c::AbstractVector,
         )
-            vector_transport_along!(M, Y.$vfield, p.$pfield, X.$vfield, c)
+            ManifoldsBase.vector_transport_along!(M, Y.$vfield, p.$pfield, X.$vfield, c)
             return Y
         end
 
@@ -209,7 +209,7 @@ macro default_manifold_fallbacks(TM, TP, TV, pfield::Symbol, vfield::Symbol)
         end
 
         function ManifoldsBase.zero_vector!(M::$TM, X::$TV, p::$TP)
-            zero_vector!(M, X.$vfield, p.$pfield)
+            ManifoldsBase.zero_vector!(M, X.$vfield, p.$pfield)
             return X
         end
     end
diff --git a/test/ManifoldsBaseTestUtils.jl b/test/ManifoldsBaseTestUtils.jl
new file mode 100644
index 00000000..fcdc9130
--- /dev/null
+++ b/test/ManifoldsBaseTestUtils.jl
@@ -0,0 +1,589 @@
+"""
+    ManifoldsBaseTestUtils
+
+A small module to collect common definitions and functions used in (several) tests of
+`ManifoldsBase.jl`.
+
+* A `TestSphere`
+*
+"""
+module ManifoldsBaseTestUtils
+
+using ManifoldsBase, LinearAlgebra, Random
+using ManifoldsBase: ℝ, ℂ, AbstractManifold, DefaultManifold, EuclideanMetric
+using ManifoldsBase: AbstractNumbers, RealNumbers, ComplexNumbers
+using ManifoldsBase:
+    @manifold_element_forwards, @manifold_vector_forwards, @default_manifold_fallbacks
+import Base: +, *, -
+#
+#
+# minimal implementation of the sphere – to test a few more involved Riemannian functions
+struct TestSphere{N,𝔽} <: AbstractManifold{𝔽} end
+TestSphere(N::Int, 𝔽 = ℝ) = TestSphere{N,𝔽}()
+
+function ManifoldsBase.change_metric!(
+    M::TestSphere,
+    Y,
+    ::ManifoldsBase.EuclideanMetric,
+    p,
+    X,
+)
+    return copyto!(M, Y, p, X)
+end
+function ManifoldsBase.change_representer!(
+    M::TestSphere,
+    Y,
+    ::ManifoldsBase.EuclideanMetric,
+    p,
+    X,
+)
+    return copyto!(M, Y, p, X)
+end
+function ManifoldsBase.check_point(M::TestSphere, p; kwargs...)
+    if !isapprox(norm(p), 1.0; kwargs...)
+        return DomainError(
+            norm(p),
+            "The point $(p) does not lie on the $(M) since its norm is not 1.",
+        )
+    end
+    return nothing
+end
+function ManifoldsBase.check_vector(M::TestSphere, p, X; kwargs...)
+    if !isapprox(abs(real(dot(p, X))), 0.0; kwargs...)
+        return DomainError(
+            abs(dot(p, X)),
+            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal in the embedding.",
+        )
+    end
+    return nothing
+end
+function ManifoldsBase.exp!(M::TestSphere, q, p, X)
+    return exp!(M, q, p, X, one(number_eltype(X)))
+end
+function ManifoldsBase.exp!(::TestSphere, q, p, X, t::Number)
+    θ = abs(t) * norm(X)
+    if θ == 0
+        copyto!(q, p)
+    else
+        X_scale = t * sin(θ) / θ
+        q .= p .* cos(θ) .+ X .* X_scale
+    end
+    return q
+end
+function ManifoldsBase.get_basis_diagonalizing(
+    M::TestSphere{n},
+    p,
+    B::DiagonalizingOrthonormalBasis{ℝ},
+) where {n}
+    A = zeros(n + 1, n + 1)
+    A[1, :] = transpose(p)
+    A[2, :] = transpose(B.frame_direction)
+    V = nullspace(A)
+    κ = ones(n)
+    if !iszero(B.frame_direction)
+        # if we have a nonzero direction for the geodesic, add it and it gets curvature zero from the tensor
+        V = hcat(B.frame_direction / norm(M, p, B.frame_direction), V)
+        κ[1] = 0 # no curvature along the geodesic direction, if x!=y
+    end
+    T = typeof(similar(B.frame_direction))
+    Ξ = [convert(T, V[:, i]) for i in 1:n]
+    return CachedBasis(B, κ, Ξ)
+end
+@inline function nzsign(z, absz = abs(z))
+    psignz = z / absz
+    return ifelse(iszero(absz), one(psignz), psignz)
+end
+function ManifoldsBase.get_coordinates_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers)
+    n = manifold_dimension(M)
+    p1 = p[1]
+    cosθ = abs(p1)
+    λ = nzsign(p1, cosθ)
+    pend, Xend = view(p, 2:(n + 1)), view(X, 2:(n + 1))
+    factor = λ * X[1] / (1 + cosθ)
+    Y .= Xend .- pend .* factor
+    return Y
+end
+function ManifoldsBase.get_vector_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers)
+    n = manifold_dimension(M)
+    p1 = p[1]
+    cosθ = abs(p1)
+    λ = nzsign(p1, cosθ)
+    pend = view(p, 2:(n + 1))
+    pX = dot(pend, X)
+    factor = pX / (1 + cosθ)
+    Y[1] = -λ * pX
+    Y[2:(n + 1)] .= X .- pend .* factor
+    return Y
+end
+ManifoldsBase.injectivity_radius(::TestSphere) = π
+ManifoldsBase.inner(::TestSphere, p, X, Y) = dot(X, Y)
+function ManifoldsBase.inverse_retract_project!(M::TestSphere, X, p, q)
+    X .= q .- p
+    project!(M, X, p, X)
+    return X
+end
+ManifoldsBase.is_flat(M::TestSphere) = manifold_dimension(M) == 1
+function ManifoldsBase.log!(::TestSphere, X, p, q)
+    cosθ = clamp(real(dot(p, q)), -1, 1)
+    X .= q .- p .* cosθ
+    nrmX = norm(X)
+    if nrmX > 0
+        X .*= acos(cosθ) / nrmX
+    end
+    return X
+end
+ManifoldsBase.manifold_dimension(::TestSphere{N}) where {N} = N
+function ManifoldsBase.parallel_transport_to!(::TestSphere, Y, p, X, q)
+    m = p .+ q
+    mnorm2 = real(dot(m, m))
+    factor = 2 * real(dot(X, q)) / mnorm2
+    Y .= X .- m .* factor
+    return Y
+end
+function ManifoldsBase.project!(::TestSphere, q, p)
+    q .= p ./ norm(p)
+    return q
+end
+function ManifoldsBase.project!(::TestSphere, Y, p, X)
+    Y .= X .- p .* real(dot(p, X))
+    return Y
+end
+function Random.rand!(M::TestSphere, pX; vector_at = nothing, σ = one(eltype(pX)))
+    return rand!(Random.default_rng(), M, pX; vector_at = vector_at, σ = σ)
+end
+function Random.rand!(
+    rng::AbstractRNG,
+    M::TestSphere,
+    pX;
+    vector_at = nothing,
+    σ = one(eltype(pX)),
+)
+    if vector_at === nothing
+        project!(M, pX, randn(rng, eltype(pX), representation_size(M)))
+    else
+        n = σ * randn(rng, eltype(pX), size(pX)) # Gaussian in embedding
+        project!(M, pX, vector_at, n) #project to TpM (keeps Gaussianness)
+    end
+    return pX
+end
+ManifoldsBase.representation_size(::TestSphere{N}) where {N} = (N + 1,)
+function ManifoldsBase.retract_project!(M::TestSphere, q, p, X, t::Number)
+    q .= p .+ t .* X
+    project!(M, q, q)
+    return q
+end
+function ManifoldsBase.riemann_tensor!(M::TestSphere, Xresult, p, X, Y, Z)
+    innerZX = inner(M, p, Z, X)
+    innerZY = inner(M, p, Z, Y)
+    Xresult .= innerZY .* X .- innerZX .* Y
+    return Xresult
+end
+function ManifoldsBase.sectional_curvature_max(M::TestSphere)
+    return ifelse(manifold_dimension(M) == 1, 0.0, 1.0)
+end
+function ManifoldsBase.sectional_curvature_min(M::TestSphere)
+    return ifelse(manifold_dimension(M) == 1, 0.0, 1.0)
+end
+
+# from Absil, Mahony, Trumpf, 2013 https://sites.uclouvain.be/absil/2013-01/Weingarten_07PA_techrep.pdf
+function ManifoldsBase.Weingarten!(::TestSphere, Y, p, X, V)
+    return Y .= -X * p' * V
+end
+
+#
+#
+# minimal implementation of SPD (for its nontrivial metric conversion)
+# ---
+struct TestSPD <: AbstractManifold{ℝ}
+    n::Int
+end
+function ManifoldsBase.change_metric!(::TestSPD, Y, ::EuclideanMetric, p, X)
+    Y .= p * X
+    return Y
+end
+function ManifoldsBase.change_representer!(::TestSPD, Y, ::EuclideanMetric, p, X)
+    Y .= p * X * p
+    return Y
+end
+function ManifoldsBase.inner(::TestSPD, p, X, Y)
+    F = cholesky(Symmetric(convert(AbstractMatrix, p)))
+    return dot((F \ Symmetric(X)), (Symmetric(Y) / F))
+end
+function spd_sqrt_and_sqrt_inv(p::AbstractMatrix)
+    e = eigen(Symmetric(p))
+    U = e.vectors
+    S = max.(e.values, floatmin(eltype(e.values)))
+    Ssqrt = Diagonal(sqrt.(S))
+    SsqrtInv = Diagonal(1 ./ sqrt.(S))
+    return (Symmetric(U * Ssqrt * transpose(U)), Symmetric(U * SsqrtInv * transpose(U)))
+end
+function ManifoldsBase.log!(::TestSPD, X, p, q)
+    (p_sqrt, p_sqrt_inv) = spd_sqrt_and_sqrt_inv(p)
+    T = Symmetric(p_sqrt_inv * convert(AbstractMatrix, q) * p_sqrt_inv)
+    e2 = eigen(T)
+    Se = Diagonal(log.(max.(e2.values, eps())))
+    pUe = p_sqrt * e2.vectors
+    return mul!(X, pUe, Se * transpose(pUe))
+end
+ManifoldsBase.representation_size(M::TestSPD) = (M.n, M.n)
+
+#
+#
+# A simple Manifold based on a projection onto a subspace
+struct ProjManifold <: AbstractManifold{ℝ} end
+ManifoldsBase.inner(::ProjManifold, p, X, Y) = dot(X, Y)
+ManifoldsBase.project!(::ProjManifold, Y, p, X) = (Y .= X .- dot(p, X) .* p)
+ManifoldsBase.representation_size(::ProjManifold) = (2, 3)
+ManifoldsBase.manifold_dimension(::ProjManifold) = 5
+ManifoldsBase.get_vector_orthonormal(::ProjManifold, p, X, N) = reverse(X)
+#
+#
+# A second simple Manifold based on a projection onto a subspace
+struct ProjectionTestManifold <: AbstractManifold{ℝ} end
+ManifoldsBase.inner(::ProjectionTestManifold, ::Any, X, Y) = dot(X, Y)
+function ManifoldsBase.project!(::ProjectionTestManifold, Y, p, X)
+    Y .= X .- dot(p, X) .* p
+    Y[end] = 0
+    return Y
+end
+ManifoldsBase.manifold_dimension(::ProjectionTestManifold) = 100
+
+#
+# Thre Non-Things to check for the correct errors in case functions are not implemented
+struct NonManifold <: AbstractManifold{ℝ} end
+struct NonBasis <: ManifoldsBase.AbstractBasis{ℝ,TangentSpaceType} end
+struct NonMPoint <: AbstractManifoldPoint end
+struct NonTVector <: TVector end
+struct NonCoTVector <: CoTVector end
+struct NotImplementedRetraction <: AbstractRetractionMethod end
+struct NotImplementedInverseRetraction <: AbstractInverseRetractionMethod end
+*(t::Float64, X::NonTVector) = X
+
+struct NonBroadcastBasisThing{T}
+    v::T
+end
++(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing) = NonBroadcastBasisThing(a.v + b.v)
+*(α, a::NonBroadcastBasisThing) = NonBroadcastBasisThing(α * a.v)
+-(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing) = NonBroadcastBasisThing(a.v - b.v)
+
+function ManifoldsBase.isapprox(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing)
+    return isapprox(a.v, b.v)
+end
+
+function ManifoldsBase.number_eltype(a::NonBroadcastBasisThing)
+    return typeof(reduce(+, one(number_eltype(eti)) for eti in a.v))
+end
+
+ManifoldsBase.allocate(a::NonBroadcastBasisThing) = NonBroadcastBasisThing(allocate(a.v))
+function ManifoldsBase.allocate(a::NonBroadcastBasisThing, ::Type{T}) where {T}
+    return NonBroadcastBasisThing(allocate(a.v, T))
+end
+function ManifoldsBase.allocate(::NonBroadcastBasisThing, ::Type{T}, s::Integer) where {T}
+    return Vector{T}(undef, s)
+end
+
+function Base.copyto!(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing)
+    copyto!(a.v, b.v)
+    return a
+end
+
+function ManifoldsBase.log!(
+    ::DefaultManifold,
+    v::NonBroadcastBasisThing,
+    x::NonBroadcastBasisThing,
+    y::NonBroadcastBasisThing,
+)
+    return copyto!(v, y - x)
+end
+
+function ManifoldsBase.exp!(
+    ::DefaultManifold,
+    y::NonBroadcastBasisThing,
+    x::NonBroadcastBasisThing,
+    v::NonBroadcastBasisThing,
+)
+    return copyto!(y, x + v)
+end
+
+function ManifoldsBase.get_basis_orthonormal(
+    ::DefaultManifold,
+    p::NonBroadcastBasisThing,
+    𝔽::RealNumbers,
+)
+    return CachedBasis(
+        DefaultOrthonormalBasis(𝔽),
+        [
+            NonBroadcastBasisThing(ManifoldsBase._euclidean_basis_vector(p.v, i)) for
+            i in eachindex(p.v)
+        ],
+    )
+end
+function ManifoldsBase.get_basis_orthogonal(
+    ::DefaultManifold,
+    p::NonBroadcastBasisThing,
+    𝔽::RealNumbers,
+)
+    return CachedBasis(
+        DefaultOrthogonalBasis(𝔽),
+        [
+            NonBroadcastBasisThing(ManifoldsBase._euclidean_basis_vector(p.v, i)) for
+            i in eachindex(p.v)
+        ],
+    )
+end
+function ManifoldsBase.get_basis_default(
+    ::DefaultManifold,
+    p::NonBroadcastBasisThing,
+    N::ManifoldsBase.RealNumbers,
+)
+    return CachedBasis(
+        DefaultBasis(N),
+        [
+            NonBroadcastBasisThing(ManifoldsBase._euclidean_basis_vector(p.v, i)) for
+            i in eachindex(p.v)
+        ],
+    )
+end
+
+function ManifoldsBase.get_coordinates_orthonormal!(
+    M::DefaultManifold,
+    Y,
+    ::NonBroadcastBasisThing,
+    X::NonBroadcastBasisThing,
+    ::RealNumbers,
+)
+    copyto!(Y, reshape(X.v, manifold_dimension(M)))
+    return Y
+end
+
+function ManifoldsBase.get_vector_orthonormal!(
+    M::DefaultManifold,
+    Y::NonBroadcastBasisThing,
+    ::NonBroadcastBasisThing,
+    X,
+    ::RealNumbers,
+)
+    copyto!(Y.v, reshape(X, representation_size(M)))
+    return Y
+end
+
+function ManifoldsBase.inner(
+    ::DefaultManifold,
+    ::NonBroadcastBasisThing,
+    X::NonBroadcastBasisThing,
+    Y::NonBroadcastBasisThing,
+)
+    return dot(X.v, Y.v)
+end
+
+ManifoldsBase._get_vector_cache_broadcast(::NonBroadcastBasisThing) = Val(false)
+
+DiagonalizingBasisProxy() = DiagonalizingOrthonormalBasis([1.0, 0.0, 0.0])
+
+#
+#
+# Vector Space types
+struct TestVectorSpaceType <: ManifoldsBase.VectorSpaceType end
+
+struct TestFiberType <: ManifoldsBase.FiberType end
+
+function ManifoldsBase.fiber_dimension(M::AbstractManifold, ::TestFiberType)
+    return 2 * manifold_dimension(M)
+end
+
+function ManifoldsBase.vector_space_dimension(M::AbstractManifold, ::TestVectorSpaceType)
+    return 2 * manifold_dimension(M)
+end
+
+#
+#
+# DefaultManifold with a few artificiual retrations
+
+struct CustomDefinedRetraction <: ManifoldsBase.AbstractRetractionMethod end
+struct CustomUndefinedRetraction <: ManifoldsBase.AbstractRetractionMethod end
+struct CustomDefinedKeywordRetraction <: ManifoldsBase.AbstractRetractionMethod end
+struct CustomDefinedKeywordInverseRetraction <:
+       ManifoldsBase.AbstractInverseRetractionMethod end
+struct CustomDefinedInverseRetraction <: ManifoldsBase.AbstractInverseRetractionMethod end
+
+struct DefaultPoint{T} <: AbstractManifoldPoint
+    value::T
+end
+Base.convert(::Type{DefaultPoint{T}}, v::T) where {T} = DefaultPoint(v)
+Base.eltype(p::DefaultPoint) = eltype(p.value)
+Base.length(p::DefaultPoint) = length(p.value)
+struct DefaultTVector{T} <: TVector
+    value::T
+end
+Base.convert(::Type{DefaultTVector{T}}, ::DefaultPoint, v::T) where {T} = DefaultTVector(v)
+Base.eltype(X::DefaultTVector) = eltype(X.value)
+function Base.fill!(X::DefaultTVector, x)
+    fill!(X.value, x)
+    return X
+end
+function ManifoldsBase.allocate_result_type(
+    ::DefaultManifold,
+    ::typeof(log),
+    ::Tuple{DefaultPoint,DefaultPoint},
+)
+    return DefaultTVector
+end
+function ManifoldsBase.allocate_result_type(
+    ::DefaultManifold,
+    ::typeof(inverse_retract),
+    ::Tuple{DefaultPoint,DefaultPoint},
+)
+    return DefaultTVector
+end
+
+ManifoldsBase.@manifold_element_forwards DefaultPoint value
+ManifoldsBase.@manifold_vector_forwards DefaultTVector value
+ManifoldsBase.@default_manifold_fallbacks ManifoldsBase.DefaultManifold DefaultPoint DefaultTVector value value
+
+function ManifoldsBase._injectivity_radius(::DefaultManifold, ::CustomDefinedRetraction)
+    return 10.0
+end
+function ManifoldsBase._retract!(
+    M::DefaultManifold,
+    q,
+    p,
+    X,
+    t::Number,
+    ::CustomDefinedKeywordRetraction;
+    kwargs...,
+)
+    return retract_custom_kw!(M, q, p, X, t; kwargs...)
+end
+function retract_custom_kw!(
+    ::DefaultManifold,
+    q::DefaultPoint,
+    p::DefaultPoint,
+    X::DefaultTVector,
+    t::Number;
+    scale = 2.0,
+)
+    q.value .= scale .* p.value .+ t .* X.value
+    return q
+end
+function ManifoldsBase._inverse_retract!(
+    M::DefaultManifold,
+    X,
+    p,
+    q,
+    ::CustomDefinedKeywordInverseRetraction;
+    kwargs...,
+)
+    return inverse_retract_custom_kw!(M, X, p, q; kwargs...)
+end
+function inverse_retract_custom_kw!(
+    ::DefaultManifold,
+    X::DefaultTVector,
+    p::DefaultPoint,
+    q::DefaultPoint;
+    scale = 2.0,
+)
+    X.value .= q.value - scale * p.value
+    return X
+end
+
+function ManifoldsBase.retract!(
+    ::DefaultManifold,
+    q::DefaultPoint,
+    p::DefaultPoint,
+    X::DefaultTVector,
+    t::Number,
+    ::CustomDefinedRetraction,
+)
+    q.value .= 2 .* p.value .+ t * X.value
+    return q
+end
+
+function ManifoldsBase.inverse_retract!(
+    ::DefaultManifold,
+    X::DefaultTVector,
+    p::DefaultPoint,
+    q::DefaultPoint,
+    ::CustomDefinedInverseRetraction,
+)
+    X.value .= q.value .- 2 .* p.value
+    return X
+end
+
+struct MatrixVectorTransport{T} <: AbstractVector{T}
+    m::Matrix{T}
+end
+# dummy retractions, inverse retracions for fallback tests - mutating should be enough
+ManifoldsBase.retract_polar!(::DefaultManifold, q, p, X, t::Number) = (q .= p .+ t .* X)
+ManifoldsBase.retract_project!(::DefaultManifold, q, p, X, t::Number) = (q .= p .+ t .* X)
+ManifoldsBase.retract_qr!(::DefaultManifold, q, p, X, t::Number) = (q .= p .+ t .* X)
+function ManifoldsBase.retract_exp_ode!(
+    ::DefaultManifold,
+    q,
+    p,
+    X,
+    t::Number,
+    m::AbstractRetractionMethod,
+    B::ManifoldsBase.AbstractBasis,
+)
+    return (q .= p .+ t .* X)
+end
+
+function ManifoldsBase.retract_pade!(
+    ::DefaultManifold,
+    q,
+    p,
+    X,
+    t::Number,
+    m::PadeRetraction,
+)
+    return (q .= p .+ t .* X)
+end
+function ManifoldsBase.retract_sasaki!(
+    ::DefaultManifold,
+    q,
+    p,
+    X,
+    t::Number,
+    ::SasakiRetraction,
+)
+    return (q .= p .+ t .* X)
+end
+ManifoldsBase.retract_softmax!(::DefaultManifold, q, p, X, t::Number) = (q .= p .+ t .* X)
+ManifoldsBase.get_embedding(M::DefaultManifold) = M # dummy embedding
+ManifoldsBase.inverse_retract_polar!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
+ManifoldsBase.inverse_retract_project!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
+ManifoldsBase.inverse_retract_qr!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
+ManifoldsBase.inverse_retract_softmax!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
+function ManifoldsBase.inverse_retract_nlsolve!(
+    ::DefaultManifold,
+    Y,
+    p,
+    q,
+    m::NLSolveInverseRetraction,
+)
+    return (Y .= q .- p)
+end
+function ManifoldsBase.vector_transport_along_project!(
+    ::DefaultManifold,
+    Y,
+    p,
+    X,
+    c::AbstractVector,
+)
+    return (Y .= X)
+end
+Base.getindex(x::MatrixVectorTransport, i) = x.m[:, i]
+Base.size(x::MatrixVectorTransport) = (size(x.m, 2),)
+
+export CustomDefinedInverseRetraction, CustomDefinedKeywordInverseRetraction
+export CustomDefinedKeywordRetraction, CustomDefinedRetraction, CustomUndefinedRetraction
+export DefaultPoint, DefaultTVector
+export DiagonalizingBasisProxy
+export MatrixVectorTransport
+export NonManifold, NonBasis, NonBroadcastBasisThing
+export NonMPoint, NonTVector, NonCoTVector
+export NotImplementedRetraction, NotImplementedInverseRetraction
+export ProjManifold, ProjectionTestManifold
+export TestSphere, TestSPD
+export TestVectorSpaceType, TestFiberType
+end
diff --git a/test/bases.jl b/test/bases.jl
index efbaa935..49e2de2b 100644
--- a/test/bases.jl
+++ b/test/bases.jl
@@ -4,153 +4,10 @@ using ManifoldsBase: DefaultManifold, ℝ, ℂ, RealNumbers, ComplexNumbers
 using ManifoldsBase: CotangentSpaceType, TangentSpaceType
 using ManifoldsBase: FVector
 using Test
-import Base: +, -, *, copyto!, isapprox
-import ManifoldsBase:
-    allocate,
-    get_vector_orthonormal!,
-    get_coordinates_orthonormal!,
-    get_basis_orthogonal,
-    get_basis_orthonormal
-
-
-struct ProjManifold <: AbstractManifold{ℝ} end
-
-ManifoldsBase.inner(::ProjManifold, p, X, Y) = dot(X, Y)
-ManifoldsBase.project!(::ProjManifold, Y, p, X) = (Y .= X .- dot(p, X) .* p)
-ManifoldsBase.representation_size(::ProjManifold) = (2, 3)
-ManifoldsBase.manifold_dimension(::ProjManifold) = 5
-ManifoldsBase.get_vector_orthonormal(::ProjManifold, p, X, N) = reverse(X)
-
-struct ProjectionTestManifold <: AbstractManifold{ℝ} end
-
-ManifoldsBase.inner(::ProjectionTestManifold, ::Any, X, Y) = dot(X, Y)
-function ManifoldsBase.project!(::ProjectionTestManifold, Y, p, X)
-    Y .= X .- dot(p, X) .* p
-    Y[end] = 0
-    return Y
-end
-ManifoldsBase.manifold_dimension(::ProjectionTestManifold) = 100
-
-struct NonManifold <: AbstractManifold{ℝ} end
-struct NonBasis <: ManifoldsBase.AbstractBasis{ℝ,TangentSpaceType} end
-
-struct NonBroadcastBasisThing{T}
-    v::T
-end
-
-+(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing) = NonBroadcastBasisThing(a.v + b.v)
-*(α, a::NonBroadcastBasisThing) = NonBroadcastBasisThing(α * a.v)
--(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing) = NonBroadcastBasisThing(a.v - b.v)
-
-isapprox(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing) = isapprox(a.v, b.v)
-
-function ManifoldsBase.number_eltype(a::NonBroadcastBasisThing)
-    return typeof(reduce(+, one(number_eltype(eti)) for eti in a.v))
-end
-
-allocate(a::NonBroadcastBasisThing) = NonBroadcastBasisThing(allocate(a.v))
-function allocate(a::NonBroadcastBasisThing, ::Type{T}) where {T}
-    return NonBroadcastBasisThing(allocate(a.v, T))
-end
-allocate(::NonBroadcastBasisThing, ::Type{T}, s::Integer) where {T} = Vector{T}(undef, s)
-
-function copyto!(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing)
-    copyto!(a.v, b.v)
-    return a
-end
-
-function ManifoldsBase.log!(
-    ::DefaultManifold,
-    v::NonBroadcastBasisThing,
-    x::NonBroadcastBasisThing,
-    y::NonBroadcastBasisThing,
-)
-    return copyto!(v, y - x)
-end
-
-function ManifoldsBase.exp!(
-    ::DefaultManifold,
-    y::NonBroadcastBasisThing,
-    x::NonBroadcastBasisThing,
-    v::NonBroadcastBasisThing,
-)
-    return copyto!(y, x + v)
-end
-
-function ManifoldsBase.get_basis_orthonormal(
-    ::DefaultManifold,
-    p::NonBroadcastBasisThing,
-    𝔽::RealNumbers,
-)
-    return CachedBasis(
-        DefaultOrthonormalBasis(𝔽),
-        [
-            NonBroadcastBasisThing(ManifoldsBase._euclidean_basis_vector(p.v, i)) for
-            i in eachindex(p.v)
-        ],
-    )
-end
-function ManifoldsBase.get_basis_orthogonal(
-    ::DefaultManifold,
-    p::NonBroadcastBasisThing,
-    𝔽::RealNumbers,
-)
-    return CachedBasis(
-        DefaultOrthogonalBasis(𝔽),
-        [
-            NonBroadcastBasisThing(ManifoldsBase._euclidean_basis_vector(p.v, i)) for
-            i in eachindex(p.v)
-        ],
-    )
-end
-function ManifoldsBase.get_basis_default(
-    ::DefaultManifold,
-    p::NonBroadcastBasisThing,
-    N::ManifoldsBase.RealNumbers,
-)
-    return CachedBasis(
-        DefaultBasis(N),
-        [
-            NonBroadcastBasisThing(ManifoldsBase._euclidean_basis_vector(p.v, i)) for
-            i in eachindex(p.v)
-        ],
-    )
-end
-
-function ManifoldsBase.get_coordinates_orthonormal!(
-    M::DefaultManifold,
-    Y,
-    ::NonBroadcastBasisThing,
-    X::NonBroadcastBasisThing,
-    ::RealNumbers,
-)
-    copyto!(Y, reshape(X.v, manifold_dimension(M)))
-    return Y
-end
-
-function ManifoldsBase.get_vector_orthonormal!(
-    M::DefaultManifold,
-    Y::NonBroadcastBasisThing,
-    ::NonBroadcastBasisThing,
-    X,
-    ::RealNumbers,
-)
-    copyto!(Y.v, reshape(X, representation_size(M)))
-    return Y
-end
-
-function ManifoldsBase.inner(
-    ::DefaultManifold,
-    ::NonBroadcastBasisThing,
-    X::NonBroadcastBasisThing,
-    Y::NonBroadcastBasisThing,
-)
-    return dot(X.v, Y.v)
-end
-
-ManifoldsBase._get_vector_cache_broadcast(::NonBroadcastBasisThing) = Val(false)
 
-DiagonalizingBasisProxy() = DiagonalizingOrthonormalBasis([1.0, 0.0, 0.0])
+s = @__DIR__
+!(s in LOAD_PATH) && (push!(LOAD_PATH, s))
+using ManifoldsBaseTestUtils
 
 @testset "Bases" begin
     @testset "Projected and arbitrary orthonormal basis" begin
diff --git a/test/default_manifold.jl b/test/default_manifold.jl
index 268c721d..ef1c0c86 100644
--- a/test/default_manifold.jl
+++ b/test/default_manifold.jl
@@ -1,22 +1,5 @@
 using ManifoldsBase
-using ManifoldsBase:
-    @manifold_element_forwards, @manifold_vector_forwards, @default_manifold_fallbacks
-using ManifoldsBase: DefaultManifold, AbstractNumbers, RealNumbers, ComplexNumbers
-import ManifoldsBase:
-    allocate_result_type,
-    check_point,
-    distance,
-    embed!,
-    exp!,
-    inner,
-    isapprox,
-    log!,
-    number_eltype,
-    parallel_transport_to!,
-    retract!,
-    inverse_retract!
-
-import Base: angle, convert
+
 using DoubleFloats
 using ForwardDiff
 using LinearAlgebra
@@ -25,192 +8,9 @@ using ReverseDiff
 using StaticArrays
 using Test
 
-struct CustomDefinedRetraction <: ManifoldsBase.AbstractRetractionMethod end
-struct CustomUndefinedRetraction <: ManifoldsBase.AbstractRetractionMethod end
-struct CustomDefinedKeywordRetraction <: ManifoldsBase.AbstractRetractionMethod end
-struct CustomDefinedKeywordInverseRetraction <:
-       ManifoldsBase.AbstractInverseRetractionMethod end
-struct CustomDefinedInverseRetraction <: ManifoldsBase.AbstractInverseRetractionMethod end
-
-struct DefaultPoint{T} <: AbstractManifoldPoint
-    value::T
-end
-DefaultPoint(v::T) where {T} = DefaultPoint{T}(v)
-convert(::Type{DefaultPoint{T}}, v::T) where {T} = DefaultPoint(v)
-Base.size(p::DefaultPoint) = size(p.value)
-Base.eltype(p::DefaultPoint) = eltype(p.value)
-Base.length(p::DefaultPoint) = length(p.value)
-function Base.copyto!(q::DefaultPoint, p::DefaultPoint)
-    copyto!(q.value, p.value)
-    return q
-end
-struct DefaultTVector{T} <: TVector
-    value::T
-end
-DefaultTVector(v::T) where {T} = DefaultTVector{T}(v)
-convert(::Type{DefaultTVector{T}}, ::DefaultPoint, v::T) where {T} = DefaultTVector(v)
-Base.size(X::DefaultTVector) = size(X.value)
-Base.eltype(X::DefaultTVector) = eltype(X.value)
-function Base.fill!(X::DefaultTVector, x)
-    fill!(X.value, x)
-    return X
-end
-function ManifoldsBase.allocate_result_type(
-    ::DefaultManifold,
-    ::typeof(log),
-    ::Tuple{DefaultPoint,DefaultPoint},
-)
-    return DefaultTVector
-end
-function ManifoldsBase.allocate_result_type(
-    ::DefaultManifold,
-    ::typeof(inverse_retract),
-    ::Tuple{DefaultPoint,DefaultPoint},
-)
-    return DefaultTVector
-end
-
-ManifoldsBase.@manifold_element_forwards DefaultPoint value
-ManifoldsBase.@manifold_vector_forwards DefaultTVector value
-ManifoldsBase.@default_manifold_fallbacks ManifoldsBase.DefaultManifold DefaultPoint DefaultTVector value value
-
-function ManifoldsBase._injectivity_radius(::DefaultManifold, ::CustomDefinedRetraction)
-    return 10.0
-end
-function ManifoldsBase._retract!(
-    M::DefaultManifold,
-    q,
-    p,
-    X,
-    t::Number,
-    ::CustomDefinedKeywordRetraction;
-    kwargs...,
-)
-    return retract_custom_kw!(M, q, p, X, t; kwargs...)
-end
-function retract_custom_kw!(
-    ::DefaultManifold,
-    q::DefaultPoint,
-    p::DefaultPoint,
-    X::DefaultTVector,
-    t::Number;
-    scale = 2.0,
-)
-    q.value .= scale .* p.value .+ t .* X.value
-    return q
-end
-function ManifoldsBase._inverse_retract!(
-    M::DefaultManifold,
-    X,
-    p,
-    q,
-    ::CustomDefinedKeywordInverseRetraction;
-    kwargs...,
-)
-    return inverse_retract_custom_kw!(M, X, p, q; kwargs...)
-end
-function inverse_retract_custom_kw!(
-    ::DefaultManifold,
-    X::DefaultTVector,
-    p::DefaultPoint,
-    q::DefaultPoint;
-    scale = 2.0,
-)
-    X.value .= q.value - scale * p.value
-    return X
-end
-
-function ManifoldsBase.retract!(
-    ::DefaultManifold,
-    q::DefaultPoint,
-    p::DefaultPoint,
-    X::DefaultTVector,
-    t::Number,
-    ::CustomDefinedRetraction,
-)
-    q.value .= 2 .* p.value .+ t * X.value
-    return q
-end
-
-function ManifoldsBase.inverse_retract!(
-    ::DefaultManifold,
-    X::DefaultTVector,
-    p::DefaultPoint,
-    q::DefaultPoint,
-    ::CustomDefinedInverseRetraction,
-)
-    X.value .= q.value .- 2 .* p.value
-    return X
-end
-
-struct MatrixVectorTransport{T} <: AbstractVector{T}
-    m::Matrix{T}
-end
-# dummy retractions, inverse retracions for fallback tests - mutating should be enough
-ManifoldsBase.retract_polar!(::DefaultManifold, q, p, X, t::Number) = (q .= p .+ t .* X)
-ManifoldsBase.retract_project!(::DefaultManifold, q, p, X, t::Number) = (q .= p .+ t .* X)
-ManifoldsBase.retract_qr!(::DefaultManifold, q, p, X, t::Number) = (q .= p .+ t .* X)
-function ManifoldsBase.retract_exp_ode!(
-    ::DefaultManifold,
-    q,
-    p,
-    X,
-    t::Number,
-    m::AbstractRetractionMethod,
-    B::ManifoldsBase.AbstractBasis,
-)
-    return (q .= p .+ t .* X)
-end
-
-function ManifoldsBase.retract_pade!(
-    ::DefaultManifold,
-    q,
-    p,
-    X,
-    t::Number,
-    m::PadeRetraction,
-)
-    return (q .= p .+ t .* X)
-end
-function ManifoldsBase.retract_sasaki!(
-    ::DefaultManifold,
-    q,
-    p,
-    X,
-    t::Number,
-    ::SasakiRetraction,
-)
-    return (q .= p .+ t .* X)
-end
-ManifoldsBase.retract_softmax!(::DefaultManifold, q, p, X, t::Number) = (q .= p .+ t .* X)
-ManifoldsBase.get_embedding(M::DefaultManifold) = M # dummy embedding
-ManifoldsBase.inverse_retract_polar!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
-ManifoldsBase.inverse_retract_project!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
-ManifoldsBase.inverse_retract_qr!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
-ManifoldsBase.inverse_retract_softmax!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
-function ManifoldsBase.inverse_retract_nlsolve!(
-    ::DefaultManifold,
-    Y,
-    p,
-    q,
-    m::NLSolveInverseRetraction,
-)
-    return (Y .= q .- p)
-end
-function ManifoldsBase.vector_transport_along_project!(
-    ::DefaultManifold,
-    Y,
-    p,
-    X,
-    c::AbstractVector,
-)
-    return (Y .= X)
-end
-
-
-Base.getindex(x::MatrixVectorTransport, i) = x.m[:, i]
-
-Base.size(x::MatrixVectorTransport) = (size(x.m, 2),)
+s = @__DIR__
+!(s in LOAD_PATH) && (push!(LOAD_PATH, s))
+using ManifoldsBaseTestUtils
 
 @testset "Testing Default (Euclidean)" begin
     M = ManifoldsBase.DefaultManifold(3)
diff --git a/test/empty_manifold.jl b/test/empty_manifold.jl
index c29d6973..6bc437b9 100644
--- a/test/empty_manifold.jl
+++ b/test/empty_manifold.jl
@@ -1,14 +1,9 @@
-using ManifoldsBase
+using ManifoldsBase, Test
+
+
 
 using Test
-import Base: *
-struct NonManifold <: AbstractManifold{ManifoldsBase.ℝ} end
-struct NonMPoint <: AbstractManifoldPoint end
-struct NonTVector <: TVector end
-struct NonCoTVector <: CoTVector end
-struct NotImplementedRetraction <: AbstractRetractionMethod end
-struct NotImplementedInverseRetraction <: AbstractInverseRetractionMethod end
-*(t::Float64, X::NonTVector) = X
+
 @testset "AbstractManifold with empty implementation" begin
     M = NonManifold()
     p = NonMPoint()
diff --git a/test/fibers.jl b/test/fibers.jl
index e64bf1f9..046efc08 100644
--- a/test/fibers.jl
+++ b/test/fibers.jl
@@ -1,20 +1,10 @@
 using RecursiveArrayTools, ManifoldsBase, Test
 using Random
 using ManifoldsBase: DefaultManifold, VectorSpaceType, ℝ, Fiber
-struct TestVectorSpaceType <: VectorSpaceType end
-
-struct TestFiberType <: ManifoldsBase.FiberType end
-
-function ManifoldsBase.fiber_dimension(M::AbstractManifold, ::TestFiberType)
-    return 2 * manifold_dimension(M)
-end
-
-function ManifoldsBase.vector_space_dimension(M::AbstractManifold, ::TestVectorSpaceType)
-    return 2 * manifold_dimension(M)
-end
-
-include("test_manifolds.jl")
 
+s = @__DIR__
+!(s in LOAD_PATH) && (push!(LOAD_PATH, s))
+using ManifoldsBaseTestUtils
 
 @testset "vector space fibers" begin
     M = DefaultManifold(3)
diff --git a/test/manifold_fallbacks.jl b/test/manifold_fallbacks.jl
index 4507ac61..ed85dcb2 100644
--- a/test/manifold_fallbacks.jl
+++ b/test/manifold_fallbacks.jl
@@ -1,7 +1,9 @@
 using Test
 using ManifoldsBase
 
-struct NonManifold <: AbstractManifold{ManifoldsBase.ℝ} end
+s = @__DIR__
+!(s in LOAD_PATH) && (push!(LOAD_PATH, s))
+using ManifoldsBaseTestUtils
 
 @testset "NotImplemented Errors" begin
     M = NonManifold()
diff --git a/test/numerical_checks.jl b/test/numerical_checks.jl
new file mode 100644
index 00000000..02fa47ba
--- /dev/null
+++ b/test/numerical_checks.jl
@@ -0,0 +1,172 @@
+using ManifoldsBase, Plots, Statistics, Test
+# don't show plots actually
+default(; show = false, reuse = true)
+
+s = @__DIR__
+!(s in LOAD_PATH) && (push!(LOAD_PATH, s))
+using ManifoldsBaseTestUtils
+
+@testset "Numerical Check functions" begin
+    @testset "Test retract checks" begin
+        M = TestSphere(10)
+        q = zeros(11)
+        q[1] = 1.0
+        p = zeros(11)
+        p[1:4] .= 1 / sqrt(4)
+        X = log(M, p, q)
+
+        @test check_retraction(M, ProjectionRetraction(), p, X; limits = (-2.5, 0.0))
+        # One call with generating a plot
+        check_retraction(M, ProjectionRetraction(), p, X; limits = (-2.5, 0.0), plot = true)
+
+        # ProjectionRetraction only works <= 1 well in stepsize
+        @test_throws ErrorException check_retraction(
+            M,
+            ProjectionRetraction(),
+            p,
+            X;
+            limits = (-2.5, 2.0),
+            error = :error,
+        )
+        @test !check_retraction(M, ProjectionRetraction(), p, X; limits = (-2.5, 2.0))
+
+        #test window size error
+        @test_throws ErrorException ManifoldsBase.find_best_slope_window(
+            zeros(2),
+            zeros(2),
+            20,
+        )
+        @test_throws ErrorException ManifoldsBase.find_best_slope_window(
+            zeros(2),
+            zeros(2),
+            [2, 20],
+        )
+        @test check_retraction(M, ExponentialRetraction(), p, X; exactness_tol = 1e-7)
+        check_retraction(
+            M,
+            ExponentialRetraction(),
+            p,
+            X;
+            plot = true,
+            exactness_tol = 1e-7,
+        )
+    end
+    @testset "Test inverse_retract checks" begin
+        M = TestSphere(10)
+        q = zeros(11)
+        q[1] = 1.0
+        p = zeros(11)
+        p[1:4] .= 1 / sqrt(4)
+        X = log(M, p, q)
+
+        @test check_inverse_retraction(
+            M,
+            ProjectionInverseRetraction(),
+            p,
+            X;
+            limits = (-2.5, 0.0),
+        )
+        # One call with generating a plot
+        check_inverse_retraction(
+            M,
+            ProjectionInverseRetraction(),
+            p,
+            X;
+            limits = (-2.5, 0.0),
+            plot = true,
+        )
+
+        # ProjectionRetraction only works <= 1 well in stepsize
+        @test_throws ErrorException check_inverse_retraction(
+            M,
+            ProjectionInverseRetraction(),
+            p,
+            X;
+            limits = (-2.5, 2.0), # yields a bit too long tangents
+            error = :error,
+        )
+        @test !check_inverse_retraction(
+            M,
+            ProjectionInverseRetraction(),
+            p,
+            X;
+            limits = (-2.5, 2.0),
+        )
+        # Check exatness case
+        @test check_inverse_retraction(
+            M,
+            LogarithmicInverseRetraction(),
+            p,
+            X;
+            exactness_tol = 1e-7,
+        )
+        check_inverse_retraction(
+            M,
+            LogarithmicInverseRetraction(),
+            p,
+            X;
+            plot = true,
+            exactness_tol = 1e-7,
+        )
+    end
+    @testset "Test vector_transport_to checks" begin
+        M = TestSphere(10)
+        q = zeros(11)
+        q[1] = 1.0
+        p = zeros(11)
+        p[1:4] .= 1 / sqrt(4)
+        X = log(M, p, q)
+        r = zeros(11)
+        r[3] = 1.0
+        Y = log(M, p, r)
+
+        @test check_vector_transport(
+            M,
+            ProjectionTransport(),
+            p,
+            X,
+            Y;
+            second_order = false,
+        )
+        # One call with generating a plot
+        check_vector_transport(
+            M,
+            ProjectionTransport(),
+            p,
+            X,
+            Y;
+            second_order = false,
+            plot = true,
+        )
+
+        # ProjectionRetraction only works <= 1 well in stepsize
+        @test_throws ErrorException check_vector_transport(
+            M,
+            ProjectionTransport(),
+            p,
+            X,
+            Y;
+            limits = (-2.5, 2.0), # yields a bit too long tangents
+            error = :error,
+        )
+        @test !check_vector_transport(
+            M,
+            ProjectionTransport(),
+            p,
+            X,
+            Y;
+            limits = (-2.5, 2.0),
+        )
+        # Check exactness case
+        @test check_vector_transport(M, ParallelTransport(), p, X, Y; exactness_tol = 1e-7)
+        check_vector_transport(
+            M,
+            ParallelTransport(),
+            p,
+            X,
+            Y;
+            plot = true,
+            exactness_tol = 1e-7,
+        )
+    end
+end
diff --git a/test/power.jl b/test/power.jl
index 68f699d2..aed5c441 100644
--- a/test/power.jl
+++ b/test/power.jl
@@ -6,7 +6,9 @@ using StaticArrays
 using LinearAlgebra
 using Random
 
-include("test_manifolds.jl")
+s = @__DIR__
+!(s in LOAD_PATH) && (push!(LOAD_PATH, s))
+using ManifoldsBaseTestUtils
 
 struct TestVectorSpaceType <: VectorSpaceType end
 
@@ -316,7 +318,7 @@ struct TestArrayRepresentation <: AbstractPowerRepresentation end
     @testset "metric conversion" begin
         M = TestSPD(3)
         N = PowerManifold(M, NestedPowerRepresentation(), 2)
-        e = EuclideanMetric()
+        e = ManifoldsBase.EuclideanMetric()
         p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1]
         q = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1]
         P = [p, q]
diff --git a/test/product_manifold.jl b/test/product_manifold.jl
index 9af2cd75..fb516c8d 100644
--- a/test/product_manifold.jl
+++ b/test/product_manifold.jl
@@ -7,7 +7,9 @@ using LinearAlgebra
 using Random
 using RecursiveArrayTools
 
-include("test_manifolds.jl")
+s = @__DIR__
+!(s in LOAD_PATH) && (push!(LOAD_PATH, s))
+using ManifoldsBaseTestUtils
 
 @testset "Product manifold" begin
     M1 = TestSphere(2)
diff --git a/test/runtests.jl b/test/runtests.jl
index 49f5099f..33137e04 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -38,4 +38,5 @@ using ManifoldsBase
     include("metric.jl")
     include("product_manifold.jl")
     include("fibers.jl")
+    include("numerical_checks.jl")
 end
diff --git a/test/test_manifolds.jl b/test/test_manifolds.jl
deleted file mode 100644
index 4ae5540f..00000000
--- a/test/test_manifolds.jl
+++ /dev/null
@@ -1,233 +0,0 @@
-using ManifoldsBase: ℝ, ℂ, DefaultManifold, RealNumbers, EuclideanMetric
-using LinearAlgebra, Random
-
-# minimal sphere implementation for testing more complicated manifolds
-struct TestSphere{N,𝔽} <: AbstractManifold{𝔽} end
-TestSphere(N::Int, 𝔽 = ℝ) = TestSphere{N,𝔽}()
-
-ManifoldsBase.representation_size(::TestSphere{N}) where {N} = (N + 1,)
-
-function ManifoldsBase.change_representer!(
-    M::TestSphere,
-    Y,
-    ::ManifoldsBase.EuclideanMetric,
-    p,
-    X,
-)
-    return copyto!(M, Y, p, X)
-end
-
-function ManifoldsBase.change_metric!(
-    M::TestSphere,
-    Y,
-    ::ManifoldsBase.EuclideanMetric,
-    p,
-    X,
-)
-    return copyto!(M, Y, p, X)
-end
-
-function ManifoldsBase.check_point(M::TestSphere, p; kwargs...)
-    if !isapprox(norm(p), 1.0; kwargs...)
-        return DomainError(
-            norm(p),
-            "The point $(p) does not lie on the $(M) since its norm is not 1.",
-        )
-    end
-    return nothing
-end
-
-function ManifoldsBase.check_vector(M::TestSphere, p, X; kwargs...)
-    if !isapprox(abs(real(dot(p, X))), 0.0; kwargs...)
-        return DomainError(
-            abs(dot(p, X)),
-            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal in the embedding.",
-        )
-    end
-    return nothing
-end
-
-function ManifoldsBase.exp!(M::TestSphere, q, p, X)
-    return exp!(M, q, p, X, one(number_eltype(X)))
-end
-function ManifoldsBase.exp!(::TestSphere, q, p, X, t::Number)
-    θ = abs(t) * norm(X)
-    if θ == 0
-        copyto!(q, p)
-    else
-        X_scale = t * sin(θ) / θ
-        q .= p .* cos(θ) .+ X .* X_scale
-    end
-    return q
-end
-
-function ManifoldsBase.get_basis_diagonalizing(
-    M::TestSphere{n},
-    p,
-    B::DiagonalizingOrthonormalBasis{ℝ},
-) where {n}
-    A = zeros(n + 1, n + 1)
-    A[1, :] = transpose(p)
-    A[2, :] = transpose(B.frame_direction)
-    V = nullspace(A)
-    κ = ones(n)
-    if !iszero(B.frame_direction)
-        # if we have a nonzero direction for the geodesic, add it and it gets curvature zero from the tensor
-        V = hcat(B.frame_direction / norm(M, p, B.frame_direction), V)
-        κ[1] = 0 # no curvature along the geodesic direction, if x!=y
-    end
-    T = typeof(similar(B.frame_direction))
-    Ξ = [convert(T, V[:, i]) for i in 1:n]
-    return CachedBasis(B, κ, Ξ)
-end
-
-@inline function nzsign(z, absz = abs(z))
-    psignz = z / absz
-    return ifelse(iszero(absz), one(psignz), psignz)
-end
-
-
-function ManifoldsBase.get_coordinates_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers)
-    n = manifold_dimension(M)
-    p1 = p[1]
-    cosθ = abs(p1)
-    λ = nzsign(p1, cosθ)
-    pend, Xend = view(p, 2:(n + 1)), view(X, 2:(n + 1))
-    factor = λ * X[1] / (1 + cosθ)
-    Y .= Xend .- pend .* factor
-    return Y
-end
-
-function ManifoldsBase.get_vector_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers)
-    n = manifold_dimension(M)
-    p1 = p[1]
-    cosθ = abs(p1)
-    λ = nzsign(p1, cosθ)
-    pend = view(p, 2:(n + 1))
-    pX = dot(pend, X)
-    factor = pX / (1 + cosθ)
-    Y[1] = -λ * pX
-    Y[2:(n + 1)] .= X .- pend .* factor
-    return Y
-end
-
-ManifoldsBase.inner(::TestSphere, p, X, Y) = dot(X, Y)
-
-ManifoldsBase.injectivity_radius(::TestSphere) = π
-
-function ManifoldsBase.inverse_retract_project!(M::TestSphere, X, p, q)
-    X .= q .- p
-    project!(M, X, p, X)
-    return X
-end
-
-ManifoldsBase.is_flat(M::TestSphere) = manifold_dimension(M) == 1
-
-function ManifoldsBase.log!(::TestSphere, X, p, q)
-    cosθ = clamp(real(dot(p, q)), -1, 1)
-    X .= q .- p .* cosθ
-    nrmX = norm(X)
-    if nrmX > 0
-        X .*= acos(cosθ) / nrmX
-    end
-    return X
-end
-
-ManifoldsBase.manifold_dimension(::TestSphere{N}) where {N} = N
-
-function ManifoldsBase.parallel_transport_to!(::TestSphere, Y, p, X, q)
-    m = p .+ q
-    mnorm2 = real(dot(m, m))
-    factor = 2 * real(dot(X, q)) / mnorm2
-    Y .= X .- m .* factor
-    return Y
-end
-
-function ManifoldsBase.project!(::TestSphere, q, p)
-    q .= p ./ norm(p)
-    return q
-end
-
-function ManifoldsBase.project!(::TestSphere, Y, p, X)
-    Y .= X .- p .* real(dot(p, X))
-    return Y
-end
-
-function Random.rand!(M::TestSphere, pX; vector_at = nothing, σ = one(eltype(pX)))
-    return rand!(Random.default_rng(), M, pX; vector_at = vector_at, σ = σ)
-end
-function Random.rand!(
-    rng::AbstractRNG,
-    M::TestSphere,
-    pX;
-    vector_at = nothing,
-    σ = one(eltype(pX)),
-)
-    if vector_at === nothing
-        project!(M, pX, randn(rng, eltype(pX), representation_size(M)))
-    else
-        n = σ * randn(rng, eltype(pX), size(pX)) # Gaussian in embedding
-        project!(M, pX, vector_at, n) #project to TpM (keeps Gaussianness)
-    end
-    return pX
-end
-
-function ManifoldsBase.riemann_tensor!(M::TestSphere, Xresult, p, X, Y, Z)
-    innerZX = inner(M, p, Z, X)
-    innerZY = inner(M, p, Z, Y)
-    Xresult .= innerZY .* X .- innerZX .* Y
-    return Xresult
-end
-
-function ManifoldsBase.sectional_curvature_max(M::TestSphere)
-    return ifelse(manifold_dimension(M) == 1, 0.0, 1.0)
-end
-function ManifoldsBase.sectional_curvature_min(M::TestSphere)
-    return ifelse(manifold_dimension(M) == 1, 0.0, 1.0)
-end
-
-# from Absil, Mahony, Trumpf, 2013 https://sites.uclouvain.be/absil/2013-01/Weingarten_07PA_techrep.pdf
-function ManifoldsBase.Weingarten!(::TestSphere, Y, p, X, V)
-    return Y .= -X * p' * V
-end
-
-# minimal implementation of SPD (for its nontrivial metric conversion)
-
-struct TestSPD <: AbstractManifold{ℝ}
-    n::Int
-end
-
-function ManifoldsBase.change_representer!(::TestSPD, Y, ::EuclideanMetric, p, X)
-    Y .= p * X * p
-    return Y
-end
-
-function ManifoldsBase.change_metric!(::TestSPD, Y, ::EuclideanMetric, p, X)
-    Y .= p * X
-    return Y
-end
-
-function ManifoldsBase.inner(::TestSPD, p, X, Y)
-    F = cholesky(Symmetric(convert(AbstractMatrix, p)))
-    return dot((F \ Symmetric(X)), (Symmetric(Y) / F))
-end
-
-function spd_sqrt_and_sqrt_inv(p::AbstractMatrix)
-    e = eigen(Symmetric(p))
-    U = e.vectors
-    S = max.(e.values, floatmin(eltype(e.values)))
-    Ssqrt = Diagonal(sqrt.(S))
-    SsqrtInv = Diagonal(1 ./ sqrt.(S))
-    return (Symmetric(U * Ssqrt * transpose(U)), Symmetric(U * SsqrtInv * transpose(U)))
-end
-
-function ManifoldsBase.log!(::TestSPD, X, p, q)
-    (p_sqrt, p_sqrt_inv) = spd_sqrt_and_sqrt_inv(p)
-    T = Symmetric(p_sqrt_inv * convert(AbstractMatrix, q) * p_sqrt_inv)
-    e2 = eigen(T)
-    Se = Diagonal(log.(max.(e2.values, eps())))
-    pUe = p_sqrt * e2.vectors
-    return mul!(X, pUe, Se * transpose(pUe))
-end
-
-ManifoldsBase.representation_size(M::TestSPD) = (M.n, M.n)
diff --git a/test/test_sphere.jl b/test/test_sphere.jl
index 20f2a53c..58b6f7c1 100644
--- a/test/test_sphere.jl
+++ b/test/test_sphere.jl
@@ -3,7 +3,10 @@ using ManifoldsBase
 using ManifoldsBase: ℝ, ℂ, DefaultManifold, RealNumbers
 using Test
 
-include("test_manifolds.jl")
+s = @__DIR__
+!(s in LOAD_PATH) && (push!(LOAD_PATH, s))
+
+using ManifoldsBaseTestUtils
 
 @testset "TestSphere" begin
     @testset "ShootingInverseRetraction" begin