diff --git a/CMakeLists.txt b/CMakeLists.txt
index f7a6565..9f1837d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -22,10 +22,6 @@ endif()
 # Build
 # =====
 
-set(XSIMDALGO_HEADERS
-    ${XSIMDALGO_INCLUDE_DIR}/xsimd_algorithm/algorithms.hpp
-)
-
 add_library(xsimd-algorithm INTERFACE)
 
 target_include_directories(xsimd-algorithm INTERFACE
@@ -37,12 +33,19 @@ target_compile_features(xsimd-algorithm INTERFACE cxx_std_20)
 target_link_libraries(xsimd-algorithm INTERFACE xsimd)
 
 OPTION(BUILD_TESTS "xsimd-algorithm test suite" OFF)
+OPTION(BUILD_BENCHMARK "xsimd-algorithm benchmark suite" OFF)
+
+add_subdirectory(test-utils)
 
 if(BUILD_TESTS)
     enable_testing()
     add_subdirectory(test)
 endif()
 
+if(BUILD_BENCHMARK)
+    add_subdirectory(benchmark)
+endif()
+
 # Installation
 # ============
 
diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt
new file mode 100644
index 0000000..5b70a13
--- /dev/null
+++ b/benchmark/CMakeLists.txt
@@ -0,0 +1,26 @@
+############################################################################
+# Copyright (c) xsimd-algorithm contributors                               #
+#                                                                          #
+# Distributed under the terms of the BSD 3-Clause License.                 #
+#                                                                          #
+# The full license is in the file LICENSE, distributed with this software. #
+############################################################################
+
+cmake_minimum_required(VERSION 3.8)
+
+project(xsimd-algorithm-benchmark)
+
+if (CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR)
+    find_package(xsimd-algorithm REQUIRED CONFIG)
+endif ()
+
+find_package(benchmark REQUIRED)
+
+set(XSIMD_ALGORITHM_BENCHMARKS
+    main.cpp
+    bench_math.cpp
+)
+
+add_executable(benchmark_xsimd_algorithm ${XSIMD_ALGORITHM_BENCHMARKS})
+target_link_libraries(benchmark_xsimd_algorithm
+    PRIVATE xsimd-algorithm xsimd::test-utils benchmark::benchmark)
diff --git a/benchmark/bench_math.cpp b/benchmark/bench_math.cpp
new file mode 100644
index 0000000..35694c6
--- /dev/null
+++ b/benchmark/bench_math.cpp
@@ -0,0 +1,124 @@
+/****************************************************************************
+ * Copyright (c) xsimd-algorithm contributors                               *
+ *                                                                          *
+ * Distributed under the terms of the BSD 3-Clause License.                 *
+ *                                                                          *
+ * The full license is in the file LICENSE, distributed with this software. *
+ ****************************************************************************/
+
+
+#include <cstddef>
+#include <format>
+#include <string_view>
+#include <type_traits>
+#include <vector>
+
+#include <benchmark/benchmark.h>
+
+#include "xsimd_test_utils/math_data.hpp"
+#include "xsimd_test_utils/utils.hpp"
+
+using xsimd::builder::alignment;
+
+namespace
+{
+    template <typename Op, typename Alloc, typename Apply>
+    void bench_unary(benchmark::State& state, Apply apply)
+    {
+        using value_type = typename Op::value_type;
+
+        auto const size = static_cast<std::size_t>(state.range(0));
+        auto [input, output] = Op::template make_input_output<Alloc>(size);
+
+        for (auto _ : state)
+        {
+            apply(xsimd::test::as_span(input), xsimd::test::as_span(output));
+            benchmark::DoNotOptimize(output.data());
+            benchmark::ClobberMemory();
+        }
+
+        state.SetItemsProcessed(static_cast<std::int64_t>(state.iterations() * size));
+        state.SetBytesProcessed(
+            static_cast<std::int64_t>(state.iterations() * size * 2 * sizeof(value_type)));
+    }
+
+    /// Sizes spanning L1-resident to memory-bound, each with and without a scalar tail.
+    template <typename T>
+    std::vector<std::int64_t> bench_sizes()
+    {
+        constexpr auto batch_size = static_cast<std::int64_t>(xsimd::batch<T>::size);
+
+        auto sizes = std::vector<std::int64_t> {};
+        for (std::int64_t size : { 64, 1024, 65536, 1 << 21 })
+        {
+            auto const whole = size - (size % batch_size);
+            sizes.push_back(whole);
+            sizes.push_back(whole + batch_size / 2 + 1);
+        }
+        return sizes;
+    }
+
+    template <typename T>
+    constexpr auto type_name()
+    {
+        if constexpr (std::is_same_v<T, float>)
+        {
+            return "f32";
+        }
+        else if constexpr (std::is_same_v<T, double>)
+        {
+            return "f64";
+        }
+    }
+
+    template <typename Op, typename Alloc, alignment aligned = {}>
+    void bench_simd(benchmark::State& state)
+    {
+        bench_unary<Op, Alloc>(
+            state,
+            [](auto in, auto out) { Op::template apply_range_simd<aligned>(in, out); });
+    }
+
+    template <typename Op, typename Alloc>
+    void bench_scalar(benchmark::State& state)
+    {
+        bench_unary<Op, Alloc>(state, [](auto in, auto out) { Op::apply_range_scalar(in, out); });
+    }
+
+    template <typename Op, typename Bench>
+    void register_bench(std::string_view variant, Bench bench_fn)
+    {
+        using value_type = typename Op::value_type;
+
+        auto* bench = benchmark::RegisterBenchmark(
+            std::format("{}/{}/{}", Op::name, type_name<value_type>(), variant),
+            bench_fn);
+        for (auto const size : bench_sizes<value_type>())
+        {
+            bench->Arg(size);
+        }
+    }
+
+    template <typename Op>
+    void register_benches()
+    {
+        using value_type = typename Op::value_type;
+        using aligned_alloc = typename xsimd::test::aligned_vector<value_type>::allocator_type;
+        using unaligned_alloc = typename xsimd::test::unaligned_vector<value_type>::allocator_type;
+
+        register_bench<Op>("scalar/aligned", bench_scalar<Op, aligned_alloc>);
+        register_bench<Op>("simd/aligned", bench_simd<Op, aligned_alloc, { .start_aligned = true }>);
+        register_bench<Op>("simd/unaligned", bench_simd<Op, unaligned_alloc>);
+    }
+
+    bool const registered = []
+    {
+        register_benches<xsimd::test::sqrt_op<float>>();
+        register_benches<xsimd::test::sqrt_op<double>>();
+        register_benches<xsimd::test::abs_op<float>>();
+        register_benches<xsimd::test::abs_op<double>>();
+        register_benches<xsimd::test::exp_op<float>>();
+        register_benches<xsimd::test::exp_op<double>>();
+        return true;
+    }();
+}
diff --git a/benchmark/main.cpp b/benchmark/main.cpp
new file mode 100644
index 0000000..68e814d
--- /dev/null
+++ b/benchmark/main.cpp
@@ -0,0 +1,11 @@
+/****************************************************************************
+ * Copyright (c) xsimd-algorithm contributors                               *
+ *                                                                          *
+ * Distributed under the terms of the BSD 3-Clause License.                 *
+ *                                                                          *
+ * The full license is in the file LICENSE, distributed with this software. *
+ ****************************************************************************/
+
+#include <benchmark/benchmark.h>
+
+BENCHMARK_MAIN();
diff --git a/environment-dev.yml b/environment-dev.yml
index 767ec13..bdbfb7a 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -5,4 +5,5 @@ dependencies:
 - cmake
 - xsimd=14.3.0
 - doctest
-- ninja
\ No newline at end of file
+- benchmark
+- ninja
diff --git a/include/xsimd_algorithm/builder.hpp b/include/xsimd_algorithm/builder.hpp
new file mode 100644
index 0000000..d8f45aa
--- /dev/null
+++ b/include/xsimd_algorithm/builder.hpp
@@ -0,0 +1,318 @@
+/****************************************************************************
+ * Copyright (c) xsimd-algorithm contributors                               *
+ *                                                                          *
+ * Distributed under the terms of the BSD 3-Clause License.                 *
+ *                                                                          *
+ * The full license is in the file LICENSE, distributed with this software. *
+ ****************************************************************************/
+
+#ifndef XSIMD_ALGORITHM_BUILDER_HPP
+#define XSIMD_ALGORITHM_BUILDER_HPP
+
+#include <algorithm>
+#include <array>
+#include <bit>
+#include <cassert>
+#include <cstddef>
+#include <cstdint>
+#include <cstring>
+#include <span>
+#include <utility>
+
+#include <xsimd/xsimd.hpp>
+
+#include "./macros.hpp"
+
+namespace xsimd::builder
+{
+    struct alignment
+    {
+        bool start_aligned = false;
+        bool end_aligned = false;
+    };
+
+    /// Return the pointer before the input with the given alignment or itself if aligned.
+    template <typename T>
+    XSIMD_INLINE auto prev_aligned(T* ptr, std::size_t alignment) -> T*
+    {
+        assert(std::has_single_bit(alignment));
+        auto const address = reinterpret_cast<std::size_t>(ptr);
+        return reinterpret_cast<T*>(address & ~(alignment - 1));
+    }
+
+    /// Return the pointer after the input with the given alignment or itself if aligned.
+    template <typename T>
+    XSIMD_INLINE auto next_aligned(T* ptr, std::size_t alignment) -> T*
+    {
+        assert(std::has_single_bit(alignment));
+        auto const address = reinterpret_cast<std::size_t>(ptr);
+        return reinterpret_cast<T*>((address + alignment - 1) & ~(alignment - 1));
+    }
+
+    template <typename T>
+    XSIMD_INLINE auto bytes_to_next_aligned(T* ptr, std::size_t alignment) -> std::size_t
+    {
+        assert(std::has_single_bit(alignment));
+        auto const address = reinterpret_cast<std::uintptr_t>(ptr);
+        return (alignment - (address & (alignment - 1))) & (alignment - 1);
+    }
+
+    /// Check if two spans are aliasing each others (overlapping).
+    template <typename T, typename U>
+    XSIMD_INLINE auto are_aliased(std::span<T> lhs, std::span<U> rhs) -> bool
+    {
+        // Comparing pointers from unrelated objects is unspecified, integers are not.
+        auto const lhs_begin = reinterpret_cast<std::uintptr_t>(lhs.data());
+        auto const rhs_begin = reinterpret_cast<std::uintptr_t>(rhs.data());
+        return (lhs_begin < rhs_begin + rhs.size_bytes()) && (rhs_begin < lhs_begin + lhs.size_bytes());
+    }
+
+    struct unary_options
+    {
+        std::size_t unroll_factor = 4;
+        bool pure = false;
+    };
+
+    /// Number of batches of T spanning as many elements as one batch of the widest of T and U.
+    ///
+    /// Pairing that many batches on each side lets both sides advance by the same number of
+    /// elements, so a mapping stays elementwise regardless of the respective lane counts.
+    template <typename T, typename U>
+    inline constexpr std::size_t batch_arity = sizeof(T) / std::min(sizeof(T), sizeof(U));
+
+    /// Load batch wrapper with an alignment as template parameter.
+    template <typename T, typename A, bool aligned>
+    XSIMD_INLINE xsimd::batch<T, A> load_batch(T const* ptr)
+    {
+        if constexpr (aligned)
+        {
+            return xsimd::batch<T, A>::load_aligned(ptr);
+        }
+        else
+        {
+            return xsimd::batch<T, A>::load_unaligned(ptr);
+        }
+    }
+
+    /// Store batch wrapper with an alignment as template parameter.
+    template <typename T, typename A, bool aligned>
+    XSIMD_INLINE void store_batch(xsimd::batch<T, A> x, T* ptr)
+    {
+        if constexpr (aligned)
+        {
+            x.store_aligned(ptr);
+        }
+        else
+        {
+            x.store_unaligned(ptr);
+        }
+    }
+
+    /// Load an array of batches.
+    template <std::size_t N, typename A, bool aligned, typename T>
+    XSIMD_INLINE auto load_batches(T const* ptr) -> std::array<xsimd::batch<T, A>, N>
+    {
+        std::array<xsimd::batch<T, A>, N> x;
+        for (std::size_t i = 0; i < N; ++i)
+        {
+            x[i] = load_batch<T, A, aligned>(ptr + i * xsimd::batch<T, A>::size);
+        }
+        return x;
+    }
+
+    /// Store an array of batches.
+    template <std::size_t N, typename A, bool aligned, typename T>
+    XSIMD_INLINE void store_batches(std::array<xsimd::batch<T, A>, N> const& x, T* ptr)
+    {
+        for (std::size_t i = 0; i < N; ++i)
+        {
+            store_batch<T, A, aligned>(x[i], ptr + i * xsimd::batch<T, A>::size);
+        }
+    }
+
+    namespace internal
+    {
+        template <typename T>
+        inline constexpr bool is_array = false;
+
+        template <typename T, std::size_t N>
+        inline constexpr bool is_array<std::array<T, N>> = true;
+
+        /// If an array contains only one element, return it.
+        template <typename B, std::size_t N>
+        XSIMD_INLINE auto const& unwrap_array(std::array<B, N> const& x)
+        {
+            if constexpr (N == 1)
+            {
+                return x[0];
+            }
+            else
+            {
+                return x;
+            }
+        }
+
+        /// Wrap user function to handle ``xsimd::batch`` as 1D array.
+        ///
+        /// Transform 1D input array as batch from alogrithm functions to batch for to
+        /// the user function, and user batch result as 1D arrays for the algorithm
+        /// functions.
+        template <typename Func>
+        XSIMD_INLINE auto wrap_params_as_1d_arrays(Func&& func)
+        {
+            return [func = std::forward<Func>(func)](auto const&... x)
+            {
+                auto res = func(internal::unwrap_array(x)...);
+                if constexpr (internal::is_array<decltype(res)>)
+                {
+                    return res;
+                }
+                else
+                {
+                    return std::array { res };
+                }
+            };
+        }
+    }
+
+    /// Map fewer elements than a full step through a scratch buffer.
+    template <
+        typename Arch = xsimd::default_arch,
+        typename T, typename U, typename Func>
+    XSIMD_INLINE void map_unary_batch(
+        T const* XSIMD_RESTRICT begin,
+        T const* XSIMD_RESTRICT end,
+        U* XSIMD_RESTRICT out,
+        Func&& func)
+    {
+        constexpr std::size_t in_arity = batch_arity<T, U>;
+        constexpr std::size_t out_arity = batch_arity<U, T>;
+        constexpr std::size_t step = in_arity * xsimd::batch<T, Arch>::size;
+        static_assert(step == out_arity * xsimd::batch<U, Arch>::size);
+        auto mapper = internal::wrap_params_as_1d_arrays(std::forward<Func>(func));
+
+        assert(begin <= end);
+        assert(static_cast<std::size_t>(end - begin) <= step);
+
+        if (begin == end) [[unlikely]]
+        {
+            return;
+        }
+
+        alignas(Arch::alignment()) std::array<T, step> input_buffer {};
+        alignas(Arch::alignment()) std::array<U, step> output_buffer;
+
+        const std::size_t count = static_cast<std::size_t>(end - begin);
+        std::memcpy(input_buffer.data(), begin, count * sizeof(T));
+        store_batches<out_arity, Arch, true>(
+            mapper(load_batches<in_arity, Arch, true>(input_buffer.data())),
+            output_buffer.data());
+        std::memcpy(out, output_buffer.data(), count * sizeof(U));
+    }
+
+    /// Apply func elementwise over in, writing as many elements to out.
+    ///
+    /// Func maps a std::array<batch<T, Arch>, batch_arity<T, U>> to a
+    /// std::array<batch<U, Arch>, batch_arity<U, T>>, both spanning the same element count.
+    /// When arity is one, a callback over plain batches is accepted as well.
+    template <
+        alignment align = {},
+        unary_options opts = {},
+        typename Arch = xsimd::default_arch,
+        typename T, typename U, typename Func>
+    XSIMD_INLINE void map_unary(std::span<T const> in, std::span<U> out, Func&& func)
+    {
+        constexpr std::size_t in_arity = batch_arity<T, U>;
+        constexpr std::size_t out_arity = batch_arity<U, T>;
+        // Elements consumed and produced by a single call to func.
+        constexpr std::size_t step = in_arity * xsimd::batch<T, Arch>::size;
+        static_assert(step == out_arity * xsimd::batch<U, Arch>::size);
+        auto mapper = internal::wrap_params_as_1d_arrays(std::forward<Func>(func));
+
+        // Input and output may not have the same alignment so it may be impossible
+        // to get both aligned. We align preferably the output (more expensive
+        // unaligned stores) or otherwise the input.
+        constexpr bool align_output = sizeof(U) >= sizeof(T);
+        constexpr bool load_is_aligned = align.start_aligned || !align_output;
+        constexpr bool store_is_aligned = align.start_aligned || align_output;
+
+        assert(in.size() == out.size());
+        assert(!are_aliased(in, out));
+
+        if (in.empty()) [[unlikely]]
+        {
+            return;
+        }
+
+        auto ot = out.data();
+        auto it = in.data();
+        auto const iend = in.data() + in.size();
+
+        if constexpr (!align.start_aligned)
+        {
+            // The span may be too short to reach the next alignment boundary.
+            const auto head = std::min(
+                align_output ? bytes_to_next_aligned(ot, Arch::alignment()) / sizeof(U)
+                             : bytes_to_next_aligned(it, Arch::alignment()) / sizeof(T),
+                in.size());
+
+            if (opts.pure && (head != 0) && (in.size() >= step))
+            {
+                // Recompute the head as a full step, the body overwrites the excess.
+                store_batches<out_arity, Arch, false>(
+                    mapper(load_batches<in_arity, Arch, false>(it)),
+                    ot);
+            }
+            else
+            {
+                map_unary_batch<Arch>(it, it + head, ot, func);
+            }
+            it += head;
+            ot += head;
+        }
+
+        // Unrolled loop processing multiple steps at a time
+        while (static_cast<std::size_t>(iend - it) >= opts.unroll_factor * step)
+        {
+            std::array<std::array<xsimd::batch<T, Arch>, in_arity>, opts.unroll_factor> x;
+            for (std::size_t u = 0; u < opts.unroll_factor; ++u)
+            {
+                x[u] = load_batches<in_arity, Arch, load_is_aligned>(it + u * step);
+            }
+            for (std::size_t u = 0; u < opts.unroll_factor; ++u)
+            {
+                store_batches<out_arity, Arch, store_is_aligned>(mapper(x[u]), ot + u * step);
+            }
+
+            it += opts.unroll_factor * step;
+            ot += opts.unroll_factor * step;
+        }
+
+        while (static_cast<std::size_t>(iend - it) >= step)
+        {
+            auto x = load_batches<in_arity, Arch, load_is_aligned>(it);
+            store_batches<out_arity, Arch, store_is_aligned>(mapper(x), ot);
+            it += step;
+            ot += step;
+        }
+
+        // Unlikely to be skipped, meant for users that know they allocate
+        // a multiple of the batch size, such as in a local buffer
+        if constexpr (!align.end_aligned)
+        {
+            auto const oend = out.data() + out.size();
+            if (opts.pure && (it != iend) && (in.size() >= step)) [[likely]]
+            {
+                // Recompute overlapping data, this time starting from the end.
+                auto x = load_batches<in_arity, Arch, false>(iend - step);
+                store_batches<out_arity, Arch, false>(mapper(x), oend - step);
+            }
+            else
+            {
+                map_unary_batch<Arch>(it, iend, ot, func);
+            }
+        }
+    }
+}
+
+#endif
diff --git a/include/xsimd_algorithm/macros.hpp b/include/xsimd_algorithm/macros.hpp
new file mode 100644
index 0000000..341f898
--- /dev/null
+++ b/include/xsimd_algorithm/macros.hpp
@@ -0,0 +1,20 @@
+/****************************************************************************
+ * Copyright (c) xsimd-algorithm contributors                               *
+ *                                                                          *
+ * Distributed under the terms of the BSD 3-Clause License.                 *
+ *                                                                          *
+ * The full license is in the file LICENSE, distributed with this software. *
+ ****************************************************************************/
+
+#ifndef XSIMD_ALGORITHM_MACRO_HPP
+#define XSIMD_ALGORITHM_MACRO_HPP
+
+#if defined(_MSC_VER) && !defined(__clang__)
+#define XSIMD_RESTRICT __restrict
+#elif defined(__GNUC__) || defined(__clang__)
+#define XSIMD_RESTRICT __restrict__
+#else
+#define XSIMD_RESTRICT
+#endif
+
+#endif
diff --git a/include/xsimd_algorithm/math.hpp b/include/xsimd_algorithm/math.hpp
new file mode 100644
index 0000000..139016a
--- /dev/null
+++ b/include/xsimd_algorithm/math.hpp
@@ -0,0 +1,55 @@
+/****************************************************************************
+ * Copyright (c) xsimd-algorithm contributors                               *
+ *                                                                          *
+ * Distributed under the terms of the BSD 3-Clause License.                 *
+ *                                                                          *
+ * The full license is in the file LICENSE, distributed with this software. *
+ ****************************************************************************/
+
+#ifndef XSIMD_ALGORITHM_MATH_HPP
+#define XSIMD_ALGORITHM_MATH_HPP
+
+#include <xsimd/xsimd.hpp>
+
+#include "./builder.hpp"
+
+namespace xsimd::algo
+{
+    template <
+        xsimd::builder::alignment align = {},
+        typename Arch = xsimd::default_arch,
+        typename T>
+    void sqrt(std::span<T const> in, std::span<T> out)
+    {
+        constexpr builder::unary_options opts = { .unroll_factor = 4, .pure = true };
+        return xsimd::builder::map_unary<align, opts, Arch>(
+            in, out, [](auto x)
+            { return sqrt(x); });
+    }
+
+    template <
+        xsimd::builder::alignment align = {},
+        typename Arch = xsimd::default_arch,
+        typename T>
+    void abs(std::span<T const> in, std::span<T> out)
+    {
+        constexpr builder::unary_options opts = { .unroll_factor = 4, .pure = true };
+        return xsimd::builder::map_unary<align, opts, Arch>(
+            in, out, [](auto x)
+            { return abs(x); });
+    }
+
+    template <
+        xsimd::builder::alignment align = {},
+        typename Arch = xsimd::default_arch,
+        typename T>
+    void exp(std::span<T const> in, std::span<T> out)
+    {
+        constexpr builder::unary_options opts = { .unroll_factor = 4, .pure = true };
+        return xsimd::builder::map_unary<align, opts, Arch>(
+            in, out, [](auto x)
+            { return exp(x); });
+    }
+}
+
+#endif
diff --git a/test-utils/CMakeLists.txt b/test-utils/CMakeLists.txt
new file mode 100644
index 0000000..7f832e4
--- /dev/null
+++ b/test-utils/CMakeLists.txt
@@ -0,0 +1,23 @@
+############################################################################
+# Copyright (c) xsimd-algorithm contributors                               #
+#                                                                          #
+# Distributed under the terms of the BSD 3-Clause License.                 #
+#                                                                          #
+# The full license is in the file LICENSE, distributed with this software. #
+############################################################################
+
+cmake_minimum_required(VERSION 3.8)
+
+project(xsimd-algorithm-test-utils)
+
+if (CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR)
+    find_package(xsimd-algorithm REQUIRED CONFIG)
+endif ()
+
+add_library(xsimd-algorithm-test-utils INTERFACE)
+add_library(xsimd::test-utils ALIAS xsimd-algorithm-test-utils)
+
+target_include_directories(xsimd-algorithm-test-utils INTERFACE
+    ${CMAKE_CURRENT_SOURCE_DIR}/include)
+
+target_link_libraries(xsimd-algorithm-test-utils INTERFACE xsimd-algorithm)
diff --git a/test-utils/include/xsimd_test_utils/math_data.hpp b/test-utils/include/xsimd_test_utils/math_data.hpp
new file mode 100644
index 0000000..736fb7e
--- /dev/null
+++ b/test-utils/include/xsimd_test_utils/math_data.hpp
@@ -0,0 +1,128 @@
+/****************************************************************************
+ * Copyright (c) xsimd-algorithm contributors                               *
+ *                                                                          *
+ * Distributed under the terms of the BSD 3-Clause License.                 *
+ *                                                                          *
+ * The full license is in the file LICENSE, distributed with this software. *
+ ****************************************************************************/
+
+#ifndef XSIMD_ALGORITHM_TEST_UTILS_MATH_DATA_HPP
+#define XSIMD_ALGORITHM_TEST_UTILS_MATH_DATA_HPP
+
+#include <cmath>
+#include <cstddef>
+#include <span>
+#include <utility>
+#include <vector>
+
+#include "xsimd_algorithm/builder.hpp"
+#include "xsimd_algorithm/math.hpp"
+
+#include "xsimd_test_utils/utils.hpp"
+
+namespace xsimd::test
+{
+    /// Derives the scalar range application from the element-wise Derived::apply.
+    template <typename Derived, typename T>
+    struct unary_op
+    {
+        using value_type = T;
+
+        static void apply_range_scalar(std::span<T const> in, std::span<T> out)
+        {
+            for (std::size_t i = 0; i < in.size(); ++i)
+            {
+                out[i] = Derived::apply(in[i]);
+            }
+        }
+
+        template <typename Alloc>
+        static std::pair<std::vector<T, Alloc>, std::vector<T, Alloc>> make_input_output(std::size_t size)
+        {
+            auto input = Derived::template make_input<Alloc>(size);
+            auto output = std::vector<T, Alloc>(input.size());
+            return { std::move(input), std::move(output) };
+        }
+    };
+
+    /*******************
+     *  Test fixtures  *
+     *******************/
+
+    template <typename T>
+    struct sqrt_op : unary_op<sqrt_op<T>, T>
+    {
+        static constexpr auto name = "sqrt";
+
+        static T apply(T x)
+        {
+            return std::sqrt(x);
+        }
+
+        template <xsimd::builder::alignment aligned = {}>
+        static void apply_range_simd(std::span<T const> in, std::span<T> out)
+        {
+            xsimd::algo::sqrt<aligned>(in, out);
+        }
+
+        template <typename Alloc>
+        static std::vector<T, Alloc> make_input(std::size_t size)
+        {
+            return make_arange<T, Alloc>(size);
+        }
+    };
+
+    template <typename T>
+    struct abs_op : unary_op<abs_op<T>, T>
+    {
+        static constexpr auto name = "abs";
+
+        static T apply(T x)
+        {
+            return std::abs(x);
+        }
+
+        template <xsimd::builder::alignment aligned = {}>
+        static void apply_range_simd(std::span<T const> in, std::span<T> out)
+        {
+            xsimd::algo::abs<aligned>(in, out);
+        }
+
+        template <typename Alloc>
+        static std::vector<T, Alloc> make_input(std::size_t size)
+        {
+            return make_arange<T, Alloc>(size, -static_cast<T>(size) / 2);
+        }
+    };
+
+    template <typename T>
+    struct exp_op : unary_op<exp_op<T>, T>
+    {
+        static constexpr auto name = "exp";
+
+        static T apply(T x)
+        {
+            return std::exp(x);
+        }
+
+        template <xsimd::builder::alignment aligned = {}>
+        static void apply_range_simd(std::span<T const> in, std::span<T> out)
+        {
+            xsimd::algo::exp<aligned>(in, out);
+        }
+
+        template <typename Alloc>
+        static std::vector<T, Alloc> make_input(std::size_t size)
+        {
+            // exp overflows past a small range, so wrap the values back into [-10, 10).
+            auto input = make_arange<T, Alloc>(size);
+            for (auto& x : input)
+            {
+                x = std::fmod(x, T { 20 }) - T { 10 };
+            }
+            return input;
+        }
+    };
+}
+
+#endif
diff --git a/test-utils/include/xsimd_test_utils/math_ops.hpp b/test-utils/include/xsimd_test_utils/math_ops.hpp
new file mode 100644
index 0000000..497df13
--- /dev/null
+++ b/test-utils/include/xsimd_test_utils/math_ops.hpp
@@ -0,0 +1,75 @@
+/****************************************************************************
+ * Copyright (c) xsimd-algorithm contributors                               *
+ *                                                                          *
+ * Distributed under the terms of the BSD 3-Clause License.                 *
+ *                                                                          *
+ * The full license is in the file LICENSE, distributed with this software. *
+ ****************************************************************************/
+
+#ifndef XSIMD_ALGORITHM_TEST_UTILS_MATH_OPS_HPP
+#define XSIMD_ALGORITHM_TEST_UTILS_MATH_OPS_HPP
+
+#include <cmath>
+#include <cstddef>
+#include <span>
+#include <vector>
+
+#include <xsimd_algorithm/builder.hpp>
+#include <xsimd_algorithm/math.hpp>
+
+#include "xsimd_test_utils/utils.hpp"
+
+namespace xsimd::test
+{
+    template <typename T>
+    struct sqrt_op
+    {
+        using value_type = T;
+
+        static constexpr auto name = "sqrt";
+
+        template <xsimd::builder::alignment aligned = {}>
+        static void apply(std::span<T const> in, std::span<T> out)
+        {
+            xsimd::algo::sqrt<aligned>(in, out);
+        }
+
+        static T scalar(T x)
+        {
+            return std::sqrt(x);
+        }
+
+        template <typename Alloc>
+        static std::vector<T, Alloc> input(std::size_t size)
+        {
+            return xsimd::test::make_arange<T, Alloc>(size);
+        }
+    };
+
+    template <typename T>
+    struct abs_op
+    {
+        using value_type = T;
+
+        static constexpr auto name = "abs";
+
+        template <xsimd::builder::alignment aligned = {}>
+        static void apply(std::span<T const> in, std::span<T> out)
+        {
+            xsimd::algo::abs<aligned>(in, out);
+        }
+
+        static T scalar(T x)
+        {
+            return std::abs(x);
+        }
+
+        template <typename Alloc>
+        static std::vector<T, Alloc> input(std::size_t size)
+        {
+            return xsimd::test::make_arange<T, Alloc>(size, -static_cast<T>(size) / 2);
+        }
+    };
+}
+
+#endif
diff --git a/test-utils/include/xsimd_test_utils/utils.hpp b/test-utils/include/xsimd_test_utils/utils.hpp
new file mode 100644
index 0000000..692885b
--- /dev/null
+++ b/test-utils/include/xsimd_test_utils/utils.hpp
@@ -0,0 +1,81 @@
+/****************************************************************************
+ * Copyright (c) xsimd-algorithm contributors                               *
+ *                                                                          *
+ * Distributed under the terms of the BSD 3-Clause License.                 *
+ *                                                                          *
+ * The full license is in the file LICENSE, distributed with this software. *
+ ****************************************************************************/
+
+#ifndef XSIMD_ALGORITHM_TEST_UTILS_UTILS_HPP
+#define XSIMD_ALGORITHM_TEST_UTILS_UTILS_HPP
+
+#include <cstddef>
+#include <numeric>
+#include <span>
+#include <vector>
+
+#include <xsimd/xsimd.hpp>
+
+namespace xsimd::test
+{
+    template <typename T, typename A = xsimd::default_arch>
+    using aligned_vector = std::vector<T, xsimd::aligned_allocator<T, A::alignment()>>;
+
+    /// An allocator returning memory guaranteed not to be aligned on @p Align.
+    ///
+    /// Over-allocates by @p Offset elements and shifts the returned pointer, so that the data
+    /// starts @p Offset * sizeof(T) bytes past an aligned address.
+    template <typename T, std::size_t Align = xsimd::default_arch::alignment(), std::size_t Offset = 1>
+    struct unaligned_allocator : private xsimd::aligned_allocator<T, Align>
+    {
+        static_assert((Offset * sizeof(T)) % Align != 0, "Shifted pointer would still be aligned");
+
+        using base_type = xsimd::aligned_allocator<T, Align>;
+        using value_type = T;
+
+        // The non-type Align parameter defeats the default allocator_traits rebind.
+        template <typename U>
+        struct rebind
+        {
+            using other = unaligned_allocator<U, Align, Offset>;
+        };
+
+        unaligned_allocator() = default;
+
+        template <typename U>
+        unaligned_allocator(unaligned_allocator<U, Align, Offset> const&)
+        {
+        }
+
+        T* allocate(std::size_t n) { return base_type::allocate(n + Offset) + Offset; }
+
+        void deallocate(T* p, std::size_t n) { base_type::deallocate(p - Offset, n + Offset); }
+
+        friend bool operator==(unaligned_allocator const&, unaligned_allocator const&) { return true; }
+    };
+
+    template <typename T, typename A = xsimd::default_arch>
+    using unaligned_vector = std::vector<T, unaligned_allocator<T, A::alignment()>>;
+
+    template <typename T, typename A>
+    std::span<T const> as_span(std::vector<T, A> const& v)
+    {
+        return std::span<T const> { v.data(), v.size() };
+    }
+
+    template <typename T, typename A>
+    std::span<T> as_span(std::vector<T, A>& v)
+    {
+        return std::span<T> { v.data(), v.size() };
+    }
+
+    template <typename T, typename Alloc = typename aligned_vector<T>::allocator_type>
+    std::vector<T, Alloc> make_arange(std::size_t size, T start = T { 0 })
+    {
+        std::vector<T, Alloc> data(size);
+        std::iota(data.begin(), data.end(), start);
+        return data;
+    }
+}
+
+#endif
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 2285b23..a6d7e87 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -48,13 +48,15 @@ endif()
 set(XSIMD_ALGORITHM_TESTS
     main.cpp
     test_arange.cpp
+    test_builder.cpp
     test_iterator.cpp
+    test_math.cpp
     test_reduce.cpp
     test_transform.cpp
 )
 
-add_executable(test_xsimd_algorithm ${XSIMD_ALGORITHM_TESTS})# ${XSIMD_ALGORITHM_HEADERS})
-target_link_libraries(test_xsimd_algorithm PRIVATE xsimd-algorithm)
+add_executable(test_xsimd_algorithm ${XSIMD_ALGORITHM_TESTS})
+target_link_libraries(test_xsimd_algorithm PRIVATE xsimd-algorithm xsimd::test-utils)
 
 option(DOWNLOAD_DOCTEST OFF)
 find_package(doctest QUIET)
diff --git a/test/test_builder.cpp b/test/test_builder.cpp
new file mode 100644
index 0000000..5088ebd
--- /dev/null
+++ b/test/test_builder.cpp
@@ -0,0 +1,81 @@
+/****************************************************************************
+ * Copyright (c) xsimd-algorithm contributors                               *
+ *                                                                          *
+ * Distributed under the terms of the BSD 3-Clause License.                 *
+ *                                                                          *
+ * The full license is in the file LICENSE, distributed with this software. *
+ ****************************************************************************/
+
+#include <array>
+#include <cstddef>
+#include <cstdint>
+
+#include <doctest/doctest.h>
+#include <xsimd_test_utils/utils.hpp>
+
+#include "xsimd_algorithm/builder.hpp"
+
+/// Map unary test where one input batch pairs with two output batches.
+TEST_CASE("map_unary int32 to int64")
+{
+    using input_type = std::int32_t;
+    using output_type = std::int64_t;
+
+    // Not a multiple of the batch size, to exercise the tail.
+    static constexpr std::size_t size = 94;
+
+    const auto input = xsimd::test::make_arange<input_type>(size);
+    auto output = xsimd::test::aligned_vector<output_type>(size);
+
+    const auto func = [](auto const& x)
+    { return xsimd::widen(x + input_type { 1 }); };
+
+    xsimd::builder::map_unary(xsimd::test::as_span(input), xsimd::test::as_span(output), func);
+
+    for (std::size_t i = 0; i < size; ++i)
+    {
+        CAPTURE(i);
+        CHECK(output[i] == static_cast<output_type>(input[i] + 1));
+    }
+}
+
+/// Map unary test where two input batches pair with one output batch.
+TEST_CASE("map_unary int64 to int32")
+{
+    using input_type = std::int64_t;
+    using output_type = std::int32_t;
+    using input_batch = xsimd::batch<input_type>;
+    using output_batch = xsimd::batch<output_type>;
+
+    // Not a multiple of the batch size, to exercise the tail.
+    static constexpr std::size_t size = 94;
+
+    const auto input = xsimd::test::make_arange<input_type>(size);
+    auto output = xsimd::test::aligned_vector<output_type>(size);
+
+    // xsimd has no narrowing counterpart to widen, so truncate by keeping the low
+    // half of each lane, those of the first batch followed by those of the second.
+    struct low_halves
+    {
+        static constexpr unsigned get(unsigned i, unsigned n)
+        {
+            return (i < n / 2) ? 2 * i : n + 2 * (i - n / 2);
+        }
+    };
+
+    const auto func = [](std::array<input_batch, 2> const& x) -> output_batch
+    {
+        return xsimd::shuffle(
+            xsimd::bitwise_cast<output_type>(x[0] + input_type { 1 }),
+            xsimd::bitwise_cast<output_type>(x[1] + input_type { 1 }),
+            xsimd::make_batch_constant<std::uint32_t, low_halves>());
+    };
+
+    xsimd::builder::map_unary(xsimd::test::as_span(input), xsimd::test::as_span(output), func);
+
+    for (std::size_t i = 0; i < size; ++i)
+    {
+        CAPTURE(i);
+        CHECK(output[i] == static_cast<output_type>(input[i] + 1));
+    }
+}
diff --git a/test/test_math.cpp b/test/test_math.cpp
new file mode 100644
index 0000000..fa995b8
--- /dev/null
+++ b/test/test_math.cpp
@@ -0,0 +1,64 @@
+/****************************************************************************
+ * Copyright (c) xsimd-algorithm contributors                               *
+ *                                                                          *
+ * Distributed under the terms of the BSD 3-Clause License.                 *
+ *                                                                          *
+ * The full license is in the file LICENSE, distributed with this software. *
+ ****************************************************************************/
+
+#include <cstddef>
+
+#include <doctest/doctest.h>
+
+#include <xsimd_test_utils/math_data.hpp>
+#include <xsimd_test_utils/utils.hpp>
+
+namespace
+{
+    template <typename Op, typename Alloc, xsimd::builder::alignment aligned = {}>
+    void check_unary_math()
+    {
+        // Not a multiple of the batch size, to exercise the tail.
+        constexpr std::size_t test_size = 94;
+
+        auto [input, output] = Op::template make_input_output<Alloc>(test_size);
+
+        Op::template apply_range_simd<aligned>(xsimd::test::as_span(input), xsimd::test::as_span(output));
+
+        for (std::size_t i = 0; i < input.size(); ++i)
+        {
+            CAPTURE(i);
+            CHECK(output[i] == doctest::Approx(Op::apply(input[i])));
+        }
+    }
+}
+
+TEST_CASE_TEMPLATE(
+    "unary math",
+    Op,
+    xsimd::test::sqrt_op<float>,
+    xsimd::test::sqrt_op<double>,
+    xsimd::test::abs_op<float>,
+    xsimd::test::abs_op<double>,
+    xsimd::test::exp_op<float>,
+    xsimd::test::exp_op<double>)
+{
+    using value_type = typename Op::value_type;
+    using aligned_allocator = typename xsimd::test::aligned_vector<value_type>::allocator_type;
+    using unaligned_allocator = typename xsimd::test::unaligned_vector<value_type>::allocator_type;
+
+    SUBCASE("aligned without header")
+    {
+        check_unary_math<Op, aligned_allocator, { .start_aligned = true }>();
+    }
+
+    SUBCASE("aligned with header")
+    {
+        check_unary_math<Op, aligned_allocator>();
+    }
+
+    SUBCASE("unaligned with header")
+    {
+        check_unary_math<Op, unaligned_allocator>();
+    }
+}