diff --git a/NEWS.md b/NEWS.md
index ba31e4d..3e583cc 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,6 @@
+# boostmath 1.5.0
+  * Functions (excluding distributions and polynomials) now support vectorised inputs
+
 # boostmath 1.4.0
   * Added `logit` and `logistic_sigmoid` functions
   * Added `log_pi` constant
diff --git a/src/boostmath/as_cpp.hpp b/src/boostmath/as_cpp.hpp
new file mode 100644
index 0000000..f812ed2
--- /dev/null
+++ b/src/boostmath/as_cpp.hpp
@@ -0,0 +1,125 @@
+#ifndef BOOSTMATH_AS_CPP_HPP
+#define BOOSTMATH_AS_CPP_HPP
+
+#include <cpp11.hpp>
+#include <complex>
+#include "./type_traits.hpp"
+
+namespace boostmath {
+  namespace internal {
+    template <typename T, std::size_t... Is, std::enable_if_t<std::negation<is_integral_scalar_type<T>>::value>* = nullptr>
+    T as_cpp_arithmetic_container_impl(SEXP x, std::index_sequence<Is...>) {
+      return {REAL_ELT(x, Is)...};
+    }
+
+    template <typename T, std::size_t... Is, std::enable_if_t<is_integral_scalar_type<T>::value>* = nullptr>
+    T as_cpp_arithmetic_container_impl(SEXP x, std::index_sequence<Is...>) {
+      return {INTEGER_ELT(x, Is)...};
+    }
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::complex<double>>::value>* = nullptr>
+  inline std::complex<double> as_cpp(SEXP x) {
+    Rcomplex r_complex = COMPLEX_ELT(x, 0);
+    return {r_complex.r, r_complex.i};
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<std::complex<double>>>::value>* = nullptr>
+  inline std::vector<std::complex<double>> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::vector<std::complex<double>> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      Rcomplex r_complex = COMPLEX_ELT(x, i);
+      result[i] = {r_complex.r, r_complex.i};
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<uint64_t>>::value>* = nullptr>
+  inline std::vector<uint64_t> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::vector<uint64_t> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      result[i] = static_cast<uint64_t>(Rf_isInteger(x) ? INTEGER_ELT(x, i) : REAL_ELT(x, i));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<int64_t>>::value>* = nullptr>
+  inline std::vector<int64_t> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::vector<int64_t> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      result[i] = static_cast<int64_t>(Rf_isInteger(x) ? INTEGER_ELT(x, i) : REAL_ELT(x, i));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<int>>::value>* = nullptr>
+  inline std::vector<int> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::vector<int> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      result[i] = static_cast<int>(Rf_isInteger(x) ? INTEGER_ELT(x, i) : REAL_ELT(x, i));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<double>>::value>* = nullptr>
+  inline std::vector<double> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::vector<double> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      result[i] = static_cast<double>(Rf_isInteger(x) ? INTEGER_ELT(x, i) : REAL_ELT(x, i));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::list<double>>::value>* = nullptr>
+  inline std::list<double> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::list<double> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      result.push_back(static_cast<double>(Rf_isInteger(x) ? INTEGER_ELT(x, i) : REAL_ELT(x, i)));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<std::vector<double>>>::value>* = nullptr>
+  inline std::vector<std::vector<double>> as_cpp(SEXP x) {
+    if (Rf_isNewList(x)) {
+      const R_xlen_t n = Rf_xlength(x);
+      std::vector<std::vector<double>> result(n);
+      for (R_xlen_t i = 0; i < n; ++i) {
+        SEXP xi = VECTOR_ELT(x, i);
+        result[i] = as_cpp<std::vector<double>>(xi);
+      }
+      return result;
+    }
+    return {as_cpp<std::vector<double>>(x)};
+  }
+
+  template <typename T, std::enable_if_t<is_arithmetic_container<T>::value>* = nullptr>
+  inline T as_cpp(SEXP x) {
+    return internal::as_cpp_arithmetic_container_impl<T>(x, std::make_index_sequence<std::tuple_size<T>::value>{});
+  }
+
+  template <typename T, std::enable_if_t<is_vector_of_arithmetic_containers<T>::value>* = nullptr>
+  inline T as_cpp(SEXP x) {
+    using array_type = typename T::value_type;
+    R_xlen_t n = Rf_xlength(x);
+    T result;
+    result.reserve(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      SEXP xi = VECTOR_ELT(x, i);
+      result.push_back(as_cpp<array_type>(xi));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<is_cpp11<T>::value>* = nullptr>
+  inline T as_cpp(SEXP x) {
+    return cpp11::as_cpp<T>(x);
+  }
+}
+#endif
diff --git a/src/boostmath/as_sexp.hpp b/src/boostmath/as_sexp.hpp
new file mode 100644
index 0000000..5fcb2ef
--- /dev/null
+++ b/src/boostmath/as_sexp.hpp
@@ -0,0 +1,152 @@
+#ifndef BOOSTMATH_AS_SEXP_HPP
+#define BOOSTMATH_AS_SEXP_HPP
+
+#include <cpp11.hpp>
+#include <complex>
+#include "./type_traits.hpp"
+#include "Rinternals.h"
+
+namespace boostmath {
+  namespace internal {
+    template <typename T, std::size_t... Is, std::enable_if_t<std::negation<is_integral_scalar_type<T>>::value>* = nullptr>
+    auto as_sexp_arithmetic_container_impl(const T& x, std::index_sequence<Is...>) {
+      SEXP data = cpp11::safe[Rf_allocVector](REALSXP, sizeof...(Is));
+      double* ptr = REAL(data);
+      ((ptr[Is] = std::get<Is>(x)), ...);
+      return data;
+    }
+
+    template <typename T, std::size_t... Is,  std::enable_if_t<is_integral_scalar_type<T>::value>* = nullptr>
+    auto as_sexp_arithmetic_container_impl(const T& x, std::index_sequence<Is...>) {
+      SEXP data = cpp11::safe[Rf_allocVector](INTSXP, sizeof...(Is));
+      int* ptr = INTEGER(data);
+      ((ptr[Is] = std::get<Is>(x)), ...);
+      return data;
+    }
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::complex<double>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    SEXP data = cpp11::safe[Rf_allocVector](CPLXSXP, 1);
+    Rcomplex* r_complex = COMPLEX(data);
+    r_complex->r = x.real();
+    r_complex->i = x.imag();
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<std::complex<double>>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](CPLXSXP, n);
+    Rcomplex* r_complex = COMPLEX(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      r_complex[i].r = x[i].real();
+      r_complex[i].i = x[i].imag();
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::list<double>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](REALSXP, n);
+    double* ptr = REAL(data);
+    size_t i = 0;
+    for (const auto& value : x) {
+      ptr[i++] = value;
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<std::list<double>>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](VECSXP, n);
+    PROTECT(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      SET_VECTOR_ELT(data, i, as_sexp(x[i]));
+    }
+    UNPROTECT(1);
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<uint64_t>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](INTSXP, n);
+    int* ptr = INTEGER(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      ptr[i] = static_cast<int>(x[i]);
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<int64_t>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](INTSXP, n);
+    int* ptr = INTEGER(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      ptr[i] = static_cast<int>(x[i]);
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<int>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](INTSXP, n);
+    int* ptr = INTEGER(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      ptr[i] = static_cast<int>(x[i]);
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<double>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](REALSXP, n);
+    double* ptr = REAL(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      ptr[i] = static_cast<double>(x[i]);
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<std::vector<double>>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](VECSXP, n);
+    PROTECT(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      SET_VECTOR_ELT(data, i, as_sexp(x[i]));
+    }
+    UNPROTECT(1);
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<is_arithmetic_container<T>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    return internal::as_sexp_arithmetic_container_impl<T>(x, std::make_index_sequence<std::tuple_size<T>::value>{});
+  }
+
+  template <typename T, std::enable_if_t<is_vector_of_arithmetic_containers<T>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](VECSXP, n);
+    PROTECT(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      SET_VECTOR_ELT(data, i, as_sexp(x[i]));
+    }
+    UNPROTECT(1);
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<is_cpp11<T>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    return cpp11::as_sexp(x);
+  }
+}
+
+#endif // BOOSTMATH_SEXP_HPP
diff --git a/src/boostmath/boost_function.hpp b/src/boostmath/boost_function.hpp
new file mode 100644
index 0000000..fe9a7ac
--- /dev/null
+++ b/src/boostmath/boost_function.hpp
@@ -0,0 +1,58 @@
+#ifndef BOOSTMATH_UTILS_HPP
+#define BOOSTMATH_UTILS_HPP
+
+#include <tuple>
+#include "as_sexp.hpp"
+#include "scalar_container.hpp"
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+namespace boostmath {
+  namespace internal {
+    template <typename TypesTuple, typename ArgsTupleT, std::size_t... Is>
+    auto make_containers_tuple_impl(ArgsTupleT&& args, std::index_sequence<Is...>) {
+      return std::make_tuple(boostmath::scalar_container<std::tuple_element_t<Is, TypesTuple>>(std::get<Is>(args))...);
+    }
+  }
+
+  template <typename TypesTupleT, typename ArgsTupleT>
+  auto make_containers_tuple(ArgsTupleT&& args) {
+    constexpr auto tuple_index_seq = std::make_index_sequence<std::tuple_size_v<std::decay_t<ArgsTupleT>>>{};
+    return internal::make_containers_tuple_impl<TypesTupleT>(std::forward<ArgsTupleT>(args), tuple_index_seq);
+  }
+
+  template <typename F, typename... ArgsT>
+  struct invoke_rtn {
+    using type = std::invoke_result_t<F, ArgsT...>;
+  };
+
+  template <typename F, typename... ArgsT>
+  struct invoke_rtn<F, std::tuple<ArgsT...>> {
+    using type = std::invoke_result_t<F, ArgsT...>;
+  };
+
+  template <typename F, typename... ArgsT>
+  using invoke_rtn_t = typename invoke_rtn<F, ArgsT...>::type;
+
+  template <typename TypesTuple, typename F, typename ArgsTupleT>
+  SEXP boost_function(F&& f, ArgsTupleT&& args) {
+    const auto containers_tuple = make_containers_tuple<TypesTuple>(std::forward<ArgsTupleT>(args));
+
+    using containers_tuple_t = decltype(containers_tuple);
+    size_t max_n = std::apply([](auto&&... conts) { return std::max({conts.size()...}); },
+                            std::forward<containers_tuple_t>(containers_tuple));
+
+    std::vector<invoke_rtn_t<F, TypesTuple>> results(max_n);
+    for (size_t i = 0; i < max_n; ++i) {
+      results[i] = std::apply([&](auto&&... conts){ return f(conts[i]...); },
+                              std::forward<containers_tuple_t>(containers_tuple));
+    }
+    if (max_n == 1) {
+      return boostmath::as_sexp(results[0]);
+    }
+    return boostmath::as_sexp(results);
+  }
+}
+
+#endif // BOOSTMATH_UTILS_HPP
diff --git a/src/boostmath/macros.hpp b/src/boostmath/macros.hpp
index 566e099..f867650 100644
--- a/src/boostmath/macros.hpp
+++ b/src/boostmath/macros.hpp
@@ -2,7 +2,9 @@
 #define BOOSTMATH_MACROS_HPP
 
 #include <cpp11/declarations.hpp>
-#include "sexp.hpp"
+#include "as_sexp.hpp"
+#include "as_cpp.hpp"
+#include "boost_function.hpp"
 
 #define NOARG_BOOST_MEMBER(name) \
   extern "C" SEXP name##_() { \
@@ -21,166 +23,161 @@
 #define UNARY_BOOST_FUNCTION(name, arg_type) \
   extern "C" SEXP name##_(SEXP x_) { \
     BEGIN_CPP11 \
-    const arg_type x = boostmath::as_cpp<arg_type>(x_); \
-    return boostmath::as_sexp(boost::math::name<double>(x)); \
+    using args_type = std::tuple<arg_type>; \
+    const auto f = [](auto&&... args) { return boost::math::name<double>(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_)); \
     END_CPP11 \
   }
 
 #define UNARY_BOOST_FUNCTION_NOTYPE(name, arg_type) \
   extern "C" SEXP name##_(SEXP x_) { \
     BEGIN_CPP11 \
-    const arg_type x = boostmath::as_cpp<arg_type>(x_); \
-    return boostmath::as_sexp(boost::math::name(x)); \
+    using args_type = std::tuple<arg_type>; \
+    const auto f = [](auto&&... args) { return boost::math::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_)); \
     END_CPP11 \
   }
 
 #define UNARY_BOOST_FUNCTION_NAMESPACE(namespace, name, arg_type) \
   extern "C" SEXP name##_(SEXP x_) { \
     BEGIN_CPP11 \
-    arg_type x = boostmath::as_cpp<arg_type>(x_); \
-    return boostmath::as_sexp(boost::math::namespace::name(x)); \
+    using args_type = std::tuple<arg_type>; \
+    const auto f = [](auto&&... args) { return boost::math::namespace::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_)); \
     END_CPP11 \
   }
 
 #define UNARY_BOOST_FUNCTION_SUFFIX(name, suffix, arg_type) \
   extern "C" SEXP name##_##suffix(SEXP x_) { \
     BEGIN_CPP11 \
-    const arg_type x = boostmath::as_cpp<arg_type>(x_); \
-    return boostmath::as_sexp(boost::math::name(x)); \
+    using args_type = std::tuple<arg_type>; \
+    const auto f = [](auto&&... args) { return boost::math::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_)); \
     END_CPP11 \
   }
 
 #define BINARY_BOOST_FUNCTION(name, arg1_type, arg2_type) \
   extern "C" SEXP name##_(SEXP x_, SEXP y_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    return boostmath::as_sexp(boost::math::name(x, y)); \
+    using args_type = std::tuple<arg1_type, arg2_type>; \
+    const auto f = [](auto&&... args) { return boost::math::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_)); \
     END_CPP11 \
   }
 
 #define BINARY_BOOST_FUNCTION_TYPE(name, template_type, arg1_type, arg2_type) \
   extern "C" SEXP name##_(SEXP x_, SEXP y_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    return boostmath::as_sexp(boost::math::name<template_type>(x, y)); \
+    using args_type = std::tuple<arg1_type, arg2_type>; \
+    const auto f = [](auto&&... args) { return boost::math::name<template_type>(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_)); \
     END_CPP11 \
   }
 
 #define BINARY_BOOST_FUNCTION_NAMESPACE(namespace, name, arg1_type, arg2_type) \
   extern "C" SEXP name##_(SEXP x_, SEXP y_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    return boostmath::as_sexp(boost::math::namespace::name(x, y)); \
+    using args_type = std::tuple<arg1_type, arg2_type>; \
+    const auto f = [](auto&&... args) { return boost::math::namespace::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_)); \
     END_CPP11 \
   }
 
 #define BINARY_BOOST_FUNCTION_VECTOR_SUFFIX(name, suffix, arg1_type, arg2_type) \
   extern "C" SEXP name##_##suffix(SEXP x_, SEXP y_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    std::vector<double> results; \
-    boost::math::name<double>(x, y, std::back_inserter(results)); \
-    return boostmath::as_sexp(results); \
+    using args_type = std::tuple<arg1_type, arg2_type>; \
+    const auto f = [](auto&&... args) { \
+      std::vector<double> results; \
+      boost::math::name<double>(args..., std::back_inserter(results)); \
+      return results; \
+    }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_)); \
     END_CPP11 \
   }
 
 #define BINARY_BOOST_FUNCTION_SUFFIX(name, suffix, arg1_type, arg2_type) \
   extern "C" SEXP name##_##suffix(SEXP x_, SEXP y_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    return boostmath::as_sexp(boost::math::name(x, y)); \
+    using args_type = std::tuple<arg1_type, arg2_type>; \
+    const auto f = [](auto&&... args) { return boost::math::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_)); \
     END_CPP11 \
   }
 
 #define TERNARY_BOOST_FUNCTION(name, arg1_type, arg2_type, arg3_type) \
   extern "C" SEXP name##_(SEXP x_, SEXP y_, SEXP z_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    const arg3_type z = boostmath::as_cpp<arg3_type>(z_); \
-    return boostmath::as_sexp(boost::math::name(x, y, z)); \
+    using args_type = std::tuple<arg1_type, arg2_type, arg3_type>; \
+    const auto f = [](auto&&... args) { return boost::math::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_, z_)); \
     END_CPP11 \
   }
 
 #define TERNARY_BOOST_FUNCTION_NAMESPACE(namespace, name, arg1_type, arg2_type, arg3_type) \
   extern "C" SEXP name##_(SEXP x_, SEXP y_, SEXP z_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    const arg3_type z = boostmath::as_cpp<arg3_type>(z_); \
-    return boostmath::as_sexp(boost::math::namespace::name(x, y, z)); \
+    using args_type = std::tuple<arg1_type, arg2_type, arg3_type>; \
+    const auto f = [](auto&&... args) { return boost::math::namespace::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_, z_)); \
     END_CPP11 \
   }
 
 #define TERNARY_BOOST_FUNCTION_VECTOR_SUFFIX(name, suffix, arg1_type, arg2_type, arg3_type) \
   extern "C" SEXP name##_##suffix(SEXP x_, SEXP y_, SEXP z_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    const arg3_type z = boostmath::as_cpp<arg3_type>(z_); \
-    std::vector<double> results; \
-    boost::math::name<double>(x, y, z, std::back_inserter(results)); \
-    return boostmath::as_sexp(results); \
+    using args_type = std::tuple<arg1_type, arg2_type, arg3_type>; \
+    const auto f = [](auto&&... args) { \
+      std::vector<double> results; \
+      boost::math::name<double>(args..., std::back_inserter(results)); \
+      return results; \
+    }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_, z_)); \
     END_CPP11 \
   }
 
 #define TERNARY_BOOST_FUNCTION_SUFFIX(name, suffix, arg1_type, arg2_type, arg3_type) \
   extern "C" SEXP name##_##suffix(SEXP x_, SEXP y_, SEXP z_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    const arg3_type z = boostmath::as_cpp<arg3_type>(z_); \
-    return boostmath::as_sexp(boost::math::name(x, y, z)); \
+    using args_type = std::tuple<arg1_type, arg2_type, arg3_type>; \
+    const auto f = [](auto&&... args) { return boost::math::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_, z_)); \
     END_CPP11 \
   }
 
 #define QUARTERNARY_BOOST_FUNCTION(name, arg1_type, arg2_type, arg3_type, arg4_type) \
   extern "C" SEXP name##_(SEXP x_, SEXP y_, SEXP z_, SEXP w_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    const arg3_type z = boostmath::as_cpp<arg3_type>(z_); \
-    const arg4_type w = boostmath::as_cpp<arg4_type>(w_); \
-    return boostmath::as_sexp(boost::math::name(x, y, z, w)); \
+    using args_type = std::tuple<arg1_type, arg2_type, arg3_type, arg4_type>; \
+    const auto f = [](auto&&... args) { return boost::math::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_, z_, w_)); \
     END_CPP11 \
   }
 
 #define QUARTERNARY_BOOST_FUNCTION_NAMESPACE(namespace, name, arg1_type, arg2_type, arg3_type, arg4_type) \
   extern "C" SEXP name##_(SEXP x_, SEXP y_, SEXP z_, SEXP w_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    const arg3_type z = boostmath::as_cpp<arg3_type>(z_); \
-    const arg4_type w = boostmath::as_cpp<arg4_type>(w_); \
-    return boostmath::as_sexp(boost::math::namespace::name(x, y, z, w)); \
+    using args_type = std::tuple<arg1_type, arg2_type, arg3_type, arg4_type>; \
+    const auto f = [](auto&&... args) { return boost::math::namespace::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_, z_, w_)); \
     END_CPP11 \
   }
 
 #define QUARTERNARY_BOOST_FUNCTION_NAMESPACE_SUFFIX(namespace, name, suffix, arg1_type, arg2_type, arg3_type, arg4_type) \
   extern "C" SEXP name##_##suffix(SEXP x_, SEXP y_, SEXP z_, SEXP w_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    const arg3_type z = boostmath::as_cpp<arg3_type>(z_); \
-    const arg4_type w = boostmath::as_cpp<arg4_type>(w_); \
-    return boostmath::as_sexp(boost::math::namespace::name(x, y, z, w)); \
+    using args_type = std::tuple<arg1_type, arg2_type, arg3_type, arg4_type>; \
+    const auto f = [](auto&&... args) { return boost::math::namespace::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_, z_, w_)); \
     END_CPP11 \
   }
 
 #define PENTA_BOOST_FUNCTION_NAMESPACE(namespace, name, arg1_type, arg2_type, arg3_type, arg4_type, arg5_type) \
   extern "C" SEXP name##_(SEXP x_, SEXP y_, SEXP z_, SEXP w_, SEXP v_) { \
     BEGIN_CPP11 \
-    const arg1_type x = boostmath::as_cpp<arg1_type>(x_); \
-    const arg2_type y = boostmath::as_cpp<arg2_type>(y_); \
-    const arg3_type z = boostmath::as_cpp<arg3_type>(z_); \
-    const arg4_type w = boostmath::as_cpp<arg4_type>(w_); \
-    const arg4_type v = boostmath::as_cpp<arg4_type>(v_); \
-    return boostmath::as_sexp(boost::math::namespace::name(x, y, z, w, v)); \
+    using args_type = std::tuple<arg1_type, arg2_type, arg3_type, arg4_type, arg5_type>; \
+    const auto f = [](auto&&... args) { return boost::math::namespace::name(args...); }; \
+    return boostmath::boost_function<args_type>(f, std::forward_as_tuple(x_, y_, z_, w_, v_)); \
     END_CPP11 \
   }
 
diff --git a/src/boostmath/scalar_container.hpp b/src/boostmath/scalar_container.hpp
new file mode 100644
index 0000000..b964ce1
--- /dev/null
+++ b/src/boostmath/scalar_container.hpp
@@ -0,0 +1,27 @@
+#ifndef BOOSTMATH_SCALAR_CONTAINER_HPP
+#define BOOSTMATH_SCALAR_CONTAINER_HPP
+
+#include <cstddef>
+#include <vector>
+#include "as_cpp.hpp"
+
+namespace boostmath {
+  template <typename T>
+  struct scalar_container {
+    std::vector<T> value_;
+    size_t size_;
+    bool is_scalar_;
+
+    scalar_container() = default;
+    scalar_container(SEXP x) :
+      value_(boostmath::as_cpp<std::vector<T>>(x)),
+      size_(value_.size()),
+      is_scalar_(size_ == 1) {}
+
+    T operator[](size_t i) const { return is_scalar_ ? value_[0] : value_[i]; }
+    size_t size() const { return size_; }
+
+  };
+}
+
+#endif
diff --git a/src/boostmath/sexp.hpp b/src/boostmath/sexp.hpp
deleted file mode 100644
index 8f41f88..0000000
--- a/src/boostmath/sexp.hpp
+++ /dev/null
@@ -1,97 +0,0 @@
-#ifndef BOOSTMATH_SEXP_HPP
-#define BOOSTMATH_SEXP_HPP
-
-#include <cpp11.hpp>
-#include <complex>
-#include "./type_traits.hpp"
-
-namespace boostmath {
-  namespace internal {
-    template <typename T, std::size_t... Is, std::enable_if_t<std::negation<is_integral_scalar_type<T>>::value>* = nullptr>
-    T as_cpp_arithmetic_container_impl(SEXP x, std::index_sequence<Is...>) {
-      return {REAL_ELT(x, Is)...};
-    }
-
-    template <typename T, std::size_t... Is, std::enable_if_t<is_integral_scalar_type<T>::value>* = nullptr>
-    T as_cpp_arithmetic_container_impl(SEXP x, std::index_sequence<Is...>) {
-      return {INTEGER_ELT(x, Is)...};
-    }
-
-    template <typename T, std::size_t... Is, std::enable_if_t<std::negation<is_integral_scalar_type<T>>::value>* = nullptr>
-    auto as_sexp_arithmetic_container_impl(const T& x, std::index_sequence<Is...>) {
-      SEXP data = cpp11::safe[Rf_allocVector](REALSXP, sizeof...(Is));
-      double* ptr = REAL(data);
-      ((ptr[Is] = std::get<Is>(x)), ...);
-      return data;
-    }
-
-    template <typename T, std::size_t... Is,  std::enable_if_t<is_integral_scalar_type<T>::value>* = nullptr>
-    auto as_sexp_arithmetic_container_impl(const T& x, std::index_sequence<Is...>) {
-      SEXP data = cpp11::safe[Rf_allocVector](INTSXP, sizeof...(Is));
-      int* ptr = INTEGER(data);
-      ((ptr[Is] = std::get<Is>(x)), ...);
-      return data;
-    }
-  }
-
-  template <typename T, std::enable_if_t<std::is_same<T, std::complex<double>>::value>* = nullptr>
-  inline std::complex<double> as_cpp(SEXP x) {
-    Rcomplex r_complex = COMPLEX_ELT(x, 0);
-    return {r_complex.r, r_complex.i};
-  }
-
-  template <typename T, std::enable_if_t<is_arithmetic_container<T>::value>* = nullptr>
-  inline T as_cpp(SEXP x) {
-    return internal::as_cpp_arithmetic_container_impl<T>(x, std::make_index_sequence<std::tuple_size<T>::value>{});
-  }
-
-  template <typename T, std::enable_if_t<is_vector_of_arithmetic_containers<T>::value>* = nullptr>
-  inline T as_cpp(SEXP x) {
-    using array_type = typename T::value_type;
-    R_xlen_t n = Rf_xlength(x);
-    T result;
-    result.reserve(n);
-    for (R_xlen_t i = 0; i < n; ++i) {
-      SEXP xi = VECTOR_ELT(x, i);
-      result.push_back(as_cpp<array_type>(xi));
-    }
-    return result;
-  }
-
-  template <typename T, std::enable_if_t<is_cpp11<T>::value>* = nullptr>
-  inline T as_cpp(SEXP x) {
-    return cpp11::as_cpp<T>(x);
-  }
-
-  template <typename T, std::enable_if_t<std::is_same<T, std::complex<double>>::value>* = nullptr>
-  inline SEXP as_sexp(const T& x) {
-    SEXP data = cpp11::safe[Rf_allocVector](CPLXSXP, 1);
-    Rcomplex* r_complex = COMPLEX(data);
-    r_complex->r = x.real();
-    r_complex->i = x.imag();
-    return data;
-  }
-
-  template <typename T, std::enable_if_t<is_arithmetic_container<T>::value>* = nullptr>
-  inline SEXP as_sexp(const T& x) {
-    return internal::as_sexp_arithmetic_container_impl<T>(x, std::make_index_sequence<std::tuple_size<T>::value>{});
-  }
-
-  template <typename T, std::enable_if_t<is_vector_of_arithmetic_containers<T>::value>* = nullptr>
-  inline SEXP as_sexp(const T& x) {
-    R_xlen_t n = x.size();
-    SEXP data = cpp11::safe[Rf_allocVector](VECSXP, n);
-    for (R_xlen_t i = 0; i < n; ++i) {
-      SET_VECTOR_ELT(data, i, as_sexp(x[i]));
-    }
-    return data;
-  }
-
-  template <typename T, std::enable_if_t<is_cpp11<T>::value>* = nullptr>
-  inline SEXP as_sexp(const T& x) {
-    return cpp11::as_sexp(x);
-  }
-
-}
-
-#endif // BOOSTMATH_SEXP_HPP
diff --git a/src/boostmath/type_traits.hpp b/src/boostmath/type_traits.hpp
index cda274b..ce9c488 100644
--- a/src/boostmath/type_traits.hpp
+++ b/src/boostmath/type_traits.hpp
@@ -1,9 +1,12 @@
 #ifndef BOOSTMATH_TYPE_TRAITS_HPP
 #define BOOSTMATH_TYPE_TRAITS_HPP
 
+#include <cstdint>
+#include <sys/types.h>
 #include <type_traits>
 #include <tuple>
 #include <vector>
+#include <list>
 
 namespace boostmath {
   template <typename... Targs>
@@ -26,18 +29,33 @@ namespace boostmath {
   struct is_arithmetic_pair<std::pair<T1, T2>> : std::conjunction<std::is_arithmetic<T1>, std::is_arithmetic<T2>> {};
 
   template <typename T>
-  struct is_vector_of_arithmetic_containers : std::false_type {};
+  struct is_arithmetic_list : std::false_type {};
 
-  template <typename T, std::size_t N>
-  struct is_vector_of_arithmetic_containers<std::vector<std::array<T, N>>> : std::is_arithmetic<T> {};
+  template <typename T>
+  struct is_arithmetic_list<std::list<T>> : std::is_arithmetic<T> {};
 
   template <typename T>
   using is_arithmetic_container = std::disjunction<is_arithmetic_tuple<T>, is_arithmetic_array<T>, is_arithmetic_pair<T>>;
 
+  template <typename T>
+  struct is_vector_of_arithmetic_containers : std::false_type {};
+
+  template <typename T>
+  struct is_vector_of_arithmetic_containers<std::vector<T>> : is_arithmetic_container<T> {};
+
+
   template <typename T>
   using is_cpp11 = std::negation<
     std::disjunction<
       std::is_same<T, std::complex<double>>,
+      std::is_same<T, std::vector<std::complex<double>>>,
+      std::is_same<T, std::vector<double>>,
+      std::is_same<T, std::list<double>>,
+      std::is_same<T, std::vector<std::list<double>>>,
+      std::is_same<T, std::vector<std::vector<double>>>,
+      std::is_same<T, std::vector<int>>,
+      std::is_same<T, std::vector<uint64_t>>,
+      std::is_same<T, std::vector<int64_t>>,
       is_arithmetic_container<T>,
       is_vector_of_arithmetic_containers<T>,
       std::is_void<T>
diff --git a/src/quadrature_and_differentiation/numerical_differentiation.cpp b/src/quadrature_and_differentiation/numerical_differentiation.cpp
index 6fb1a97..d864156 100644
--- a/src/quadrature_and_differentiation/numerical_differentiation.cpp
+++ b/src/quadrature_and_differentiation/numerical_differentiation.cpp
@@ -1,7 +1,8 @@
 #include <cpp11.hpp>
 #include <cpp11/declarations.hpp>
 #include <boost/math/differentiation/finite_difference.hpp>
-#include "../boostmath/sexp.hpp"
+#include "../boostmath/as_sexp.hpp"
+#include "../boostmath/as_cpp.hpp"
 
 extern "C" {
   SEXP finite_difference_derivative_(SEXP f_, SEXP x_, SEXP order_) {
diff --git a/src/rootfinding_and_minimisation/with_derivatives.cpp b/src/rootfinding_and_minimisation/with_derivatives.cpp
index ddcac38..2684d09 100644
--- a/src/rootfinding_and_minimisation/with_derivatives.cpp
+++ b/src/rootfinding_and_minimisation/with_derivatives.cpp
@@ -1,6 +1,6 @@
 #include <cpp11/declarations.hpp>
 #include <boost/math/tools/roots.hpp>
-#include "../boostmath/sexp.hpp"
+#include "../boostmath/as_cpp.hpp"
 
 extern "C" {
   SEXP newton_raphson_iterate_(SEXP f_, SEXP guess_, SEXP min_, SEXP max_, SEXP digits_, SEXP max_iter_) {
@@ -65,4 +65,4 @@ extern "C" {
     return cpp11::as_sexp(result);
     END_CPP11
   }
-}
\ No newline at end of file
+}
diff --git a/src/special_functions/logistic_functions.cpp b/src/special_functions/logistic_functions.cpp
index 65ddeec..e50865d 100644
--- a/src/special_functions/logistic_functions.cpp
+++ b/src/special_functions/logistic_functions.cpp
@@ -3,4 +3,4 @@
 #include "../boostmath/macros.hpp"
 
 UNARY_BOOST_FUNCTION(logit, double)
-UNARY_BOOST_FUNCTION(logistic_sigmoid, double)
\ No newline at end of file
+UNARY_BOOST_FUNCTION(logistic_sigmoid, double)