Update Eigen to commit:489dbbc6516f6729fa8e5c8bddc9ea446d5ef032 CHANGELOG ========= 489dbbc65 - make fixed_size matrices conform to std::is_standard_layout 283d871a3 - Add missing EIGEN_DEVICE_FUNCTION decorations. 0d366f653 - Vectorize erfc(x) for double and improve erfc(x) for float. 8adf43640 - more avx predux_any bc424f617 - add missing avx predux_any functions e52ac76ca - use EIGEN_CPLUSPLUS instead of checking cpp version 122be167c - Revert "make fixed-size objects trivially move assignable" d49021212 - Tensor Roll / Circular Shift / Rotate bb73be8a2 - make fixed-size objects trivially move assignable PiperOrigin-RevId: 696264736 Change-Id: I84aad26ecbe5fd7afee89965e0d6c6f05aca48e6

diff --git a/Eigen/src/Core/DenseStorage.h b/Eigen/src/Core/DenseStorage.h
index 90bc94e..d62586c 100644
--- a/Eigen/src/Core/DenseStorage.h
+++ b/Eigen/src/Core/DenseStorage.h

@@ -30,10 +30,10 @@
 #if defined(EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT)
 #define EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(Alignment)
 #else
-#define EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(Alignment)                                              \
-  eigen_assert((internal::is_constant_evaluated() || (std::uintptr_t(array) % Alignment == 0)) && \
-               "this assertion is explained here: "                                               \
-               "http://eigen.tuxfamily.org/dox-devel/group__TopicUnalignedArrayAssert.html"       \
+#define EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(Alignment)                                        \
+  eigen_assert((is_constant_evaluated() || (std::uintptr_t(array) % Alignment == 0)) &&     \
+               "this assertion is explained here: "                                         \
+               "http://eigen.tuxfamily.org/dox-devel/group__TopicUnalignedArrayAssert.html" \
                " **** READ THIS WEB PAGE !!! ****");
 #endif
 
@@ -48,6 +48,7 @@
  * Static array. If the MatrixOrArrayOptions require auto-alignment, the array will be automatically aligned:
  * to 16 bytes boundary if the total size is a multiple of 16 bytes.
  */
+
 template <typename T, int Size, int MatrixOrArrayOptions,
           int Alignment = (MatrixOrArrayOptions & DontAlign) ? 0 : compute_default_alignment<T, Size>::value>
 struct plain_array {
@@ -78,92 +79,454 @@
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr plain_array() = default;
 };
 
-// this class is intended to be inherited by DenseStorage to take advantage of empty base optimization
-template <int Rows, int Cols>
-struct DenseStorageIndices {
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices() = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(const DenseStorageIndices&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(DenseStorageIndices&&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices& operator=(const DenseStorageIndices&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices& operator=(DenseStorageIndices&&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(Index /*rows*/, Index /*cols*/) {}
+template <typename T, int Size, int Options, int Alignment>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap_plain_array(plain_array<T, Size, Options, Alignment>& a,
+                                                                      plain_array<T, Size, Options, Alignment>& b,
+                                                                      Index a_size, Index b_size) {
+  Index common_size = numext::mini(a_size, b_size);
+  std::swap_ranges(a.array, a.array + common_size, b.array);
+  if (a_size > b_size)
+    smart_copy(a.array + common_size, a.array + a_size, b.array + common_size);
+  else if (b_size > a_size)
+    smart_copy(b.array + common_size, b.array + b_size, a.array + common_size);
+}
+
+template <typename T, int Size, int Rows, int Cols, int Options>
+class DenseStorage_impl {
+  plain_array<T, Size, Options> m_data;
+
+ public:
+#ifndef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl&) = default;
+#else
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN(Index size = Size)
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl& other) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN(Index size = Size)
+    smart_copy(other.m_data.array, other.m_data.array + Size, m_data.array);
+  }
+#endif
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(Index /*size*/, Index /*rows*/, Index /*cols*/) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(const DenseStorage_impl&) = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage_impl& other) {
+    numext::swap(m_data, other.m_data);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index /*size*/, Index /*rows*/,
+                                                                          Index /*cols*/) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index /*size*/, Index /*rows*/, Index /*cols*/) {}
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return Rows; }
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return Cols; }
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return Rows * Cols; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void set(Index /*rows*/, Index /*cols*/) {}
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void swap(DenseStorageIndices& /*other*/) noexcept {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return m_data.array; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return m_data.array; }
 };
-template <int Rows>
-struct DenseStorageIndices<Rows, Dynamic> {
-  Index m_cols;
+template <typename T, int Size, int Cols, int Options>
+class DenseStorage_impl<T, Size, Dynamic, Cols, Options> {
+  plain_array<T, Size, Options> m_data;
+  Index m_rows = 0;
 
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices() : m_cols(0) {}
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(const DenseStorageIndices&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(DenseStorageIndices&&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices& operator=(const DenseStorageIndices&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices& operator=(DenseStorageIndices&&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(Index /*rows*/, Index cols) : m_cols(cols) {}
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return Rows; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return m_cols; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return Rows * m_cols; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void set(Index /*rows*/, Index cols) { m_cols = cols; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void swap(DenseStorageIndices& other) noexcept {
-    numext::swap(m_cols, other.m_cols);
+ public:
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl& other)
+      : m_rows(other.m_rows) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN(Index size = other.size())
+    smart_copy(other.m_data.array, other.m_data.array + other.size(), m_data.array);
   }
-};
-template <int Cols>
-struct DenseStorageIndices<Dynamic, Cols> {
-  Index m_rows;
-
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices() : m_rows(0) {}
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(const DenseStorageIndices&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(DenseStorageIndices&&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices& operator=(const DenseStorageIndices&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices& operator=(DenseStorageIndices&&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(Index rows, Index /*cols*/) : m_rows(rows) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(Index size, Index rows, Index /*cols*/)
+      : m_rows(rows) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
+    EIGEN_UNUSED_VARIABLE(size)
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(const DenseStorage_impl& other) {
+    smart_copy(other.m_data.array, other.m_data.array + other.size(), m_data.array);
+    m_rows = other.m_rows;
+    return *this;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage_impl& other) {
+    swap_plain_array(m_data, other.m_data, size(), other.size());
+    numext::swap(m_rows, other.m_rows);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index /*size*/, Index rows, Index /*cols*/) {
+    m_rows = rows;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index /*size*/, Index rows, Index /*cols*/) {
+    m_rows = rows;
+  }
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return m_rows; }
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return Cols; }
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return m_rows * Cols; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void set(Index rows, Index /*cols*/) { m_rows = rows; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void swap(DenseStorageIndices& other) noexcept {
-    numext::swap(m_rows, other.m_rows);
-  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return m_data.array; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return m_data.array; }
 };
-template <>
-struct DenseStorageIndices<Dynamic, Dynamic> {
-  Index m_rows;
-  Index m_cols;
+template <typename T, int Size, int Rows, int Options>
+class DenseStorage_impl<T, Size, Rows, Dynamic, Options> {
+  plain_array<T, Size, Options> m_data;
+  Index m_cols = 0;
 
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices() : m_rows(0), m_cols(0) {}
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(const DenseStorageIndices&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(DenseStorageIndices&&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices& operator=(const DenseStorageIndices&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices& operator=(DenseStorageIndices&&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorageIndices(Index rows, Index cols)
-      : m_rows(rows), m_cols(cols) {}
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return m_rows; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return m_cols; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return m_rows * m_cols; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void set(Index rows, Index cols) {
-    m_rows = rows;
+ public:
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl& other)
+      : m_cols(other.m_cols) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN(Index size = other.size())
+    smart_copy(other.m_data.array, other.m_data.array + other.size(), m_data.array);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(Index size, Index /*rows*/, Index cols)
+      : m_cols(cols) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
+    EIGEN_UNUSED_VARIABLE(size)
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(const DenseStorage_impl& other) {
+    smart_copy(other.m_data.array, other.m_data.array + other.size(), m_data.array);
+    m_cols = other.m_cols;
+    return *this;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage_impl& other) {
+    swap_plain_array(m_data, other.m_data, size(), other.size());
+    numext::swap(m_cols, other.m_cols);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index /*size*/, Index /*rows*/, Index cols) {
     m_cols = cols;
   }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void swap(DenseStorageIndices& other) noexcept {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index /*size*/, Index /*rows*/, Index cols) {
+    m_cols = cols;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return Rows; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return Rows * m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return m_data.array; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return m_data.array; }
+};
+template <typename T, int Size, int Options>
+class DenseStorage_impl<T, Size, Dynamic, Dynamic, Options> {
+  plain_array<T, Size, Options> m_data;
+  Index m_rows = 0;
+  Index m_cols = 0;
+
+ public:
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl& other)
+      : m_rows(other.m_rows), m_cols(other.m_cols) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN(Index size = other.size())
+    smart_copy(other.m_data.array, other.m_data.array + other.size(), m_data.array);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(Index size, Index rows, Index cols)
+      : m_rows(rows), m_cols(cols) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
+    EIGEN_UNUSED_VARIABLE(size)
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(const DenseStorage_impl& other) {
+    smart_copy(other.m_data.array, other.m_data.array + other.size(), m_data.array);
+    m_rows = other.m_rows;
+    m_cols = other.m_cols;
+    return *this;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage_impl& other) {
+    swap_plain_array(m_data, other.m_data, size(), other.size());
     numext::swap(m_rows, other.m_rows);
     numext::swap(m_cols, other.m_cols);
   }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index /*size*/, Index rows, Index cols) {
+    m_rows = rows;
+    m_cols = cols;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index /*size*/, Index rows, Index cols) {
+    m_rows = rows;
+    m_cols = cols;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return m_rows; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return m_rows * m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return m_data.array; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return m_data.array; }
 };
-
-template <int Size, int Rows, int Cols>
-struct use_trivial_ctors {
-  static constexpr bool value = (Size >= 0) && (Rows >= 0) && (Cols >= 0) && (Size == Rows * Cols);
+// null matrix variants
+template <typename T, int Rows, int Cols, int Options>
+class DenseStorage_impl<T, 0, Rows, Cols, Options> {
+ public:
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl&) = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(Index /*size*/, Index /*rows*/, Index /*cols*/) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(const DenseStorage_impl&) = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage_impl&) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index /*size*/, Index /*rows*/,
+                                                                          Index /*cols*/) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index /*size*/, Index /*rows*/, Index /*cols*/) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return Rows; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return Cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return Rows * Cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return nullptr; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return nullptr; }
 };
+template <typename T, int Cols, int Options>
+class DenseStorage_impl<T, 0, Dynamic, Cols, Options> {
+  Index m_rows = 0;
 
+ public:
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl&) = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(Index /*size*/, Index rows, Index /*cols*/)
+      : m_rows(rows) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(const DenseStorage_impl&) = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage_impl& other) noexcept {
+    numext::swap(m_rows, other.m_rows);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index /*size*/, Index rows, Index /*cols*/) {
+    m_rows = rows;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index /*size*/, Index rows, Index /*cols*/) {
+    m_rows = rows;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return m_rows; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return Cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return m_rows * Cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return nullptr; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return nullptr; }
+};
+template <typename T, int Rows, int Options>
+class DenseStorage_impl<T, 0, Rows, Dynamic, Options> {
+  Index m_cols = 0;
+
+ public:
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl&) = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(Index /*size*/, Index /*rows*/, Index cols)
+      : m_cols(cols) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(const DenseStorage_impl&) = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage_impl& other) noexcept {
+    numext::swap(m_cols, other.m_cols);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index /*size*/, Index /*rows*/, Index cols) {
+    m_cols = cols;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index /*size*/, Index /*rows*/, Index cols) {
+    m_cols = cols;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return Rows; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return Rows * m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return nullptr; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return nullptr; }
+};
+template <typename T, int Options>
+class DenseStorage_impl<T, 0, Dynamic, Dynamic, Options> {
+  Index m_rows = 0;
+  Index m_cols = 0;
+
+ public:
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl&) = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(Index /*size*/, Index rows, Index cols)
+      : m_rows(rows), m_cols(cols) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(const DenseStorage_impl&) = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage_impl& other) noexcept {
+    numext::swap(m_rows, other.m_rows);
+    numext::swap(m_cols, other.m_cols);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index /*size*/, Index rows, Index cols) {
+    m_rows = rows;
+    m_cols = cols;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index /*size*/, Index rows, Index cols) {
+    m_rows = rows;
+    m_cols = cols;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return m_rows; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return m_rows * m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return nullptr; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return nullptr; }
+};
+// fixed-size matrix with dynamic memory allocation not currently supported
+template <typename T, int Rows, int Cols, int Options>
+class DenseStorage_impl<T, Dynamic, Rows, Cols, Options> {};
+// dynamic-sized variants
+template <typename T, int Cols, int Options>
+class DenseStorage_impl<T, Dynamic, Dynamic, Cols, Options> {
+  static constexpr bool Align = (Options & DontAlign) == 0;
+  T* m_data = nullptr;
+  Index m_rows = 0;
+
+ public:
+  static constexpr int Size = Dynamic;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl& other)
+      : m_data(conditional_aligned_new_auto<T, Align>(other.size())), m_rows(other.m_rows) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN(Index size = other.size())
+    smart_copy(other.m_data, other.m_data + other.size(), m_data);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(Index size, Index rows, Index /*cols*/)
+      : m_data(conditional_aligned_new_auto<T, Align>(size)), m_rows(rows) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(DenseStorage_impl&& other) noexcept
+      : m_data(other.m_data), m_rows(other.m_rows) {
+    other.m_data = nullptr;
+    other.m_rows = 0;
+  }
+  EIGEN_DEVICE_FUNC ~DenseStorage_impl() { conditional_aligned_delete_auto<T, Align>(m_data, size()); }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(const DenseStorage_impl& other) {
+    resize(other.size(), other.rows(), other.cols());
+    smart_copy(other.m_data, other.m_data + other.size(), m_data);
+    return *this;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(DenseStorage_impl&& other) noexcept {
+    this->swap(other);
+    return *this;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage_impl& other) noexcept {
+    numext::swap(m_data, other.m_data);
+    numext::swap(m_rows, other.m_rows);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index size, Index rows, Index /*cols*/) {
+    m_data = conditional_aligned_realloc_new_auto<T, Align>(m_data, size, this->size());
+    m_rows = rows;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index size, Index rows, Index /*cols*/) {
+    Index oldSize = this->size();
+    if (oldSize != size) {
+      conditional_aligned_delete_auto<T, Align>(m_data, oldSize);
+      m_data = conditional_aligned_new_auto<T, Align>(size);
+      EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
+    }
+    m_rows = rows;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return m_rows; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return Cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return m_rows * Cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return m_data; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return m_data; }
+};
+template <typename T, int Rows, int Options>
+class DenseStorage_impl<T, Dynamic, Rows, Dynamic, Options> {
+  static constexpr bool Align = (Options & DontAlign) == 0;
+  T* m_data = nullptr;
+  Index m_cols = 0;
+
+ public:
+  static constexpr int Size = Dynamic;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl& other)
+      : m_data(conditional_aligned_new_auto<T, Align>(other.size())), m_cols(other.m_cols) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN(Index size = other.size())
+    smart_copy(other.m_data, other.m_data + other.size(), m_data);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(Index size, Index /*rows*/, Index cols)
+      : m_data(conditional_aligned_new_auto<T, Align>(size)), m_cols(cols) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(DenseStorage_impl&& other) noexcept
+      : m_data(other.m_data), m_cols(other.m_cols) {
+    other.m_data = nullptr;
+    other.m_cols = 0;
+  }
+  EIGEN_DEVICE_FUNC ~DenseStorage_impl() { conditional_aligned_delete_auto<T, Align>(m_data, size()); }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(const DenseStorage_impl& other) {
+    resize(other.size(), other.rows(), other.cols());
+    smart_copy(other.m_data, other.m_data + other.size(), m_data);
+    return *this;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(DenseStorage_impl&& other) noexcept {
+    this->swap(other);
+    return *this;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage_impl& other) noexcept {
+    numext::swap(m_data, other.m_data);
+    numext::swap(m_cols, other.m_cols);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index size, Index /*rows*/, Index cols) {
+    m_data = conditional_aligned_realloc_new_auto<T, Align>(m_data, size, this->size());
+    m_cols = cols;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index size, Index /*rows*/, Index cols) {
+    Index oldSize = this->size();
+    if (oldSize != size) {
+      conditional_aligned_delete_auto<T, Align>(m_data, oldSize);
+      m_data = conditional_aligned_new_auto<T, Align>(size);
+      EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
+    }
+    m_cols = cols;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return Rows; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return Rows * m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return m_data; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return m_data; }
+};
+template <typename T, int Options>
+class DenseStorage_impl<T, Dynamic, Dynamic, Dynamic, Options> {
+  static constexpr bool Align = (Options & DontAlign) == 0;
+  T* m_data = nullptr;
+  Index m_rows = 0;
+  Index m_cols = 0;
+
+ public:
+  static constexpr int Size = Dynamic;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(const DenseStorage_impl& other)
+      : m_data(conditional_aligned_new_auto<T, Align>(other.size())), m_rows(other.m_rows), m_cols(other.m_cols) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN(Index size = other.size())
+    smart_copy(other.m_data, other.m_data + other.size(), m_data);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(Index size, Index rows, Index cols)
+      : m_data(conditional_aligned_new_auto<T, Align>(size)), m_rows(rows), m_cols(cols) {
+    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl(DenseStorage_impl&& other) noexcept
+      : m_data(other.m_data), m_rows(other.m_rows), m_cols(other.m_cols) {
+    other.m_data = nullptr;
+    other.m_rows = 0;
+    other.m_cols = 0;
+  }
+  EIGEN_DEVICE_FUNC ~DenseStorage_impl() { conditional_aligned_delete_auto<T, Align>(m_data, size()); }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(const DenseStorage_impl& other) {
+    resize(other.size(), other.rows(), other.cols());
+    smart_copy(other.m_data, other.m_data + other.size(), m_data);
+    return *this;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage_impl& operator=(DenseStorage_impl&& other) noexcept {
+    this->swap(other);
+    return *this;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage_impl& other) noexcept {
+    numext::swap(m_data, other.m_data);
+    numext::swap(m_rows, other.m_rows);
+    numext::swap(m_cols, other.m_cols);
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index size, Index rows, Index cols) {
+    m_data = conditional_aligned_realloc_new_auto<T, Align>(m_data, size, this->size());
+    m_rows = rows;
+    m_cols = cols;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index size, Index rows, Index cols) {
+    Index oldSize = this->size();
+    if (oldSize != size) {
+      conditional_aligned_delete_auto<T, Align>(m_data, oldSize);
+      m_data = conditional_aligned_new_auto<T, Align>(size);
+      EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
+    }
+    m_rows = rows;
+    m_cols = cols;
+  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index rows() const { return m_rows; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index cols() const { return m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Index size() const { return m_rows * m_cols; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return m_data; }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return m_data; }
+};
+template <typename T, int Size, int Rows, int Cols>
+struct use_default_move {
+  static constexpr bool DynamicObject = Size == Dynamic;
+  static constexpr bool TrivialObject =
+      (!NumTraits<T>::RequireInitialization) && (Rows >= 0) && (Cols >= 0) && (Size == Rows * Cols);
+  static constexpr bool value = DynamicObject || TrivialObject;
+};
 }  // end namespace internal
 
 /** \internal
  *
- * \class DenseStorage
+ * \class DenseStorage_impl
  * \ingroup Core_Module
  *
  * \brief Stores the data of a matrix
@@ -174,201 +537,42 @@
  * \sa Matrix
  */
 template <typename T, int Size, int Rows, int Cols, int Options,
-          bool Trivial = internal::use_trivial_ctors<Size, Rows, Cols>::value>
-class DenseStorage;
-
-// fixed-size storage with fixed dimensions
-template <typename T, int Size, int Rows, int Cols, int Options>
-class DenseStorage<T, Size, Rows, Cols, Options, true> : internal::DenseStorageIndices<Rows, Cols> {
-  using Base = internal::DenseStorageIndices<Rows, Cols>;
-
-  internal::plain_array<T, Size, Options> m_data;
+          bool Trivial = internal::use_default_move<T, Size, Rows, Cols>::value>
+class DenseStorage : public internal::DenseStorage_impl<T, Size, Rows, Cols, Options> {
+  using Base = internal::DenseStorage_impl<T, Size, Rows, Cols, Options>;
 
  public:
-  using Base::cols;
-  using Base::rows;
-#ifndef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE DenseStorage() = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(const DenseStorage&) = default;
-#else
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE DenseStorage() { EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN(Index size = Size) }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(const DenseStorage& other)
-      : Base(other), m_data(other.m_data) {
-    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN(Index size = Size)
-  }
-#endif
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(DenseStorage&&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage& operator=(const DenseStorage&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage& operator=(DenseStorage&&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(Index size, Index rows, Index cols) : Base(rows, cols) {
-    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
-    EIGEN_UNUSED_VARIABLE(size);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage& other) {
-    numext::swap(m_data, other.m_data);
-    Base::swap(other);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index /*size*/, Index rows, Index cols) {
-    Base::set(rows, cols);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index /*size*/, Index rows, Index cols) {
-    Base::set(rows, cols);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return m_data.array; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return m_data.array; }
-};
-// fixed-size storage with dynamic dimensions
-template <typename T, int Size, int Rows, int Cols, int Options>
-class DenseStorage<T, Size, Rows, Cols, Options, false> : internal::DenseStorageIndices<Rows, Cols> {
-  using Base = internal::DenseStorageIndices<Rows, Cols>;
-
-  internal::plain_array<T, Size, Options> m_data;
-
- public:
-  using Base::cols;
-  using Base::rows;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE DenseStorage() = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(const DenseStorage& other) : Base(other), m_data() {
-    Index size = other.size();
-    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
-    internal::smart_copy(other.m_data.array, other.m_data.array + size, m_data.array);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(DenseStorage&& other) : Base(other), m_data() {
-    Index size = other.size();
-    internal::smart_move(other.m_data.array, other.m_data.array + size, m_data.array);
-    other.resize(Size, 0, 0);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage& operator=(const DenseStorage& other) {
-    Base::set(other.rows(), other.cols());
-    Index size = other.size();
-    internal::smart_copy(other.m_data.array, other.m_data.array + size, m_data.array);
-    return *this;
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage& operator=(DenseStorage&& other) {
-    Base::set(other.rows(), other.cols());
-    Index size = other.size();
-    internal::smart_move(other.m_data.array, other.m_data.array + size, m_data.array);
-    other.resize(Size, 0, 0);
-    return *this;
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(Index size, Index rows, Index cols) : Base(rows, cols) {
-    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
-    EIGEN_UNUSED_VARIABLE(size);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage& other) {
-    Index thisSize = this->size();
-    Index otherSize = other.size();
-    Index commonSize = numext::mini(thisSize, otherSize);
-    std::swap_ranges(m_data.array, m_data.array + commonSize, other.m_data.array);
-    if (thisSize > otherSize)
-      internal::smart_move(m_data.array + commonSize, m_data.array + thisSize, other.m_data.array + commonSize);
-    else if (otherSize > thisSize)
-      internal::smart_move(other.m_data.array + commonSize, other.m_data.array + otherSize, m_data.array + commonSize);
-    Base::swap(other);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index /*size*/, Index rows, Index cols) {
-    Base::set(rows, cols);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index /*size*/, Index rows, Index cols) {
-    Base::set(rows, cols);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return m_data.array; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return m_data.array; }
-};
-// null matrix specialization
-template <typename T, int Rows, int Cols, int Options>
-class DenseStorage<T, 0, Rows, Cols, Options, true> : internal::DenseStorageIndices<Rows, Cols> {
-  using Base = internal::DenseStorageIndices<Rows, Cols>;
-
- public:
-  using Base::cols;
-  using Base::rows;
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage() = default;
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(const DenseStorage&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(DenseStorage&&) = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(Index size, Index rows, Index cols)
+      : Base(size, rows, cols) {}
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage& operator=(const DenseStorage&) = default;
+  // if DenseStorage meets the requirements of use_default_move, then use the move construction and move assignment
+  // operation defined in DenseStorage_impl, or the compiler-generated version if none is defined
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(DenseStorage&&) = default;
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage& operator=(DenseStorage&&) = default;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(Index /*size*/, Index rows, Index cols)
-      : Base(rows, cols) {}
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage& other) noexcept { Base::swap(other); }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index /*size*/, Index rows, Index cols) {
-    Base::set(rows, cols);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index /*size*/, Index rows, Index cols) {
-    Base::set(rows, cols);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return nullptr; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return nullptr; }
 };
-// dynamic matrix specialization
-template <typename T, int Rows, int Cols, int Options>
-class DenseStorage<T, Dynamic, Rows, Cols, Options, false> : internal::DenseStorageIndices<Rows, Cols> {
-  using Base = internal::DenseStorageIndices<Rows, Cols>;
-  static constexpr int Size = Dynamic;
-  static constexpr bool Align = (Options & DontAlign) == 0;
-
-  T* m_data;
+template <typename T, int Size, int Rows, int Cols, int Options>
+class DenseStorage<T, Size, Rows, Cols, Options, false>
+    : public internal::DenseStorage_impl<T, Size, Rows, Cols, Options> {
+  using Base = internal::DenseStorage_impl<T, Size, Rows, Cols, Options>;
 
  public:
-  using Base::cols;
-  using Base::rows;
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage() : m_data(nullptr) {}
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(const DenseStorage& other)
-      : Base(other), m_data(internal::conditional_aligned_new_auto<T, Align>(other.size())) {
-    Index size = other.size();
-    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
-    internal::smart_copy(other.m_data, other.m_data + size, m_data);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(DenseStorage&& other) noexcept
-      : Base(other), m_data(other.m_data) {
-    other.set(0, 0);
-    other.m_data = nullptr;
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage& operator=(const DenseStorage& other) {
-    Base::set(other.rows(), other.cols());
-    Index size = other.size();
-    m_data = internal::conditional_aligned_new_auto<T, Align>(size);
-    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
-    internal::smart_copy(other.m_data, other.m_data + size, m_data);
-    return *this;
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage& operator=(DenseStorage&& other) noexcept {
-    this->swap(other);
-    return *this;
-  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage() = default;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(const DenseStorage&) = default;
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(Index size, Index rows, Index cols)
-      : Base(rows, cols), m_data(internal::conditional_aligned_new_auto<T, Align>(size)) {
-    EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
+      : Base(size, rows, cols) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage& operator=(const DenseStorage&) = default;
+  // if DenseStorage does not meet the requirements of use_default_move, then defer to the copy construction and copy
+  // assignment behavior
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage(DenseStorage&& other)
+      : DenseStorage(static_cast<const DenseStorage&>(other)) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr DenseStorage& operator=(DenseStorage&& other) {
+    *this = other;
+    return *this;
   }
-  EIGEN_DEVICE_FUNC ~DenseStorage() {
-    Index size = this->size();
-    internal::conditional_aligned_delete_auto<T, Align>(m_data, size);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void swap(DenseStorage& other) noexcept {
-    numext::swap(m_data, other.m_data);
-    Base::swap(other);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void conservativeResize(Index size, Index rows, Index cols) {
-    Index oldSize = this->size();
-    m_data = internal::conditional_aligned_realloc_new_auto<T, Align>(m_data, size, oldSize);
-    Base::set(rows, cols);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr void resize(Index size, Index rows, Index cols) {
-    Index oldSize = this->size();
-    if (size != oldSize) {
-      internal::conditional_aligned_delete_auto<T, Align>(m_data, oldSize);
-      if (size > 0)  // >0 and not simply !=0 to let the compiler knows that size cannot be negative
-      {
-        m_data = internal::conditional_aligned_new_auto<T, Align>(size);
-        EIGEN_INTERNAL_DENSE_STORAGE_CTOR_PLUGIN({})
-      } else
-        m_data = nullptr;
-    }
-    Base::set(rows, cols);
-  }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr T* data() { return m_data; }
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr const T* data() const { return m_data; }
 };
+
 }  // end namespace Eigen
 
 #endif  // EIGEN_MATRIX_H

diff --git a/Eigen/src/Core/StlIterators.h b/Eigen/src/Core/StlIterators.h
index bb897f8..25d4575 100644
--- a/Eigen/src/Core/StlIterators.h
+++ b/Eigen/src/Core/StlIterators.h

@@ -325,7 +325,7 @@
  public:
   typedef Index difference_type;
   typedef typename XprType::Scalar value_type;
-#if __cplusplus >= 202002L
+#if EIGEN_CPLUSPLUS >= 202002L
   typedef std::conditional_t<XprType::InnerStrideAtCompileTime == 1, std::contiguous_iterator_tag,
                              std::random_access_iterator_tag>
       iterator_category;

diff --git a/Eigen/src/Core/arch/AVX/PacketMath.h b/Eigen/src/Core/arch/AVX/PacketMath.h
index 77fd7e6..1980e92 100644
--- a/Eigen/src/Core/arch/AVX/PacketMath.h
+++ b/Eigen/src/Core/arch/AVX/PacketMath.h

@@ -124,6 +124,7 @@
     HasRsqrt = 1,
     HasTanh = EIGEN_FAST_MATH,
     HasErf = EIGEN_FAST_MATH,
+    HasErfc = EIGEN_FAST_MATH,
     HasBlend = 1
   };
 };
@@ -144,6 +145,7 @@
 #endif
     HasTanh = EIGEN_FAST_MATH,
     HasLog = 1,
+    HasErfc = 1,
     HasExp = 1,
     HasSqrt = 1,
     HasRsqrt = 1,
@@ -659,6 +661,16 @@
   __m128i r = _mm_add_epi64(_mm256_castsi256_si128(a), _mm256_extractf128_si256(a, 1));
   return numext::bit_cast<uint64_t>(_mm_extract_epi64_0(r) + _mm_extract_epi64_1(r));
 }
+
+template <>
+EIGEN_STRONG_INLINE bool predux_any(const Packet4l& a) {
+  return _mm256_movemask_pd(_mm256_castsi256_pd(a)) != 0;
+}
+template <>
+EIGEN_STRONG_INLINE bool predux_any(const Packet4ul& a) {
+  return _mm256_movemask_pd(_mm256_castsi256_pd(a)) != 0;
+}
+
 #define MM256_SHUFFLE_EPI64(A, B, M) _mm256_shuffle_pd(_mm256_castsi256_pd(A), _mm256_castsi256_pd(B), M)
 EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet4l, 4>& kernel) {
   __m256d T0 = MM256_SHUFFLE_EPI64(kernel.packet[0], kernel.packet[1], 15);
@@ -2002,6 +2014,11 @@
 }
 
 template <>
+EIGEN_STRONG_INLINE bool predux_any(const Packet4d& x) {
+  return _mm256_movemask_pd(x) != 0;
+}
+
+template <>
 EIGEN_STRONG_INLINE bool predux_any(const Packet8i& x) {
   return _mm256_movemask_ps(_mm256_castsi256_ps(x)) != 0;
 }
@@ -2010,6 +2027,15 @@
   return _mm256_movemask_ps(_mm256_castsi256_ps(x)) != 0;
 }
 
+template <>
+EIGEN_STRONG_INLINE bool predux_any(const Packet8h& x) {
+  return _mm_movemask_epi8(x) != 0;
+}
+template <>
+EIGEN_STRONG_INLINE bool predux_any(const Packet8bf& x) {
+  return _mm_movemask_epi8(x) != 0;
+}
+
 EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet8f, 8>& kernel) {
   __m256 T0 = _mm256_unpacklo_ps(kernel.packet[0], kernel.packet[1]);
   __m256 T1 = _mm256_unpackhi_ps(kernel.packet[0], kernel.packet[1]);

diff --git a/Eigen/src/Core/arch/AVX512/PacketMath.h b/Eigen/src/Core/arch/AVX512/PacketMath.h
index 0a167c8..78d17d5 100644
--- a/Eigen/src/Core/arch/AVX512/PacketMath.h
+++ b/Eigen/src/Core/arch/AVX512/PacketMath.h

@@ -133,6 +133,7 @@
     HasReciprocal = EIGEN_FAST_MATH,
     HasTanh = EIGEN_FAST_MATH,
     HasErf = EIGEN_FAST_MATH,
+    HasErfc = EIGEN_FAST_MATH,
     HasCmp = 1,
     HasDiv = 1
   };
@@ -154,6 +155,7 @@
     HasExp = 1,
     HasATan = 1,
     HasTanh = EIGEN_FAST_MATH,
+    HasErfc = EIGEN_FAST_MATH,
     HasATanh = 1,
     HasCmp = 1,
     HasDiv = 1
@@ -1640,16 +1642,23 @@
 }
 
 template <>
-EIGEN_STRONG_INLINE bool predux_any(const Packet16f& x) {
-  Packet16i xi = _mm512_castps_si512(x);
-  __mmask16 tmp = _mm512_test_epi32_mask(xi, xi);
-  return !_mm512_kortestz(tmp, tmp);
+EIGEN_STRONG_INLINE bool predux_any(const Packet16f& a) {
+  return _mm512_reduce_or_epi32(_mm512_castps_si512(a)) != 0;
 }
 
 template <>
-EIGEN_STRONG_INLINE bool predux_any(const Packet16i& x) {
-  __mmask16 tmp = _mm512_test_epi32_mask(x, x);
-  return !_mm512_kortestz(tmp, tmp);
+EIGEN_STRONG_INLINE bool predux_any(const Packet16i& a) {
+  return _mm512_reduce_or_epi32(a) != 0;
+}
+
+template <>
+EIGEN_STRONG_INLINE bool predux_any(const Packet8d& a) {
+  return _mm512_reduce_or_epi64(_mm512_castpd_si512(a)) != 0;
+}
+
+template <>
+EIGEN_STRONG_INLINE bool predux_any(const Packet8l& a) {
+  return _mm512_reduce_or_epi64(a) != 0;
 }
 
 #define PACK_OUTPUT(OUTPUT, INPUT, INDEX, STRIDE) \

diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h
index 7401c0b..da26cd4 100644
--- a/Eigen/src/Core/arch/AltiVec/PacketMath.h
+++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h

@@ -193,6 +193,7 @@
 #endif
     HasTanh = EIGEN_FAST_MATH,
     HasErf = EIGEN_FAST_MATH,
+    HasErfc = EIGEN_FAST_MATH,
 #else
     HasSqrt = 0,
     HasRsqrt = 0,
@@ -3182,6 +3183,7 @@
     HasSin = EIGEN_FAST_MATH,
     HasCos = EIGEN_FAST_MATH,
     HasTanh = EIGEN_FAST_MATH,
+    HasErfc = EIGEN_FAST_MATH,
     HasATanh = 1,
     HasATan = 0,
     HasLog = 0,

diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
index a184931..4e441b4 100644
--- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
+++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h

@@ -1634,7 +1634,7 @@
 // This function splits x into the nearest integer n and fractional part r,
 // such that x = n + r holds exactly.
 template <typename Packet>
-EIGEN_STRONG_INLINE void absolute_split(const Packet& x, Packet& n, Packet& r) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void absolute_split(const Packet& x, Packet& n, Packet& r) {
   n = pround(x);
   r = psub(x, n);
 }
@@ -1642,7 +1642,7 @@
 // This function computes the sum {s, r}, such that x + y = s_hi + s_lo
 // holds exactly, and s_hi = fl(x+y), if |x| >= |y|.
 template <typename Packet>
-EIGEN_STRONG_INLINE void fast_twosum(const Packet& x, const Packet& y, Packet& s_hi, Packet& s_lo) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void fast_twosum(const Packet& x, const Packet& y, Packet& s_hi, Packet& s_lo) {
   s_hi = padd(x, y);
   const Packet t = psub(s_hi, x);
   s_lo = psub(y, t);
@@ -1654,11 +1654,18 @@
 // {p_hi, p_lo} such that x * y = p_hi + p_lo holds exactly and
 // p_hi = fl(x * y).
 template <typename Packet>
-EIGEN_STRONG_INLINE void twoprod(const Packet& x, const Packet& y, Packet& p_hi, Packet& p_lo) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x, const Packet& y, Packet& p_hi, Packet& p_lo) {
   p_hi = pmul(x, y);
   p_lo = pmsub(x, y, p_hi);
 }
 
+// A version of twoprod that takes x, y, and fl(x*y) as input and returns the p_lo such that
+// x * y = xy + p_lo holds exactly.
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet twoprod_low(const Packet& x, const Packet& y, const Packet& xy) {
+  return pmsub(x, y, xy);
+}
+
 #else
 
 // This function implements the Veltkamp splitting. Given a floating point
@@ -1667,7 +1674,7 @@
 // This is Algorithm 3 from Jean-Michel Muller, "Elementary Functions",
 // 3rd edition, Birkh\"auser, 2016.
 template <typename Packet>
-EIGEN_STRONG_INLINE void veltkamp_splitting(const Packet& x, Packet& x_hi, Packet& x_lo) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void veltkamp_splitting(const Packet& x, Packet& x_hi, Packet& x_lo) {
   typedef typename unpacket_traits<Packet>::type Scalar;
   EIGEN_CONSTEXPR int shift = (NumTraits<Scalar>::digits() + 1) / 2;
   const Scalar shift_scale = Scalar(uint64_t(1) << shift);  // Scalar constructor not necessarily constexpr.
@@ -1682,7 +1689,7 @@
 // {p_hi, p_lo} such that x * y = p_hi + p_lo holds exactly and
 // p_hi = fl(x * y).
 template <typename Packet>
-EIGEN_STRONG_INLINE void twoprod(const Packet& x, const Packet& y, Packet& p_hi, Packet& p_lo) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x, const Packet& y, Packet& p_hi, Packet& p_lo) {
   Packet x_hi, x_lo, y_hi, y_lo;
   veltkamp_splitting(x, x_hi, x_lo);
   veltkamp_splitting(y, y_hi, y_lo);
@@ -1694,6 +1701,21 @@
   p_lo = pmadd(x_lo, y_lo, p_lo);
 }
 
+// A version of twoprod that takes x, y, and fl(x*y) as input and returns the p_lo such that
+// x * y = xy + p_lo holds exactly.
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet twoprod_low(const Packet& x, const Packet& y, const Packet& xy) {
+  Packet x_hi, x_lo, y_hi, y_lo;
+  veltkamp_splitting(x, x_hi, x_lo);
+  veltkamp_splitting(y, y_hi, y_lo);
+
+  Packet p_lo = pmadd(x_hi, y_hi, pnegate(xy));
+  p_lo = pmadd(x_hi, y_lo, p_lo);
+  p_lo = pmadd(x_lo, y_hi, p_lo);
+  p_lo = pmadd(x_lo, y_lo, p_lo);
+  return p_lo;
+}
+
 #endif  // EIGEN_VECTORIZE_FMA
 
 // This function implements Dekker's algorithm for the addition
@@ -1703,8 +1725,8 @@
 // This is Algorithm 5 from Jean-Michel Muller, "Elementary Functions",
 // 3rd edition, Birkh\"auser, 2016.
 template <typename Packet>
-EIGEN_STRONG_INLINE void twosum(const Packet& x_hi, const Packet& x_lo, const Packet& y_hi, const Packet& y_lo,
-                                Packet& s_hi, Packet& s_lo) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twosum(const Packet& x_hi, const Packet& x_lo, const Packet& y_hi,
+                                                  const Packet& y_lo, Packet& s_hi, Packet& s_lo) {
   const Packet x_greater_mask = pcmp_lt(pabs(y_hi), pabs(x_hi));
   Packet r_hi_1, r_lo_1;
   fast_twosum(x_hi, y_hi, r_hi_1, r_lo_1);
@@ -1722,8 +1744,8 @@
 // This is a version of twosum for double word numbers,
 // which assumes that |x_hi| >= |y_hi|.
 template <typename Packet>
-EIGEN_STRONG_INLINE void fast_twosum(const Packet& x_hi, const Packet& x_lo, const Packet& y_hi, const Packet& y_lo,
-                                     Packet& s_hi, Packet& s_lo) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void fast_twosum(const Packet& x_hi, const Packet& x_lo, const Packet& y_hi,
+                                                       const Packet& y_lo, Packet& s_hi, Packet& s_lo) {
   Packet r_hi, r_lo;
   fast_twosum(x_hi, y_hi, r_hi, r_lo);
   const Packet s = padd(padd(y_lo, r_lo), x_lo);
@@ -1734,8 +1756,8 @@
 // double word number {y_hi, y_lo} number, with the assumption
 // that |x| >= |y_hi|.
 template <typename Packet>
-EIGEN_STRONG_INLINE void fast_twosum(const Packet& x, const Packet& y_hi, const Packet& y_lo, Packet& s_hi,
-                                     Packet& s_lo) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void fast_twosum(const Packet& x, const Packet& y_hi, const Packet& y_lo,
+                                                       Packet& s_hi, Packet& s_lo) {
   Packet r_hi, r_lo;
   fast_twosum(x, y_hi, r_hi, r_lo);
   const Packet s = padd(y_lo, r_lo);
@@ -1751,7 +1773,8 @@
 // This is Algorithm 7 from Jean-Michel Muller, "Elementary Functions",
 // 3rd edition, Birkh\"auser, 2016.
 template <typename Packet>
-EIGEN_STRONG_INLINE void twoprod(const Packet& x_hi, const Packet& x_lo, const Packet& y, Packet& p_hi, Packet& p_lo) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x_hi, const Packet& x_lo, const Packet& y,
+                                                   Packet& p_hi, Packet& p_lo) {
   Packet c_hi, c_lo1;
   twoprod(x_hi, y, c_hi, c_lo1);
   const Packet c_lo2 = pmul(x_lo, y);
@@ -1768,8 +1791,8 @@
 // of less than 2*2^{-2p}, where p is the number of significand bit
 // in the floating point type.
 template <typename Packet>
-EIGEN_STRONG_INLINE void twoprod(const Packet& x_hi, const Packet& x_lo, const Packet& y_hi, const Packet& y_lo,
-                                 Packet& p_hi, Packet& p_lo) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x_hi, const Packet& x_lo, const Packet& y_hi,
+                                                   const Packet& y_lo, Packet& p_hi, Packet& p_lo) {
   Packet p_hi_hi, p_hi_lo;
   twoprod(x_hi, x_lo, y_hi, p_hi_hi, p_hi_lo);
   Packet p_lo_hi, p_lo_lo;
@@ -1782,7 +1805,8 @@
 // for basic building blocks of double-word arithmetic", Joldes, Muller, & Popescu,
 // 2017. https://hal.archives-ouvertes.fr/hal-01351529
 template <typename Packet>
-void doubleword_div_fp(const Packet& x_hi, const Packet& x_lo, const Packet& y, Packet& z_hi, Packet& z_lo) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void doubleword_div_fp(const Packet& x_hi, const Packet& x_lo, const Packet& y,
+                                                             Packet& z_hi, Packet& z_lo) {
   const Packet t_hi = pdiv(x_hi, y);
   Packet pi_hi, pi_lo;
   twoprod(t_hi, y, pi_hi, pi_lo);
@@ -1797,7 +1821,7 @@
 template <typename Scalar>
 struct accurate_log2 {
   template <typename Packet>
-  EIGEN_STRONG_INLINE void operator()(const Packet& x, Packet& log2_x_hi, Packet& log2_x_lo) {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(const Packet& x, Packet& log2_x_hi, Packet& log2_x_lo) {
     log2_x_hi = plog2(x);
     log2_x_lo = pzero(x);
   }
@@ -1812,7 +1836,7 @@
 template <>
 struct accurate_log2<float> {
   template <typename Packet>
-  EIGEN_STRONG_INLINE void operator()(const Packet& z, Packet& log2_x_hi, Packet& log2_x_lo) {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(const Packet& z, Packet& log2_x_hi, Packet& log2_x_lo) {
     // The function log(1+x)/x is approximated in the interval
     // [1/sqrt(2)-1;sqrt(2)-1] by a degree 10 polynomial of the form
     //  Q(x) = (C0 + x * (C1 + x * (C2 + x * (C3 + x * P(x))))),
@@ -1892,7 +1916,7 @@
 template <>
 struct accurate_log2<double> {
   template <typename Packet>
-  EIGEN_STRONG_INLINE void operator()(const Packet& x, Packet& log2_x_hi, Packet& log2_x_lo) {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(const Packet& x, Packet& log2_x_hi, Packet& log2_x_lo) {
     // We use a transformation of variables:
     //    r = c * (x-1) / (x+1),
     // such that
@@ -1978,7 +2002,7 @@
 template <typename Scalar>
 struct fast_accurate_exp2 {
   template <typename Packet>
-  EIGEN_STRONG_INLINE Packet operator()(const Packet& x) {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet operator()(const Packet& x) {
     return generic_exp2(x);
   }
 };
@@ -1990,7 +2014,7 @@
 template <>
 struct fast_accurate_exp2<float> {
   template <typename Packet>
-  EIGEN_STRONG_INLINE Packet operator()(const Packet& x) {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet operator()(const Packet& x) {
     // This function approximates exp2(x) by a degree 6 polynomial of the form
     // Q(x) = 1 + x * (C + x * P(x)), where the degree 4 polynomial P(x) is evaluated in
     // single precision, and the remaining steps are evaluated with extra precision using
@@ -2047,7 +2071,7 @@
 template <>
 struct fast_accurate_exp2<double> {
   template <typename Packet>
-  EIGEN_STRONG_INLINE Packet operator()(const Packet& x) {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet operator()(const Packet& x) {
     // This function approximates exp2(x) by a degree 10 polynomial of the form
     // Q(x) = 1 + x * (C + x * P(x)), where the degree 8 polynomial P(x) is evaluated in
     // single precision, and the remaining steps are evaluated with extra precision using
@@ -2113,7 +2137,7 @@
 // TODO(rmlarsen): We should probably add this as a packet up 'ppow', to make it
 // easier to specialize or turn off for specific types and/or backends.x
 template <typename Packet>
-EIGEN_STRONG_INLINE Packet generic_pow_impl(const Packet& x, const Packet& y) {
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet generic_pow_impl(const Packet& x, const Packet& y) {
   typedef typename unpacket_traits<Packet>::type Scalar;
   // Split x into exponent e_x and mantissa m_x.
   Packet e_x;

diff --git a/Eigen/src/Core/arch/NEON/PacketMath.h b/Eigen/src/Core/arch/NEON/PacketMath.h
index 414f0f5..2f401fd 100644
--- a/Eigen/src/Core/arch/NEON/PacketMath.h
+++ b/Eigen/src/Core/arch/NEON/PacketMath.h

@@ -209,6 +209,7 @@
     HasRsqrt = 1,
     HasTanh = EIGEN_FAST_MATH,
     HasErf = EIGEN_FAST_MATH,
+    HasErfc = EIGEN_FAST_MATH,
     HasBessel = 0,  // Issues with accuracy.
     HasNdtri = 0
   };
@@ -5140,7 +5141,8 @@
     HasSqrt = 1,
     HasRsqrt = 1,
     HasTanh = EIGEN_FAST_MATH,
-    HasErf = 0
+    HasErf = 0,
+    HasErfc = EIGEN_FAST_MATH
   };
 };
 

diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h
index 37f6048..c749763 100644
--- a/Eigen/src/Core/arch/SSE/PacketMath.h
+++ b/Eigen/src/Core/arch/SSE/PacketMath.h

@@ -197,6 +197,7 @@
     HasRsqrt = 1,
     HasTanh = EIGEN_FAST_MATH,
     HasErf = EIGEN_FAST_MATH,
+    HasErfc = EIGEN_FAST_MATH,
     HasBlend = 1,
     HasSign = 0  // The manually vectorized version is slightly slower for SSE.
   };
@@ -216,6 +217,7 @@
     HasCos = EIGEN_FAST_MATH,
     HasTanh = EIGEN_FAST_MATH,
     HasLog = 1,
+    HasErfc = EIGEN_FAST_MATH,
     HasExp = 1,
     HasSqrt = 1,
     HasRsqrt = 1,

diff --git a/Eigen/src/Core/arch/SVE/PacketMath.h b/Eigen/src/Core/arch/SVE/PacketMath.h
index 51bbfe0..952d756 100644
--- a/Eigen/src/Core/arch/SVE/PacketMath.h
+++ b/Eigen/src/Core/arch/SVE/PacketMath.h

@@ -360,7 +360,8 @@
     HasExp = 1,
     HasSqrt = 1,
     HasTanh = EIGEN_FAST_MATH,
-    HasErf = EIGEN_FAST_MATH
+    HasErf = EIGEN_FAST_MATH,
+    HasErfc = EIGEN_FAST_MATH
   };
 };
 

diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h
index d59071f..2acdd9d 100644
--- a/Eigen/src/Core/util/Memory.h
+++ b/Eigen/src/Core/util/Memory.h

@@ -147,7 +147,7 @@
   check_that_malloc_is_allowed();
   EIGEN_USING_STD(malloc)
   void* original = malloc(size + alignment);
-  if (original == 0) return 0;
+  if (original == nullptr) return nullptr;
   uint8_t offset = static_cast<uint8_t>(alignment - (reinterpret_cast<std::size_t>(original) & (alignment - 1)));
   void* aligned = static_cast<void*>(static_cast<uint8_t*>(original) + offset);
   *(static_cast<uint8_t*>(aligned) - 1) = offset;
@@ -391,7 +391,8 @@
 
 template <typename T>
 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE void check_size_for_overflow(std::size_t size) {
-  if (size > std::size_t(-1) / sizeof(T)) throw_std_bad_alloc();
+  constexpr std::size_t max_elements = PTRDIFF_MAX / sizeof(T);
+  if (size > max_elements) throw_std_bad_alloc();
 }
 
 /** \internal Allocates \a size objects of type T. The returned pointer is guaranteed to have 16 bytes alignment.
@@ -473,7 +474,7 @@
 
 template <typename T, bool Align>
 EIGEN_DEVICE_FUNC inline T* conditional_aligned_new_auto(std::size_t size) {
-  if (size == 0) return 0;  // short-cut. Also fixes Bug 884
+  if (size == 0) return nullptr;  // short-cut. Also fixes Bug 884
   check_size_for_overflow<T>(size);
   T* result = static_cast<T*>(conditional_aligned_malloc<Align>(sizeof(T) * size));
   if (NumTraits<T>::RequireInitialization) {

diff --git a/test/MovableScalar.h b/test/MovableScalar.h
index c8bf546..6a046de 100644
--- a/test/MovableScalar.h
+++ b/test/MovableScalar.h

@@ -10,24 +10,42 @@
 #ifndef EIGEN_MISC_MOVABLE_SCALAR_H
 #define EIGEN_MISC_MOVABLE_SCALAR_H
 
-#include <vector>
-
 namespace Eigen {
-template <typename Scalar, typename Base = std::vector<Scalar>>
-struct MovableScalar : public Base {
-  MovableScalar() = default;
-  ~MovableScalar() = default;
-  MovableScalar(const MovableScalar&) = default;
-  MovableScalar(MovableScalar&& other) = default;
-  MovableScalar& operator=(const MovableScalar&) = default;
-  MovableScalar& operator=(MovableScalar&& other) = default;
-  MovableScalar(Scalar scalar) : Base(100, scalar) {}
+template <typename Scalar>
+struct MovableScalar {
+  MovableScalar() : m_data(new Scalar) {}
+  ~MovableScalar() { delete m_data; }
+  MovableScalar(const MovableScalar& other) : m_data(new Scalar) { set(other.get()); }
+  MovableScalar(MovableScalar&& other) noexcept : m_data(other.m_data) { other.m_data = nullptr; }
+  MovableScalar& operator=(const MovableScalar& other) {
+    set(other.get());
+    return *this;
+  }
+  MovableScalar& operator=(MovableScalar&& other) noexcept {
+    m_data = other.m_data;
+    other.m_data = nullptr;
+    return *this;
+  }
+  MovableScalar(const Scalar& scalar) : m_data(new Scalar) { set(scalar); }
 
-  operator Scalar() const { return this->size() > 0 ? this->back() : Scalar(); }
+  operator Scalar() const { return get(); }
+
+ private:
+  void set(const Scalar& value) {
+    eigen_assert(m_data != nullptr);
+    // suppress compiler warnings
+    if (m_data != nullptr) *m_data = value;
+  }
+  Scalar get() const {
+    eigen_assert(m_data != nullptr);
+    // suppress compiler warnings
+    return m_data == nullptr ? Scalar() : *m_data;
+  }
+  Scalar* m_data = nullptr;
 };
 
 template <typename Scalar>
-struct NumTraits<MovableScalar<Scalar>> : GenericNumTraits<Scalar> {
+struct NumTraits<MovableScalar<Scalar>> : NumTraits<Scalar> {
   enum { RequireInitialization = 1 };
 };
 

diff --git a/test/dense_storage.cpp b/test/dense_storage.cpp
index 5d0083f..0f6ab64 100644
--- a/test/dense_storage.cpp
+++ b/test/dense_storage.cpp

@@ -11,29 +11,38 @@
 
 #include "main.h"
 #include "AnnoyingScalar.h"
+#include "MovableScalar.h"
 #include "SafeScalar.h"
 
 #include <Eigen/Core>
 
 using DenseStorageD3x3 = Eigen::DenseStorage<double, 9, 3, 3, 0>;
-static_assert(std::is_trivially_move_constructible<DenseStorageD3x3>::value,
-              "DenseStorage not trivially_move_constructible");
-static_assert(std::is_trivially_move_assignable<DenseStorageD3x3>::value, "DenseStorage not trivially_move_assignable");
 #if !defined(EIGEN_DENSE_STORAGE_CTOR_PLUGIN)
 static_assert(std::is_trivially_copy_constructible<DenseStorageD3x3>::value,
               "DenseStorage not trivially_copy_constructible");
+static_assert(std::is_trivially_move_constructible<DenseStorageD3x3>::value,
+              "DenseStorage not trivially_move_constructible");
 static_assert(std::is_trivially_copy_assignable<DenseStorageD3x3>::value, "DenseStorage not trivially_copy_assignable");
-static_assert(std::is_trivially_copyable<DenseStorageD3x3>::value, "DenseStorage not trivially_copyable");
+static_assert(std::is_trivially_move_assignable<DenseStorageD3x3>::value, "DenseStorage not trivially_move_assignable");
 #endif
-
+// all plain object types conform to standard layout
+static_assert(std::is_standard_layout<Matrix4f>::value, "Matrix4f not standard_layout");
+static_assert(std::is_standard_layout<Array4f>::value, "Array4f not standard_layout");
+static_assert(std::is_standard_layout<VectorXf>::value, "VectorXf not standard_layout");
+static_assert(std::is_standard_layout<ArrayXf>::value, "ArrayXf not standard_layout");
+static_assert(std::is_standard_layout<MatrixXf>::value, "MatrixXf not standard_layout");
+static_assert(std::is_standard_layout<ArrayXXf>::value, "ArrayXXf not standard_layout");
+// all fixed-size, fixed-dimension plain object types are trivially default constructible
+static_assert(std::is_trivially_default_constructible<Matrix4f>::value, "Matrix4f not trivially_default_constructible");
+static_assert(std::is_trivially_default_constructible<Array4f>::value, "Array4f not trivially_default_constructible");
+// all fixed-size, fixed-dimension plain object types are trivially move constructible
 static_assert(std::is_trivially_move_constructible<Matrix4f>::value, "Matrix4f not trivially_move_constructible");
 static_assert(std::is_trivially_move_constructible<Array4f>::value, "Array4f not trivially_move_constructible");
 #if !defined(EIGEN_DENSE_STORAGE_CTOR_PLUGIN)
+// all fixed-size, fixed-dimension plain object types are trivially copy constructible
 static_assert(std::is_trivially_copy_constructible<Matrix4f>::value, "Matrix4f not trivially_copy_constructible");
 static_assert(std::is_trivially_copy_constructible<Array4f>::value, "Array4f not trivially_copy_constructible");
 #endif
-static_assert(std::is_trivially_default_constructible<Matrix4f>::value, "Matrix4f not trivially_default_constructible");
-static_assert(std::is_trivially_default_constructible<Array4f>::value, "Array4f not trivially_default_constructible");
 
 template <typename T, int Size, int Rows, int Cols>
 void dense_storage_copy(int rows, int cols) {
@@ -42,7 +51,7 @@
   const int size = rows * cols;
   DenseStorageType reference(size, rows, cols);
   T* raw_reference = reference.data();
-  for (int i = 0; i < size; ++i) raw_reference[i] = static_cast<T>(i);
+  for (int i = 0; i < size; ++i) raw_reference[i] = internal::random<T>();
 
   DenseStorageType copied_reference(reference);
   const T* raw_copied_reference = copied_reference.data();
@@ -56,7 +65,7 @@
   const int size = rows * cols;
   DenseStorageType reference(size, rows, cols);
   T* raw_reference = reference.data();
-  for (int i = 0; i < size; ++i) raw_reference[i] = static_cast<T>(i);
+  for (int i = 0; i < size; ++i) raw_reference[i] = internal::random<T>();
 
   DenseStorageType copied_reference;
   copied_reference = reference;
@@ -65,30 +74,25 @@
 }
 
 template <typename T, int Size, int Rows, int Cols>
-void dense_storage_swap(int rows0, int cols0, int rows1, int cols1) {
+void dense_storage_swap(int rowsa, int colsa, int rowsb, int colsb) {
   typedef DenseStorage<T, Size, Rows, Cols, 0> DenseStorageType;
 
-  const int size0 = rows0 * cols0;
-  DenseStorageType a(size0, rows0, cols0);
-  for (int i = 0; i < size0; ++i) {
-    a.data()[i] = static_cast<T>(i);
-  }
+  const int sizea = rowsa * colsa;
+  ArrayX<T> referencea(sizea);
+  referencea.setRandom();
+  DenseStorageType a(sizea, rowsa, colsa);
+  for (int i = 0; i < sizea; ++i) a.data()[i] = referencea(i);
 
-  const int size1 = rows1 * cols1;
-  DenseStorageType b(size1, rows1, cols1);
-  for (int i = 0; i < size1; ++i) {
-    b.data()[i] = static_cast<T>(-i);
-  }
+  const int sizeb = rowsb * colsb;
+  ArrayX<T> referenceb(sizeb);
+  referenceb.setRandom();
+  DenseStorageType b(sizeb, rowsb, colsb);
+  for (int i = 0; i < sizeb; ++i) b.data()[i] = referenceb(i);
 
   a.swap(b);
 
-  for (int i = 0; i < size0; ++i) {
-    VERIFY_IS_EQUAL(b.data()[i], static_cast<T>(i));
-  }
-
-  for (int i = 0; i < size1; ++i) {
-    VERIFY_IS_EQUAL(a.data()[i], static_cast<T>(-i));
-  }
+  for (int i = 0; i < sizea; i++) VERIFY_IS_EQUAL(b.data()[i], referencea(i));
+  for (int i = 0; i < sizeb; i++) VERIFY_IS_EQUAL(a.data()[i], referenceb(i));
 }
 
 template <typename T, int Size, std::size_t Alignment>
@@ -104,12 +108,12 @@
   };
   VERIFY_IS_EQUAL(std::alignment_of<Nested1>::value, Alignment);
 
-  VERIFY_IS_EQUAL((std::alignment_of<internal::plain_array<T, Size, AutoAlign, Alignment> >::value), Alignment);
+  VERIFY_IS_EQUAL((std::alignment_of<internal::plain_array<T, Size, AutoAlign, Alignment>>::value), Alignment);
 
   const std::size_t default_alignment = internal::compute_default_alignment<T, Size>::value;
   if (default_alignment > 0) {
-    VERIFY_IS_EQUAL((std::alignment_of<DenseStorage<T, Size, 1, 1, AutoAlign> >::value), default_alignment);
-    VERIFY_IS_EQUAL((std::alignment_of<Matrix<T, Size, 1, AutoAlign> >::value), default_alignment);
+    VERIFY_IS_EQUAL((std::alignment_of<DenseStorage<T, Size, 1, 1, AutoAlign>>::value), default_alignment);
+    VERIFY_IS_EQUAL((std::alignment_of<Matrix<T, Size, 1, AutoAlign>>::value), default_alignment);
     struct Nested2 {
       Matrix<T, Size, 1, AutoAlign> mat;
     };
@@ -185,11 +189,90 @@
   dense_storage_alignment<T, 16, 64>();
 }
 
+template <typename PlainType>
+void plaintype_tests() {
+  constexpr int RowsAtCompileTime = PlainType::RowsAtCompileTime;
+  constexpr int ColsAtCompileTime = PlainType::ColsAtCompileTime;
+  constexpr int MaxRowsAtCompileTime = PlainType::MaxRowsAtCompileTime;
+  constexpr int MaxColsAtCompileTime = PlainType::MaxColsAtCompileTime;
+  const Index expectedDefaultRows = RowsAtCompileTime == Dynamic ? 0 : RowsAtCompileTime;
+  const Index expectedDefaultCols = ColsAtCompileTime == Dynamic ? 0 : ColsAtCompileTime;
+  const Index minRows = RowsAtCompileTime == Dynamic ? 0 : RowsAtCompileTime;
+  const Index minCols = ColsAtCompileTime == Dynamic ? 0 : ColsAtCompileTime;
+  const Index maxRows = MaxRowsAtCompileTime == Dynamic ? 100 : MaxRowsAtCompileTime;
+  const Index maxCols = MaxColsAtCompileTime == Dynamic ? 100 : MaxColsAtCompileTime;
+  const Index rows = internal::random<Index>(minRows, maxRows);
+  const Index cols = internal::random<Index>(minCols, maxCols);
+  // default construction
+  PlainType m0;
+  VERIFY_IS_EQUAL(m0.rows(), expectedDefaultRows);
+  VERIFY_IS_EQUAL(m0.cols(), expectedDefaultCols);
+  m0.resize(rows, cols);
+  m0.setRandom();
+  // copy construction
+  PlainType m1(m0);
+  VERIFY_IS_EQUAL(m1.rows(), m0.rows());
+  VERIFY_IS_EQUAL(m1.cols(), m0.cols());
+  VERIFY_IS_CWISE_EQUAL(m1, m0);
+  // move construction
+  PlainType m2(std::move(m1));
+  VERIFY_IS_EQUAL(m2.rows(), m0.rows());
+  VERIFY_IS_EQUAL(m2.cols(), m0.cols());
+  VERIFY_IS_CWISE_EQUAL(m2, m0);
+  // check that object is usable after move construction
+  m1.resize(minRows, minCols);
+  m1.setRandom();
+  // copy assignment
+  m1 = m0;
+  VERIFY_IS_EQUAL(m1.rows(), m0.rows());
+  VERIFY_IS_EQUAL(m1.cols(), m0.cols());
+  VERIFY_IS_CWISE_EQUAL(m1, m0);
+  // move assignment
+  m2.resize(minRows, minCols);
+  m2.setRandom();
+  m2 = std::move(m1);
+  VERIFY_IS_EQUAL(m2.rows(), m0.rows());
+  VERIFY_IS_EQUAL(m2.cols(), m0.cols());
+  VERIFY_IS_CWISE_EQUAL(m2, m0);
+  // check that object is usable after move assignment
+  m1.resize(minRows, minCols);
+  m1.setRandom();
+  m1 = m2;
+  VERIFY_IS_EQUAL(m1.rows(), m0.rows());
+  VERIFY_IS_EQUAL(m1.cols(), m0.cols());
+  VERIFY_IS_CWISE_EQUAL(m1, m0);
+}
+
 EIGEN_DECLARE_TEST(dense_storage) {
   dense_storage_tests<int>();
   dense_storage_tests<float>();
-  dense_storage_tests<SafeScalar<float> >();
+  dense_storage_tests<SafeScalar<float>>();
+  dense_storage_tests<MovableScalar<float>>();
   dense_storage_tests<AnnoyingScalar>();
+  for (int i = 0; i < g_repeat; i++) {
+    plaintype_tests<Matrix<float, 0, 0, ColMajor>>();
+    plaintype_tests<Matrix<float, Dynamic, Dynamic, ColMajor, 0, 0>>();
+
+    plaintype_tests<Matrix<float, 16, 16, ColMajor>>();
+    plaintype_tests<Matrix<float, 16, Dynamic, ColMajor>>();
+    plaintype_tests<Matrix<float, Dynamic, Dynamic, ColMajor>>();
+    plaintype_tests<Matrix<float, Dynamic, Dynamic, ColMajor, 16, 16>>();
+
+    plaintype_tests<Matrix<SafeScalar<float>, 16, 16, ColMajor>>();
+    plaintype_tests<Matrix<SafeScalar<float>, 16, Dynamic, ColMajor>>();
+    plaintype_tests<Matrix<SafeScalar<float>, Dynamic, Dynamic, ColMajor>>();
+    plaintype_tests<Matrix<SafeScalar<float>, Dynamic, Dynamic, ColMajor, 16, 16>>();
+
+    plaintype_tests<Matrix<MovableScalar<float>, 16, 16, ColMajor>>();
+    plaintype_tests<Matrix<MovableScalar<float>, 16, Dynamic, ColMajor>>();
+    plaintype_tests<Matrix<MovableScalar<float>, Dynamic, Dynamic, ColMajor>>();
+    plaintype_tests<Matrix<MovableScalar<float>, Dynamic, Dynamic, ColMajor, 16, 16>>();
+
+    plaintype_tests<Matrix<AnnoyingScalar, 16, 16, ColMajor>>();
+    plaintype_tests<Matrix<AnnoyingScalar, 16, Dynamic, ColMajor>>();
+    plaintype_tests<Matrix<AnnoyingScalar, Dynamic, Dynamic, ColMajor>>();
+    plaintype_tests<Matrix<AnnoyingScalar, Dynamic, Dynamic, ColMajor, 16, 16>>();
+  }
 }
 
 #undef EIGEN_TESTING_PLAINOBJECT_CTOR
\ No newline at end of file

diff --git a/unsupported/Eigen/CXX11/Tensor b/unsupported/Eigen/CXX11/Tensor
index 1b6cc7e..290a0c0 100644
--- a/unsupported/Eigen/CXX11/Tensor
+++ b/unsupported/Eigen/CXX11/Tensor

@@ -109,6 +109,7 @@
 #include "src/Tensor/TensorMorphing.h"
 #include "src/Tensor/TensorPadding.h"
 #include "src/Tensor/TensorReverse.h"
+#include "src/Tensor/TensorRoll.h"
 #include "src/Tensor/TensorShuffling.h"
 #include "src/Tensor/TensorStriding.h"
 #include "src/Tensor/TensorCustomOp.h"

diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h b/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h
index 2c2c781..fc3f3b7 100644
--- a/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h
+++ b/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h

@@ -946,6 +946,11 @@
     reverse(const ReverseDimensions& rev) const {
       return TensorReverseOp<const ReverseDimensions, const Derived>(derived(), rev);
     }
+    template <typename Rolls> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+    const TensorRollOp<const Rolls, const Derived>
+    roll(const Rolls& rolls) const {
+      return TensorRollOp<const Rolls, const Derived>(derived(), rolls);
+    }
     template <typename PaddingDimensions> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
     const TensorPaddingOp<const PaddingDimensions, const Derived>
     pad(const PaddingDimensions& padding) const {
@@ -1166,6 +1171,17 @@
       return TensorReverseOp<const ReverseDimensions, Derived>(derived(), rev);
     }
 
+    template <typename Rolls> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+    const TensorRollOp<const Rolls, const Derived>
+    roll(const Rolls& roll) const {
+      return TensorRollOp<const Rolls, const Derived>(derived(), roll);
+    }
+    template <typename Rolls> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+    TensorRollOp<const Rolls, Derived>
+    roll(const Rolls& roll) {
+      return TensorRollOp<const Rolls, Derived>(derived(), roll);
+    }
+
     template <typename Shuffle> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
     const TensorShufflingOp<const Shuffle, const Derived>
     shuffle(const Shuffle& shfl) const {

diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorForwardDeclarations.h b/unsupported/Eigen/CXX11/src/Tensor/TensorForwardDeclarations.h
index 3bc3a5b..49c20a4 100644
--- a/unsupported/Eigen/CXX11/src/Tensor/TensorForwardDeclarations.h
+++ b/unsupported/Eigen/CXX11/src/Tensor/TensorForwardDeclarations.h

@@ -111,6 +111,8 @@
 class TensorSlicingOp;
 template <typename ReverseDimensions, typename XprType>
 class TensorReverseOp;
+template <typename Rolls, typename XprType>
+class TensorRollOp;
 template <typename PaddingDimensions, typename XprType>
 class TensorPaddingOp;
 template <typename Shuffle, typename XprType>

diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorRoll.h b/unsupported/Eigen/CXX11/src/Tensor/TensorRoll.h
new file mode 100644
index 0000000..d5b203a
--- /dev/null
+++ b/unsupported/Eigen/CXX11/src/Tensor/TensorRoll.h

@@ -0,0 +1,361 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2024 Tobias Wood tobias@spinicist.org.uk
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_CXX11_TENSOR_TENSOR_ROLL_H
+#define EIGEN_CXX11_TENSOR_TENSOR_ROLL_H
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+namespace Eigen {
+
+/** \class TensorRoll
+ * \ingroup CXX11_Tensor_Module
+ *
+ * \brief Tensor roll (circular shift) elements class.
+ *
+ */
+namespace internal {
+template <typename RollDimensions, typename XprType>
+struct traits<TensorRollOp<RollDimensions, XprType> > : public traits<XprType> {
+  typedef typename XprType::Scalar Scalar;
+  typedef traits<XprType> XprTraits;
+  typedef typename XprTraits::StorageKind StorageKind;
+  typedef typename XprTraits::Index Index;
+  typedef typename XprType::Nested Nested;
+  typedef std::remove_reference_t<Nested> Nested_;
+  static constexpr int NumDimensions = XprTraits::NumDimensions;
+  static constexpr int Layout = XprTraits::Layout;
+  typedef typename XprTraits::PointerType PointerType;
+};
+
+template <typename RollDimensions, typename XprType>
+struct eval<TensorRollOp<RollDimensions, XprType>, Eigen::Dense> {
+  typedef const TensorRollOp<RollDimensions, XprType>& type;
+};
+
+template <typename RollDimensions, typename XprType>
+struct nested<TensorRollOp<RollDimensions, XprType>, 1, typename eval<TensorRollOp<RollDimensions, XprType> >::type> {
+  typedef TensorRollOp<RollDimensions, XprType> type;
+};
+
+}  // end namespace internal
+
+template <typename RollDimensions, typename XprType>
+class TensorRollOp : public TensorBase<TensorRollOp<RollDimensions, XprType>, WriteAccessors> {
+ public:
+  typedef TensorBase<TensorRollOp<RollDimensions, XprType>, WriteAccessors> Base;
+  typedef typename Eigen::internal::traits<TensorRollOp>::Scalar Scalar;
+  typedef typename Eigen::NumTraits<Scalar>::Real RealScalar;
+  typedef typename XprType::CoeffReturnType CoeffReturnType;
+  typedef typename Eigen::internal::nested<TensorRollOp>::type Nested;
+  typedef typename Eigen::internal::traits<TensorRollOp>::StorageKind StorageKind;
+  typedef typename Eigen::internal::traits<TensorRollOp>::Index Index;
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorRollOp(const XprType& expr, const RollDimensions& roll_dims)
+      : m_xpr(expr), m_roll_dims(roll_dims) {}
+
+  EIGEN_DEVICE_FUNC const RollDimensions& roll() const { return m_roll_dims; }
+
+  EIGEN_DEVICE_FUNC const internal::remove_all_t<typename XprType::Nested>& expression() const { return m_xpr; }
+
+  EIGEN_TENSOR_INHERIT_ASSIGNMENT_OPERATORS(TensorRollOp)
+
+ protected:
+  typename XprType::Nested m_xpr;
+  const RollDimensions m_roll_dims;
+};
+
+// Eval as rvalue
+template <typename RollDimensions, typename ArgType, typename Device>
+struct TensorEvaluator<const TensorRollOp<RollDimensions, ArgType>, Device> {
+  typedef TensorRollOp<RollDimensions, ArgType> XprType;
+  typedef typename XprType::Index Index;
+  static constexpr int NumDims = internal::array_size<RollDimensions>::value;
+  typedef DSizes<Index, NumDims> Dimensions;
+  typedef typename XprType::Scalar Scalar;
+  typedef typename XprType::CoeffReturnType CoeffReturnType;
+  typedef typename PacketType<CoeffReturnType, Device>::type PacketReturnType;
+  static constexpr int PacketSize = PacketType<CoeffReturnType, Device>::size;
+  typedef StorageMemory<CoeffReturnType, Device> Storage;
+  typedef typename Storage::Type EvaluatorPointerType;
+
+  static constexpr int Layout = TensorEvaluator<ArgType, Device>::Layout;
+  enum {
+    IsAligned = false,
+    PacketAccess = TensorEvaluator<ArgType, Device>::PacketAccess,
+    BlockAccess = NumDims > 0,
+    PreferBlockAccess = true,
+    CoordAccess = false,  // to be implemented
+    RawAccess = false
+  };
+
+  typedef internal::TensorIntDivisor<Index> IndexDivisor;
+
+  //===- Tensor block evaluation strategy (see TensorBlock.h) -------------===//
+  using TensorBlockDesc = internal::TensorBlockDescriptor<NumDims, Index>;
+  using TensorBlockScratch = internal::TensorBlockScratchAllocator<Device>;
+  using ArgTensorBlock = typename TensorEvaluator<const ArgType, Device>::TensorBlock;
+  using TensorBlock = typename internal::TensorMaterializedBlock<CoeffReturnType, NumDims, Layout, Index>;
+  //===--------------------------------------------------------------------===//
+
+  EIGEN_STRONG_INLINE TensorEvaluator(const XprType& op, const Device& device)
+      : m_impl(op.expression(), device), m_rolls(op.roll()), m_device(device) {
+    EIGEN_STATIC_ASSERT((NumDims > 0), Must_Have_At_Least_One_Dimension_To_Roll);
+
+    // Compute strides
+    m_dimensions = m_impl.dimensions();
+    if (static_cast<int>(Layout) == static_cast<int>(ColMajor)) {
+      m_strides[0] = 1;
+      for (int i = 1; i < NumDims; ++i) {
+        m_strides[i] = m_strides[i - 1] * m_dimensions[i - 1];
+        if (m_strides[i] > 0) m_fast_strides[i] = IndexDivisor(m_strides[i]);
+      }
+    } else {
+      m_strides[NumDims - 1] = 1;
+      for (int i = NumDims - 2; i >= 0; --i) {
+        m_strides[i] = m_strides[i + 1] * m_dimensions[i + 1];
+        if (m_strides[i] > 0) m_fast_strides[i] = IndexDivisor(m_strides[i]);
+      }
+    }
+  }
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Dimensions& dimensions() const { return m_dimensions; }
+
+  EIGEN_STRONG_INLINE bool evalSubExprsIfNeeded(EvaluatorPointerType) {
+    m_impl.evalSubExprsIfNeeded(nullptr);
+    return true;
+  }
+
+#ifdef EIGEN_USE_THREADS
+  template <typename EvalSubExprsCallback>
+  EIGEN_STRONG_INLINE void evalSubExprsIfNeededAsync(EvaluatorPointerType, EvalSubExprsCallback done) {
+    m_impl.evalSubExprsIfNeededAsync(nullptr, [done](bool) { done(true); });
+  }
+#endif  // EIGEN_USE_THREADS
+
+  EIGEN_STRONG_INLINE void cleanup() { m_impl.cleanup(); }
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index roll(Index const i, Index const r, Index const n) const {
+    auto const tmp = (i + r) % n;
+    if (tmp < 0) {
+      return tmp + n;
+    } else {
+      return tmp;
+    }
+  }
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE array<Index, NumDims> rollCoords(array<Index, NumDims> const& coords) const {
+    array<Index, NumDims> rolledCoords;
+    for (int id = 0; id < NumDims; id++) {
+      eigen_assert(coords[id] < m_dimensions[id]);
+      rolledCoords[id] = roll(coords[id], m_rolls[id], m_dimensions[id]);
+    }
+    return rolledCoords;
+  }
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index rollIndex(Index index) const {
+    eigen_assert(index < dimensions().TotalSize());
+    Index rolledIndex = 0;
+    if (static_cast<int>(Layout) == static_cast<int>(ColMajor)) {
+      EIGEN_UNROLL_LOOP
+      for (int i = NumDims - 1; i > 0; --i) {
+        Index idx = index / m_fast_strides[i];
+        index -= idx * m_strides[i];
+        rolledIndex += roll(idx, m_rolls[i], m_dimensions[i]) * m_strides[i];
+      }
+      rolledIndex += roll(index, m_rolls[0], m_dimensions[0]);
+    } else {
+      EIGEN_UNROLL_LOOP
+      for (int i = 0; i < NumDims - 1; ++i) {
+        Index idx = index / m_fast_strides[i];
+        index -= idx * m_strides[i];
+        rolledIndex += roll(idx, m_rolls[i], m_dimensions[i]) * m_strides[i];
+      }
+      rolledIndex += roll(index, m_rolls[NumDims - 1], m_dimensions[NumDims - 1]);
+    }
+    return rolledIndex;
+  }
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index index) const {
+    return m_impl.coeff(rollIndex(index));
+  }
+
+  template <int LoadMode>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PacketReturnType packet(Index index) const {
+    eigen_assert(index + PacketSize - 1 < dimensions().TotalSize());
+    EIGEN_ALIGN_MAX std::remove_const_t<CoeffReturnType> values[PacketSize];
+    EIGEN_UNROLL_LOOP
+    for (int i = 0; i < PacketSize; ++i) {
+      values[i] = coeff(index + i);
+    }
+    PacketReturnType rslt = internal::pload<PacketReturnType>(values);
+    return rslt;
+  }
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE internal::TensorBlockResourceRequirements getResourceRequirements() const {
+    const size_t target_size = m_device.lastLevelCacheSize();
+    return internal::TensorBlockResourceRequirements::skewed<Scalar>(target_size).addCostPerCoeff({0, 0, 24});
+  }
+
+  struct BlockIteratorState {
+    Index stride;
+    Index span;
+    Index size;
+    Index count;
+  };
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorBlock block(TensorBlockDesc& desc, TensorBlockScratch& scratch,
+                                                          bool /*root_of_expr_ast*/ = false) const {
+    static const bool is_col_major = static_cast<int>(Layout) == static_cast<int>(ColMajor);
+
+    // Compute spatial coordinates for the first block element.
+    array<Index, NumDims> coords;
+    extract_coordinates(desc.offset(), coords);
+    array<Index, NumDims> initial_coords = coords;
+    Index offset = 0;  // Offset in the output block buffer.
+
+    // Initialize output block iterator state. Dimension in this array are
+    // always in inner_most -> outer_most order (col major layout).
+    array<BlockIteratorState, NumDims> it;
+    for (int i = 0; i < NumDims; ++i) {
+      const int dim = is_col_major ? i : NumDims - 1 - i;
+      it[i].size = desc.dimension(dim);
+      it[i].stride = i == 0 ? 1 : (it[i - 1].size * it[i - 1].stride);
+      it[i].span = it[i].stride * (it[i].size - 1);
+      it[i].count = 0;
+    }
+    eigen_assert(it[0].stride == 1);
+
+    // Prepare storage for the materialized generator result.
+    const typename TensorBlock::Storage block_storage = TensorBlock::prepareStorage(desc, scratch);
+    CoeffReturnType* block_buffer = block_storage.data();
+
+    static const int inner_dim = is_col_major ? 0 : NumDims - 1;
+    const Index inner_dim_size = it[0].size;
+
+    while (it[NumDims - 1].count < it[NumDims - 1].size) {
+      Index i = 0;
+      for (; i < inner_dim_size; ++i) {
+        auto const rolled = rollCoords(coords);
+        auto const index = is_col_major ? m_dimensions.IndexOfColMajor(rolled) : m_dimensions.IndexOfRowMajor(rolled);
+        *(block_buffer + offset + i) = m_impl.coeff(index);
+        coords[inner_dim]++;
+      }
+      coords[inner_dim] = initial_coords[inner_dim];
+
+      if (NumDims == 1) break;  // For the 1d tensor we need to generate only one inner-most dimension.
+
+      // Update offset.
+      for (i = 1; i < NumDims; ++i) {
+        if (++it[i].count < it[i].size) {
+          offset += it[i].stride;
+          coords[is_col_major ? i : NumDims - 1 - i]++;
+          break;
+        }
+        if (i != NumDims - 1) it[i].count = 0;
+        coords[is_col_major ? i : NumDims - 1 - i] = initial_coords[is_col_major ? i : NumDims - 1 - i];
+        offset -= it[i].span;
+      }
+    }
+
+    return block_storage.AsTensorMaterializedBlock();
+  }
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorOpCost costPerCoeff(bool vectorized) const {
+    double compute_cost = NumDims * (2 * TensorOpCost::AddCost<Index>() + 2 * TensorOpCost::MulCost<Index>() +
+                                     TensorOpCost::DivCost<Index>());
+    for (int i = 0; i < NumDims; ++i) {
+      compute_cost += 2 * TensorOpCost::AddCost<Index>();
+    }
+    return m_impl.costPerCoeff(vectorized) + TensorOpCost(0, 0, compute_cost, false /* vectorized */, PacketSize);
+  }
+
+  EIGEN_DEVICE_FUNC typename Storage::Type data() const { return nullptr; }
+
+ protected:
+  Dimensions m_dimensions;
+  array<Index, NumDims> m_strides;
+  array<IndexDivisor, NumDims> m_fast_strides;
+  TensorEvaluator<ArgType, Device> m_impl;
+  RollDimensions m_rolls;
+  const Device EIGEN_DEVICE_REF m_device;
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void extract_coordinates(Index index, array<Index, NumDims>& coords) const {
+    if (static_cast<int>(Layout) == static_cast<int>(ColMajor)) {
+      for (int i = NumDims - 1; i > 0; --i) {
+        const Index idx = index / m_fast_strides[i];
+        index -= idx * m_strides[i];
+        coords[i] = idx;
+      }
+      coords[0] = index;
+    } else {
+      for (int i = 0; i < NumDims - 1; ++i) {
+        const Index idx = index / m_fast_strides[i];
+        index -= idx * m_strides[i];
+        coords[i] = idx;
+      }
+      coords[NumDims - 1] = index;
+    }
+  }
+
+ private:
+};
+
+// Eval as lvalue
+
+template <typename RollDimensions, typename ArgType, typename Device>
+struct TensorEvaluator<TensorRollOp<RollDimensions, ArgType>, Device>
+    : public TensorEvaluator<const TensorRollOp<RollDimensions, ArgType>, Device> {
+  typedef TensorEvaluator<const TensorRollOp<RollDimensions, ArgType>, Device> Base;
+  typedef TensorRollOp<RollDimensions, ArgType> XprType;
+  typedef typename XprType::Index Index;
+  static constexpr int NumDims = internal::array_size<RollDimensions>::value;
+  typedef DSizes<Index, NumDims> Dimensions;
+
+  static constexpr int Layout = TensorEvaluator<ArgType, Device>::Layout;
+  enum {
+    IsAligned = false,
+    PacketAccess = TensorEvaluator<ArgType, Device>::PacketAccess,
+    BlockAccess = false,
+    PreferBlockAccess = false,
+    CoordAccess = false,
+    RawAccess = false
+  };
+  EIGEN_STRONG_INLINE TensorEvaluator(const XprType& op, const Device& device) : Base(op, device) {}
+
+  typedef typename XprType::Scalar Scalar;
+  typedef typename XprType::CoeffReturnType CoeffReturnType;
+  typedef typename PacketType<CoeffReturnType, Device>::type PacketReturnType;
+  static constexpr int PacketSize = PacketType<CoeffReturnType, Device>::size;
+
+  //===- Tensor block evaluation strategy (see TensorBlock.h) -------------===//
+  typedef internal::TensorBlockNotImplemented TensorBlock;
+  //===--------------------------------------------------------------------===//
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Dimensions& dimensions() const { return this->m_dimensions; }
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar& coeffRef(Index index) const {
+    return this->m_impl.coeffRef(this->rollIndex(index));
+  }
+
+  template <int StoreMode>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void writePacket(Index index, const PacketReturnType& x) const {
+    eigen_assert(index + PacketSize - 1 < dimensions().TotalSize());
+    EIGEN_ALIGN_MAX CoeffReturnType values[PacketSize];
+    internal::pstore<CoeffReturnType, PacketReturnType>(values, x);
+    EIGEN_UNROLL_LOOP
+    for (int i = 0; i < PacketSize; ++i) {
+      this->coeffRef(index + i) = values[i];
+    }
+  }
+};
+
+}  // end namespace Eigen
+
+#endif  // EIGEN_CXX11_TENSOR_TENSOR_ROLL_H

diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
index 5169f1c..0b266f9 100644
--- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
+++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h

@@ -345,11 +345,17 @@
 /***************************************************************************
  * Implementation of erfc, requires C++11/C99                               *
  ****************************************************************************/
+template <typename Scalar>
+struct generic_fast_erfc {
+  template <typename T>
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T run(const T& x_in);
+};
+
+template <>
 template <typename T>
-EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_fast_erfc_float(const T& x) {
-  const T x_abs = pmin(pabs(x), pset1<T>(10.0f));
-  const T one = pset1<T>(1.0f);
-  const T x_abs_gt_one_mask = pcmp_lt(one, x_abs);
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_fast_erfc<float>::run(const T& x_in) {
+  constexpr float kClamp = 11.0f;
+  const T x = pmin(pmax(x_in, pset1<T>(-kClamp)), pset1<T>(kClamp));
 
   // erfc(x) = 1 + x * S(x^2), |x| <= 1.
   //
@@ -360,10 +366,12 @@
                              2.67075151205062866210937500000e-02, -1.12800106406211853027343750000e-01,
                              3.76122951507568359375000000000e-01, -1.12837910652160644531250000000e+00};
   const T x2 = pmul(x, x);
+  const T one = pset1<T>(1.0);
   const T erfc_small = pmadd(x, ppolevl<T, 5>::run(x2, alpha), one);
 
   // Return early if we don't need the more expensive approximation for any
   // entry in a.
+  const T x_abs_gt_one_mask = pcmp_lt(one, x2);
   if (!predux_any(x_abs_gt_one_mask)) return erfc_small;
 
   // erfc(x) = exp(-x^2) * 1/x * P(1/x^2) / Q(1/x^2), 1 < x < 9.
@@ -377,23 +385,103 @@
   constexpr float delta[] = {1.7251677811145782470703125e-02f, 3.9137163758277893066406250e-01f,
                              1.0000000000000000000000000e+00f, 6.2173241376876831054687500e-01f,
                              9.5662862062454223632812500e-02f};
-  const T z = pexp(pnegate(x2));
+  const T x2_lo = twoprod_low(x, x, x2);
+  // Here we use that
+  //   exp(-x^2) = exp(-(x2+x2_lo)^2) ~= exp(-x2)*exp(-x2_lo) ~= exp(-x2)*(1-x2_lo)
+  // since x2_lo < kClamp * eps << 1 in the region we care about. This trick reduces the max error
+  // from 34 ulps to below 5 ulps.
+  const T exp2_hi = pexp(pnegate(x2));
+  const T z = pnmadd(exp2_hi, x2_lo, exp2_hi);
   const T q2 = preciprocal(x2);
   const T num = ppolevl<T, 3>::run(q2, gamma);
-  const T denom = pmul(x_abs, ppolevl<T, 4>::run(q2, delta));
+  const T denom = pmul(x, ppolevl<T, 4>::run(q2, delta));
   const T r = pdiv(num, denom);
-  // If x < -1 then use erfc(x) = 2 - erfc(|x|).
-  const T x_negative = pcmp_lt(x, pset1<T>(0.0f));
-  const T erfc_large = pselect(x_negative, pnmadd(z, r, pset1<T>(2.0f)), pmul(z, r));
-
+  const T maybe_two = pand(pcmp_lt(x, pset1<T>(0.0)), pset1<T>(2.0));
+  const T erfc_large = pmadd(z, r, maybe_two);
   return pselect(x_abs_gt_one_mask, erfc_large, erfc_small);
 }
 
-template <typename Scalar>
-struct erfc_impl {
-  EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), THIS_TYPE_IS_NOT_SUPPORTED)
+template <>
+template <typename T>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_fast_erfc<double>::run(const T& x_in) {
+  // Clamp x to [-27:27] beyond which erfc(x) is either two or zero (below the underflow threshold).
+  // This avoids having to deal with twoprod(x,x) producing NaN for sufficiently large x.
+  constexpr double kClamp = 28.0;
+  const T x = pmin(pmax(x_in, pset1<T>(-kClamp)), pset1<T>(kClamp));
 
-  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Scalar run(const Scalar) { return Scalar(0); }
+  // erfc(x) = 1 + x * S(x^2) / T(x^2), |x| <= 1.
+  //
+  // Coefficients for S and T generated with Rminimax command:
+  //  ./ratapprox --function="erfc(x)-1" --dom='[-1,1]' --type=[9,10]
+  //  --num="odd" --numF="[D]" --den="even" --denF="[D]" --log --dispCoeff="dec"
+  constexpr double alpha[] = {-1.9493725660006057018823477644531294572516344487667083740234375e-04,
+                              -1.8272566210022942682217328425053892715368419885635375976562500e-03,
+                              -4.5303363351690106863856044583371840417385101318359375000000000e-02,
+                              -1.4215015503619179981775744181504705920815467834472656250000000e-01,
+                              -1.1283791670955125585606992899556644260883331298828125000000000e+00};
+  constexpr double beta[] = {2.0294484101083099089526257108317963684385176748037338256835938e-05,
+                             6.8117805899186819641732970609382391558028757572174072265625000e-04,
+                             1.0582026056098614921752165685120417037978768348693847656250000e-02,
+                             9.3252603143757495374188692949246615171432495117187500000000000e-02,
+                             4.5931062818368939559832142549566924571990966796875000000000000e-01,
+                             1.0};
+  const T x2 = pmul(x, x);
+  const T num_small = ppolevl<T, 4>::run(x2, alpha);
+  const T denom_small = ppolevl<T, 5>::run(x2, beta);
+  const T one = pset1<T>(1.0);
+  const T erfc_small = pmadd(x, pdiv(num_small, denom_small), one);
+
+  // Return early if we don't need the more expensive approximation for any
+  // entry in a.
+  const T x_abs_gt_one_mask = pcmp_lt(one, x2);
+  if (!predux_any(x_abs_gt_one_mask)) return erfc_small;
+
+  // erfc(x) = exp(-x^2) * 1/x * P(x) / Q(x), 1 < x < 27.
+  //
+  // Coefficients for P and Q generated with Rminimax command:
+  //  ./ratapprox --function="erfc(1/sqrt(x))*exp(1/x)/sqrt(x)"  --dom='[0.0013717,1]' --type=[9,9] --numF="[D]"
+  //  --denF="[D]" --log --dispCoeff="dec"
+  constexpr double gamma[] = {1.5252844933226974316088642158462107545346952974796295166015625e-04,
+                              1.0909912393738931124520519233556115068495273590087890625000000e-02,
+                              1.0628604636755033252537572252549580298364162445068359375000000e-01,
+                              3.3492472973137982217295416376146022230386734008789062500000000e-01,
+                              4.5065776215933289750026347064704168587923049926757812500000000e-01,
+                              2.9433039130294824659017649537418037652969360351562500000000000e-01,
+                              9.8792676360600226170838311645638896152377128601074218750000000e-02,
+                              1.7095935395503719655962981960328761488199234008789062500000000e-02,
+                              1.4249109729504577659398023570247460156679153442382812500000000e-03,
+                              4.4567378313647954771875570045835956989321857690811157226562500e-05};
+  constexpr double delta[] = {2.041985103115789845773520028160419315099716186523437500000000e-03,
+                              5.316030659946043707142493417450168635696172714233398437500000e-02,
+                              3.426242193784684864077405563875799998641014099121093750000000e-01,
+                              8.565637124308049799026321124983951449394226074218750000000000e-01,
+                              1.000000000000000000000000000000000000000000000000000000000000e+00,
+                              5.968805280570776972126623149961233139038085937500000000000000e-01,
+                              1.890922854723317836356244470152887515723705291748046875000000e-01,
+                              3.152505418656005586885981983868987299501895904541015625000000e-02,
+                              2.565085751861882583380047861965067568235099315643310546875000e-03,
+                              7.899362131678837697403017248376499992446042597293853759765625e-05};
+
+  const T x2_lo = twoprod_low(x, x, x2);
+  // Here we use that
+  //   exp(-x^2) = exp(-(x2+x2_lo)^2) ~= exp(-x2)*exp(-x2_lo) ~= exp(-x2)*(1-x2_lo)
+  // since x2_lo < kClamp *eps << 1 in the region we care about. This trick reduces the max error
+  // from 258 ulps to below 7 ulps.
+  const T exp2_hi = pexp(pnegate(x2));
+  const T z = pnmadd(exp2_hi, x2_lo, exp2_hi);
+  const T q2 = preciprocal(x2);
+  const T num_large = ppolevl<T, 9>::run(q2, gamma);
+  const T denom_large = pmul(x, ppolevl<T, 9>::run(q2, delta));
+  const T r = pdiv(num_large, denom_large);
+  const T maybe_two = pand(pcmp_lt(x, pset1<T>(0.0)), pset1<T>(2.0));
+  const T erfc_large = pmadd(z, r, maybe_two);
+  return pselect(x_abs_gt_one_mask, erfc_large, erfc_small);
+}
+
+template <typename T>
+struct erfc_impl {
+  typedef typename unpacket_traits<T>::type Scalar;
+  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { return generic_fast_erfc<Scalar>::run(x); }
 };
 
 template <typename Scalar>
@@ -408,7 +496,7 @@
 #if defined(SYCL_DEVICE_ONLY)
     return cl::sycl::erfc(x);
 #else
-    return generic_fast_erfc_float(x);
+    return generic_fast_erfc<float>::run(x);
 #endif
   }
 };
@@ -419,7 +507,7 @@
 #if defined(SYCL_DEVICE_ONLY)
     return cl::sycl::erfc(x);
 #else
-    return ::erfc(x);
+    return generic_fast_erfc<double>::run(x);
 #endif
   }
 };

diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index 8af6130..c270458 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt

@@ -198,8 +198,10 @@
 ei_add_test(cxx11_tensor_padding)
 ei_add_test(cxx11_tensor_patch)
 ei_add_test(cxx11_tensor_random)
+ei_add_test(cxx11_tensor_reverse)
 ei_add_test(cxx11_tensor_reduction)
 ei_add_test(cxx11_tensor_ref)
+ei_add_test(cxx11_tensor_roll)
 ei_add_test(cxx11_tensor_roundings)
 ei_add_test(cxx11_tensor_scan)
 ei_add_test(cxx11_tensor_shuffling)

diff --git a/unsupported/test/cxx11_tensor_roll.cpp b/unsupported/test/cxx11_tensor_roll.cpp
new file mode 100644
index 0000000..59f5efe
--- /dev/null
+++ b/unsupported/test/cxx11_tensor_roll.cpp

@@ -0,0 +1,156 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2024 Tobias Wood tobias@spinicist.org.uk
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "main.h"
+
+#include <Eigen/CXX11/Tensor>
+
+using Eigen::array;
+using Eigen::Tensor;
+
+template <int DataLayout>
+static void test_simple_roll() {
+  Tensor<float, 4, DataLayout> tensor(2, 3, 5, 7);
+  tensor.setRandom();
+
+  array<Index, 4> dim_roll;
+  dim_roll[0] = 0;
+  dim_roll[1] = 1;
+  dim_roll[2] = 4;
+  dim_roll[3] = 8;
+
+  Tensor<float, 4, DataLayout> rolled_tensor;
+  rolled_tensor = tensor.roll(dim_roll);
+
+  VERIFY_IS_EQUAL(rolled_tensor.dimension(0), 2);
+  VERIFY_IS_EQUAL(rolled_tensor.dimension(1), 3);
+  VERIFY_IS_EQUAL(rolled_tensor.dimension(2), 5);
+  VERIFY_IS_EQUAL(rolled_tensor.dimension(3), 7);
+
+  for (int i = 0; i < 2; ++i) {
+    for (int j = 0; j < 3; ++j) {
+      for (int k = 0; k < 5; ++k) {
+        for (int l = 0; l < 7; ++l) {
+          VERIFY_IS_EQUAL(tensor(i, (j + 1) % 3, (k + 4) % 5, (l + 8) % 7), rolled_tensor(i, j, k, l));
+        }
+      }
+    }
+  }
+
+  dim_roll[0] = -3;
+  dim_roll[1] = -2;
+  dim_roll[2] = -1;
+  dim_roll[3] = 0;
+
+  rolled_tensor = tensor.roll(dim_roll);
+
+  VERIFY_IS_EQUAL(rolled_tensor.dimension(0), 2);
+  VERIFY_IS_EQUAL(rolled_tensor.dimension(1), 3);
+  VERIFY_IS_EQUAL(rolled_tensor.dimension(2), 5);
+  VERIFY_IS_EQUAL(rolled_tensor.dimension(3), 7);
+
+  for (int i = 0; i < 2; ++i) {
+    for (int j = 0; j < 3; ++j) {
+      for (int k = 0; k < 5; ++k) {
+        for (int l = 0; l < 7; ++l) {
+          VERIFY_IS_EQUAL(tensor((i + 1) % 2, (j + 1) % 3, (k + 4) % 5, l), rolled_tensor(i, j, k, l));
+        }
+      }
+    }
+  }
+}
+
+template <int DataLayout>
+static void test_expr_roll(bool LValue) {
+  Tensor<float, 4, DataLayout> tensor(2, 3, 5, 7);
+  tensor.setRandom();
+
+  array<bool, 4> dim_roll;
+  dim_roll[0] = 2;
+  dim_roll[1] = 1;
+  dim_roll[2] = 0;
+  dim_roll[3] = 3;
+
+  Tensor<float, 4, DataLayout> expected(tensor.dimensions());
+  if (LValue) {
+    expected.roll(dim_roll) = tensor;
+  } else {
+    expected = tensor.roll(dim_roll);
+  }
+
+  Tensor<float, 4, DataLayout> result(tensor.dimensions());
+
+  array<ptrdiff_t, 4> src_slice_dim;
+  src_slice_dim[0] = tensor.dimension(0);
+  src_slice_dim[1] = tensor.dimension(1);
+  src_slice_dim[2] = 1;
+  src_slice_dim[3] = tensor.dimension(3);
+  array<ptrdiff_t, 4> src_slice_start;
+  src_slice_start[0] = 0;
+  src_slice_start[1] = 0;
+  src_slice_start[2] = 0;
+  src_slice_start[3] = 0;
+  array<ptrdiff_t, 4> dst_slice_dim = src_slice_dim;
+  array<ptrdiff_t, 4> dst_slice_start = src_slice_start;
+
+  for (int i = 0; i < tensor.dimension(2); ++i) {
+    if (LValue) {
+      result.slice(dst_slice_start, dst_slice_dim).roll(dim_roll) = tensor.slice(src_slice_start, src_slice_dim);
+    } else {
+      result.slice(dst_slice_start, dst_slice_dim) = tensor.slice(src_slice_start, src_slice_dim).roll(dim_roll);
+    }
+    src_slice_start[2] += 1;
+    dst_slice_start[2] += 1;
+  }
+
+  VERIFY_IS_EQUAL(result.dimension(0), tensor.dimension(0));
+  VERIFY_IS_EQUAL(result.dimension(1), tensor.dimension(1));
+  VERIFY_IS_EQUAL(result.dimension(2), tensor.dimension(2));
+  VERIFY_IS_EQUAL(result.dimension(3), tensor.dimension(3));
+
+  for (int i = 0; i < expected.dimension(0); ++i) {
+    for (int j = 0; j < expected.dimension(1); ++j) {
+      for (int k = 0; k < expected.dimension(2); ++k) {
+        for (int l = 0; l < expected.dimension(3); ++l) {
+          VERIFY_IS_EQUAL(result(i, j, k, l), expected(i, j, k, l));
+        }
+      }
+    }
+  }
+
+  dst_slice_start[2] = 0;
+  result.setRandom();
+  for (int i = 0; i < tensor.dimension(2); ++i) {
+    if (LValue) {
+      result.slice(dst_slice_start, dst_slice_dim).roll(dim_roll) = tensor.slice(dst_slice_start, dst_slice_dim);
+    } else {
+      result.slice(dst_slice_start, dst_slice_dim) = tensor.roll(dim_roll).slice(dst_slice_start, dst_slice_dim);
+    }
+    dst_slice_start[2] += 1;
+  }
+
+  for (int i = 0; i < expected.dimension(0); ++i) {
+    for (int j = 0; j < expected.dimension(1); ++j) {
+      for (int k = 0; k < expected.dimension(2); ++k) {
+        for (int l = 0; l < expected.dimension(3); ++l) {
+          VERIFY_IS_EQUAL(result(i, j, k, l), expected(i, j, k, l));
+        }
+      }
+    }
+  }
+}
+
+EIGEN_DECLARE_TEST(cxx11_tensor_roll) {
+  CALL_SUBTEST(test_simple_roll<ColMajor>());
+  CALL_SUBTEST(test_simple_roll<RowMajor>());
+  CALL_SUBTEST(test_expr_roll<ColMajor>(true));
+  CALL_SUBTEST(test_expr_roll<RowMajor>(true));
+  CALL_SUBTEST(test_expr_roll<ColMajor>(false));
+  CALL_SUBTEST(test_expr_roll<RowMajor>(false));
+}

diff --git a/unsupported/test/special_packetmath.cpp b/unsupported/test/special_packetmath.cpp
index 1ec5c3d..044ce67 100644
--- a/unsupported/test/special_packetmath.cpp
+++ b/unsupported/test/special_packetmath.cpp

@@ -117,6 +117,8 @@
 #if EIGEN_HAS_C99_MATH
   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLGamma, std::lgamma, internal::plgamma);
   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErf, std::erf, internal::perf);
+  // FIXME(rmlarsen): This test occasionally fails due to difference in tiny subnormal results
+  // near the underflow boundary. I am not sure which version is correct.
   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErfc, std::erfc, internal::perfc);
 #endif
 }