diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml
index 5df542f..29013c7 100644
--- a/.github/workflows/linux.yml
+++ b/.github/workflows/linux.yml
@@ -17,7 +17,6 @@ jobs:
           - { compiler: 'gcc',   version: '13',  flags: 'enable_xtl_complex' }
           - { compiler: 'gcc',   version: '14',  flags: 'avx' }
           - { compiler: 'gcc',   version: '13', flags: 'avx512' }
-          - { compiler: 'gcc',   version: '12', flags: 'i386' }
           - { compiler: 'gcc',   version: '13', flags: 'avx512pf' }
           - { compiler: 'gcc',   version: '13', flags: 'avx512vbmi' }
           - { compiler: 'gcc',   version: '14', flags: 'avx512vbmi2' }
@@ -34,10 +33,6 @@ jobs:
         GCC_VERSION=${{ matrix.sys.version }}
         sudo apt-get update
         sudo apt-get --no-install-suggests --no-install-recommends install g++-$GCC_VERSION
-        sudo dpkg --add-architecture i386
-        sudo add-apt-repository ppa:ubuntu-toolchain-r/test
-        sudo apt-get update
-        sudo apt-get --no-install-suggests --no-install-recommends install gcc-$GCC_VERSION-multilib g++-$GCC_VERSION-multilib linux-libc-dev:i386
         CC=gcc-$GCC_VERSION
         echo "CC=$CC" >> $GITHUB_ENV
         CXX=g++-$GCC_VERSION
@@ -95,9 +90,6 @@ jobs:
         if [[ '${{ matrix.sys.flags }}' == 'avx512vnni' ]]; then
           CMAKE_EXTRA_ARGS="$CMAKE_EXTRA_ARGS -DTARGET_ARCH=knm"
         fi
-        if [[ '${{ matrix.sys.flags }}' == 'i386' ]]; then
-          CXX_FLAGS="$CXX_FLAGS -m32"
-        fi
         if [[ '${{ matrix.sys.flags }}' == 'force_no_instr_set' ]]; then
           :
         else
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..e214e36
--- /dev/null
+++ b/benchmark/CMakeLists.txt
@@ -0,0 +1,27 @@
+############################################################################
+# 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
+    bench_map.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_map.cpp b/benchmark/bench_map.cpp
new file mode 100644
index 0000000..562c53a
--- /dev/null
+++ b/benchmark/bench_map.cpp
@@ -0,0 +1,158 @@
+/****************************************************************************
+ * 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 <cstdint>
+
+#include <benchmark/benchmark.h>
+
+#include "map_binary_utils.hpp"
+#include "map_unary_utils.hpp"
+#include "xsimd_algorithm/map.hpp"
+#include "xsimd_test_utils/map_binary_data.hpp"
+#include "xsimd_test_utils/map_unary_data.hpp"
+#include "xsimd_test_utils/utils.hpp"
+
+namespace
+{
+    using xsimd::bench::bench_binary_scalar;
+    using xsimd::bench::bench_map_binary;
+    using xsimd::bench::bench_map_unary;
+    using xsimd::bench::bench_scalar;
+    using xsimd::bench::bench_transform;
+    using xsimd::bench::register_binary_bench;
+    using xsimd::bench::register_bench;
+    using xsimd::alignment_options;
+    using xsimd::map_options;
+
+    template <typename Op>
+    void register_benches()
+    {
+        using input_t = typename Op::input_t;
+        using arch = xsimd::default_arch;
+        using aligned_alloc = typename xsimd::test::aligned_vector<input_t, arch>::allocator_type;
+        using unaligned_alloc = typename xsimd::test::unaligned_vector<input_t, arch>::allocator_type;
+
+        constexpr auto noalign = alignment_options {};
+        constexpr auto noopts = map_options { .unroll_factor = 1, .pure = false };
+
+        register_bench<Op, arch>("aligned/scalar", bench_scalar<Op, aligned_alloc>);
+        register_bench<Op, arch>("aligned/simd/transform", bench_transform<Op, aligned_alloc, arch>);
+        register_bench<Op, arch>(
+            "aligned/simd/map",
+            bench_map_unary<Op, aligned_alloc, arch, noalign, noopts>);
+        register_bench<Op, arch>(
+            "aligned/simd/map:pure",
+            bench_map_unary<
+                Op, aligned_alloc, arch,
+                noalign, map_options { .unroll_factor = 1, .pure = true }>);
+        register_bench<Op, arch>(
+            "aligned/simd/map:unroll4",
+            bench_map_unary<
+                Op, aligned_alloc, arch, noalign, map_options { .unroll_factor = 4 }>);
+        register_bench<Op, arch>(
+            "aligned/simd/map:noheader",
+            bench_map_unary<
+                Op, aligned_alloc, arch,
+                alignment_options { .start_aligned = true }, noopts>);
+        register_bench<Op, arch>(
+            "aligned/simd/map:noheader+pure+unroll4",
+            bench_map_unary<
+                Op, aligned_alloc, arch,
+                alignment_options { .start_aligned = true }, map_options { .unroll_factor = 4, .pure = true }>);
+
+        register_bench<Op, arch>("unaligned/scalar", bench_scalar<Op, unaligned_alloc>);
+        register_bench<Op, arch>("unaligned/simd/transform", bench_transform<Op, unaligned_alloc, arch>);
+        register_bench<Op, arch>(
+            "unaligned/simd/map",
+            bench_map_unary<Op, unaligned_alloc, arch, noalign, noopts>);
+        register_bench<Op, arch>(
+            "unaligned/simd/map:pure",
+            bench_map_unary<
+                Op, unaligned_alloc, arch,
+                noalign, map_options { .unroll_factor = 1, .pure = true }>);
+        register_bench<Op, arch>(
+            "unaligned/simd/map:unroll4",
+            bench_map_unary<
+                Op, unaligned_alloc, arch, noalign, map_options { .unroll_factor = 4 }>);
+        register_bench<Op, arch>(
+            "unaligned/simd/map:pure+unroll4",
+            bench_map_unary<
+                Op, unaligned_alloc, arch,
+                noalign, map_options { .unroll_factor = 4, .pure = true }>);
+    }
+
+    template <typename Op>
+    void register_binary_benches()
+    {
+        using lhs_t = typename Op::lhs_t;
+        using arch = xsimd::default_arch;
+        using aligned_alloc = typename xsimd::test::aligned_vector<lhs_t, arch>::allocator_type;
+        using unaligned_alloc = typename xsimd::test::unaligned_vector<lhs_t, arch>::allocator_type;
+
+        constexpr auto noalign = alignment_options {};
+        constexpr auto noopts = map_options { .unroll_factor = 1, .pure = false };
+
+        register_binary_bench<Op, arch>("aligned/scalar", bench_binary_scalar<Op, aligned_alloc>);
+        register_binary_bench<Op, arch>(
+            "aligned/simd/map",
+            bench_map_binary<Op, aligned_alloc, arch, noalign, noopts>);
+        register_binary_bench<Op, arch>(
+            "aligned/simd/map:pure",
+            bench_map_binary<
+                Op, aligned_alloc, arch,
+                noalign, map_options { .unroll_factor = 1, .pure = true }>);
+        register_binary_bench<Op, arch>(
+            "aligned/simd/map:unroll4",
+            bench_map_binary<
+                Op, aligned_alloc, arch, noalign, map_options { .unroll_factor = 4 }>);
+        register_binary_bench<Op, arch>(
+            "aligned/simd/map:noheader",
+            bench_map_binary<
+                Op, aligned_alloc, arch,
+                alignment_options { .start_aligned = true }, noopts>);
+        register_binary_bench<Op, arch>(
+            "aligned/simd/map:noheader+pure+unroll4",
+            bench_map_binary<
+                Op, aligned_alloc, arch,
+                alignment_options { .start_aligned = true }, map_options { .unroll_factor = 4, .pure = true }>);
+
+        register_binary_bench<Op, arch>("unaligned/scalar", bench_binary_scalar<Op, unaligned_alloc>);
+        register_binary_bench<Op, arch>(
+            "unaligned/simd/map",
+            bench_map_binary<Op, unaligned_alloc, arch, noalign, noopts>);
+        register_binary_bench<Op, arch>(
+            "unaligned/simd/map:pure",
+            bench_map_binary<
+                Op, unaligned_alloc, arch,
+                noalign, map_options { .unroll_factor = 1, .pure = true }>);
+        register_binary_bench<Op, arch>(
+            "unaligned/simd/map:unroll4",
+            bench_map_binary<
+                Op, unaligned_alloc, arch, noalign, map_options { .unroll_factor = 4 }>);
+        register_binary_bench<Op, arch>(
+            "unaligned/simd/map:pure+unroll4",
+            bench_map_binary<
+                Op, unaligned_alloc, arch,
+                noalign, map_options { .unroll_factor = 4, .pure = true }>);
+    }
+
+    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>>();
+        register_binary_benches<xsimd::test::add_op<std::int32_t>>();
+        register_binary_benches<xsimd::test::multiply_op<float>>();
+        register_binary_benches<xsimd::test::mixed_multiply_op>();
+        return true;
+    }();
+}
diff --git a/benchmark/bench_math.cpp b/benchmark/bench_math.cpp
new file mode 100644
index 0000000..b5d4d21
--- /dev/null
+++ b/benchmark/bench_math.cpp
@@ -0,0 +1,84 @@
+/****************************************************************************
+ * 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 <benchmark/benchmark.h>
+
+#include "map_binary_utils.hpp"
+#include "map_unary_utils.hpp"
+#include "xsimd_test_utils/map_binary_data.hpp"
+#include "xsimd_test_utils/map_unary_data.hpp"
+#include "xsimd_test_utils/utils.hpp"
+
+namespace
+{
+    using xsimd::bench::bench_binary_scalar;
+    using xsimd::bench::bench_map_binary;
+    using xsimd::bench::bench_map_unary;
+    using xsimd::bench::bench_scalar;
+    using xsimd::bench::register_binary_bench;
+    using xsimd::bench::register_bench;
+    using xsimd::alignment_options;
+
+    /// Register math benchmarks.
+    ///
+    /// To avoid an explosion of benchmarks, we only add a simple aligned benchmark.
+    /// This will let us know the performance of the xsimd wrappers.
+    /// See bench_map for benchmarks on the different flavor of mapping, alignment,
+    /// headers and trailers.
+    ///
+    /// This benchmark aims to test raw the performance of intrinsic, unrelated to
+    /// how they are iterated on (alignment, memory etc). To do so, they aim to stay
+    /// in L1 cache.
+    template <typename Op>
+    void register_benches()
+    {
+        using input_t = typename Op::input_t;
+        using arch = xsimd::default_arch;
+        using aligned_alloc = typename xsimd::test::aligned_vector<input_t, arch>::allocator_type;
+
+        register_bench<Op, arch>(
+            "hot/scalar", bench_scalar<Op, aligned_alloc>, /* sizes = */ { 1024 });
+        register_bench<Op, arch>(
+            "hot/simd",
+            bench_map_unary<Op, aligned_alloc, arch, alignment_options { .start_aligned = true, .end_aligned = true }>,
+            /* sizes = */ { 1024 });
+    }
+
+    template <typename Op>
+    void register_binary_benches()
+    {
+        using lhs_t = typename Op::lhs_t;
+        using arch = xsimd::default_arch;
+        using aligned_alloc = typename xsimd::test::aligned_vector<lhs_t, arch>::allocator_type;
+
+        register_binary_bench<Op, arch>(
+            "hot/scalar", bench_binary_scalar<Op, aligned_alloc>, /* sizes = */ { 1024 });
+        register_binary_bench<Op, arch>(
+            "hot/simd",
+            bench_map_binary<Op, aligned_alloc, arch, alignment_options { .start_aligned = true, .end_aligned = true }>,
+            /* sizes = */ { 1024 });
+    }
+
+    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>>();
+        register_benches<xsimd::test::widen_op<std::int8_t>>();
+        register_benches<xsimd::test::widen_op<std::int16_t>>();
+        register_benches<xsimd::test::widen_op<std::int32_t>>();
+        register_binary_benches<xsimd::test::add_op<std::int32_t>>();
+        register_binary_benches<xsimd::test::multiply_op<float>>();
+        return true;
+    }();
+}
diff --git a/benchmark/bench_utils.hpp b/benchmark/bench_utils.hpp
new file mode 100644
index 0000000..f7b3686
--- /dev/null
+++ b/benchmark/bench_utils.hpp
@@ -0,0 +1,91 @@
+/****************************************************************************
+ * 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_BENCHMARK_BENCH_UTILS_HPP
+#define XSIMD_ALGORITHM_BENCHMARK_BENCH_UTILS_HPP
+
+#include <cstdint>
+#include <type_traits>
+#include <vector>
+
+#include <xsimd/xsimd.hpp>
+
+namespace xsimd::bench
+{
+    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()
+    {
+        // Avoid failing compiler check such as std::is_same<std::uint64_t, long long> due
+        // to it not being an alias.
+        constexpr bool is_int = std::is_integral_v<T> && !std::is_same_v<T, bool>;
+        constexpr bool is_sint = is_int && std::is_signed_v<T>;
+        constexpr bool is_uint = is_int && std::is_unsigned_v<T>;
+
+        if constexpr (std::is_same_v<T, bool>)
+        {
+            return "bool";
+        }
+        else if constexpr (is_sint && sizeof(T) == 1)
+        {
+            return "i8";
+        }
+        else if constexpr (is_uint && sizeof(T) == 1)
+        {
+            return "u8";
+        }
+        else if constexpr (is_sint && sizeof(T) == 2)
+        {
+            return "i16";
+        }
+        else if constexpr (is_uint && sizeof(T) == 2)
+        {
+            return "u16";
+        }
+        else if constexpr (is_sint && sizeof(T) == 4)
+        {
+            return "i32";
+        }
+        else if constexpr (is_uint && sizeof(T) == 4)
+        {
+            return "u32";
+        }
+        else if constexpr (is_sint && sizeof(T) == 8)
+        {
+            return "i64";
+        }
+        else if constexpr (is_uint && sizeof(T) == 8)
+        {
+            return "u64";
+        }
+        else if constexpr (std::is_same_v<T, float>)
+        {
+            return "f32";
+        }
+        else if constexpr (std::is_same_v<T, double>)
+        {
+            return "f64";
+        }
+    }
+}
+
+#endif
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/benchmark/map_unary_utils.hpp b/benchmark/map_unary_utils.hpp
new file mode 100644
index 0000000..68b54c4
--- /dev/null
+++ b/benchmark/map_unary_utils.hpp
@@ -0,0 +1,92 @@
+/****************************************************************************
+ * 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 <cstdint>
+#include <format>
+#include <string_view>
+#include <vector>
+
+#include <benchmark/benchmark.h>
+
+#include "bench_utils.hpp"
+#include "xsimd_algorithm/map.hpp"
+#include "xsimd_test_utils/utils.hpp"
+
+namespace xsimd::bench
+{
+    using xsimd::alignment_options;
+    using xsimd::map_options;
+
+    template <typename Op, typename Alloc, typename Apply>
+    void bench_unary(benchmark::State& state, Apply apply)
+    {
+        using input_t = typename Op::input_t;
+
+        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(input_t)));
+    }
+
+    template <
+        typename Op,
+        typename Alloc,
+        typename Arch,
+        alignment_options aligned = alignment_options {},
+        map_options opts = map_options {}>
+    void bench_map_unary(benchmark::State& state)
+    {
+        bench_unary<Op, Alloc>(
+            state,
+            [](auto in, auto out)
+            { Op::template range_apply_map_unary<aligned, opts, Arch>(in, out); });
+    }
+
+    template <typename Op, typename Alloc, typename Arch>
+    void bench_transform(benchmark::State& state)
+    {
+        bench_unary<Op, Alloc>(
+            state,
+            [](auto in, auto out)
+            { Op::template range_apply_transform<Arch>(in, out); });
+    }
+
+    template <typename Op, typename Alloc>
+    void bench_scalar(benchmark::State& state)
+    {
+        bench_unary<Op, Alloc>(state, [](auto in, auto out)
+                               { Op::range_apply_scalar(in, out); });
+    }
+
+    template <typename Op, typename Arch, typename Bench>
+    void register_bench(
+        std::string_view variant,
+        Bench bench_fn,
+        std::vector<std::int64_t> const& sizes = xsimd::bench::bench_sizes<typename Op::input_t>())
+    {
+        using input_t = typename Op::input_t;
+
+        auto* bench = benchmark::RegisterBenchmark(
+            std::format("{}/{}/{}/{}", Arch::name(), Op::name, xsimd::bench::type_name<input_t>(), variant),
+            bench_fn);
+        for (auto const size : sizes)
+        {
+            bench->Arg(size);
+        }
+    }
+}
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/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/map.hpp b/include/xsimd_algorithm/map.hpp
new file mode 100644
index 0000000..12e9ed9
--- /dev/null
+++ b/include/xsimd_algorithm/map.hpp
@@ -0,0 +1,499 @@
+/****************************************************************************
+ * 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_MAP_HPP
+#define XSIMD_ALGORITHM_MAP_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
+{
+    /// 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, std::size_t et, typename U, std::size_t eu>
+    XSIMD_INLINE auto are_aliased(std::span<T, et> lhs, std::span<U, eu> 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());
+    }
+
+    /// Copy fewer than 2 * k elements without a call to memcpy.
+    ///
+    /// For count in [k, 2k), two overlapping copies of k elements cover the range; since k is a
+    /// compile-time constant, each memcpy compiles down to a few fixed-size loads and stores.
+    template <std::size_t k, typename T>
+    XSIMD_INLINE void copy_small(T* XSIMD_RESTRICT dst, T const* XSIMD_RESTRICT src, std::size_t count)
+    {
+        assert(std::has_single_bit(k));
+        assert(count < 2 * k);
+        if (count >= k)
+        {
+            std::memcpy(dst, src, k * sizeof(T));
+            std::memcpy(dst + count - k, src + count - k, k * sizeof(T));
+        }
+        else if (k > 1)
+        {
+            copy_small<k / 2>(dst, src, count);
+        }
+    }
+
+    /// 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);
+        }
+    }
+
+    struct alignment_options
+    {
+        bool start_aligned = false;
+        bool end_aligned = false;
+    };
+
+    struct map_options
+    {
+        std::size_t unroll_factor = 4;
+        bool pure = false;
+    };
+
+    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 algorithm 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 };
+                }
+            };
+        }
+
+        template <map_options opts, alignment_options align, typename A, typename Out, typename... In>
+        struct map_helper
+        {
+            static constexpr std::size_t n_input = sizeof...(In);
+            static constexpr auto input_elem_size = std::array { sizeof(In)... };
+
+            static constexpr std::size_t min_elem_size = std::min({ sizeof(Out), sizeof(In)... });
+            static constexpr std::size_t max_elem_size = std::max({ sizeof(Out), sizeof(In)... });
+
+            /// Inputs and output may not have the same alignment so it may be impossible
+            /// to get all aligned. We align preferably the output (more expensive unaligned
+            /// stores) or otherwise one of the input.
+            static constexpr bool align_output = sizeof(Out) == max_elem_size;
+
+            static constexpr std::array<bool, n_input> get_align_inputs()
+            {
+                std::array<bool, n_input> out {};
+                bool found = false;
+                for (std::size_t k = 0; k < out.size(); ++k)
+                {
+                    if (found || align_output)
+                    {
+                        out[k] = false;
+                    }
+                    else
+                    {
+                        found = input_elem_size[k] == max_elem_size;
+                        out[k] = found;
+                    }
+                }
+                return out;
+            }
+
+            static constexpr std::array<bool, n_input> align_inputs = get_align_inputs();
+
+            static constexpr bool output_is_aligned = align_output || align.start_aligned;
+
+            static constexpr std::array<bool, n_input> get_input_is_aligned()
+            {
+                std::array<bool, n_input> out = align_inputs;
+                for (bool& b : out)
+                {
+                    b = b || align.start_aligned;
+                }
+                return out;
+            }
+
+            static constexpr std::array<bool, n_input> input_is_aligned = get_input_is_aligned();
+
+            /// Number of batches of T spanning as many elements as one batch of the widest element.
+            ///
+            /// Pairing that many batches on all 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>
+            static constexpr std::size_t batch_arity()
+            {
+                return sizeof(T) / min_elem_size;
+            }
+
+            template <typename T>
+            using batch_array = std::array<xsimd::batch<T, A>, batch_arity<T>()>;
+
+            static constexpr std::size_t chunk_size = batch_arity<Out>() * xsimd::batch<Out, A>::size;
+
+            /// Load an array of batches.
+            template <bool aligned, typename T>
+            XSIMD_INLINE static auto load_batches(T const* ptr) -> batch_array<T>
+            {
+                batch_array<T> x;
+                for (std::size_t i = 0; i < x.size(); ++i)
+                {
+                    x[i] = load_batch<T, A, aligned>(ptr + i * xsimd::batch<T, A>::size);
+                }
+                return x;
+            }
+
+            /// Store an array of batches.
+            template <bool aligned, typename T>
+            XSIMD_INLINE static void store_batches(batch_array<T> const& x, T* ptr)
+            {
+                for (std::size_t i = 0; i < x.size(); ++i)
+                {
+                    store_batch<T, A, aligned>(x[i], ptr + i * xsimd::batch<T, A>::size);
+                }
+            }
+
+            /// Map a single unaligned chunk.
+            template <typename Func>
+            XSIMD_INLINE static void map_chunk_unaligned(
+                In const* XSIMD_RESTRICT... in,
+                Out* XSIMD_RESTRICT out,
+                Func&& func)
+            {
+                store_batches<false>(func(load_batches<false>(in)...), out);
+            }
+
+            /// Map multiple chunks with a compile-time unrolled loop.
+            template <
+                std::array<bool, n_input> load_aligned,
+                bool store_aligned,
+                typename Func,
+                std::size_t... step>
+            XSIMD_INLINE static void map_unrolled(
+                In const* XSIMD_RESTRICT... in,
+                Out* XSIMD_RESTRICT out,
+                Func&& func,
+                std::index_sequence<step...>)
+            {
+                static constexpr std::size_t factor = sizeof...(step);
+
+                constexpr auto read = []<class T, std::size_t... a>(T const* ptr, std::index_sequence<a...>)
+                {
+                    std::array<batch_array<T>, factor> x;
+                    auto load_one = [&](std::size_t s)
+                    {
+                        // Expand the parameter pack a for each alignment.
+                        ((x[s] = load_batches<load_aligned[a]>(ptr + s * chunk_size)), ...);
+                    };
+                    // Expand the step parameter pack: repeat the unrolled operation.
+                    (load_one(step), ...);
+                    return x;
+                };
+
+                constexpr auto map = [](auto* out, auto&& f, auto const&... x)
+                {
+                    auto map_one = [&](std::size_t s)
+                    {
+                        // Expand the parameter pack a for each input.
+                        store_batches<store_aligned>(f(x[s]...), out + s * chunk_size);
+                    };
+                    // Expand the step parameter pack: repeat the unrolled operation.
+                    (map_one(step), ...);
+                };
+
+                map(out, std::forward<Func>(func), read(in, std::index_sequence_for<In...> {})...);
+            }
+
+            /// Map chunks in a loop with given unrolling factor.
+            ///
+            /// Return number of elements mapped.
+            template <
+                std::array<bool, n_input> load_aligned,
+                bool store_aligned,
+                std::size_t unroll_factor,
+                typename Func>
+            XSIMD_INLINE static auto map_loop(
+                In const* XSIMD_RESTRICT... in,
+                Out* XSIMD_RESTRICT out,
+                std::size_t count,
+                Func&& func) -> std::size_t
+            {
+                constexpr auto steps = std::make_index_sequence<unroll_factor>();
+                constexpr std::size_t total_step_size = unroll_factor * chunk_size;
+
+                std::size_t remaining = count;
+                while (remaining >= total_step_size)
+                {
+                    map_unrolled<load_aligned, store_aligned>(in..., out, func, steps);
+                    ((in += total_step_size), ...);
+                    out += total_step_size;
+                    remaining -= total_step_size;
+                }
+                return count - remaining;
+            }
+
+            /// Map fewer elements than a full chunk through a scratch buffer.
+            template <typename Func>
+            XSIMD_INLINE static void map_chunk_partial(
+                In const* XSIMD_RESTRICT... begin,
+                Out* XSIMD_RESTRICT out,
+                std::size_t count,
+                Func&& func)
+            {
+                assert(count <= chunk_size);
+                if (count == 0) [[unlikely]]
+                {
+                    return;
+                }
+
+                constexpr auto read = []<typename T>(T const* in, std::size_t cnt)
+                {
+                    alignas(A::alignment()) std::array<T, chunk_size> in_buffer = {};
+                    copy_small<chunk_size>(in_buffer.data(), in, cnt);
+                    return load_batches<true>(in_buffer.data());
+                };
+
+                alignas(A::alignment()) std::array<Out, chunk_size> out_buffer;
+                store_batches<true>(func(read(begin, count)...), out_buffer.data());
+                copy_small<chunk_size>(out, out_buffer.data(), count);
+            }
+
+            /// Given some pointers, return the number of element to process until desired alignment.
+            ///
+            /// The desired alignment is given via the compile-time parameters.
+            /// Only one can be true.
+            template <std::array<bool, n_input> align_in, bool align_out>
+            XSIMD_INLINE static auto elems_to_alignment(In const*... in, Out* out) -> std::size_t
+            {
+                if constexpr (align_out)
+                {
+                    return bytes_to_next_aligned(out, A::alignment()) / sizeof(Out);
+                }
+                else
+                {
+                    constexpr auto iter = std::find(align_in.begin(), align_in.end(), true);
+                    static_assert(iter < align_in.end());
+                    constexpr auto idx = iter - align_in.begin();
+                    auto to_align = std::array { in... }[idx];
+                    return bytes_to_next_aligned(to_align, A::alignment()) / input_elem_size[idx];
+                }
+            }
+
+            template <typename Func>
+            XSIMD_INLINE static void map_n(
+                In const* XSIMD_RESTRICT... in,
+                Out* XSIMD_RESTRICT out,
+                std::size_t count,
+                Func&& func)
+            {
+                const auto advance = [&](std::size_t n)
+                {
+                    ((in += n), ...);
+                    out += n;
+                    count -= n;
+                };
+
+                if (count == 0) [[unlikely]]
+                {
+                    return;
+                }
+
+                if constexpr (!align.start_aligned)
+                {
+                    // The span may be too short to reach the next alignment boundary.
+                    const auto to_alignment = elems_to_alignment<align_inputs, align_output>(in..., out);
+                    const auto head = std::min(to_alignment, count);
+
+                    if (opts.pure && (head != 0) && (count >= chunk_size))
+                    {
+                        // Recompute the head as a full step, the body overwrites the excess.
+                        map_chunk_unaligned(in..., out, func);
+                    }
+                    else
+                    {
+                        map_chunk_partial(in..., out, head, func);
+                    }
+                    advance(head);
+                }
+
+                // Unrolled loop processing multiple chunks at a time.
+                auto processed = map_loop<input_is_aligned, output_is_aligned, opts.unroll_factor>(
+                    in..., out, count, func);
+                advance(processed);
+
+                // Regular simd loop one chunk at a time.
+                processed = map_loop<input_is_aligned, output_is_aligned, 1>(in..., out, count, func);
+                advance(processed);
+
+                // 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)
+                {
+                    if (opts.pure && (count != 0) && (count >= chunk_size)) [[likely]]
+                    {
+                        // Recompute overlapping data, this time starting from the end.
+                        map_chunk_unaligned((in + count - chunk_size)..., out + count - chunk_size, func);
+                    }
+                    else
+                    {
+                        map_chunk_partial(in..., out, count, func);
+                    }
+                }
+            }
+        };
+    }
+
+    /// Apply func elementwise over in, writing as many elements to out.
+    ///
+    /// Func maps as many input as given to the function.
+    /// If input as of different sizes, then the larger ones must be passed as an array of
+    /// as many batches as the factor to the smallest element size, so that the function
+    /// processes a fixed amount of elements.
+    template <
+        alignment_options align = alignment_options {},
+        map_options opts = map_options {},
+        typename Arch = xsimd::default_arch,
+        typename Func,
+        typename Out,
+        typename... In>
+    XSIMD_INLINE void map_n(Func&& func, Out&& out, In&&... in)
+    {
+        using H = internal::map_helper<
+            opts,
+            align,
+            Arch,
+            typename std::remove_cvref_t<Out>::value_type,
+            typename std::remove_cvref_t<In>::value_type...>;
+
+        assert((... && (in.size() == out.size())));
+        assert((... && !are_aliased(std::span(in.data(), in.size()), std::span(out.data(), out.size()))));
+
+        auto mapper = internal::wrap_params_as_1d_arrays(std::forward<Func>(func));
+        return H::map_n(in.data()..., out.data(), out.size(), mapper);
+    }
+
+    template <
+        alignment_options align = alignment_options {},
+        map_options opts = map_options {},
+        typename Arch = xsimd::default_arch,
+        typename Func,
+        typename Out,
+        typename In>
+    XSIMD_INLINE void map_unary(In&& in, Out&& out, Func&& func)
+    {
+        return map_n<align, opts, Arch>(func, out, in);
+    }
+
+    template <
+        alignment_options align = alignment_options {},
+        map_options opts = map_options {},
+        typename Arch = xsimd::default_arch,
+        typename Func,
+        typename Out,
+        typename Lhs,
+        typename Rhs>
+    XSIMD_INLINE void map_binary(Lhs&& lhs, Rhs&& rhs, Out&& out, Func&& func)
+    {
+        return map_n<align, opts, Arch>(func, out, lhs, rhs);
+    }
+}
+
+#endif
diff --git a/include/xsimd_algorithm/math.hpp b/include/xsimd_algorithm/math.hpp
new file mode 100644
index 0000000..514b0c4
--- /dev/null
+++ b/include/xsimd_algorithm/math.hpp
@@ -0,0 +1,79 @@
+/****************************************************************************
+ * 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 "./map.hpp"
+
+namespace xsimd::algo
+{
+    template <
+        xsimd::alignment_options align = xsimd::alignment_options{},
+        typename Arch = xsimd::default_arch,
+        typename T>
+    void sqrt(std::span<T const> in, std::span<T> out)
+    {
+        constexpr map_options opts = { .unroll_factor = 4, .pure = true };
+        return xsimd::map_unary<align, opts, Arch>(
+            in, out, [](auto x)
+            { return sqrt(x); });
+    }
+
+    template <
+        xsimd::alignment_options align = xsimd::alignment_options{},
+        typename Arch = xsimd::default_arch,
+        typename T>
+    void abs(std::span<T const> in, std::span<T> out)
+    {
+        constexpr map_options opts = { .unroll_factor = 4, .pure = true };
+        return xsimd::map_unary<align, opts, Arch>(
+            in, out, [](auto x)
+            { return abs(x); });
+    }
+
+    template <
+        xsimd::alignment_options align = xsimd::alignment_options{},
+        typename Arch = xsimd::default_arch,
+        typename T>
+    void exp(std::span<T const> in, std::span<T> out)
+    {
+        constexpr map_options opts = { .unroll_factor = 4, .pure = true };
+        return xsimd::map_unary<align, opts, Arch>(
+            in, out, [](auto x)
+            { return exp(x); });
+    }
+
+    template <
+        xsimd::alignment_options align = xsimd::alignment_options{},
+        typename Arch = xsimd::default_arch,
+        typename T>
+    void add(std::span<T const> lhs, std::span<T const> rhs, std::span<T> out)
+    {
+        constexpr map_options opts = { .unroll_factor = 4, .pure = true };
+        return xsimd::map_binary<align, opts, Arch>(
+            lhs, rhs, out, [](auto x, auto y)
+            { return x + y; });
+    }
+
+    template <
+        xsimd::alignment_options align = xsimd::alignment_options{},
+        typename Arch = xsimd::default_arch,
+        typename T>
+    void multiply(std::span<T const> lhs, std::span<T const> rhs, std::span<T> out)
+    {
+        constexpr map_options opts = { .unroll_factor = 4, .pure = true };
+        return xsimd::map_binary<align, opts, Arch>(
+            lhs, rhs, out, [](auto x, auto y)
+            { return x * y; });
+    }
+}
+
+#endif
diff --git a/include/xsimd_algorithm/stl/transform.hpp b/include/xsimd_algorithm/stl/transform.hpp
index d757ff8..2c53738 100644
--- a/include/xsimd_algorithm/stl/transform.hpp
+++ b/include/xsimd_algorithm/stl/transform.hpp
@@ -16,7 +16,7 @@
 #include <iterator>
 #include <type_traits>
 
-#include "xsimd/xsimd.hpp"
+#include <xsimd/xsimd.hpp>
 
 namespace xsimd
 {
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/map_binary_data.hpp b/test-utils/include/xsimd_test_utils/map_binary_data.hpp
new file mode 100644
index 0000000..4df143e
--- /dev/null
+++ b/test-utils/include/xsimd_test_utils/map_binary_data.hpp
@@ -0,0 +1,122 @@
+/****************************************************************************
+ * 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_MAP_BINARY_DATA_HPP
+#define XSIMD_ALGORITHM_TEST_UTILS_MAP_BINARY_DATA_HPP
+
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <memory>
+#include <span>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+#include "xsimd_algorithm/map.hpp"
+
+#include "xsimd_test_utils/utils.hpp"
+
+namespace xsimd::test
+{
+    using xsimd::alignment_options;
+    using xsimd::map_options;
+
+    /// Derives the scalar range application from the element-wise Derived::apply.
+    template <typename Derived, typename Lhs, typename Rhs = Lhs, typename Out = Lhs>
+    struct binary_op
+    {
+        using lhs_t = Lhs;
+        using rhs_t = Rhs;
+        using output_t = Out;
+
+        template <typename Alloc, typename T>
+        using rebind_allocator = typename std::allocator_traits<Alloc>::template rebind_alloc<T>;
+
+        static void range_apply_scalar(std::span<lhs_t const> lhs, std::span<rhs_t const> rhs, std::span<output_t> out)
+        {
+            for (std::size_t i = 0; i < lhs.size(); ++i)
+            {
+                out[i] = Derived::apply_scalar(lhs[i], rhs[i]);
+            }
+        }
+
+        template <
+            alignment_options aligned = alignment_options {},
+            map_options opts = map_options {},
+            typename Arch = xsimd::default_arch>
+        static void range_apply_map_binary(std::span<lhs_t const> lhs, std::span<rhs_t const> rhs, std::span<output_t> out)
+        {
+            return xsimd::map_binary<aligned, opts, Arch>(
+                lhs, rhs, out, [](auto x, auto y)
+                { return Derived::apply_batch(x, y); });
+        }
+
+        template <typename Alloc>
+        static auto make_input_output(std::size_t size)
+            -> std::tuple<
+                std::vector<lhs_t, rebind_allocator<Alloc, lhs_t>>,
+                std::vector<rhs_t, rebind_allocator<Alloc, rhs_t>>,
+                std::vector<output_t, rebind_allocator<Alloc, output_t>>>
+        {
+            auto lhs = make_arange<lhs_t, rebind_allocator<Alloc, lhs_t>>(size, -static_cast<lhs_t>(size / 2));
+            auto rhs = make_arange<rhs_t, rebind_allocator<Alloc, rhs_t>>(size, static_cast<rhs_t>(3));
+            auto output = std::vector<output_t, rebind_allocator<Alloc, output_t>>(size);
+            return { std::move(lhs), std::move(rhs), std::move(output) };
+        }
+    };
+
+    /*******************
+     *  Test fixtures  *
+     *******************/
+
+    template <typename T>
+    struct add_op : binary_op<add_op<T>, T>
+    {
+        static constexpr auto name = "add";
+        static constexpr bool pure = true;
+
+        static auto apply_scalar(T x, T y) { return static_cast<T>(x + y); }
+
+        template <typename A>
+        static auto apply_batch(xsimd::batch<T, A> x, xsimd::batch<T, A> y) { return x + y; }
+    };
+
+    template <typename T>
+    struct multiply_op : binary_op<multiply_op<T>, T>
+    {
+        static constexpr auto name = "multiply";
+        static constexpr bool pure = true;
+
+        static auto apply_scalar(T x, T y) { return static_cast<T>(x * y); }
+
+        template <typename A>
+        static auto apply_batch(xsimd::batch<T, A> x, xsimd::batch<T, A> y) { return x * y; }
+    };
+
+    /// Multiply an int32 by a double, one int32 batch with two double batches.
+    struct mixed_multiply_op : binary_op<mixed_multiply_op, std::int32_t, double, double>
+    {
+        static constexpr auto name = "mixed_multiply";
+        static constexpr bool pure = true;
+
+        static auto apply_scalar(std::int32_t x, double y) { return static_cast<double>(x) * y; }
+
+        template <typename A>
+        static auto apply_batch(xsimd::batch<std::int32_t, A> x, std::array<xsimd::batch<double, A>, 2> const& y)
+        {
+            auto const wide = xsimd::widen(x);
+            return std::array {
+                xsimd::batch_cast<double>(wide[0]) * y[0],
+                xsimd::batch_cast<double>(wide[1]) * y[1],
+            };
+        }
+    };
+}
+
+#endif
diff --git a/test-utils/include/xsimd_test_utils/map_unary_data.hpp b/test-utils/include/xsimd_test_utils/map_unary_data.hpp
new file mode 100644
index 0000000..e1e3fd9
--- /dev/null
+++ b/test-utils/include/xsimd_test_utils/map_unary_data.hpp
@@ -0,0 +1,171 @@
+/****************************************************************************
+ * 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 <memory>
+#include <span>
+#include <utility>
+#include <vector>
+
+#include "xsimd_algorithm/map.hpp"
+#include "xsimd_algorithm/stl/transform.hpp"
+
+#include "xsimd_test_utils/utils.hpp"
+
+namespace xsimd::test
+{
+    using xsimd::alignment_options;
+    using xsimd::map_options;
+
+    /// Derives the scalar range application from the element-wise Derived::apply.
+    template <typename Derived, typename In, typename Out = In>
+    struct unary_op
+    {
+        using input_t = In;
+        using output_t = Out;
+
+        /// Allocator of output_t matching an allocator of input_t.
+        template <typename Alloc>
+        using output_allocator = typename std::allocator_traits<Alloc>::template rebind_alloc<output_t>;
+
+        static void range_apply_scalar(std::span<input_t const> in, std::span<output_t> out)
+        {
+            for (std::size_t i = 0; i < in.size(); ++i)
+            {
+                out[i] = Derived::apply_scalar(in[i]);
+            }
+        }
+
+        template <
+            alignment_options aligned = alignment_options {},
+            map_options opts = map_options {},
+            typename Arch = xsimd::default_arch>
+        static void range_apply_map_unary(std::span<input_t const> in, std::span<output_t> out)
+        {
+            return xsimd::map_unary<aligned, opts, Arch>(
+                in, out, [](auto x)
+                { return Derived::apply_batch(x); });
+        }
+
+        template <typename Arch = xsimd::default_arch>
+        inline static void range_apply_transform(std::span<input_t const> in, std::span<output_t> out)
+        {
+            return xsimd::transform<Arch>(
+                in.data(), in.data() + in.size(), out.data(),
+                []<typename T>(T x)
+                {
+                    if constexpr (xsimd::is_batch<T>::value)
+                    {
+                        return Derived::apply_batch(x);
+                    }
+                    else
+                    {
+                        return Derived::apply_scalar(x);
+                    }
+                });
+        }
+
+        template <typename Alloc>
+        static auto make_input_output(std::size_t size)
+            -> std::pair<std::vector<input_t, Alloc>, std::vector<output_t, output_allocator<Alloc>>>
+        {
+            auto input = Derived::template make_input<Alloc>(size);
+            auto output = std::vector<output_t, output_allocator<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 constexpr bool pure = true;
+
+        static auto apply_scalar(T x) { return std::sqrt(x); }
+
+        template <typename A>
+        static auto apply_batch(xsimd::batch<T, A> x) { return xsimd::sqrt(x); }
+
+        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 constexpr bool pure = true;
+
+        static auto apply_scalar(T x) { return std::abs(x); }
+
+        template <typename A>
+        static auto apply_batch(xsimd::batch<T, A> x) { return xsimd::abs(x); }
+
+        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 constexpr bool pure = true;
+
+        static auto apply_scalar(T x) { return std::exp(x); }
+
+        template <typename A>
+        static auto apply_batch(xsimd::batch<T, A> x) { return xsimd::exp(x); }
+
+        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;
+        }
+    };
+
+    /// Sign-extend to the type with twice as many bytes, one input batch to two output batches.
+    template <typename T>
+    struct widen_op : unary_op<widen_op<T>, T, xsimd::widen_t<T>>
+    {
+        static constexpr auto name = "widen";
+        static constexpr bool pure = true;
+
+        static auto apply_scalar(T x) { return static_cast<xsimd::widen_t<T>>(x); }
+
+        template <typename A>
+        static auto apply_batch(xsimd::batch<T, A> x) { return xsimd::widen(x); }
+
+        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));
+        }
+    };
+}
+
+#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..4414e98
--- /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/map.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::alignment_options aligned = xsimd::alignment_options{}>
+        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::alignment_options aligned = xsimd::alignment_options{}>
+        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..7aba53b 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -48,13 +48,15 @@ endif()
 set(XSIMD_ALGORITHM_TESTS
     main.cpp
     test_arange.cpp
+    test_map.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_map.cpp b/test/test_map.cpp
new file mode 100644
index 0000000..d3f7d05
--- /dev/null
+++ b/test/test_map.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/map.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::map_unary(input, 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::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..26c1cf9
--- /dev/null
+++ b/test/test_math.cpp
@@ -0,0 +1,126 @@
+/****************************************************************************
+ * 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 <type_traits>
+
+#include <doctest/doctest.h>
+
+#include <xsimd_test_utils/map_binary_data.hpp>
+#include <xsimd_test_utils/map_unary_data.hpp>
+#include <xsimd_test_utils/utils.hpp>
+
+namespace
+{
+    template <typename Op, typename Alloc, xsimd::alignment_options aligned = xsimd::alignment_options {}>
+    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 range_apply_map_unary<aligned>(xsimd::test::as_span(input), xsimd::test::as_span(output));
+
+        for (std::size_t i = 0; i < input.size(); ++i)
+        {
+            CAPTURE(i);
+            if constexpr (std::is_floating_point_v<typename Op::output_t>)
+            {
+                CHECK(output[i] == doctest::Approx(Op::apply_scalar(input[i])));
+            }
+            else
+            {
+                CHECK(output[i] == Op::apply_scalar(input[i]));
+            }
+        }
+    }
+
+    template <typename Op, typename Alloc, xsimd::alignment_options aligned = xsimd::alignment_options {}>
+    void check_binary_math()
+    {
+        // Not a multiple of the batch size, to exercise the tail.
+        constexpr std::size_t test_size = 94;
+
+        auto [lhs, rhs, output] = Op::template make_input_output<Alloc>(test_size);
+
+        Op::template range_apply_map_binary<aligned>(
+            xsimd::test::as_span(lhs), xsimd::test::as_span(rhs), xsimd::test::as_span(output));
+
+        for (std::size_t i = 0; i < lhs.size(); ++i)
+        {
+            CAPTURE(i);
+            if constexpr (std::is_floating_point_v<typename Op::output_t>)
+            {
+                CHECK(output[i] == doctest::Approx(Op::apply_scalar(lhs[i], rhs[i])));
+            }
+            else
+            {
+                CHECK(output[i] == Op::apply_scalar(lhs[i], rhs[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>,
+    xsimd::test::widen_op<std::int32_t>)
+{
+    using input_t = typename Op::input_t;
+    using aligned_allocator = typename xsimd::test::aligned_vector<input_t>::allocator_type;
+    using unaligned_allocator = typename xsimd::test::unaligned_vector<input_t>::allocator_type;
+
+    SUBCASE("aligned without header")
+    {
+        check_unary_math<Op, aligned_allocator, xsimd::alignment_options { .start_aligned = true }>();
+    }
+
+    SUBCASE("aligned with header")
+    {
+        check_unary_math<Op, aligned_allocator>();
+    }
+
+    SUBCASE("unaligned with header")
+    {
+        check_unary_math<Op, unaligned_allocator>();
+    }
+}
+
+TEST_CASE_TEMPLATE(
+    "binary math",
+    Op,
+    xsimd::test::add_op<std::int32_t>,
+    xsimd::test::multiply_op<float>,
+    xsimd::test::mixed_multiply_op)
+{
+    using lhs_t = typename Op::lhs_t;
+    using aligned_allocator = typename xsimd::test::aligned_vector<lhs_t>::allocator_type;
+    using unaligned_allocator = typename xsimd::test::unaligned_vector<lhs_t>::allocator_type;
+
+    SUBCASE("aligned without header")
+    {
+        check_binary_math<Op, aligned_allocator, xsimd::alignment_options { .start_aligned = true }>();
+    }
+
+    SUBCASE("aligned with header")
+    {
+        check_binary_math<Op, aligned_allocator>();
+    }
+
+    SUBCASE("unaligned with header")
+    {
+        check_binary_math<Op, unaligned_allocator>();
+    }
+}