Integrate Nonius benchmark into Catch2

Changes done to Nonius:
* Moved things into "Catch::Benchmark" namespace
* Benchmarks were integrated with `TEST_CASE`/`SECTION`/`GENERATE` macros
* Removed Nonius's parameters for benchmarks, Generators should be used instead
* Added relevant methods to the reporter interface (default-implemented, to avoid
breaking existing 3rd party reporters)
* Async processing is guarded with `_REENTRANT` macro for GCC/Clang, used by default
on MSVC
* Added a macro `CATCH_CONFIG_DISABLE_BENCHMARKING` that removes all traces of
benchmarking from Catch

include/internal/benchmark/catch_benchmark.hpp

@@ -0,0 +1,122 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+ // Benchmark
+#ifndef TWOBLUECUBES_CATCH_BENCHMARK_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_BENCHMARK_HPP_INCLUDED
+
+#include "../catch_config.hpp"
+#include "../catch_context.h"
+#include "../catch_interfaces_reporter.h"
+#include "../catch_test_registry.h"
+
+#include "catch_chronometer.hpp"
+#include "catch_clock.hpp"
+#include "catch_environment.hpp"
+#include "catch_execution_plan.hpp"
+#include "detail/catch_estimate_clock.hpp"
+#include "detail/catch_complete_invoke.hpp"
+#include "detail/catch_analyse.hpp"
+#include "detail/catch_benchmark_function.hpp"
+#include "detail/catch_run_for_at_least.hpp"
+
+#include <algorithm>
+#include <functional>
+#include <string>
+#include <vector>
+#include <cmath>
+
+namespace Catch {
+    namespace Benchmark {
+        struct Benchmark {
+            Benchmark(std::string &&name)
+                : name(std::move(name)) {}
+
+            template <class FUN>
+            Benchmark(std::string &&name, FUN &&func)
+                : fun(std::move(func)), name(std::move(name)) {}
+
+            template <typename Clock>
+            ExecutionPlan<FloatDuration<Clock>> prepare(const IConfig &cfg, Environment<FloatDuration<Clock>> env) const {
+                auto min_time = env.clock_resolution.mean * Detail::minimum_ticks;
+                auto run_time = std::max(min_time, std::chrono::duration_cast<decltype(min_time)>(Detail::warmup_time));
+                auto&& test = Detail::run_for_at_least<Clock>(std::chrono::duration_cast<ClockDuration<Clock>>(run_time), 1, fun);
+                int new_iters = static_cast<int>(std::ceil(min_time * test.iterations / test.elapsed));
+                return { new_iters, test.elapsed / test.iterations * new_iters * cfg.benchmarkSamples(), fun, std::chrono::duration_cast<FloatDuration<Clock>>(Detail::warmup_time), Detail::warmup_iterations };
+            }
+
+            template <typename Clock = default_clock>
+            void run() {
+                IConfigPtr cfg = getCurrentContext().getConfig();
+
+                auto env = Detail::measure_environment<Clock>();
+
+                getResultCapture().benchmarkPreparing(name);
+                CATCH_TRY{
+                    auto plan = user_code([&] {
+                        return prepare<Clock>(*cfg, env);
+                    });
+
+                    BenchmarkInfo info {
+                        name,
+                        plan.estimated_duration.count(),
+                        plan.iterations_per_sample,
+                        cfg->benchmarkSamples(),
+                        cfg->benchmarkResamples(),
+                        env.clock_resolution.mean.count(),
+                        env.clock_cost.mean.count()
+                    };
+
+                    getResultCapture().benchmarkStarting(info);
+
+                    auto samples = user_code([&] {
+                        return plan.template run<Clock>(*cfg, env);
+                    });
+
+                    auto analysis = Detail::analyse(*cfg, env, samples.begin(), samples.end());
+                    BenchmarkStats<std::chrono::duration<double, std::nano>> stats{ info, analysis.samples, analysis.mean, analysis.standard_deviation, analysis.outliers, analysis.outlier_variance };
+                    getResultCapture().benchmarkEnded(stats);
+
+                } CATCH_CATCH_ALL{
+                    if (translateActiveException() != Detail::benchmarkErrorMsg) // benchmark errors have been reported, otherwise rethrow.
+                        std::rethrow_exception(std::current_exception());
+                }
+            }
+
+            // sets lambda to be used in fun *and* executes benchmark!
+            template <typename Fun,
+                typename std::enable_if<!Detail::is_related<Fun, Benchmark>::value, int>::type = 0>
+                Benchmark & operator=(Fun func) {
+                fun = Detail::BenchmarkFunction(func);
+                run();
+                return *this;
+            }
+
+            explicit operator bool() {
+                return true;
+            }
+
+        private:
+            Detail::BenchmarkFunction fun;
+            std::string name;
+        };
+    }
+} // namespace Catch
+
+#define INTERNAL_CATCH_GET_1_ARG(arg1, arg2, ...) arg1
+#define INTERNAL_CATCH_GET_2_ARG(arg1, arg2, ...) arg2
+
+#define INTERNAL_CATCH_BENCHMARK(BenchmarkName, name, benchmarkIndex)\
+    if( Catch::Benchmark::Benchmark BenchmarkName{name} ) \
+        BenchmarkName = [&](int benchmarkIndex)
+
+#define INTERNAL_CATCH_BENCHMARK_ADVANCED(BenchmarkName, name)\
+    if( Catch::Benchmark::Benchmark BenchmarkName{name} ) \
+        BenchmarkName = [&]
+
+#endif // TWOBLUECUBES_CATCH_BENCHMARK_HPP_INCLUDED

include/internal/benchmark/catch_chronometer.hpp

@@ -0,0 +1,71 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// User-facing chronometer
+
+#ifndef TWOBLUECUBES_CATCH_CHRONOMETER_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_CHRONOMETER_HPP_INCLUDED
+
+#include "catch_clock.hpp"
+#include "catch_optimizer.hpp"
+#include "detail/catch_complete_invoke.hpp"
+#include "../catch_meta.hpp"
+
+namespace Catch {
+    namespace Benchmark {
+        namespace Detail {
+            struct ChronometerConcept {
+                virtual void start() = 0;
+                virtual void finish() = 0;
+                virtual ~ChronometerConcept() = default;
+            };
+            template <typename Clock>
+            struct ChronometerModel final : public ChronometerConcept {
+                void start() override { started = Clock::now(); }
+                void finish() override { finished = Clock::now(); }
+
+                ClockDuration<Clock> elapsed() const { return finished - started; }
+
+                TimePoint<Clock> started;
+                TimePoint<Clock> finished;
+            };
+        } // namespace Detail
+
+        struct Chronometer {
+        public:
+            template <typename Fun>
+            void measure(Fun&& fun) { measure(std::forward<Fun>(fun), is_callable<Fun(int)>()); }
+
+            int runs() const { return k; }
+
+            Chronometer(Detail::ChronometerConcept& meter, int k)
+                : impl(&meter)
+                , k(k) {}
+
+        private:
+            template <typename Fun>
+            void measure(Fun&& fun, std::false_type) {
+                measure([&fun](int) { return fun(); }, std::true_type());
+            }
+
+            template <typename Fun>
+            void measure(Fun&& fun, std::true_type) {
+                Detail::optimizer_barrier();
+                impl->start();
+                for (int i = 0; i < k; ++i) invoke_deoptimized(fun, i);
+                impl->finish();
+                Detail::optimizer_barrier();
+            }
+
+            Detail::ChronometerConcept* impl;
+            int k;
+        };
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_CHRONOMETER_HPP_INCLUDED

include/internal/benchmark/catch_clock.hpp

@@ -0,0 +1,46 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// Clocks
+
+#ifndef TWOBLUECUBES_CATCH_CLOCK_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_CLOCK_HPP_INCLUDED
+
+#include <chrono>
+#include <ratio>
+
+namespace Catch {
+    namespace Benchmark {
+        template <unsigned Num, unsigned Den = 1>
+        using ratio = std::ratio<Num, Den>;
+        using milli = ratio<1, 1000>;
+        using micro = ratio<1, 1000000>;
+        using nano = ratio<1, 1000000000>;
+
+        template <typename Clock>
+        using ClockDuration = typename Clock::duration;
+        template <typename Clock>
+        using FloatDuration = std::chrono::duration<double, typename Clock::period>;
+
+        template <typename Clock>
+        using TimePoint = typename Clock::time_point;
+
+        using default_clock = std::chrono::high_resolution_clock;
+
+        template <typename Clock>
+        struct now {
+            TimePoint<Clock> operator()() const {
+                return Clock::now();
+            }
+        };
+
+        using fp_seconds = std::chrono::duration<double, ratio<1>>;
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_CLOCK_HPP_INCLUDED

include/internal/benchmark/catch_constructor.hpp

@@ -0,0 +1,73 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// Constructor and destructor helpers
+
+#ifndef TWOBLUECUBES_CATCH_CONSTRUCTOR_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_CONSTRUCTOR_HPP_INCLUDED
+
+#include <type_traits>
+
+namespace Catch {
+    namespace Detail {
+        template <typename T, bool Destruct>
+        struct ObjectStorage
+        {
+            using TStorage = typename std::aligned_storage<sizeof(T), std::alignment_of<T>::value>::type;
+
+            ObjectStorage() : data() {}
+
+            ObjectStorage(const ObjectStorage& other)
+            {
+                new(&data) T(other.stored_object());
+            }
+
+            ObjectStorage(ObjectStorage&& other)
+            {
+                new(&data) T(std::move(other.stored_object()));
+            }
+
+            ~ObjectStorage() { destruct_on_exit<T>(); }
+
+            template <typename... Args>
+            void construct(Args&&... args)
+            {
+                new (&data) T(std::forward<Args>(args)...);
+            }
+
+            template <bool AllowManualDestruction = !Destruct>
+            typename std::enable_if<AllowManualDestruction>::type destruct()
+            {
+                stored_object().~T();
+            }
+
+        private:
+            // If this is a constructor benchmark, destruct the underlying object
+            template <typename U>
+            void destruct_on_exit(typename std::enable_if<Destruct, U>::type* = 0) { destruct<true>(); }
+            // Otherwise, don't
+            template <typename U>
+            void destruct_on_exit(typename std::enable_if<!Destruct, U>::type* = 0) { }
+
+            T& stored_object()
+            {
+                return *static_cast<T*>(static_cast<void*>(&data));
+            }
+
+            TStorage data;
+        };
+    }
+
+    template <typename T>
+    using storage_for = Detail::ObjectStorage<T, true>;
+
+    template <typename T>
+    using destructable_object = Detail::ObjectStorage<T, false>;
+}
+
+#endif // TWOBLUECUBES_CATCH_CONSTRUCTOR_HPP_INCLUDED

include/internal/benchmark/catch_environment.hpp

@@ -0,0 +1,38 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// Environment information
+
+#ifndef TWOBLUECUBES_CATCH_ENVIRONMENT_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_ENVIRONMENT_HPP_INCLUDED
+
+#include "catch_clock.hpp"
+#include "catch_outlier_classification.hpp"
+
+namespace Catch {
+    namespace Benchmark {
+        template <typename Duration>
+        struct EnvironmentEstimate {
+            Duration mean;
+            OutlierClassification outliers;
+
+            template <typename Duration2>
+            operator EnvironmentEstimate<Duration2>() const {
+                return { mean, outliers };
+            }
+        };
+        template <typename Clock>
+        struct Environment {
+            using clock_type = Clock;
+            EnvironmentEstimate<FloatDuration<Clock>> clock_resolution;
+            EnvironmentEstimate<FloatDuration<Clock>> clock_cost;
+        };
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_ENVIRONMENT_HPP_INCLUDED

include/internal/benchmark/catch_estimate.hpp

@@ -0,0 +1,31 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+ // Statistics estimates
+
+#ifndef TWOBLUECUBES_CATCH_ESTIMATE_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_ESTIMATE_HPP_INCLUDED
+
+namespace Catch {
+    namespace Benchmark {
+        template <typename Duration>
+        struct Estimate {
+            Duration point;
+            Duration lower_bound;
+            Duration upper_bound;
+            double confidence_interval;
+
+            template <typename Duration2>
+            operator Estimate<Duration2>() const {
+                return { point, lower_bound, upper_bound, confidence_interval };
+            }
+        };
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_ESTIMATE_HPP_INCLUDED

include/internal/benchmark/catch_execution_plan.hpp

@@ -0,0 +1,58 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+ // Execution plan
+
+#ifndef TWOBLUECUBES_CATCH_EXECUTION_PLAN_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_EXECUTION_PLAN_HPP_INCLUDED
+
+#include "../catch_config.hpp"
+#include "catch_clock.hpp"
+#include "catch_environment.hpp"
+#include "detail/catch_benchmark_function.hpp"
+#include "detail/catch_repeat.hpp"
+#include "detail/catch_run_for_at_least.hpp"
+
+#include <algorithm>
+
+namespace Catch {
+    namespace Benchmark {
+        template <typename Duration>
+        struct ExecutionPlan {
+            int iterations_per_sample;
+            Duration estimated_duration;
+            Detail::BenchmarkFunction benchmark;
+            Duration warmup_time;
+            int warmup_iterations;
+
+            template <typename Duration2>
+            operator ExecutionPlan<Duration2>() const {
+                return { iterations_per_sample, estimated_duration, benchmark, warmup_time, warmup_iterations };
+            }
+
+            template <typename Clock>
+            std::vector<FloatDuration<Clock>> run(const IConfig &cfg, Environment<FloatDuration<Clock>> env) const {
+                // warmup a bit
+                Detail::run_for_at_least<Clock>(std::chrono::duration_cast<ClockDuration<Clock>>(warmup_time), warmup_iterations, Detail::repeat(now<Clock>{}));
+
+                std::vector<FloatDuration<Clock>> times;
+                times.reserve(cfg.benchmarkSamples());
+                std::generate_n(std::back_inserter(times), cfg.benchmarkSamples(), [this, env] {
+                    Detail::ChronometerModel<Clock> model;
+                    this->benchmark(Chronometer(model, iterations_per_sample));
+                    auto sample_time = model.elapsed() - env.clock_cost.mean;
+                    if (sample_time < FloatDuration<Clock>::zero()) sample_time = FloatDuration<Clock>::zero();
+                    return sample_time / iterations_per_sample;
+                });
+                return times;
+            }
+        };
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_EXECUTION_PLAN_HPP_INCLUDED

include/internal/benchmark/catch_optimizer.hpp

@@ -0,0 +1,68 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+ // Hinting the optimizer
+
+#ifndef TWOBLUECUBES_CATCH_OPTIMIZER_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_OPTIMIZER_HPP_INCLUDED
+
+#if defined(_MSC_VER)
+#   include <atomic> // atomic_thread_fence
+#endif
+
+namespace Catch {
+    namespace Benchmark {
+#if defined(__GNUC__) || defined(__clang__)
+        template <typename T>
+        inline void keep_memory(T* p) {
+            asm volatile("" : : "g"(p) : "memory");
+        }
+        inline void keep_memory() {
+            asm volatile("" : : : "memory");
+        }
+
+        namespace Detail {
+            inline void optimizer_barrier() { keep_memory(); }
+        } // namespace Detail
+#elif defined(_MSC_VER)
+
+#pragma optimize("", off)
+        template <typename T>
+        inline void keep_memory(T* p) {
+            // thanks @milleniumbug
+            *reinterpret_cast<char volatile*>(p) = *reinterpret_cast<char const volatile*>(p);
+        }
+        // TODO equivalent keep_memory()
+#pragma optimize("", on)
+
+        namespace Detail {
+            inline void optimizer_barrier() {
+                std::atomic_thread_fence(std::memory_order_seq_cst);
+            }
+        } // namespace Detail
+
+#endif
+
+        template <typename T>
+        inline void deoptimize_value(T&& x) {
+            keep_memory(&x);
+        }
+
+        template <typename Fn, typename... Args>
+        inline auto invoke_deoptimized(Fn&& fn, Args&&... args) -> typename std::enable_if<!std::is_same<void, decltype(fn(args...))>::value>::type {
+            deoptimize_value(std::forward<Fn>(fn) (std::forward<Args...>(args...)));
+        }
+
+        template <typename Fn, typename... Args>
+        inline auto invoke_deoptimized(Fn&& fn, Args&&... args) -> typename std::enable_if<std::is_same<void, decltype(fn(args...))>::value>::type {
+            std::forward<Fn>(fn) (std::forward<Args...>(args...));
+        }
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_OPTIMIZER_HPP_INCLUDED

include/internal/benchmark/catch_outlier_classification.hpp

@@ -0,0 +1,29 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// Outlier information
+#ifndef TWOBLUECUBES_CATCH_OUTLIERS_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_OUTLIERS_HPP_INCLUDED
+
+namespace Catch {
+    namespace Benchmark {
+        struct OutlierClassification {
+            int samples_seen = 0;
+            int low_severe = 0;     // more than 3 times IQR below Q1
+            int low_mild = 0;       // 1.5 to 3 times IQR below Q1
+            int high_mild = 0;      // 1.5 to 3 times IQR above Q3
+            int high_severe = 0;    // more than 3 times IQR above Q3
+
+            int total() const {
+                return low_severe + low_mild + high_mild + high_severe;
+            }
+        };
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_OUTLIERS_HPP_INCLUDED

include/internal/benchmark/catch_sample_analysis.hpp

@@ -0,0 +1,50 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// Benchmark results
+
+#ifndef TWOBLUECUBES_CATCH_BENCHMARK_RESULTS_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_BENCHMARK_RESULTS_HPP_INCLUDED
+
+#include "catch_clock.hpp"
+#include "catch_estimate.hpp"
+#include "catch_outlier_classification.hpp"
+
+#include <algorithm>
+#include <vector>
+#include <string>
+#include <iterator>
+
+namespace Catch {
+    namespace Benchmark {
+        template <typename Duration>
+        struct SampleAnalysis {
+            std::vector<Duration> samples;
+            Estimate<Duration> mean;
+            Estimate<Duration> standard_deviation;
+            OutlierClassification outliers;
+            double outlier_variance;
+
+            template <typename Duration2>
+            operator SampleAnalysis<Duration2>() const {
+                std::vector<Duration2> samples2;
+                samples2.reserve(samples.size());
+                std::transform(samples.begin(), samples.end(), std::back_inserter(samples2), [](Duration d) { return Duration2(d); });
+                return {
+                    std::move(samples2),
+                    mean,
+                    standard_deviation,
+                    outliers,
+                    outlier_variance,
+                };
+            }
+        };
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_BENCHMARK_RESULTS_HPP_INCLUDED

include/internal/benchmark/detail/catch_analyse.hpp

@@ -0,0 +1,78 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+ // Run and analyse one benchmark
+
+#ifndef TWOBLUECUBES_CATCH_DETAIL_ANALYSE_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_DETAIL_ANALYSE_HPP_INCLUDED
+
+#include "../catch_clock.hpp"
+#include "../catch_sample_analysis.hpp"
+#include "catch_stats.hpp"
+
+#include <algorithm>
+#include <iterator>
+#include <vector>
+
+namespace Catch {
+    namespace Benchmark {
+        namespace Detail {
+            template <typename Duration, typename Iterator>
+            SampleAnalysis<Duration> analyse(const IConfig &cfg, Environment<Duration>, Iterator first, Iterator last) {
+                if (!cfg.benchmarkNoAnalysis()) {
+                    std::vector<double> samples;
+                    samples.reserve(last - first);
+                    std::transform(first, last, std::back_inserter(samples), [](Duration d) { return d.count(); });
+
+                    auto analysis = Catch::Benchmark::Detail::analyse_samples(cfg.benchmarkConfidenceInterval(), cfg.benchmarkResamples(), samples.begin(), samples.end());
+                    auto outliers = Catch::Benchmark::Detail::classify_outliers(samples.begin(), samples.end());
+
+                    auto wrap_estimate = [](Estimate<double> e) {
+                        return Estimate<Duration> {
+                            Duration(e.point),
+                                Duration(e.lower_bound),
+                                Duration(e.upper_bound),
+                                e.confidence_interval,
+                        };
+                    };
+                    std::vector<Duration> samples2;
+                    samples2.reserve(samples.size());
+                    std::transform(samples.begin(), samples.end(), std::back_inserter(samples2), [](double d) { return Duration(d); });
+                    return {
+                        std::move(samples2),
+                        wrap_estimate(analysis.mean),
+                        wrap_estimate(analysis.standard_deviation),
+                        outliers,
+                        analysis.outlier_variance,
+                    };
+                } else {
+                    std::vector<Duration> samples; 
+                    samples.reserve(last - first);
+
+                    Duration mean = Duration(0);
+                    int i = 0;
+                    for (auto it = first; it < last; ++it, ++i) {
+                        samples.push_back(Duration(*it));
+                        mean += Duration(*it);
+                    }
+                    mean /= i;
+
+                    return {
+                        std::move(samples),
+                        Estimate<Duration>{mean, mean, mean, 0.0},
+                        Estimate<Duration>{Duration(0), Duration(0), Duration(0), 0.0},
+                        OutlierClassification{},
+                        0.0
+                    };
+                }
+            }
+        } // namespace Detail
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_DETAIL_ANALYSE_HPP_INCLUDED

include/internal/benchmark/detail/catch_benchmark_function.hpp

@@ -0,0 +1,105 @@
+    /*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+ // Dumb std::function implementation for consistent call overhead
+
+#ifndef TWOBLUECUBES_CATCH_DETAIL_BENCHMARK_FUNCTION_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_DETAIL_BENCHMARK_FUNCTION_HPP_INCLUDED
+
+#include "../catch_chronometer.hpp"
+#include "catch_complete_invoke.hpp"
+#include "../../catch_meta.hpp"
+
+#include <cassert>
+#include <type_traits>
+#include <utility>
+#include <memory>
+
+namespace Catch {
+    namespace Benchmark {
+        namespace Detail {
+            template <typename T>
+            using Decay = typename std::decay<T>::type;
+            template <typename T, typename U>
+            struct is_related
+                : std::is_same<Decay<T>, Decay<U>> {};
+
+            /// We need to reinvent std::function because every piece of code that might add overhead
+            /// in a measurement context needs to have consistent performance characteristics so that we
+            /// can account for it in the measurement.
+            /// Implementations of std::function with optimizations that aren't always applicable, like
+            /// small buffer optimizations, are not uncommon.
+            /// This is effectively an implementation of std::function without any such optimizations;
+            /// it may be slow, but it is consistently slow.
+            struct BenchmarkFunction {
+            private:
+                struct callable {
+                    virtual void call(Chronometer meter) const = 0;
+                    virtual callable* clone() const = 0;
+                    virtual ~callable() = default;
+                };
+                template <typename Fun>
+                struct model : public callable {
+                    model(Fun&& fun) : fun(std::move(fun)) {}
+                    model(Fun const& fun) : fun(fun) {}
+
+                    model<Fun>* clone() const override { return new model<Fun>(*this); }
+
+                    void call(Chronometer meter) const override {
+                        call(meter, is_callable<Fun(Chronometer)>());
+                    }
+                    void call(Chronometer meter, std::true_type) const {
+                        fun(meter);
+                    }
+                    void call(Chronometer meter, std::false_type) const {
+                        meter.measure(fun);
+                    }
+
+                    Fun fun;
+                };
+
+                struct do_nothing { void operator()() const {} };
+
+                template <typename T>
+                BenchmarkFunction(model<T>* c) : f(c) {}
+
+            public:
+                BenchmarkFunction()
+                    : f(new model<do_nothing>{ {} }) {}
+
+                template <typename Fun,
+                    typename std::enable_if<!is_related<Fun, BenchmarkFunction>::value, int>::type = 0>
+                    BenchmarkFunction(Fun&& fun)
+                    : f(new model<typename std::decay<Fun>::type>(std::forward<Fun>(fun))) {}
+
+                BenchmarkFunction(BenchmarkFunction&& that)
+                    : f(std::move(that.f)) {}
+
+                BenchmarkFunction(BenchmarkFunction const& that)
+                    : f(that.f->clone()) {}
+
+                BenchmarkFunction& operator=(BenchmarkFunction&& that) {
+                    f = std::move(that.f);
+                    return *this;
+                }
+
+                BenchmarkFunction& operator=(BenchmarkFunction const& that) {
+                    f.reset(that.f->clone());
+                    return *this;
+                }
+
+                void operator()(Chronometer meter) const { f->call(meter); }
+
+            private:
+                std::unique_ptr<callable> f;
+            };
+        } // namespace Detail
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_DETAIL_BENCHMARK_FUNCTION_HPP_INCLUDED

include/internal/benchmark/detail/catch_complete_invoke.hpp

@@ -0,0 +1,69 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// Invoke with a special case for void
+
+#ifndef TWOBLUECUBES_CATCH_DETAIL_COMPLETE_INVOKE_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_DETAIL_COMPLETE_INVOKE_HPP_INCLUDED
+
+#include "../../catch_enforce.h"
+
+#include <type_traits>
+#include <utility>
+
+namespace Catch {
+    namespace Benchmark {
+        namespace Detail {
+            template <typename T>
+            struct CompleteType { using type = T; };
+            template <>
+            struct CompleteType<void> { struct type {}; };
+
+            template <typename T>
+            using CompleteType_t = typename CompleteType<T>::type;
+
+            template <typename Result>
+            struct CompleteInvoker {
+                template <typename Fun, typename... Args>
+                static Result invoke(Fun&& fun, Args&&... args) {
+                    return std::forward<Fun>(fun)(std::forward<Args>(args)...);
+                }
+            };
+            template <>
+            struct CompleteInvoker<void> {
+                template <typename Fun, typename... Args>
+                static CompleteType_t<void> invoke(Fun&& fun, Args&&... args) {
+                    std::forward<Fun>(fun)(std::forward<Args>(args)...);
+                    return {};
+                }
+            };
+            template <typename Sig>
+            using ResultOf_t = typename std::result_of<Sig>::type;
+
+            // invoke and not return void :(
+            template <typename Fun, typename... Args>
+            CompleteType_t<ResultOf_t<Fun(Args...)>> complete_invoke(Fun&& fun, Args&&... args) {
+                return CompleteInvoker<ResultOf_t<Fun(Args...)>>::invoke(std::forward<Fun>(fun), std::forward<Args>(args)...);
+            }
+
+            const std::string benchmarkErrorMsg = "a benchmark failed to run successfully";
+        } // namespace Detail
+
+        template <typename Fun>
+        Detail::CompleteType_t<Detail::ResultOf_t<Fun()>> user_code(Fun&& fun) {
+            CATCH_TRY{
+                return Detail::complete_invoke(std::forward<Fun>(fun));
+            } CATCH_CATCH_ALL{
+                getResultCapture().benchmarkFailed(translateActiveException());
+                CATCH_RUNTIME_ERROR(Detail::benchmarkErrorMsg);
+            }
+        }
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_DETAIL_COMPLETE_INVOKE_HPP_INCLUDED

include/internal/benchmark/detail/catch_estimate_clock.hpp

@@ -0,0 +1,113 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+ // Environment measurement
+
+#ifndef TWOBLUECUBES_CATCH_DETAIL_ESTIMATE_CLOCK_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_DETAIL_ESTIMATE_CLOCK_HPP_INCLUDED
+
+#include "../catch_clock.hpp"
+#include "../catch_environment.hpp"
+#include "catch_stats.hpp"
+#include "catch_measure.hpp"
+#include "catch_run_for_at_least.hpp"
+#include "../catch_clock.hpp"
+
+#include <algorithm>
+#include <iterator>
+#include <tuple>
+#include <vector>
+#include <cmath>
+
+namespace Catch {
+    namespace Benchmark {
+        namespace Detail {
+            template <typename Clock>
+            std::vector<double> resolution(int k) {
+                std::vector<TimePoint<Clock>> times;
+                times.reserve(k + 1);
+                std::generate_n(std::back_inserter(times), k + 1, now<Clock>{});
+
+                std::vector<double> deltas;
+                deltas.reserve(k);
+                std::transform(std::next(times.begin()), times.end(), times.begin(),
+                    std::back_inserter(deltas),
+                    [](TimePoint<Clock> a, TimePoint<Clock> b) { return static_cast<double>((a - b).count()); });
+
+                return deltas;
+            }
+
+            const auto warmup_iterations = 10000;
+            const auto warmup_time = std::chrono::milliseconds(100);
+            const auto minimum_ticks = 1000;
+            const auto warmup_seed = 10000;
+            const auto clock_resolution_estimation_time = std::chrono::milliseconds(500);
+            const auto clock_cost_estimation_time_limit = std::chrono::seconds(1);
+            const auto clock_cost_estimation_tick_limit = 100000;
+            const auto clock_cost_estimation_time = std::chrono::milliseconds(10);
+            const auto clock_cost_estimation_iterations = 10000;
+
+            template <typename Clock>
+            int warmup() {
+                return run_for_at_least<Clock>(std::chrono::duration_cast<ClockDuration<Clock>>(warmup_time), warmup_seed, &resolution<Clock>)
+                    .iterations;
+            }
+            template <typename Clock>
+            EnvironmentEstimate<FloatDuration<Clock>> estimate_clock_resolution(int iterations) {
+                auto r = run_for_at_least<Clock>(std::chrono::duration_cast<ClockDuration<Clock>>(clock_resolution_estimation_time), iterations, &resolution<Clock>)
+                    .result;
+                return {
+                    FloatDuration<Clock>(mean(r.begin(), r.end())),
+                    classify_outliers(r.begin(), r.end()),
+                };
+            }
+            template <typename Clock>
+            EnvironmentEstimate<FloatDuration<Clock>> estimate_clock_cost(FloatDuration<Clock> resolution) {
+                auto time_limit = std::min(resolution * clock_cost_estimation_tick_limit, FloatDuration<Clock>(clock_cost_estimation_time_limit));
+                auto time_clock = [](int k) {
+                    return Detail::measure<Clock>([k] {
+                        for (int i = 0; i < k; ++i) {
+                            volatile auto ignored = Clock::now();
+                            (void)ignored;
+                        }
+                    }).elapsed;
+                };
+                time_clock(1);
+                int iters = clock_cost_estimation_iterations;
+                auto&& r = run_for_at_least<Clock>(std::chrono::duration_cast<ClockDuration<Clock>>(clock_cost_estimation_time), iters, time_clock);
+                std::vector<double> times;
+                int nsamples = static_cast<int>(std::ceil(time_limit / r.elapsed));
+                times.reserve(nsamples);
+                std::generate_n(std::back_inserter(times), nsamples, [time_clock, &r] {
+                    return static_cast<double>((time_clock(r.iterations) / r.iterations).count());
+                });
+                return {
+                    FloatDuration<Clock>(mean(times.begin(), times.end())),
+                    classify_outliers(times.begin(), times.end()),
+                };
+            }
+
+            template <typename Clock>
+            Environment<FloatDuration<Clock>> measure_environment() {
+                static Environment<FloatDuration<Clock>>* env = nullptr;
+                if (env) {
+                    return *env;
+                }
+
+                auto iters = Detail::warmup<Clock>();
+                auto resolution = Detail::estimate_clock_resolution<Clock>(iters);
+                auto cost = Detail::estimate_clock_cost<Clock>(resolution.mean);
+
+                env = new Environment<FloatDuration<Clock>>{ resolution, cost };
+                return *env;
+            }
+        } // namespace Detail
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_DETAIL_ESTIMATE_CLOCK_HPP_INCLUDED

include/internal/benchmark/detail/catch_measure.hpp

@@ -0,0 +1,35 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// Measure
+
+#ifndef TWOBLUECUBES_CATCH_DETAIL_MEASURE_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_DETAIL_MEASURE_HPP_INCLUDED
+
+#include "../catch_clock.hpp"
+#include "catch_complete_invoke.hpp"
+#include "catch_timing.hpp"
+
+#include <utility>
+
+namespace Catch {
+    namespace Benchmark {
+        namespace Detail {
+            template <typename Clock, typename Fun, typename... Args>
+            TimingOf<Clock, Fun(Args...)> measure(Fun&& fun, Args&&... args) {
+                auto start = Clock::now();
+                auto&& r = Detail::complete_invoke(fun, std::forward<Args>(args)...);
+                auto end = Clock::now();
+                auto delta = end - start;
+                return { delta, std::forward<decltype(r)>(r), 1 };
+            }
+        } // namespace Detail
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_DETAIL_MEASURE_HPP_INCLUDED

include/internal/benchmark/detail/catch_repeat.hpp

@@ -0,0 +1,37 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// repeat algorithm
+
+#ifndef TWOBLUECUBES_CATCH_DETAIL_REPEAT_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_DETAIL_REPEAT_HPP_INCLUDED
+
+#include <type_traits>
+#include <utility>
+
+namespace Catch {
+    namespace Benchmark {
+        namespace Detail {
+            template <typename Fun>
+            struct repeater {
+                void operator()(int k) const {
+                    for (int i = 0; i < k; ++i) {
+                        fun();
+                    }
+                }
+                Fun fun;
+            };
+            template <typename Fun>
+            repeater<typename std::decay<Fun>::type> repeat(Fun&& fun) {
+                return { std::forward<Fun>(fun) };
+            }
+        } // namespace Detail
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_DETAIL_REPEAT_HPP_INCLUDED

include/internal/benchmark/detail/catch_run_for_at_least.hpp

@@ -0,0 +1,65 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// Run a function for a minimum amount of time
+
+#ifndef TWOBLUECUBES_CATCH_RUN_FOR_AT_LEAST_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_RUN_FOR_AT_LEAST_HPP_INCLUDED
+
+#include "../catch_clock.hpp"
+#include "../catch_chronometer.hpp"
+#include "catch_measure.hpp"
+#include "catch_complete_invoke.hpp"
+#include "catch_timing.hpp"
+#include "../../catch_meta.hpp"
+
+#include <utility>
+#include <type_traits>
+
+namespace Catch {
+    namespace Benchmark {
+        namespace Detail {
+            template <typename Clock, typename Fun>
+            TimingOf<Clock, Fun(int)> measure_one(Fun&& fun, int iters, std::false_type) {
+                return Detail::measure<Clock>(fun, iters);
+            }
+            template <typename Clock, typename Fun>
+            TimingOf<Clock, Fun(Chronometer)> measure_one(Fun&& fun, int iters, std::true_type) {
+                Detail::ChronometerModel<Clock> meter;
+                auto&& result = Detail::complete_invoke(fun, Chronometer(meter, iters));
+
+                return { meter.elapsed(), std::move(result), iters };
+            }
+
+            template <typename Clock, typename Fun>
+            using run_for_at_least_argument_t = typename std::conditional<is_callable<Fun(Chronometer)>::value, Chronometer, int>::type;
+
+            struct optimized_away_error : std::exception {
+                const char* what() const noexcept override {
+                    return "could not measure benchmark, maybe it was optimized away";
+                }
+            };
+
+            template <typename Clock, typename Fun>
+            TimingOf<Clock, Fun(run_for_at_least_argument_t<Clock, Fun>)> run_for_at_least(ClockDuration<Clock> how_long, int seed, Fun&& fun) {
+                auto iters = seed;
+                while (iters < (1 << 30)) {
+                    auto&& Timing = measure_one<Clock>(fun, iters, is_callable<Fun(Chronometer)>());
+
+                    if (Timing.elapsed >= how_long) {
+                        return { Timing.elapsed, std::move(Timing.result), iters };
+                    }
+                    iters *= 2;
+                }
+                throw optimized_away_error{};
+            }
+        } // namespace Detail
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_RUN_FOR_AT_LEAST_HPP_INCLUDED

include/internal/benchmark/detail/catch_stats.hpp

@@ -0,0 +1,342 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// Statistical analysis tools
+
+#ifndef TWOBLUECUBES_CATCH_DETAIL_ANALYSIS_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_DETAIL_ANALYSIS_HPP_INCLUDED
+
+#include "../catch_clock.hpp"
+#include "../catch_estimate.hpp"
+#include "../catch_outlier_classification.hpp"
+
+#include <algorithm>
+#include <cassert>
+#include <functional>
+#include <iterator>
+#include <vector>
+#include <array>
+#include <random>
+#include <numeric>
+#include <tuple>
+#include <cmath>
+#include <utility>
+#include <cstddef>
+
+#ifdef CATCH_USE_ASYNC
+#include <future>
+#endif
+
+namespace Catch {
+    namespace Benchmark {
+        namespace Detail {
+            using sample = std::vector<double>;
+
+            template <typename Iterator>
+            double weighted_average_quantile(int k, int q, Iterator first, Iterator last) {
+                auto count = last - first;
+                double idx = (count - 1) * k / static_cast<double>(q);
+                int j = static_cast<int>(idx);
+                double g = idx - j;
+                std::nth_element(first, first + j, last);
+                auto xj = first[j];
+                if (g == 0) return xj;
+
+                auto xj1 = *std::min_element(first + (j + 1), last);
+                return xj + g * (xj1 - xj);
+            }
+
+            template <typename Iterator>
+            OutlierClassification classify_outliers(Iterator first, Iterator last) {
+                std::vector<double> copy(first, last);
+
+                auto q1 = weighted_average_quantile(1, 4, copy.begin(), copy.end());
+                auto q3 = weighted_average_quantile(3, 4, copy.begin(), copy.end());
+                auto iqr = q3 - q1;
+                auto los = q1 - (iqr * 3.);
+                auto lom = q1 - (iqr * 1.5);
+                auto him = q3 + (iqr * 1.5);
+                auto his = q3 + (iqr * 3.);
+
+                OutlierClassification o;
+                for (; first != last; ++first) {
+                    auto&& t = *first;
+                    if (t < los) ++o.low_severe;
+                    else if (t < lom) ++o.low_mild;
+                    else if (t > his) ++o.high_severe;
+                    else if (t > him) ++o.high_mild;
+                    ++o.samples_seen;
+                }
+                return o;
+            }
+
+            template <typename Iterator>
+            double mean(Iterator first, Iterator last) {
+                auto count = last - first;
+                double sum = std::accumulate(first, last, 0.);
+                return sum / count;
+            }
+
+            template <typename Iterator>
+            double standard_deviation(Iterator first, Iterator last) {
+                auto m = mean(first, last);
+                double variance = std::accumulate(first, last, 0., [m](double a, double b) {
+                    double diff = b - m;
+                    return a + diff * diff;
+                }) / (last - first);
+                return std::sqrt(variance);
+            }
+
+            template <typename URng, typename Iterator, typename Estimator>
+            sample resample(URng& rng, int resamples, Iterator first, Iterator last, Estimator& estimator) {
+                auto n = last - first;
+                std::uniform_int_distribution<decltype(n)> dist(0, n - 1);
+
+                sample out;
+                out.reserve(resamples);
+                std::generate_n(std::back_inserter(out), resamples, [n, first, &estimator, &dist, &rng] {
+                    std::vector<double> resampled;
+                    resampled.reserve(n);
+                    std::generate_n(std::back_inserter(resampled), n, [first, &dist, &rng] { return first[dist(rng)]; });
+                    return estimator(resampled.begin(), resampled.end());
+                });
+                std::sort(out.begin(), out.end());
+                return out;
+            }
+
+            template <typename Estimator, typename Iterator>
+            sample jackknife(Estimator&& estimator, Iterator first, Iterator last) {
+                auto n = last - first;
+                auto second = std::next(first);
+                sample results;
+                results.reserve(n);
+
+                for (auto it = first; it != last; ++it) {
+                    std::iter_swap(it, first);
+                    results.push_back(estimator(second, last));
+                }
+
+                return results;
+            }
+
+            inline double normal_cdf(double x) {
+                return std::erfc(-x / std::sqrt(2.0)) / 2.0;
+            }
+
+            inline double erf_inv(double x) {
+                // Code accompanying the article "Approximating the erfinv function" in GPU Computing Gems, Volume 2
+                double w, p;
+
+                w = -log((1.0 - x)*(1.0 + x));
+
+                if (w < 6.250000) {
+                    w = w - 3.125000;
+                    p = -3.6444120640178196996e-21;
+                    p = -1.685059138182016589e-19 + p * w;
+                    p = 1.2858480715256400167e-18 + p * w;
+                    p = 1.115787767802518096e-17 + p * w;
+                    p = -1.333171662854620906e-16 + p * w;
+                    p = 2.0972767875968561637e-17 + p * w;
+                    p = 6.6376381343583238325e-15 + p * w;
+                    p = -4.0545662729752068639e-14 + p * w;
+                    p = -8.1519341976054721522e-14 + p * w;
+                    p = 2.6335093153082322977e-12 + p * w;
+                    p = -1.2975133253453532498e-11 + p * w;
+                    p = -5.4154120542946279317e-11 + p * w;
+                    p = 1.051212273321532285e-09 + p * w;
+                    p = -4.1126339803469836976e-09 + p * w;
+                    p = -2.9070369957882005086e-08 + p * w;
+                    p = 4.2347877827932403518e-07 + p * w;
+                    p = -1.3654692000834678645e-06 + p * w;
+                    p = -1.3882523362786468719e-05 + p * w;
+                    p = 0.0001867342080340571352 + p * w;
+                    p = -0.00074070253416626697512 + p * w;
+                    p = -0.0060336708714301490533 + p * w;
+                    p = 0.24015818242558961693 + p * w;
+                    p = 1.6536545626831027356 + p * w;
+                } else if (w < 16.000000) {
+                    w = sqrt(w) - 3.250000;
+                    p = 2.2137376921775787049e-09;
+                    p = 9.0756561938885390979e-08 + p * w;
+                    p = -2.7517406297064545428e-07 + p * w;
+                    p = 1.8239629214389227755e-08 + p * w;
+                    p = 1.5027403968909827627e-06 + p * w;
+                    p = -4.013867526981545969e-06 + p * w;
+                    p = 2.9234449089955446044e-06 + p * w;
+                    p = 1.2475304481671778723e-05 + p * w;
+                    p = -4.7318229009055733981e-05 + p * w;
+                    p = 6.8284851459573175448e-05 + p * w;
+                    p = 2.4031110387097893999e-05 + p * w;
+                    p = -0.0003550375203628474796 + p * w;
+                    p = 0.00095328937973738049703 + p * w;
+                    p = -0.0016882755560235047313 + p * w;
+                    p = 0.0024914420961078508066 + p * w;
+                    p = -0.0037512085075692412107 + p * w;
+                    p = 0.005370914553590063617 + p * w;
+                    p = 1.0052589676941592334 + p * w;
+                    p = 3.0838856104922207635 + p * w;
+                } else {
+                    w = sqrt(w) - 5.000000;
+                    p = -2.7109920616438573243e-11;
+                    p = -2.5556418169965252055e-10 + p * w;
+                    p = 1.5076572693500548083e-09 + p * w;
+                    p = -3.7894654401267369937e-09 + p * w;
+                    p = 7.6157012080783393804e-09 + p * w;
+                    p = -1.4960026627149240478e-08 + p * w;
+                    p = 2.9147953450901080826e-08 + p * w;
+                    p = -6.7711997758452339498e-08 + p * w;
+                    p = 2.2900482228026654717e-07 + p * w;
+                    p = -9.9298272942317002539e-07 + p * w;
+                    p = 4.5260625972231537039e-06 + p * w;
+                    p = -1.9681778105531670567e-05 + p * w;
+                    p = 7.5995277030017761139e-05 + p * w;
+                    p = -0.00021503011930044477347 + p * w;
+                    p = -0.00013871931833623122026 + p * w;
+                    p = 1.0103004648645343977 + p * w;
+                    p = 4.8499064014085844221 + p * w;
+                }
+                return p * x;
+            }
+
+            inline double erfc_inv(double x) {
+                return erf_inv(1.0 - x);
+            }
+
+            inline double normal_quantile(double p) {
+                static const double ROOT_TWO = std::sqrt(2.0);
+
+                double result = 0.0;
+                assert(p >= 0 && p <= 1);
+                if (p < 0 || p > 1) {
+                    return result;
+                }
+
+                result = -erfc_inv(2.0 * p);
+                // result *= normal distribution standard deviation (1.0) * sqrt(2)
+                result *= /*sd * */ ROOT_TWO;
+                // result += normal disttribution mean (0)
+                return result;
+            }
+
+            template <typename Iterator, typename Estimator>
+            Estimate<double> bootstrap(double confidence_level, Iterator first, Iterator last, sample const& resample, Estimator&& estimator) {
+                auto n_samples = last - first;
+
+                double point = estimator(first, last);
+                // Degenerate case with a single sample
+                if (n_samples == 1) return { point, point, point, confidence_level };
+
+                sample jack = jackknife(estimator, first, last);
+                double jack_mean = mean(jack.begin(), jack.end());
+                double sum_squares, sum_cubes;
+                std::tie(sum_squares, sum_cubes) = std::accumulate(jack.begin(), jack.end(), std::make_pair(0., 0.), [jack_mean](std::pair<double, double> sqcb, double x) -> std::pair<double, double> {
+                    auto d = jack_mean - x;
+                    auto d2 = d * d;
+                    auto d3 = d2 * d;
+                    return { sqcb.first + d2, sqcb.second + d3 };
+                });
+
+                double accel = sum_cubes / (6 * std::pow(sum_squares, 1.5));
+                int n = static_cast<int>(resample.size());
+                double prob_n = std::count_if(resample.begin(), resample.end(), [point](double x) { return x < point; }) / (double)n;
+                // degenerate case with uniform samples
+                if (prob_n == 0) return { point, point, point, confidence_level };
+
+                double bias = normal_quantile(prob_n);
+                double z1 = normal_quantile((1. - confidence_level) / 2.);
+
+                auto cumn = [n](double x) -> int {
+                    return std::lround(normal_cdf(x) * n); };
+                auto a = [bias, accel](double b) { return bias + b / (1. - accel * b); };
+                double b1 = bias + z1;
+                double b2 = bias - z1;
+                double a1 = a(b1);
+                double a2 = a(b2);
+                auto lo = std::max(cumn(a1), 0);
+                auto hi = std::min(cumn(a2), n - 1);
+
+                return { point, resample[lo], resample[hi], confidence_level };
+            }
+
+            inline double outlier_variance(Estimate<double> mean, Estimate<double> stddev, int n) {
+                double sb = stddev.point;
+                double mn = mean.point / n;
+                double mg_min = mn / 2.;
+                double sg = std::min(mg_min / 4., sb / std::sqrt(n));
+                double sg2 = sg * sg;
+                double sb2 = sb * sb;
+
+                auto c_max = [n, mn, sb2, sg2](double x) -> double {
+                    double k = mn - x;
+                    double d = k * k;
+                    double nd = n * d;
+                    double k0 = -n * nd;
+                    double k1 = sb2 - n * sg2 + nd;
+                    double det = k1 * k1 - 4 * sg2 * k0;
+                    return (int)(-2. * k0 / (k1 + std::sqrt(det)));
+                };
+
+                auto var_out = [n, sb2, sg2](double c) {
+                    double nc = n - c;
+                    return (nc / n) * (sb2 - nc * sg2);
+                };
+
+                return std::min(var_out(1), var_out(std::min(c_max(0.), c_max(mg_min)))) / sb2;
+            }
+
+            struct bootstrap_analysis {
+                Estimate<double> mean;
+                Estimate<double> standard_deviation;
+                double outlier_variance;
+            };
+
+            template <typename Iterator>
+            bootstrap_analysis analyse_samples(double confidence_level, int n_resamples, Iterator first, Iterator last) {
+                static std::random_device entropy;
+
+                auto n = static_cast<int>(last - first); // seriously, one can't use integral types without hell in C++
+
+                auto mean = &Detail::mean<Iterator>;
+                auto stddev = &Detail::standard_deviation<Iterator>;
+
+#ifdef CATCH_USE_ASYNC
+                auto Estimate = [=](double(*f)(Iterator, Iterator)) {
+                    auto seed = entropy();
+                    return std::async(std::launch::async, [=] {
+                        std::mt19937 rng(seed);
+                        auto resampled = resample(rng, n_resamples, first, last, f);
+                        return bootstrap(confidence_level, first, last, resampled, f);
+                    });
+                };
+
+                auto mean_future = Estimate(mean);
+                auto stddev_future = Estimate(stddev);
+
+                auto mean_estimate = mean_future.get();
+                auto stddev_estimate = stddev_future.get();
+#else
+                auto Estimate = [=](double(*f)(Iterator, Iterator)) {
+                    auto seed = entropy();
+                    std::mt19937 rng(seed);
+                    auto resampled = resample(rng, n_resamples, first, last, f);
+                    return bootstrap(confidence_level, first, last, resampled, f);
+                };
+
+                auto mean_estimate = Estimate(mean);
+                auto stddev_estimate = Estimate(stddev);
+#endif // CATCH_USE_ASYNC
+
+                double outlier_variance = Detail::outlier_variance(mean_estimate, stddev_estimate, n);
+
+                return { mean_estimate, stddev_estimate, outlier_variance };
+            }
+        } // namespace Detail
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_DETAIL_ANALYSIS_HPP_INCLUDED

include/internal/benchmark/detail/catch_timing.hpp

@@ -0,0 +1,33 @@
+/*
+ *  Created by Joachim on 16/04/2019.
+ *  Adapted from donated nonius code.
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See accompanying
+ *  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// Timing
+
+#ifndef TWOBLUECUBES_CATCH_DETAIL_TIMING_HPP_INCLUDED
+#define TWOBLUECUBES_CATCH_DETAIL_TIMING_HPP_INCLUDED
+
+#include "../catch_clock.hpp"
+#include "catch_complete_invoke.hpp"
+
+#include <tuple>
+#include <type_traits>
+
+namespace Catch {
+    namespace Benchmark {
+        template <typename Duration, typename Result>
+        struct Timing {
+            Duration elapsed;
+            Result result;
+            int iterations;
+        };
+        template <typename Clock, typename Sig>
+        using TimingOf = Timing<ClockDuration<Clock>, Detail::CompleteType_t<Detail::ResultOf_t<Sig>>>;
+    } // namespace Benchmark
+} // namespace Catch
+
+#endif // TWOBLUECUBES_CATCH_DETAIL_TIMING_HPP_INCLUDED