diff --git a/util/math_reference.h b/util/math_reference.h
index 51fa1e15d..3ca887d34 100644
--- a/util/math_reference.h
+++ b/util/math_reference.h
@@ -115,6 +115,7 @@
 #endif
 
 namespace reference {
+
 /* two argument relational reference */
 template <typename T>
 auto isequal(T a, T b) {
@@ -406,7 +407,6 @@ T sub_sat(T x, T y) {
     }
   }
 }
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(sub_sat)
 
 /* upsample */
 uint16_t upsample(uint8_t h, uint8_t l);
@@ -449,23 +449,6 @@ struct upsample_t<int32_t> {
   using type = int64_t;
 };
 
-template <typename T, int N>
-sycl::vec<typename upsample_t<T>::type, N> upsample(
-    sycl::vec<T, N> a, sycl::vec<typename std::make_unsigned<T>::type, N> b) {
-  return sycl_cts::math::run_func_on_vector<typename upsample_t<T>::type, T, N>(
-      [](T x, T y) { return upsample(x, y); }, a, b);
-}
-// FIXME: AdaptiveCpp does not support marray
-#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
-template <typename T, size_t N>
-sycl::marray<typename upsample_t<T>::type, N> upsample(
-    sycl::marray<T, N> a,
-    sycl::marray<typename std::make_unsigned<T>::type, N> b) {
-  return sycl_cts::math::run_func_on_marray<typename upsample_t<T>::type, T, N>(
-      [](T x, T y) { return upsample(x, y); }, a, b);
-}
-#endif
-
 /* return number of non zero bits in x */
 template <typename T>
 T popcount(T x) {
@@ -474,49 +457,15 @@ T popcount(T x) {
     if (x & (1ull << i)) lz++;
   return lz;
 }
-MAKE_VEC_AND_MARRAY_VERSIONS(popcount)
 
 /* fast multiply add 24bits */
 sycl_cts::resultRef<int32_t> mad24(int32_t x, int32_t y, int32_t z);
 sycl_cts::resultRef<uint32_t> mad24(uint32_t x, uint32_t y, uint32_t z);
 
-template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> mad24(sycl::vec<T, N> a, sycl::vec<T, N> b,
-                                           sycl::vec<T, N> c) {
-  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
-      [](T x, T y, T z) { return mad24(x, y, z); }, a, b, c);
-}
-// FIXME: AdaptiveCpp does not support marray
-#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
-template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> mad24(sycl::marray<T, N> a,
-                                              sycl::marray<T, N> b,
-                                              sycl::marray<T, N> c) {
-  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
-      [](T x, T y, T z) { return mad24(x, y, z); }, a, b, c);
-}
-#endif
-
 /* fast multiply 24bits */
 sycl_cts::resultRef<int32_t> mul24(int32_t x, int32_t y);
 sycl_cts::resultRef<uint32_t> mul24(uint32_t x, uint32_t y);
 
-template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> mul24(sycl::vec<T, N> a,
-                                           sycl::vec<T, N> b) {
-  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
-      [](T x, T y) { return mul24(x, y); }, a, b);
-}
-// FIXME: AdaptiveCpp does not support marray
-#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
-template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> mul24(sycl::marray<T, N> a,
-                                              sycl::marray<T, N> b) {
-  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
-      [](T x, T y) { return mul24(x, y); }, a, b);
-}
-#endif
-
 // Common functions
 
 float degrees(float a);
@@ -536,18 +485,6 @@ struct higher_accuracy<float> {
   using type = double;
 };
 
-template <typename T, int N>
-struct higher_accuracy<sycl::vec<T, N>> {
-  using type = sycl::vec<typename higher_accuracy<T>::type, N>;
-};
-// FIXME: AdaptiveCpp does not support marray
-#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
-template <typename T, size_t N>
-struct higher_accuracy<sycl::marray<T, N>> {
-  using type = sycl::marray<typename higher_accuracy<T>::type, N>;
-};
-#endif
-
 template <typename T>
 T acos(T a) {
   return std::acos(static_cast<typename higher_accuracy<T>::type>(a));
@@ -661,6 +598,12 @@ using std::ilogb;
 using std::ldexp;
 using std::lgamma;
 
+template <typename T>
+T lgamma_r(T a, int *b) {
+  *b = (std::tgamma(a) > 0) ? 1 : -1;
+  return std::lgamma(a);
+}
+
 template <typename T>
 T log(T a) {
   return std::log(static_cast<typename higher_accuracy<T>::type>(a));
@@ -707,9 +650,7 @@ T minmag(T a, T b) {
 }
 
 using std::modf;
-
 float nan(unsigned int a);
-
 using std::nextafter;
 
 template <typename T>
@@ -739,7 +680,6 @@ T remquo(T x, T y, int *quo) {
 }
 
 using std::rint;
-MAKE_VEC_AND_MARRAY_VERSIONS(rint)
 
 template <typename T>
 T rootn(T a, int b) {
@@ -817,6 +757,17 @@ sycl::mfloat4 cross(sycl::mfloat4 p0, sycl::mfloat4 p1);
 sycl::mfloat3 cross(sycl::mfloat3 p0, sycl::mfloat3 p1);
 #endif
 
+template <typename T>
+T dot(T p0, T p1) {
+  return p0 * p1;
+}
+
+template <typename T>
+T normalize(T p) {
+  if (p < 0) return -1;
+  return 1;
+}
+
 
 #if SYCL_CTS_ENABLE_HALF_TESTS
 
@@ -943,38 +894,26 @@ sycl::mdouble3 cross(sycl::mdouble3 p0, sycl::mdouble3 p1);
 // non-templated conditional half / double overloads for correct resolution
 // with DPC++.
 
-template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> abs(sycl::vec<T, N> a) {
-  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
-      [](T x) { return abs(x); }, a);
-}
-#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
-template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> abs(sycl::marray<T, N> a) {
-  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
-      [](T x) { return abs(x); }, a);
-}
-#endif
+// Integer functions
 
 template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> abs_diff(sycl::vec<T, N> a,
-                                              sycl::vec<T, N> b) {
-  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
-      [](T x, T y) { return abs_diff(x, y); }, a, b);
+int any(sycl::vec<T, N> a) {
+  for (int i = 0; i < N; i++) {
+    if (any(a[i]) == 1) return true;
+  }
+  return false;
 }
 // FIXME: AdaptiveCpp does not support marray
 #ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
 template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> abs_diff(sycl::marray<T, N> a,
-                                                 sycl::marray<T, N> b) {
-  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
-      [](T x, T y) { return abs_diff(x, y); }, a, b);
+bool any(sycl::marray<T, N> a) {
+  for (size_t i = 0; i < N; i++) {
+    if (any(a[i]) == 1) return true;
+  }
+  return false;
 }
 #endif
 
-MAKE_VEC_AND_MARRAY_VERSIONS(acospi)
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(add_sat)
-
 template <typename T, int N>
 int all(sycl::vec<T, N> a) {
   for (int i = 0; i < N; i++) {
@@ -993,37 +932,65 @@ bool all(sycl::marray<T, N> a) {
 }
 #endif
 
-template <typename T, int N>
-int any(sycl::vec<T, N> a) {
+MAKE_VEC_AND_MARRAY_VERSIONS_3ARGS(bitselect)
+
+template <typename T, typename K, int N>
+sycl::vec<T, N> select(sycl::vec<T, N> a, sycl::vec<T, N> b,
+                       sycl::vec<K, N> c) {
+  sycl::vec<T, N> res;
   for (int i = 0; i < N; i++) {
-    if (any(a[i]) == 1) return true;
+    if (any(c[i]) == 1)
+      res[i] = b[i];
+    else
+      res[i] = a[i];
   }
-  return false;
+  return res;
 }
 // FIXME: AdaptiveCpp does not support marray
 #ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
 template <typename T, size_t N>
-bool any(sycl::marray<T, N> a) {
+sycl::marray<T, N> select(sycl::marray<T, N> a, sycl::marray<T, N> b,
+                          sycl::marray<bool, N> c) {
+  sycl::marray<T, N> res;
   for (size_t i = 0; i < N; i++) {
-    if (any(a[i]) == 1) return true;
+    res[i] = c[i] ? b[i] : a[i];
   }
-  return false;
+  return res;
 }
 #endif
 
-MAKE_VEC_AND_MARRAY_VERSIONS(acos)
-MAKE_VEC_AND_MARRAY_VERSIONS(acosh)
-MAKE_VEC_AND_MARRAY_VERSIONS(asin)
-MAKE_VEC_AND_MARRAY_VERSIONS(asinh)
-MAKE_VEC_AND_MARRAY_VERSIONS(atan)
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(atan2)
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(atan2pi)
-MAKE_VEC_AND_MARRAY_VERSIONS(atanh)
-MAKE_VEC_AND_MARRAY_VERSIONS(asinpi)
-MAKE_VEC_AND_MARRAY_VERSIONS(atanpi)
-MAKE_VEC_AND_MARRAY_VERSIONS_3ARGS(bitselect)
-MAKE_VEC_AND_MARRAY_VERSIONS(cbrt)
-MAKE_VEC_AND_MARRAY_VERSIONS(ceil)
+template <typename T, int N>
+sycl_cts::resultRef<sycl::vec<T, N>> abs(sycl::vec<T, N> a) {
+  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
+      [](T x) { return abs(x); }, a);
+}
+#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
+template <typename T, size_t N>
+sycl_cts::resultRef<sycl::marray<T, N>> abs(sycl::marray<T, N> a) {
+  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
+      [](T x) { return abs(x); }, a);
+}
+#endif
+
+template <typename T, int N>
+sycl_cts::resultRef<sycl::vec<T, N>> abs_diff(sycl::vec<T, N> a,
+                                              sycl::vec<T, N> b) {
+  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
+      [](T x, T y) { return abs_diff(x, y); }, a, b);
+}
+// FIXME: AdaptiveCpp does not support marray
+#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
+template <typename T, size_t N>
+sycl_cts::resultRef<sycl::marray<T, N>> abs_diff(sycl::marray<T, N> a,
+                                                 sycl::marray<T, N> b) {
+  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
+      [](T x, T y) { return abs_diff(x, y); }, a, b);
+}
+#endif
+
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(add_sat)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(hadd)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(rhadd)
 
 template <typename T, int N>
 sycl_cts::resultRef<sycl::vec<T, N>> clamp(sycl::vec<T, N> a, sycl::vec<T, N> b,
@@ -1070,67 +1037,316 @@ sycl_cts::resultRef<sycl::marray<T, N>> clamp(sycl::marray<T, N> a, T b, T c) {
 
 MAKE_VEC_AND_MARRAY_VERSIONS(clz)
 MAKE_VEC_AND_MARRAY_VERSIONS(ctz)
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(copysign)
-MAKE_VEC_AND_MARRAY_VERSIONS(cos)
-MAKE_VEC_AND_MARRAY_VERSIONS(cosh)
-MAKE_VEC_AND_MARRAY_VERSIONS(cospi)
-MAKE_VEC_AND_MARRAY_VERSIONS(degrees)
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(divide)
-MAKE_VEC_AND_MARRAY_VERSIONS(erf)
-MAKE_VEC_AND_MARRAY_VERSIONS(erfc)
-MAKE_VEC_AND_MARRAY_VERSIONS(exp)
-MAKE_VEC_AND_MARRAY_VERSIONS(exp2)
-MAKE_VEC_AND_MARRAY_VERSIONS(exp10)
-MAKE_VEC_AND_MARRAY_VERSIONS(expm1)
-MAKE_VEC_AND_MARRAY_VERSIONS(fabs)
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(fdim)
-MAKE_VEC_AND_MARRAY_VERSIONS(floor)
-MAKE_VEC_AND_MARRAY_VERSIONS_3ARGS(fma)
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(fmax)
-MAKE_VEC_AND_MARRAY_VERSIONS_WITH_SCALAR(fmax)
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(fmin)
-MAKE_VEC_AND_MARRAY_VERSIONS_WITH_SCALAR(fmin)
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(fmod)
+
+MAKE_VEC_AND_MARRAY_VERSIONS_3ARGS(mad_sat)
 
 template <typename T, int N>
-sycl::vec<T, N> fract(sycl::vec<T, N> a, sycl::vec<T, N> *b) {
-  sycl::vec<T, N> res;
-  sycl::vec<T, N> resPtr;
-  for (int i = 0; i < N; i++) {
-    T value;
-    res[i] = reference::fract(a[i], &value);
-    resPtr[i] = value;
-  }
-  *b = resPtr;
-  return res;
+sycl_cts::resultRef<sycl::vec<T, N>> max(sycl::vec<T, N> a, sycl::vec<T, N> b) {
+  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
+      [](T x, T y) { return max(x, y); }, a, b);
+}
+template <typename T, int N>
+sycl_cts::resultRef<sycl::vec<T, N>> max(sycl::vec<T, N> a, T b) {
+  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
+      [](T x, T y) { return max(x, y); }, a, b);
 }
 // FIXME: AdaptiveCpp does not support marray
 #ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
 template <typename T, size_t N>
-sycl::marray<T, N> fract(sycl::marray<T, N> a, sycl::marray<T, N> *b) {
-  sycl::marray<T, N> res;
-  sycl::marray<T, N> resPtr;
-  for (size_t i = 0; i < N; i++) {
-    T value;
-    res[i] = reference::fract(a[i], &value);
-    resPtr[i] = value;
-  }
-  *b = resPtr;
-  return res;
+sycl_cts::resultRef<sycl::marray<T, N>> max(sycl::marray<T, N> a,
+                                            sycl::marray<T, N> b) {
+  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
+      [](T x, T y) { return max(x, y); }, a, b);
+}
+template <typename T, size_t N>
+sycl_cts::resultRef<sycl::marray<T, N>> max(sycl::marray<T, N> a, T b) {
+  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
+      [](T x, T y) { return max(x, y); }, a, b);
 }
 #endif
 
 template <typename T, int N>
-sycl::vec<T, N> frexp(sycl::vec<T, N> a, sycl::vec<int, N> *b) {
-  sycl::vec<T, N> res;
-  sycl::vec<int, N> resPtr;
-  for (int i = 0; i < N; i++) {
-    int value;
-    res[i] = reference::frexp(a[i], &value);
-    resPtr[i] = value;
-  }
-  *b = resPtr;
-  return res;
+sycl_cts::resultRef<sycl::vec<T, N>> min(sycl::vec<T, N> a, sycl::vec<T, N> b) {
+  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
+      [](T x, T y) { return min(x, y); }, a, b);
+}
+template <typename T, int N>
+sycl_cts::resultRef<sycl::vec<T, N>> min(sycl::vec<T, N> a, T b) {
+  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
+      [](T x, T y) { return min(x, y); }, a, b);
+}
+// FIXME: AdaptiveCpp does not support marray
+#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
+template <typename T, size_t N>
+sycl_cts::resultRef<sycl::marray<T, N>> min(sycl::marray<T, N> a,
+                                            sycl::marray<T, N> b) {
+  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
+      [](T x, T y) { return min(x, y); }, a, b);
+}
+template <typename T, size_t N>
+sycl_cts::resultRef<sycl::marray<T, N>> min(sycl::marray<T, N> a, T b) {
+  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
+      [](T x, T y) { return min(x, y); }, a, b);
+}
+#endif
+
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(mul_hi)
+MAKE_VEC_AND_MARRAY_VERSIONS_3ARGS(mad_hi)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(rotate)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(sub_sat)
+
+template <typename T, int N>
+sycl::vec<typename upsample_t<T>::type, N> upsample(
+    sycl::vec<T, N> a, sycl::vec<typename std::make_unsigned<T>::type, N> b) {
+  return sycl_cts::math::run_func_on_vector<typename upsample_t<T>::type, T, N>(
+      [](T x, T y) { return upsample(x, y); }, a, b);
+}
+// FIXME: AdaptiveCpp does not support marray
+#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
+template <typename T, size_t N>
+sycl::marray<typename upsample_t<T>::type, N> upsample(
+    sycl::marray<T, N> a,
+    sycl::marray<typename std::make_unsigned<T>::type, N> b) {
+  return sycl_cts::math::run_func_on_marray<typename upsample_t<T>::type, T, N>(
+      [](T x, T y) { return upsample(x, y); }, a, b);
+}
+#endif
+
+MAKE_VEC_AND_MARRAY_VERSIONS(popcount)
+
+template <typename T, int N>
+sycl_cts::resultRef<sycl::vec<T, N>> mad24(sycl::vec<T, N> a, sycl::vec<T, N> b,
+                                           sycl::vec<T, N> c) {
+  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
+      [](T x, T y, T z) { return mad24(x, y, z); }, a, b, c);
+}
+// FIXME: AdaptiveCpp does not support marray
+#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
+template <typename T, size_t N>
+sycl_cts::resultRef<sycl::marray<T, N>> mad24(sycl::marray<T, N> a,
+                                              sycl::marray<T, N> b,
+                                              sycl::marray<T, N> c) {
+  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
+      [](T x, T y, T z) { return mad24(x, y, z); }, a, b, c);
+}
+#endif
+
+template <typename T, int N>
+sycl_cts::resultRef<sycl::vec<T, N>> mul24(sycl::vec<T, N> a,
+                                           sycl::vec<T, N> b) {
+  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
+      [](T x, T y) { return mul24(x, y); }, a, b);
+}
+// FIXME: AdaptiveCpp does not support marray
+#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
+template <typename T, size_t N>
+sycl_cts::resultRef<sycl::marray<T, N>> mul24(sycl::marray<T, N> a,
+                                              sycl::marray<T, N> b) {
+  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
+      [](T x, T y) { return mul24(x, y); }, a, b);
+}
+#endif
+
+MAKE_VEC_AND_MARRAY_VERSIONS(degrees)
+
+template <typename T, int N>
+sycl_cts::resultRef<sycl::vec<T, N>> mix(sycl::vec<T, N> a, sycl::vec<T, N> b,
+                                         sycl::vec<T, N> c) {
+  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
+      [](T x, T y, T z) { return mix(x, y, z); }, a, b, c);
+}
+template <typename T, int N>
+sycl_cts::resultRef<sycl::vec<T, N>> mix(sycl::vec<T, N> a, sycl::vec<T, N> b,
+                                         T c) {
+  sycl::vec<T, N> res;
+  std::map<int, bool> undefined;
+  for (int i = 0; i < N; i++) {
+    sycl_cts::resultRef<T> element = mix(a[i], b[i], c);
+    if (element.undefined.empty())
+      res[i] = element.res;
+    else
+      undefined[i] = true;
+  }
+  return sycl_cts::resultRef<sycl::vec<T, N>>(res, undefined);
+}
+// FIXME: AdaptiveCpp does not support marray
+#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
+template <typename T, size_t N>
+sycl_cts::resultRef<sycl::marray<T, N>> mix(sycl::marray<T, N> a,
+                                            sycl::marray<T, N> b,
+                                            sycl::marray<T, N> c) {
+  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
+      [](T x, T y, T z) { return mix(x, y, z); }, a, b, c);
+}
+template <typename T, size_t N>
+sycl_cts::resultRef<sycl::marray<T, N>> mix(sycl::marray<T, N> a,
+                                            sycl::marray<T, N> b, T c) {
+  sycl::marray<T, N> res;
+  std::map<int, bool> undefined;
+  for (size_t i = 0; i < N; i++) {
+    sycl_cts::resultRef<T> element = mix(a[i], b[i], c);
+    if (element.undefined.empty())
+      res[i] = element.res;
+    else
+      undefined[i] = true;
+  }
+  return sycl_cts::resultRef<sycl::marray<T, N>>(res, undefined);
+}
+#endif
+
+MAKE_VEC_AND_MARRAY_VERSIONS(radians)
+
+template <typename T, int N>
+sycl::vec<T, N> step(T a, sycl::vec<T, N> b) {
+  sycl::vec<T, N> res;
+  for (int i = 0; i < N; i++) {
+    res[i] = step(a, b[i]);
+  }
+  return res;
+}
+// FIXME: AdaptiveCpp does not support marray
+#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
+template <typename T, size_t N>
+sycl::marray<T, N> step(T a, sycl::marray<T, N> b) {
+  sycl::marray<T, N> res;
+  for (size_t i = 0; i < N; i++) {
+    res[i] = step(a, b[i]);
+  }
+  return res;
+}
+#endif
+
+template <typename T, int N>
+sycl_cts::resultRef<sycl::vec<T, N>> smoothstep(sycl::vec<T, N> a,
+                                                sycl::vec<T, N> b,
+                                                sycl::vec<T, N> c) {
+  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
+      [](T x, T y, T z) { return smoothstep(x, y, z); }, a, b, c);
+}
+template <typename T, int N>
+sycl_cts::resultRef<sycl::vec<T, N>> smoothstep(T a, T b, sycl::vec<T, N> c) {
+  sycl::vec<T, N> res;
+  std::map<int, bool> undefined;
+  for (int i = 0; i < N; i++) {
+    sycl_cts::resultRef<T> element = smoothstep(a, b, c[i]);
+    if (element.undefined.empty())
+      res[i] = element.res;
+    else
+      undefined[i] = true;
+  }
+  return sycl_cts::resultRef<sycl::vec<T, N>>(res, undefined);
+}
+// FIXME: AdaptiveCpp does not support marray
+#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
+template <typename T, size_t N>
+sycl_cts::resultRef<sycl::marray<T, N>> smoothstep(sycl::marray<T, N> a,
+                                                   sycl::marray<T, N> b,
+                                                   sycl::marray<T, N> c) {
+  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
+      [](T x, T y, T z) { return smoothstep(x, y, z); }, a, b, c);
+}
+template <typename T, size_t N>
+sycl_cts::resultRef<sycl::marray<T, N>> smoothstep(T a, T b,
+                                                   sycl::marray<T, N> c) {
+  sycl::marray<T, N> res;
+  std::map<int, bool> undefined;
+  for (size_t i = 0; i < N; i++) {
+    sycl_cts::resultRef<T> element = smoothstep(a, b, c[i]);
+    if (element.undefined.empty())
+      res[i] = element.res;
+    else
+      undefined[i] = true;
+  }
+  return sycl_cts::resultRef<sycl::marray<T, N>>(res, undefined);
+}
+#endif
+
+// Math functions
+
+template <typename T, int N>
+struct higher_accuracy<sycl::vec<T, N>> {
+  using type = sycl::vec<typename higher_accuracy<T>::type, N>;
+};
+// FIXME: AdaptiveCpp does not support marray
+#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
+template <typename T, size_t N>
+struct higher_accuracy<sycl::marray<T, N>> {
+  using type = sycl::marray<typename higher_accuracy<T>::type, N>;
+};
+#endif
+
+MAKE_VEC_AND_MARRAY_VERSIONS(acos)
+MAKE_VEC_AND_MARRAY_VERSIONS(acosh)
+MAKE_VEC_AND_MARRAY_VERSIONS(acospi)
+MAKE_VEC_AND_MARRAY_VERSIONS(asin)
+MAKE_VEC_AND_MARRAY_VERSIONS(asinh)
+MAKE_VEC_AND_MARRAY_VERSIONS(asinpi)
+MAKE_VEC_AND_MARRAY_VERSIONS(atan)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(atan2)
+MAKE_VEC_AND_MARRAY_VERSIONS(atanh)
+MAKE_VEC_AND_MARRAY_VERSIONS(atanpi)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(atan2pi)
+MAKE_VEC_AND_MARRAY_VERSIONS(cbrt)
+MAKE_VEC_AND_MARRAY_VERSIONS(ceil)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(copysign)
+MAKE_VEC_AND_MARRAY_VERSIONS(cos)
+MAKE_VEC_AND_MARRAY_VERSIONS(cosh)
+MAKE_VEC_AND_MARRAY_VERSIONS(cospi)
+MAKE_VEC_AND_MARRAY_VERSIONS(erfc)
+MAKE_VEC_AND_MARRAY_VERSIONS(erf)
+MAKE_VEC_AND_MARRAY_VERSIONS(exp)
+MAKE_VEC_AND_MARRAY_VERSIONS(exp2)
+MAKE_VEC_AND_MARRAY_VERSIONS(exp10)
+MAKE_VEC_AND_MARRAY_VERSIONS(expm1)
+MAKE_VEC_AND_MARRAY_VERSIONS(fabs)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(fdim)
+MAKE_VEC_AND_MARRAY_VERSIONS(floor)
+MAKE_VEC_AND_MARRAY_VERSIONS_3ARGS(fma)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(fmax)
+MAKE_VEC_AND_MARRAY_VERSIONS_WITH_SCALAR(fmax)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(fmin)
+MAKE_VEC_AND_MARRAY_VERSIONS_WITH_SCALAR(fmin)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(fmod)
+
+template <typename T, int N>
+sycl::vec<T, N> fract(sycl::vec<T, N> a, sycl::vec<T, N> *b) {
+  sycl::vec<T, N> res;
+  sycl::vec<T, N> resPtr;
+  for (int i = 0; i < N; i++) {
+    T value;
+    res[i] = reference::fract(a[i], &value);
+    resPtr[i] = value;
+  }
+  *b = resPtr;
+  return res;
+}
+// FIXME: AdaptiveCpp does not support marray
+#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
+template <typename T, size_t N>
+sycl::marray<T, N> fract(sycl::marray<T, N> a, sycl::marray<T, N> *b) {
+  sycl::marray<T, N> res;
+  sycl::marray<T, N> resPtr;
+  for (size_t i = 0; i < N; i++) {
+    T value;
+    res[i] = reference::fract(a[i], &value);
+    resPtr[i] = value;
+  }
+  *b = resPtr;
+  return res;
+}
+#endif
+
+template <typename T, int N>
+sycl::vec<T, N> frexp(sycl::vec<T, N> a, sycl::vec<int, N> *b) {
+  sycl::vec<T, N> res;
+  sycl::vec<int, N> resPtr;
+  for (int i = 0; i < N; i++) {
+    int value;
+    res[i] = reference::frexp(a[i], &value);
+    resPtr[i] = value;
+  }
+  *b = resPtr;
+  return res;
 }
 // FIXME: AdaptiveCpp does not support marray
 #ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
@@ -1148,7 +1364,6 @@ sycl::marray<T, N> frexp(sycl::marray<T, N> a, sycl::marray<int, N> *b) {
 }
 #endif
 
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(hadd)
 MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(hypot)
 
 template <typename T, int N>
@@ -1212,11 +1427,6 @@ sycl::marray<T, N> ldexp(sycl::marray<T, N> a, int b) {
 
 MAKE_VEC_AND_MARRAY_VERSIONS(lgamma)
 
-template <typename T>
-T lgamma_r(T a, int *b) {
-  *b = (std::tgamma(a) > 0) ? 1 : -1;
-  return std::lgamma(a);
-}
 template <typename T, int N>
 sycl::vec<T, N> lgamma_r(sycl::vec<T, N> a, sycl::vec<int, N> *b) {
   sycl::vec<T, N> res;
@@ -1250,107 +1460,9 @@ MAKE_VEC_AND_MARRAY_VERSIONS(log2)
 MAKE_VEC_AND_MARRAY_VERSIONS(log10)
 MAKE_VEC_AND_MARRAY_VERSIONS(log1p)
 MAKE_VEC_AND_MARRAY_VERSIONS(logb)
-MAKE_VEC_AND_MARRAY_VERSIONS_3ARGS(mad)
-MAKE_VEC_AND_MARRAY_VERSIONS_3ARGS(mad_hi)
-MAKE_VEC_AND_MARRAY_VERSIONS_3ARGS(mad_sat)
-
-template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> max(sycl::vec<T, N> a, sycl::vec<T, N> b) {
-  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
-      [](T x, T y) { return max(x, y); }, a, b);
-}
-template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> max(sycl::vec<T, N> a, T b) {
-  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
-      [](T x, T y) { return max(x, y); }, a, b);
-}
-// FIXME: AdaptiveCpp does not support marray
-#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
-template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> max(sycl::marray<T, N> a,
-                                            sycl::marray<T, N> b) {
-  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
-      [](T x, T y) { return max(x, y); }, a, b);
-}
-template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> max(sycl::marray<T, N> a, T b) {
-  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
-      [](T x, T y) { return max(x, y); }, a, b);
-}
-#endif
 
+MAKE_VEC_AND_MARRAY_VERSIONS_3ARGS(mad)
 MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(maxmag)
-
-template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> min(sycl::vec<T, N> a, sycl::vec<T, N> b) {
-  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
-      [](T x, T y) { return min(x, y); }, a, b);
-}
-template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> min(sycl::vec<T, N> a, T b) {
-  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
-      [](T x, T y) { return min(x, y); }, a, b);
-}
-// FIXME: AdaptiveCpp does not support marray
-#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
-template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> min(sycl::marray<T, N> a,
-                                            sycl::marray<T, N> b) {
-  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
-      [](T x, T y) { return min(x, y); }, a, b);
-}
-template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> min(sycl::marray<T, N> a, T b) {
-  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
-      [](T x, T y) { return min(x, y); }, a, b);
-}
-#endif
-
-template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> mix(sycl::vec<T, N> a, sycl::vec<T, N> b,
-                                         sycl::vec<T, N> c) {
-  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
-      [](T x, T y, T z) { return mix(x, y, z); }, a, b, c);
-}
-template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> mix(sycl::vec<T, N> a, sycl::vec<T, N> b,
-                                         T c) {
-  sycl::vec<T, N> res;
-  std::map<int, bool> undefined;
-  for (int i = 0; i < N; i++) {
-    sycl_cts::resultRef<T> element = mix(a[i], b[i], c);
-    if (element.undefined.empty())
-      res[i] = element.res;
-    else
-      undefined[i] = true;
-  }
-  return sycl_cts::resultRef<sycl::vec<T, N>>(res, undefined);
-}
-// FIXME: AdaptiveCpp does not support marray
-#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
-template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> mix(sycl::marray<T, N> a,
-                                            sycl::marray<T, N> b,
-                                            sycl::marray<T, N> c) {
-  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
-      [](T x, T y, T z) { return mix(x, y, z); }, a, b, c);
-}
-template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> mix(sycl::marray<T, N> a,
-                                            sycl::marray<T, N> b, T c) {
-  sycl::marray<T, N> res;
-  std::map<int, bool> undefined;
-  for (size_t i = 0; i < N; i++) {
-    sycl_cts::resultRef<T> element = mix(a[i], b[i], c);
-    if (element.undefined.empty())
-      res[i] = element.res;
-    else
-      undefined[i] = true;
-  }
-  return sycl_cts::resultRef<sycl::marray<T, N>>(res, undefined);
-}
-#endif
-
 MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(minmag)
 
 template <typename T, int N>
@@ -1381,8 +1493,6 @@ sycl::marray<T, N> modf(sycl::marray<T, N> a, sycl::marray<T, N> *b) {
 }
 #endif
 
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(mul_hi)
-
 template <int N>
 sycl::vec<float, N> nan(sycl::vec<unsigned int, N> a) {
   return sycl_cts::math::run_func_on_vector<float, unsigned int, N>(
@@ -1398,7 +1508,6 @@ sycl::marray<float, N> nan(sycl::marray<unsigned int, N> a) {
 #endif
 
 MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(nextafter)
-
 MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(pow)
 
 template <typename T, int N>
@@ -1437,8 +1546,6 @@ sycl_cts::resultRef<sycl::marray<T, N>> powr(sycl::marray<T, N> a,
 }
 #endif
 
-MAKE_VEC_AND_MARRAY_VERSIONS(radians)
-MAKE_VEC_AND_MARRAY_VERSIONS(recip)
 MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(remainder)
 
 template <typename T, int N>
@@ -1471,7 +1578,7 @@ sycl::marray<T, N> remquo(sycl::marray<T, N> a, sycl::marray<T, N> b,
 }
 #endif
 
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(rhadd)
+MAKE_VEC_AND_MARRAY_VERSIONS(rint)
 
 template <typename T, int N>
 sycl::vec<T, N> rootn(sycl::vec<T, N> a, sycl::vec<int, N> b) {
@@ -1493,36 +1600,8 @@ sycl::marray<T, N> rootn(sycl::marray<T, N> a, sycl::marray<int, N> b) {
 }
 #endif
 
-
-MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(rotate)
 MAKE_VEC_AND_MARRAY_VERSIONS(round)
 MAKE_VEC_AND_MARRAY_VERSIONS(rsqrt)
-
-template <typename T, typename K, int N>
-sycl::vec<T, N> select(sycl::vec<T, N> a, sycl::vec<T, N> b,
-                       sycl::vec<K, N> c) {
-  sycl::vec<T, N> res;
-  for (int i = 0; i < N; i++) {
-    if (any(c[i]) == 1)
-      res[i] = b[i];
-    else
-      res[i] = a[i];
-  }
-  return res;
-}
-// FIXME: AdaptiveCpp does not support marray
-#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
-template <typename T, size_t N>
-sycl::marray<T, N> select(sycl::marray<T, N> a, sycl::marray<T, N> b,
-                          sycl::marray<bool, N> c) {
-  sycl::marray<T, N> res;
-  for (size_t i = 0; i < N; i++) {
-    res[i] = c[i] ? b[i] : a[i];
-  }
-  return res;
-}
-#endif
-
 MAKE_VEC_AND_MARRAY_VERSIONS(sign)
 
 template <typename T, int N>
@@ -1556,89 +1635,18 @@ sycl::marray<T, N> sincos(sycl::marray<T, N> a, sycl::marray<T, N> *b) {
 MAKE_VEC_AND_MARRAY_VERSIONS(sin)
 MAKE_VEC_AND_MARRAY_VERSIONS(sinh)
 MAKE_VEC_AND_MARRAY_VERSIONS(sinpi)
-
-template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> smoothstep(sycl::vec<T, N> a,
-                                                sycl::vec<T, N> b,
-                                                sycl::vec<T, N> c) {
-  return sycl_cts::math::run_func_on_vector_result_ref<T, N>(
-      [](T x, T y, T z) { return smoothstep(x, y, z); }, a, b, c);
-}
-template <typename T, int N>
-sycl_cts::resultRef<sycl::vec<T, N>> smoothstep(T a, T b, sycl::vec<T, N> c) {
-  sycl::vec<T, N> res;
-  std::map<int, bool> undefined;
-  for (int i = 0; i < N; i++) {
-    sycl_cts::resultRef<T> element = smoothstep(a, b, c[i]);
-    if (element.undefined.empty())
-      res[i] = element.res;
-    else
-      undefined[i] = true;
-  }
-  return sycl_cts::resultRef<sycl::vec<T, N>>(res, undefined);
-}
-// FIXME: AdaptiveCpp does not support marray
-#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
-template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> smoothstep(sycl::marray<T, N> a,
-                                                   sycl::marray<T, N> b,
-                                                   sycl::marray<T, N> c) {
-  return sycl_cts::math::run_func_on_marray_result_ref<T, N>(
-      [](T x, T y, T z) { return smoothstep(x, y, z); }, a, b, c);
-}
-template <typename T, size_t N>
-sycl_cts::resultRef<sycl::marray<T, N>> smoothstep(T a, T b,
-                                                   sycl::marray<T, N> c) {
-  sycl::marray<T, N> res;
-  std::map<int, bool> undefined;
-  for (size_t i = 0; i < N; i++) {
-    sycl_cts::resultRef<T> element = smoothstep(a, b, c[i]);
-    if (element.undefined.empty())
-      res[i] = element.res;
-    else
-      undefined[i] = true;
-  }
-  return sycl_cts::resultRef<sycl::marray<T, N>>(res, undefined);
-}
-#endif
-
 MAKE_VEC_AND_MARRAY_VERSIONS(sqrt)
-
-template <typename T, int N>
-sycl::vec<T, N> step(T a, sycl::vec<T, N> b) {
-  sycl::vec<T, N> res;
-  for (int i = 0; i < N; i++) {
-    res[i] = step(a, b[i]);
-  }
-  return res;
-}
-// FIXME: AdaptiveCpp does not support marray
-#ifndef SYCL_CTS_COMPILING_WITH_ADAPTIVECPP
-template <typename T, size_t N>
-sycl::marray<T, N> step(T a, sycl::marray<T, N> b) {
-  sycl::marray<T, N> res;
-  for (size_t i = 0; i < N; i++) {
-    res[i] = step(a, b[i]);
-  }
-  return res;
-}
-#endif
-
 MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(step)
 MAKE_VEC_AND_MARRAY_VERSIONS(tan)
 MAKE_VEC_AND_MARRAY_VERSIONS(tanh)
 MAKE_VEC_AND_MARRAY_VERSIONS(tanpi)
 MAKE_VEC_AND_MARRAY_VERSIONS(tgamma)
 MAKE_VEC_AND_MARRAY_VERSIONS(trunc)
-
+MAKE_VEC_AND_MARRAY_VERSIONS(recip)
+MAKE_VEC_AND_MARRAY_VERSIONS_2ARGS(divide)
 
 // Geometric functions
 
-template <typename T>
-T dot(T p0, T p1) {
-  return p0 * p1;
-}
-
 template <typename T, int N>
 T dot(sycl::vec<T, N> a, sycl::vec<T, N> b) {
   T res = 0;
@@ -1655,6 +1663,8 @@ T dot(sycl::marray<T, N> a, sycl::marray<T, N> b) {
 }
 #endif
 
+// Generic functions over both scalars and vec / marray types - these need to be defined last
+
 template <typename T>
 auto length(T p) {
   return reference::sqrt(reference::dot(p, p));
@@ -1665,27 +1675,6 @@ auto distance(T p0, T p1) {
   return reference::length(p0 - p1);
 }
 
-template <typename T>
-T normalize(T p) {
-  if (p < 0) return -1;
-  return 1;
-}
-
-template <typename T>
-float fast_length(T p0) {
-  return reference::sqrt(fast_dot(p0));
-}
-
-template <typename T>
-float fast_distance(T p0, T p1) {
-  return reference::fast_length(p0 - p1);
-}
-
-template <typename T>
-T fast_normalize(T p0) {
-  return p0 * reference::rsqrt(fast_dot(p0));
-}
-
 template <typename T, int N>
 sycl::vec<T, N> normalize(sycl::vec<T, N> a) {
   sycl::vec<T, N> res;
@@ -1706,6 +1695,21 @@ sycl::marray<T, N> normalize(sycl::marray<T, N> a) {
 }
 #endif
 
+template <typename T>
+float fast_length(T p0) {
+  return reference::sqrt(fast_dot(p0));
+}
+
+template <typename T>
+float fast_distance(T p0, T p1) {
+  return reference::fast_length(p0 - p1);
+}
+
+template <typename T>
+T fast_normalize(T p0) {
+  return p0 * reference::rsqrt(fast_dot(p0));
+}
+
 }  // namespace reference
 
 #endif  // __SYCLCTS_UTIL_MATH_REFERENCE_H