diff --git a/CMakeLists.txt b/CMakeLists.txt
index d01c68f..0c7d022 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -144,10 +144,22 @@ endif()
 
 message(STATUS "MATH: ${MATH_BACKEND}")
 
-if (MATH_BACKEND STREQUAL "Fastor")
+if (MATH_BACKEND STREQUAL "Arma")
+    find_package(BLAS REQUIRED)
+    find_package(LAPACK REQUIRED)
+    find_package(Armadillo REQUIRED)
+    target_link_libraries(operon_operon PUBLIC BLAS::BLAS LAPACK::LAPACK ${ADMADILLO_LIBRARIES})
+    target_compile_definitions(operon_operon PUBLIC OPERON_MATH_ARMA)
+elseif (MATH_BACKEND STREQUAL "Blaze")
+    find_package(blaze REQUIRED)
+    find_package(sleef REQUIRED)
+    target_link_libraries(operon_operon INTERFACE sleef::sleef)
+    target_compile_definitions(operon_operon PUBLIC OPERON_MATH_BLAZE BLAZE_USE_SHARED_MEMORY_PARALLELIZATION=0 BLAZE_USE_SLEEF=1 EIGEN_DONT_PARALLELIZE)
+elseif (MATH_BACKEND STREQUAL "Fastor")
     find_package(Fastor REQUIRED)
-    target_link_libraries(operon_operon PUBLIC Fastor::Fastor)
-    target_compile_definitions(operon_operon PUBLIC OPERON_MATH_FASTOR)
+    find_package(sleef REQUIRED)
+    target_link_libraries(operon_operon PUBLIC Fastor::Fastor sleef::sleef)
+    target_compile_definitions(operon_operon PUBLIC OPERON_MATH_FASTOR FASTOR_USE_SLEEF_U35)
 elseif (MATH_BACKEND STREQUAL "Eve")
     target_compile_definitions(operon_operon PUBLIC OPERON_MATH_EVE)
 elseif (MATH_BACKEND STREQUAL "Vdt")
@@ -156,7 +168,7 @@ elseif (MATH_BACKEND STREQUAL "Vdt")
     target_compile_definitions(operon_operon PUBLIC OPERON_MATH_VDT)
 elseif (MATH_BACKEND STREQUAL "Eigen")
     message(STATUS "Using Eigen backend")
-    target_compile_definitions(operon_operon PUBLIC OPERON_MATH_EIGEN)
+    target_compile_definitions(operon_operon PUBLIC OPERON_MATH_EIGEN EIGEN_DONT_PARALLELIZE)
 elseif(MATH_BACKEND STREQUAL "Stl")
     target_compile_definitions(operon_operon PUBLIC OPERON_MATH_STL)
 elseif (MATH_BACKEND STREQUAL "Fast_v1")
diff --git a/flake.lock b/flake.lock
index 3ff915e..308ee1d 100644
--- a/flake.lock
+++ b/flake.lock
@@ -97,11 +97,11 @@
         ]
       },
       "locked": {
-        "lastModified": 1714463628,
-        "narHash": "sha256-58N7PhmMvArtOXurpGLa2QtDJCfE9v5YvEnz7gPf3wo=",
+        "lastModified": 1715862037,
+        "narHash": "sha256-s/D7vvGELAvk0oU4qnoAoOTIIYw6jXT3spPPkwMWxpU=",
         "owner": "foolnotion",
         "repo": "nur-pkg",
-        "rev": "eae6690e3cf1f5866ed0713b82b0d0d2b5575bcb",
+        "rev": "396182e1a0e2fc5a4f6c6137bbb45c732049574b",
         "type": "github"
       },
       "original": {
@@ -118,11 +118,11 @@
         ]
       },
       "locked": {
-        "lastModified": 1714463628,
-        "narHash": "sha256-58N7PhmMvArtOXurpGLa2QtDJCfE9v5YvEnz7gPf3wo=",
+        "lastModified": 1715862037,
+        "narHash": "sha256-s/D7vvGELAvk0oU4qnoAoOTIIYw6jXT3spPPkwMWxpU=",
         "owner": "foolnotion",
         "repo": "nur-pkg",
-        "rev": "eae6690e3cf1f5866ed0713b82b0d0d2b5575bcb",
+        "rev": "396182e1a0e2fc5a4f6c6137bbb45c732049574b",
         "type": "github"
       },
       "original": {
@@ -136,11 +136,11 @@
         "nixpkgs": "nixpkgs_2"
       },
       "locked": {
-        "lastModified": 1714463628,
-        "narHash": "sha256-58N7PhmMvArtOXurpGLa2QtDJCfE9v5YvEnz7gPf3wo=",
+        "lastModified": 1715862037,
+        "narHash": "sha256-s/D7vvGELAvk0oU4qnoAoOTIIYw6jXT3spPPkwMWxpU=",
         "owner": "foolnotion",
         "repo": "nur-pkg",
-        "rev": "eae6690e3cf1f5866ed0713b82b0d0d2b5575bcb",
+        "rev": "396182e1a0e2fc5a4f6c6137bbb45c732049574b",
         "type": "github"
       },
       "original": {
@@ -154,11 +154,11 @@
         "nixpkgs": "nixpkgs_3"
       },
       "locked": {
-        "lastModified": 1714463628,
-        "narHash": "sha256-58N7PhmMvArtOXurpGLa2QtDJCfE9v5YvEnz7gPf3wo=",
+        "lastModified": 1715862037,
+        "narHash": "sha256-s/D7vvGELAvk0oU4qnoAoOTIIYw6jXT3spPPkwMWxpU=",
         "owner": "foolnotion",
         "repo": "nur-pkg",
-        "rev": "eae6690e3cf1f5866ed0713b82b0d0d2b5575bcb",
+        "rev": "396182e1a0e2fc5a4f6c6137bbb45c732049574b",
         "type": "github"
       },
       "original": {
@@ -191,11 +191,11 @@
     },
     "nixpkgs": {
       "locked": {
-        "lastModified": 1714463833,
-        "narHash": "sha256-2QTF0WufldsBDkB7dl0WkaxN+V1svbgfRQlXVGayolQ=",
+        "lastModified": 1715861982,
+        "narHash": "sha256-UhgFvRo4o7ygyQcYV/izRHLsrW00VqL4mxbFggkBAKg=",
         "owner": "nixos",
         "repo": "nixpkgs",
-        "rev": "38cc39329a91a925aba2fb32db33fd759a5306fa",
+        "rev": "9499b5362868c001e7317d4b11a8e21ae788918d",
         "type": "github"
       },
       "original": {
@@ -207,11 +207,11 @@
     },
     "nixpkgs_2": {
       "locked": {
-        "lastModified": 1714463833,
-        "narHash": "sha256-2QTF0WufldsBDkB7dl0WkaxN+V1svbgfRQlXVGayolQ=",
+        "lastModified": 1715861982,
+        "narHash": "sha256-UhgFvRo4o7ygyQcYV/izRHLsrW00VqL4mxbFggkBAKg=",
         "owner": "NixOS",
         "repo": "nixpkgs",
-        "rev": "38cc39329a91a925aba2fb32db33fd759a5306fa",
+        "rev": "9499b5362868c001e7317d4b11a8e21ae788918d",
         "type": "github"
       },
       "original": {
@@ -223,11 +223,11 @@
     },
     "nixpkgs_3": {
       "locked": {
-        "lastModified": 1714463833,
-        "narHash": "sha256-2QTF0WufldsBDkB7dl0WkaxN+V1svbgfRQlXVGayolQ=",
+        "lastModified": 1715861982,
+        "narHash": "sha256-UhgFvRo4o7ygyQcYV/izRHLsrW00VqL4mxbFggkBAKg=",
         "owner": "NixOS",
         "repo": "nixpkgs",
-        "rev": "38cc39329a91a925aba2fb32db33fd759a5306fa",
+        "rev": "9499b5362868c001e7317d4b11a8e21ae788918d",
         "type": "github"
       },
       "original": {
diff --git a/flake.nix b/flake.nix
index e3da756..2143bb8 100644
--- a/flake.nix
+++ b/flake.nix
@@ -44,6 +44,8 @@
 
           buildInputs = (with pkgs; [
             aria-csv
+            armadillo
+            blaze
             ceres-solver
             cpp-sort
             cxxopts
@@ -63,6 +65,7 @@
             pratt-parser.packages.${system}.default
             simdutf_4 # required by scnlib
             scnlib
+            sleef
             taskflow
             unordered_dense
             vdt.packages.${system}.default
@@ -75,6 +78,7 @@
             ned14-quickcpplib
             ned14-status-code
             xad
+            xsimd
             xxHash
             zstd
           ]);
diff --git a/include/operon/interpreter/backend/arma.hpp b/include/operon/interpreter/backend/arma.hpp
new file mode 100644
index 0000000..1859559
--- /dev/null
+++ b/include/operon/interpreter/backend/arma.hpp
@@ -0,0 +1,9 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright 2019-2024 Heal Research
+
+#ifndef OPERON_BACKEND_ARMA_HPP
+#define OPERON_BACKEND_ARMA_HPP
+
+#include "arma/derivatives.hpp"
+
+#endif
\ No newline at end of file
diff --git a/include/operon/interpreter/backend/arma/derivatives.hpp b/include/operon/interpreter/backend/arma/derivatives.hpp
new file mode 100644
index 0000000..183928f
--- /dev/null
+++ b/include/operon/interpreter/backend/arma/derivatives.hpp
@@ -0,0 +1,192 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright 2019-2023 Heal Research
+
+#ifndef OPERON_BACKEND_ARMA_DERIVATIVES_HPP
+#define OPERON_BACKEND_ARMA_DERIVATIVES_HPP
+
+#include "operon/core/node.hpp"
+#include "functions.hpp"
+
+namespace Operon::Backend {
+namespace detail {
+    template<typename T>
+    inline auto IsNaN(T value) { return std::isnan(value); }
+
+    template<typename Compare>
+    struct FComp {
+        auto operator()(auto x, auto y) const {
+            using T = std::common_type_t<decltype(x), decltype(y)>;
+            if ((IsNaN(x) && IsNaN(y)) || (x == y)) {
+                return std::numeric_limits<T>::quiet_NaN();
+            }
+            if (IsNaN(x)) { return T{0}; }
+            if (IsNaN(y)) { return T{1}; }
+            return static_cast<T>(Compare{}(T{x}, T{y}));
+        }
+    };
+} // namespace detail
+
+    template<typename T, std::size_t S>
+    auto Add(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> /*primal*/, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j).fill(T{1});
+    }
+
+    template<typename T, std::size_t S>
+    auto Mul(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        Col(trace, j) = Col(primal, i) / Col(primal, j);
+    }
+
+    template<typename T, std::size_t S>
+    auto Sub(std::vector<Operon::Node> const& nodes, Backend::View<T const, S> /*primal*/, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        auto v = (nodes[i].Arity == 1 || j < i-1) ? T{-1} : T{+1};
+        Col(trace, j).fill(v);
+    }
+
+    template<typename T, std::size_t S>
+    auto Div(std::vector<Operon::Node> const& nodes, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        auto const& n = nodes[i];
+        if (n.Arity == 1) {
+            Col(trace, j) = -T{1} / (arma::square(Col(primal, j)));
+        } else {
+            Col(trace, j) = (j == i-1 ? T{1} : T{-1}) * Col(primal, i) / Col(primal, j);
+        }
+    }
+
+    template<typename T, std::size_t S>
+    auto Aq(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        if (j == i-1) {
+            Col(trace, j) = Col(primal, i) / Col(primal, j);
+        } else {
+            Col(trace, j) = -Col(primal, j) % arma::pow(Col(primal, i), T{3}) / arma::square(Col(primal, i-1));
+        }
+    }
+
+    template<typename T, std::size_t S>
+    auto Pow(std::vector<Operon::Node> const& nodes, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        if (j == i-1) {
+            auto const k = j - (nodes[j].Length + 1);
+            Col(trace, j) = Col(primal, i) % Col(primal, k) / Col(primal, j);
+        } else {
+            auto const k = i-1;
+            Col(trace, j) = Col(primal, i) % arma::log(Col(primal, k));
+        }
+    }
+
+    template<typename T, std::size_t S>
+    auto Min(std::vector<Operon::Node> const& nodes, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        auto k = j == i - 1 ? (j - nodes[j].Length - 1) : i - 1;
+        auto const* jp = Ptr(primal, j);
+        auto const* kp = Ptr(primal, k);
+        std::transform(jp, jp+S, kp, Ptr(trace, j), detail::FComp<std::less<>>{});
+    }
+
+    template<typename T, std::size_t S>
+    auto Max(std::vector<Operon::Node> const& nodes, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        auto k = j == i - 1 ? (j - nodes[j].Length - 1) : i - 1;
+        auto const* jp = Ptr(primal, j);
+        auto const* kp = Ptr(primal, k);
+        std::transform(jp, jp+S, kp, Ptr(trace, j), detail::FComp<std::greater<>>{});
+    }
+
+    template<typename T, std::size_t S>
+    auto Square(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{2} * Col(primal, j);
+    }
+
+    template<typename T, std::size_t S>
+    auto Abs(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = arma::sign(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Ceil(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = arma::ceil(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Floor(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = arma::floor(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Exp(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        Col(trace, j) = Col(primal, i);
+    }
+
+    template<typename T, std::size_t S>
+    auto Log(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1} / (Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Log1p(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1} / (T{1} + Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Logabs(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = arma::sign(Col(primal, j)) / arma::abs(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sin(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = arma::cos(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cos(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = -arma::sin(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Tan(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1} + arma::square(arma::tan(Col(primal, j)));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sinh(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = arma::cosh(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cosh(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = arma::sinh(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Tanh(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1} - arma::square(arma::tanh(Col(primal, j)));
+    }
+
+    template<typename T, std::size_t S>
+    auto Asin(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1} / (arma::sqrt(T{1} - arma::square(Col(primal, j))));
+    }
+
+    template<typename T, std::size_t S>
+    auto Acos(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = -T{1} / (arma::sqrt(((T{1} - arma::square(Col(primal, j))))));
+    }
+
+    template<typename T, std::size_t S>
+    auto Atan(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1} / (T{1} + arma::square(Col(primal, j)));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sqrt(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        Col(trace, j) = T{1} / (T{2} * Col(primal, i));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sqrtabs(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        Col(trace, j) = arma::sign(Col(primal, j)) / (T{2} * Col(primal, i));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cbrt(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        Col(trace, j) = T{1} / (T{3} * arma::square(Col(primal, i)));
+    }
+}  // namespace Operon::Backend
+
+#endif
\ No newline at end of file
diff --git a/include/operon/interpreter/backend/arma/functions.hpp b/include/operon/interpreter/backend/arma/functions.hpp
new file mode 100644
index 0000000..5ff8de9
--- /dev/null
+++ b/include/operon/interpreter/backend/arma/functions.hpp
@@ -0,0 +1,204 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright 2019-2023 Heal Research
+
+#ifndef OPERON_BACKEND_ARMA_FUNCTIONS_HPP
+#define OPERON_BACKEND_ARMA_FUNCTIONS_HPP
+
+#define ARMA_DONT_USE_WRAPPER
+#include <armadillo>
+#include "operon/interpreter/backend/backend.hpp"
+
+namespace Operon::Backend {
+    template<typename T, std::size_t S>
+    auto Map(T* ptr) {
+        return arma::Mat<Operon::Scalar>(ptr, S, 1, /*copy_aux_mem =*/false, /*strict=*/true);
+    }
+
+    template<typename T, std::size_t S>
+    auto Map(T const* ptr) {
+        using U = std::remove_const_t<T>;
+        return arma::Mat<Operon::Scalar>(const_cast<U*>(ptr), S, 1, /*copy_aux_mem =*/false, /*strict=*/true);
+    }
+
+    template<typename T, std::size_t S>
+    auto Col(Backend::View<T, S> view, std::integral auto col) {
+        return Map<T, S>(view.data_handle() + col * S);
+    }
+
+    template<typename T, std::size_t S>
+    auto Col(Backend::View<T const, S> view, std::integral auto col) {
+        return Map<T const, S>(view.data_handle() + col * S);
+    }
+
+    // utility
+    template<typename T, std::size_t S>
+    auto Fill(T* res, T value) {
+        Map<T, S>(res).fill(value);
+    }
+
+    // n-ary functions
+    template<typename T, std::size_t S>
+    auto Add(T* res, auto const*... args) {
+        Map<T, S>(res) = (Map<T const, S>(args) + ...);
+    }
+
+    template<typename T, std::size_t S>
+    auto Mul(T* res, auto const*... args) {
+        Map<T, S>(res) = (Map<T const, S>(args) % ...);
+    }
+
+    template<typename T, std::size_t S>
+    auto Sub(T* res, auto* first, auto const*... rest) {
+        static_assert(sizeof...(rest) > 0);
+        Map<T, S>(res) = Map<T const, S>(first) - (Map<T const, S>(rest) + ...);
+    }
+
+    template<typename T, std::size_t S>
+    auto Div(T* res, auto const* first, auto const*... rest) {
+        static_assert(sizeof...(rest) > 0);
+        Map<T, S>(res) = Map<T const, S>(first) / (Map<T const, S>(rest) * ...);
+    }
+
+    template<typename T, std::size_t S>
+    auto Min(T* res, auto const* first, auto const*... args) {
+        static_assert(sizeof...(args) > 0);
+        for (auto i = 0UL; i < S; ++i) {
+            res[i] = std::min({first[i], args[i]...});
+        }
+    }
+
+    template<typename T, std::size_t S>
+    auto Max(T* res, auto* first, auto const*... args) {
+        for (auto i = 0UL; i < S; ++i) {
+            res[i] = std::max({first[i], args[i]...});
+        }
+    }
+
+    // binary functions
+    template<typename T, std::size_t S>
+    auto Aq(T* res, T const* a, T const* b) {
+        Map<T, S>(res) = Map<T const, S>(a) / arma::sqrt((T{1} + arma::square(Map<T const, S>(b))));
+    }
+
+    template<typename T, std::size_t S>
+    auto Pow(T* res, T const* a, T const* b) {
+        Map<T, S>(res) = arma::pow(Map<T const, S>(a), Map<T const, S>(b));
+    }
+
+    // unary functions
+    template<typename T, std::size_t S>
+    auto Cpy(T* res, T const* arg) {
+        Map<T, S>(res) = Map<T const, S>(arg);
+    }
+
+    template<typename T, std::size_t S>
+    auto Neg(T* res, T const* arg) {
+        Map<T, S>(res) = -Map<T const, S>(arg);
+    }
+
+    template<typename T, std::size_t S>
+    auto Inv(T* res, T const* arg) {
+        Map<T, S>(res) = T{1} / Map<T const, S>(arg);
+    }
+
+    template<typename T, std::size_t S>
+    auto Abs(T* res, T const* arg) {
+        Map<T, S>(res) = arma::abs(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Ceil(T* res, T const* arg) {
+        Map<T, S>(res) = arma::ceil(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Floor(T* res, T const* arg) {
+        Map<T, S>(res) = arma::floor(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Square(T* res, T const* arg) {
+        Map<T, S>(res) = arma::square(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Exp(T* res, T const* arg) {
+        Map<T, S>(res) = arma::exp(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Log(T* res, T const* arg) {
+        Map<T, S>(res) = arma::log(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Log1p(T* res, T const* arg) {
+        Map<T, S>(res) = arma::log1p(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Logabs(T* res, T const* arg) {
+        Map<T, S>(res) = arma::log(arma::abs(Map<T const, S>(arg)));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sin(T* res, T const* arg) {
+        Map<T, S>(res) = arma::sin(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cos(T* res, T const* arg) {
+        Map<T, S>(res) = arma::cos(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Tan(T* res, T const* arg) {
+        Map<T, S>(res) = arma::tan(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Asin(T* res, T const* arg) {
+        Map<T, S>(res) = arma::asin(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Acos(T* res, T const* arg) {
+        Map<T, S>(res) = arma::acos(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Atan(T* res, T const* arg) {
+        Map<T, S>(res) = arma::atan(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sinh(T* res, T const* arg) {
+        Map<T, S>(res) = arma::sinh(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cosh(T* res, T const* arg) {
+        Map<T, S>(res) = arma::cosh(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Tanh(T* res, T const* arg) {
+        Map<T, S>(res) = arma::tanh(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sqrt(T* res, T const* arg) {
+        Map<T, S>(res) = arma::sqrt(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sqrtabs(T* res, T const* arg) {
+        Map<T, S>(res) = arma::sqrt(arma::abs(Map<T const, S>(arg)));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cbrt(T* res, T const* arg) {
+        Map<T, S>(res) = arma::cbrt(Map<T const, S>(arg));
+    }
+}  // namespace Operon::Backend
+#endif
diff --git a/include/operon/interpreter/backend/blaze.hpp b/include/operon/interpreter/backend/blaze.hpp
new file mode 100644
index 0000000..78c36c5
--- /dev/null
+++ b/include/operon/interpreter/backend/blaze.hpp
@@ -0,0 +1,9 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright 2019-2024 Heal Research
+
+#ifndef OPERON_BACKEND_BLAZE_HPP
+#define OPERON_BACKEND_BLAZE_HPP
+
+#include "blaze/derivatives.hpp"
+
+#endif
\ No newline at end of file
diff --git a/include/operon/interpreter/backend/blaze/derivatives.hpp b/include/operon/interpreter/backend/blaze/derivatives.hpp
new file mode 100644
index 0000000..0c1f591
--- /dev/null
+++ b/include/operon/interpreter/backend/blaze/derivatives.hpp
@@ -0,0 +1,194 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright 2019-2023 Heal Research
+
+#ifndef OPERON_BACKEND_BLAZE_DERIVATIVES_HPP
+#define OPERON_BACKEND_BLAZE_DERIVATIVES_HPP
+
+#include "operon/core/node.hpp"
+#include "functions.hpp"
+
+namespace Operon::Backend {
+namespace detail {
+    template<typename T>
+    inline auto IsNaN(T value) { return std::isnan(value); }
+
+    template<typename Compare>
+    struct FComp {
+        auto operator()(auto x, auto y) const {
+            using T = std::common_type_t<decltype(x), decltype(y)>;
+            if ((IsNaN(x) && IsNaN(y)) || (x == y)) {
+                return std::numeric_limits<T>::quiet_NaN();
+            }
+            if (IsNaN(x)) { return T{0}; }
+            if (IsNaN(y)) { return T{1}; }
+            return static_cast<T>(Compare{}(T{x}, T{y}));
+        }
+    };
+} // namespace detail
+
+    template<typename T, std::size_t S>
+    auto Add(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> /*primal*/, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1};
+    }
+
+    template<typename T, std::size_t S>
+    auto Mul(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        Col(trace, j) = Col(primal, i) / Col(primal, j);
+    }
+
+    template<typename T, std::size_t S>
+    auto Sub(std::vector<Operon::Node> const& nodes, Backend::View<T const, S> /*primal*/, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        auto v = (nodes[i].Arity == 1 || j < i-1) ? T{-1} : T{+1};
+        Col(trace, j) = v;
+    }
+
+    template<typename T, std::size_t S>
+    auto Div(std::vector<Operon::Node> const& nodes, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        auto const& n = nodes[i];
+        if (n.Arity == 1) {
+            Col(trace, j) = -T{1} / (Col(primal, j) * Col(primal, j));
+        } else {
+            Col(trace, j) = (j == i-1 ? T{1} : T{-1}) * Col(primal, i) / Col(primal, j);
+        }
+    }
+
+    template<typename T, std::size_t S>
+    auto Aq(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        if (j == i-1) {
+            Col(trace, j) = Col(primal, i) / Col(primal, j);
+        } else {
+            Col(trace, j) = -Col(primal, j) * blaze::pow(Col(primal, i), T{3}) / (Col(primal, i-1) * Col(primal, i-1));
+        }
+    }
+
+    template<typename T, std::size_t S>
+    auto Pow(std::vector<Operon::Node> const& nodes, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        if (j == i-1) {
+            auto const k = j - (nodes[j].Length + 1);
+            Col(trace, j) = Col(primal, i) * Col(primal, k) / Col(primal, j);
+        } else {
+            auto const k = i-1;
+            Col(trace, j) = Col(primal, i) * blaze::log(Col(primal, k));
+        }
+    }
+
+    template<typename T, std::size_t S>
+    auto Min(std::vector<Operon::Node> const& nodes, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        auto k = j == i - 1 ? (j - nodes[j].Length - 1) : i - 1;
+        auto const* a = Ptr(primal, j);
+        std::transform(a, a + S, Ptr(primal, k), Ptr(trace, j), detail::FComp<std::less<>>{});
+    }
+
+    template<typename T, std::size_t S>
+    auto Max(std::vector<Operon::Node> const& nodes, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        auto k = j == i - 1 ? (j - nodes[j].Length - 1) : i - 1;
+        auto const* a = Ptr(primal, j);
+        std::transform(a, a + S, Ptr(primal, k), Ptr(trace, j), detail::FComp<std::greater<>>{});
+    }
+
+    template<typename T, std::size_t S>
+    auto Square(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{2} * Col(primal, j);
+    }
+
+    template<typename T, std::size_t S>
+    auto Abs(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = blaze::sign(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Ceil(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = blaze::ceil(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Floor(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = blaze::floor(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Exp(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        // fmt::print("dexp 1: {} => {}\n", std::span{Ptr(trace, j), S}, std::span{Ptr(primal, i), S});
+        // Col(trace, j) = Col(primal, i);
+        std::copy_n(Ptr(primal, i), S, Ptr(trace, j));
+        // fmt::print("dexp 2: {} => {}\n", std::span{Ptr(trace, j), S}, std::span{Ptr(primal, i), S});
+    }
+
+    template<typename T, std::size_t S>
+    auto Log(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1} / (Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Log1p(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1} / (T{1} + Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Logabs(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = blaze::sign(Col(primal, j)) / blaze::abs(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sin(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = blaze::cos(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cos(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = -blaze::sin(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Tan(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        blaze::DynamicVector<T> t = blaze::tan(Col(primal, j));
+        Col(trace, j) = T{1} + t*t;
+    }
+
+    template<typename T, std::size_t S>
+    auto Sinh(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = blaze::cosh(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cosh(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = blaze::sinh(Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Tanh(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1} - blaze::pow(blaze::tanh(Col(primal, j)), T{2});
+    }
+
+    template<typename T, std::size_t S>
+    auto Asin(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1} / (blaze::sqrt(T{1} - Col(primal, j) * Col(primal, j)));
+    }
+
+    template<typename T, std::size_t S>
+    auto Acos(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = -T{1} / (blaze::sqrt(T{1} - Col(primal, j) * Col(primal, j)));
+    }
+
+    template<typename T, std::size_t S>
+    auto Atan(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto /*i*/, std::integral auto j) {
+        Col(trace, j) = T{1} / (T{1} + Col(primal, j) * Col(primal, j));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sqrt(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        Col(trace, j) = T{1} / (T{2} * Col(primal, i));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sqrtabs(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        Col(trace, j) = blaze::sign(Col(primal, j)) / (T{2} * Col(primal, i));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cbrt(std::vector<Operon::Node> const& /*nodes*/, Backend::View<T const, S> primal, Backend::View<T> trace, std::integral auto i, std::integral auto j) {
+        Col(trace, j) = T{1} / (T{3} * Col(primal, i) * Col(primal, i));
+    }
+}  // namespace Operon::Backend
+
+#endif
diff --git a/include/operon/interpreter/backend/blaze/functions.hpp b/include/operon/interpreter/backend/blaze/functions.hpp
new file mode 100644
index 0000000..496a117
--- /dev/null
+++ b/include/operon/interpreter/backend/blaze/functions.hpp
@@ -0,0 +1,215 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright 2019-2023 Heal Research
+
+#ifndef OPERON_BACKEND_BLAZE_FUNCTIONS_HPP
+#define OPERON_BACKEND_BLAZE_FUNCTIONS_HPP
+
+#include <blaze/Blaze.h>
+#include "operon/interpreter/backend/backend.hpp"
+
+namespace Operon::Backend {
+    namespace detail {
+        template<typename T>
+        using CVector = blaze::CustomVector<T, blaze::unaligned, blaze::unpadded, blaze::columnVector>;
+    } // namespace detail
+
+    template<typename T, std::size_t S>
+    inline auto Map(T* res) -> detail::CVector<T> {
+        return detail::CVector<T>(res, S);
+    }
+
+    template<typename T, std::size_t S>
+    inline auto Map(T const* res) -> detail::CVector<std::remove_const_t<T>> {
+        using U = std::remove_const_t<T>;
+        return detail::CVector<U>(const_cast<U*>(res), S);
+    }
+
+    template<typename T, std::size_t S>
+    auto Col(Backend::View<T, S> view, std::integral auto col) {
+        return Map<T, S>(view.data_handle() + col * S);
+    }
+
+    template<typename T, std::size_t S>
+    auto Col(Backend::View<T const, S> view, std::integral auto col) {
+        return Map<T, S>(view.data_handle() + col * S);
+    }
+
+    // utility
+    template<typename T, std::size_t S>
+    auto Fill(T* res, T value) {
+        std::fill_n(res, S, value);
+    }
+
+    template<typename T, std::size_t S>
+    auto Fill(T* res, int n, T value) {
+        std::fill_n(res, n, value);
+    }
+
+    // n-ary functions
+    template<typename T, std::size_t S>
+    auto Add(T* res, auto const*... args) {
+        Map<T, S>(res) = (Map<T const, S>(args) + ...);
+    }
+
+    template<typename T, std::size_t S>
+    auto Mul(T* res, auto const*... args) {
+        Map<T, S>(res) = (Map<T const, S>(args) * ...);
+    }
+
+    template<typename T, std::size_t S>
+    auto Sub(T* res, auto* first, auto const*... rest) {
+        static_assert(sizeof...(rest) > 0);
+        Map<T, S>(res) = Map<T const, S>(first) - (Map<T const, S>(rest) + ...);
+    }
+
+    template<typename T, std::size_t S>
+    auto Div(T* res, auto const* first, auto const*... rest) {
+        static_assert(sizeof...(rest) > 0);
+        Map<T, S>(res) = Map<T const, S>(first) / (Map<T const, S>(rest) * ...);
+    }
+
+    template<typename T, std::size_t S>
+    auto Min(T* res, auto const* first, auto const*... args) {
+        static_assert(sizeof...(args) > 0);
+        for (auto i = 0U; i < S; ++i) {
+            res[i] = std::min({first[i], args[i]...});
+        }
+    }
+
+    template<typename T, std::size_t S>
+    auto Max(T* res, auto* first, auto const*... args) {
+        static_assert(sizeof...(args) > 0);
+        for (auto i = 0U; i < S; ++i) {
+            res[i] = std::max({first[i], args[i]...});
+        }
+    }
+
+    // binary functions
+    template<typename T, std::size_t S>
+    auto Aq(T* res, T const* a, T const* b) {
+        blaze::DynamicVector<T> x = Map<T const, S>(b);
+        Map<T, S>(res) = Map<T const, S>(a) / blaze::sqrt(T{1} + x * x);
+    }
+
+    template<typename T, std::size_t S>
+    auto Pow(T* res, T const* a, T const* b) {
+        Map<T, S>(res) = blaze::pow(Map<T const, S>(a), Map<T const, S>(b));
+    }
+
+    // unary functions
+    template<typename T, std::size_t S>
+    auto Cpy(T* res, T const* arg) {
+        Map<T, S>(res) = Map<T const, S>(arg);
+    }
+
+    template<typename T, std::size_t S>
+    auto Neg(T* res, T const* arg) {
+        Map<T, S>(res) = -Map<T const, S>(arg);
+    }
+
+    template<typename T, std::size_t S>
+    auto Inv(T* res, T const* arg) {
+        Map<T, S>(res) = T{1} / Map<T const, S>(arg);
+    }
+
+    template<typename T, std::size_t S>
+    auto Abs(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::abs(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Ceil(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::ceil(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Floor(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::floor(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Square(T* res, T const* arg) {
+        std::transform(arg, arg+S, res, [](auto x) { return x * x; });
+    }
+
+    template<typename T, std::size_t S>
+    auto Exp(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::exp(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Log(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::log(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Log1p(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::log1p(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Logabs(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::log(blaze::abs(Map<T const, S>(arg)));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sin(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::sin(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cos(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::cos(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Tan(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::tan(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Asin(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::asin(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Acos(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::acos(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Atan(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::atan(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sinh(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::sinh(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cosh(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::cosh(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Tanh(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::tanh(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sqrt(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::sqrt(Map<T const, S>(arg));
+    }
+
+    template<typename T, std::size_t S>
+    auto Sqrtabs(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::sqrt(blaze::abs(Map<T const, S>(arg)));
+    }
+
+    template<typename T, std::size_t S>
+    auto Cbrt(T* res, T const* arg) {
+        Map<T, S>(res) = blaze::cbrt(Map<T const, S>(arg));
+    }
+} // namespace Operon::Backend
+#endif
diff --git a/include/operon/interpreter/functions.hpp b/include/operon/interpreter/functions.hpp
index 37bc9c8..99e3f60 100644
--- a/include/operon/interpreter/functions.hpp
+++ b/include/operon/interpreter/functions.hpp
@@ -9,12 +9,18 @@
 #include "operon/interpreter/backend/eigen.hpp"
 #elif defined(OPERON_MATH_EVE)
 #include "operon/interpreter/backend/eve.hpp"
+#elif defined(OPERON_MATH_ARMA)
+#include "operon/interpreter/backend/arma.hpp"
+#elif defined(OPERON_MATH_BLAZE)
+#include "operon/interpreter/backend/blaze.hpp"
 #elif defined(OPERON_MATH_FASTOR)
 #include "operon/interpreter/backend/fastor.hpp"
 #elif defined(OPERON_MATH_STL)
 #include "operon/interpreter/backend/plain.hpp"
 #elif defined(OPERON_MATH_VDT)
 #include "operon/interpreter/backend/vdt.hpp"
+#elif defined(OPERON_MATH_XTENSOR)
+#include "operon/interpreter/backend/xtensor.hpp"
 #elif defined(OPERON_MATH_FAST_V1)
 #include "operon/interpreter/backend/fast_v1.hpp"
 #elif defined(OPERON_MATH_FAST_V2)
@@ -27,7 +33,8 @@ namespace Operon {
     // utility
     template<typename T, std::size_t S>
     auto Fill(Backend::View<T, S> view, int idx, T value) {
-        Backend::Fill<T, S>(view.data_handle() + idx * S, value);
+        auto* p = view.data_handle() + idx * S;
+        std::fill_n(p, S, value);
     };
 
     // detect missing specializations