Added missing NEON pcasts, update packetmath tests.

The NEON `pcast` operators are all implemented and tested for existing
packets. This requires adding a `pcast(a,b,c,d,e,f,g,h)` for casting
between `int64_t` and `int8_t` in `GenericPacketMath.h`.

Removed incorrect `HasHalfPacket`  definition for NEON's
`Packet2l`/`Packet2ul`.

Adjustments were also made to the `packetmath` tests. These include
- minor bug fixes for cast tests (i.e. 4:1 casts, only casting for
  packets that are vectorizable)
- added 8:1 cast tests
- random number generation
  - original had uninteresting 0 to 0 casts for many casts between
    floating-point and integers, and exhibited signed overflow
    undefined behavior

Tested:
```
$ aarch64-linux-gnu-g++ -static -I./ '-DEIGEN_TEST_PART_ALL=1' test/packetmath.cpp -o packetmath
$ adb push packetmath /data/local/tmp/
$ adb shell "/data/local/tmp/packetmath"
```

4 files changed, 1884 insertions, 612 deletions

diff --git a/Eigen/src/Core/GenericPacketMath.h b/Eigen/src/Core/GenericPacketMath.h
index cb8848731..39bba577a 100644
--- a/Eigen/src/Core/GenericPacketMath.h
+++ b/Eigen/src/Core/GenericPacketMath.h
@@ -162,12 +162,17 @@ EIGEN_DEVICE_FUNC inline TgtPacket
 pcast(const SrcPacket& a, const SrcPacket& /*b*/) {
   return static_cast<TgtPacket>(a);
 }
-
 template <typename SrcPacket, typename TgtPacket>
 EIGEN_DEVICE_FUNC inline TgtPacket
 pcast(const SrcPacket& a, const SrcPacket& /*b*/, const SrcPacket& /*c*/, const SrcPacket& /*d*/) {
   return static_cast<TgtPacket>(a);
 }
+template <typename SrcPacket, typename TgtPacket>
+EIGEN_DEVICE_FUNC inline TgtPacket
+pcast(const SrcPacket& a, const SrcPacket& /*b*/, const SrcPacket& /*c*/, const SrcPacket& /*d*/,
+      const SrcPacket& /*e*/, const SrcPacket& /*f*/, const SrcPacket& /*g*/, const SrcPacket& /*h*/) {
+  return static_cast<TgtPacket>(a);
+}
 
 /** \internal \returns reinterpret_cast<Target>(a) */
 template <typename Target, typename Packet>
diff --git a/Eigen/src/Core/arch/NEON/PacketMath.h b/Eigen/src/Core/arch/NEON/PacketMath.h
index ddd43c377..340e1b1c5 100644
--- a/Eigen/src/Core/arch/NEON/PacketMath.h
+++ b/Eigen/src/Core/arch/NEON/PacketMath.h
@@ -338,7 +338,7 @@ struct packet_traits<int64_t> : default_packet_traits
     Vectorizable = 1,
     AlignedOnScalar = 1,
     size = 2,
-    HasHalfPacket = 1,
+    HasHalfPacket = 0,
 
     HasCmp       = 1,
     HasAdd       = 1,
@@ -368,7 +368,7 @@ struct packet_traits<uint64_t> : default_packet_traits
     Vectorizable = 1,
     AlignedOnScalar = 1,
     size = 2,
-    HasHalfPacket = 1,
+    HasHalfPacket = 0,
 
     HasCmp       = 1,
     HasAdd       = 1,
diff --git a/Eigen/src/Core/arch/NEON/TypeCasting.h b/Eigen/src/Core/arch/NEON/TypeCasting.h
index 68d24dc5c..80be213d2 100644
--- a/Eigen/src/Core/arch/NEON/TypeCasting.h
+++ b/Eigen/src/Core/arch/NEON/TypeCasting.h
@@ -2,6 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2018 Rasmus Munk Larsen <rmlarsen@google.com>
+// Copyright (C) 2020 Antonio Sanchez <cantonios@google.com>
 //
 // 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
@@ -14,265 +15,1397 @@ namespace Eigen {
 
 namespace internal {
 
-template<> struct type_casting_traits<float,numext::int32_t>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-template<> struct type_casting_traits<float,numext::uint32_t>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-template<> struct type_casting_traits<numext::int8_t,numext::uint8_t>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-template<> struct type_casting_traits<numext::uint8_t,numext::int8_t>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-template<> struct type_casting_traits<numext::int16_t,numext::uint16_t>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-template<> struct type_casting_traits<numext::uint16_t,numext::int16_t>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-template<> struct type_casting_traits<numext::int32_t,float>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-template<> struct type_casting_traits<numext::int32_t,numext::uint32_t>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-template<> struct type_casting_traits<numext::uint32_t,float>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-template<> struct type_casting_traits<numext::uint32_t,numext::int32_t>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-template<> struct type_casting_traits<numext::int64_t,numext::uint64_t>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-template<> struct type_casting_traits<numext::uint64_t,numext::int64_t>
-{ enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 }; };
-
-template<> EIGEN_STRONG_INLINE Packet2f pcast<Packet2i,Packet2f>(const Packet2i& a) { return vcvt_f32_s32(a); }
-template<> EIGEN_STRONG_INLINE Packet2f pcast<Packet2ui,Packet2f>(const Packet2ui& a) { return vcvt_f32_u32(a); }
-template<> EIGEN_STRONG_INLINE Packet2f pcast<Packet2l,Packet2f>(const Packet2l& a)
-{ return vcvt_f32_s32(vmovn_s64(a)); }
-template<> EIGEN_STRONG_INLINE Packet2f pcast<Packet2ul,Packet2f>(const Packet2ul& a)
-{ return vcvt_f32_u32(vmovn_u64(a)); }
-template<> EIGEN_STRONG_INLINE Packet4f pcast<Packet4c,Packet4f>(const Packet4c& a)
-{ return vcvtq_f32_s32(vmovl_s16(vget_low_s16(vmovl_s8(vreinterpret_s8_s32(vdup_n_s32(a)))))); }
-template<> EIGEN_STRONG_INLINE Packet4f pcast<Packet4uc,Packet4f>(const Packet4uc& a)
-{ return vcvtq_f32_s32(vreinterpretq_s32_u32(vmovl_u16(vget_low_u16(vmovl_u8(vreinterpret_u8_u32(vdup_n_u32(a))))))); }
-template<> EIGEN_STRONG_INLINE Packet4f pcast<Packet4s,Packet4f>(const Packet4s& a)
-{ return vcvtq_f32_s32(vmovl_s16(a)); }
-template<> EIGEN_STRONG_INLINE Packet4f pcast<Packet4us,Packet4f>(const Packet4us& a)
-{ return vcvtq_f32_s32(vreinterpretq_s32_u32(vmovl_u16(a))); }
-template<> EIGEN_STRONG_INLINE Packet4f pcast<Packet4i,Packet4f>(const Packet4i& a) { return vcvtq_f32_s32(a); }
-template<> EIGEN_STRONG_INLINE Packet4f pcast<Packet4ui,Packet4f>(const Packet4ui& a) { return vcvtq_f32_u32(a); }
-template<> EIGEN_STRONG_INLINE Packet4c pcast<Packet4f,Packet4c>(const Packet4f& a)
-{
-  const int16x4_t b = vmovn_s32(vcvtq_s32_f32(a));
-  return vget_lane_s32(vreinterpret_s32_s8(vmovn_s16(vcombine_s16(b, b))), 0);
-}
-template<> EIGEN_STRONG_INLINE Packet4c pcast<Packet4uc,Packet4c>(const Packet4uc& a)
-{ return static_cast<Packet4c>(a); }
-template<> EIGEN_STRONG_INLINE Packet4c pcast<Packet4s,Packet4c>(const Packet4s& a)
-{ return vget_lane_s32(vreinterpret_s32_s8(vmovn_s16(vcombine_s16(a, a))), 0); }
-template<> EIGEN_STRONG_INLINE Packet4c pcast<Packet4us,Packet4c>(const Packet4us& a)
-{
-  const int16x4_t b = vreinterpret_s16_u16(a);
-  return vget_lane_s32(vreinterpret_s32_s8(vmovn_s16(vcombine_s16(b, b))), 0);
-}
-template<> EIGEN_STRONG_INLINE Packet4c pcast<Packet4i,Packet4c>(const Packet4i& a)
-{
-  const int16x4_t b = vmovn_s32(a);
-  return vget_lane_s32(vreinterpret_s32_s8(vmovn_s16(vcombine_s16(b, b))), 0);
-}
-template<> EIGEN_STRONG_INLINE Packet4c pcast<Packet4ui,Packet4c>(const Packet4ui& a)
-{
-  const int16x4_t b = vmovn_s32(vreinterpretq_s32_u32(a));
-  return vget_lane_s32(vreinterpret_s32_s8(vmovn_s16(vcombine_s16(b, b))), 0);
-}
-template<> EIGEN_STRONG_INLINE Packet8c pcast<Packet8uc,Packet8c>(const Packet8uc& a) { return vreinterpret_s8_u8(a); }
-template<> EIGEN_STRONG_INLINE Packet8c pcast<Packet8s,Packet8c>(const Packet8s& a) { return vmovn_s16(a); }
-template<> EIGEN_STRONG_INLINE Packet8c pcast<Packet8us,Packet8c>(const Packet8us& a)
-{ return vreinterpret_s8_u8(vmovn_u16(a)); }
-template<> EIGEN_STRONG_INLINE Packet16c pcast<Packet16uc,Packet16c>(const Packet16uc& a)
-{ return vreinterpretq_s8_u8(a); }
-template<> EIGEN_STRONG_INLINE Packet4uc pcast<Packet4f,Packet4uc>(const Packet4f& a)
-{
-  const uint16x4_t b = vmovn_u32(vreinterpretq_u32_s32(vcvtq_s32_f32(a)));
-  return vget_lane_u32(vreinterpret_u32_u8(vmovn_u16(vcombine_u16(b, b))), 0);
-}
-template<> EIGEN_STRONG_INLINE Packet4uc pcast<Packet4i,Packet4uc>(const Packet4i& a)
-{
-  const uint16x4_t b = vmovn_u32(vreinterpretq_u32_s32(a));
-  return vget_lane_u32(vreinterpret_u32_u8(vmovn_u16(vcombine_u16(b, b))), 0);
-}
-template<> EIGEN_STRONG_INLINE Packet4uc pcast<Packet4ui,Packet4uc>(const Packet4ui& a)
-{
-  const uint16x4_t b = vmovn_u32(a);
-  return vget_lane_u32(vreinterpret_u32_u8(vmovn_u16(vcombine_u16(b, b))), 0);
-}
-template<> EIGEN_STRONG_INLINE Packet4uc pcast<Packet4c,Packet4uc>(const Packet4c& a)
-{ return static_cast<Packet4uc>(a); }
-template<> EIGEN_STRONG_INLINE Packet4uc pcast<Packet4s,Packet4uc>(const Packet4s& a)
-{
-  const uint16x4_t b = vreinterpret_u16_s16(a);
-  return vget_lane_u32(vreinterpret_u32_u8(vmovn_u16(vcombine_u16(b, b))), 0);
-}
-template<> EIGEN_STRONG_INLINE Packet4uc pcast<Packet4us,Packet4uc>(const Packet4us& a)
-{ return vget_lane_u32(vreinterpret_u32_u8(vmovn_u16(vcombine_u16(a, a))), 0); }
-template<> EIGEN_STRONG_INLINE Packet8uc pcast<Packet8c,Packet8uc>(const Packet8c& a) { return vreinterpret_u8_s8(a); }
-template<> EIGEN_STRONG_INLINE Packet8uc pcast<Packet8s,Packet8uc>(const Packet8s& a)
-{ return vreinterpret_u8_s8(vmovn_s16(a)); }
-template<> EIGEN_STRONG_INLINE Packet8uc pcast<Packet8us,Packet8uc>(const Packet8us& a) { return vmovn_u16(a); }
-template<> EIGEN_STRONG_INLINE Packet16uc pcast<Packet16c,Packet16uc>(const Packet16c& a)
-{ return vreinterpretq_u8_s8(a); }
-template<> EIGEN_STRONG_INLINE Packet4s pcast<Packet4f,Packet4s>(const Packet4f& a)
-{ return vmovn_s32(vcvtq_s32_f32(a)); }
-template<> EIGEN_STRONG_INLINE Packet4s pcast<Packet4c,Packet4s>(const Packet4c& a)
-{ return vget_low_s16(vmovl_s8(vreinterpret_s8_s32(vdup_n_s32(a)))); }
-template<> EIGEN_STRONG_INLINE Packet4s pcast<Packet4uc,Packet4s>(const Packet4uc& a)
-{ return vget_low_s16(vreinterpretq_s16_u16(vmovl_u8(vreinterpret_u8_u32(vdup_n_u32(a))))); }
-template<> EIGEN_STRONG_INLINE Packet4s pcast<Packet4us,Packet4s>(const Packet4us& a)
-{ return vreinterpret_s16_u16(a); }
-template<> EIGEN_STRONG_INLINE Packet4s pcast<Packet4i,Packet4s>(const Packet4i& a) { return vmovn_s32(a); }
-template<> EIGEN_STRONG_INLINE Packet4s pcast<Packet4ui,Packet4s>(const Packet4ui& a)
-{ return vmovn_s32(vreinterpretq_s32_u32(a)); }
-template<> EIGEN_STRONG_INLINE Packet8s pcast<Packet8uc,Packet8s>(const Packet8uc& a)
-{ return vreinterpretq_s16_u16(vmovl_u8(a)); }
-template<> EIGEN_STRONG_INLINE Packet8s pcast<Packet8c,Packet8s>(const Packet8c& a) { return vmovl_s8(a); }
-template<> EIGEN_STRONG_INLINE Packet8s pcast<Packet8us,Packet8s>(const Packet8us& a)
-{ return vreinterpretq_s16_u16(a); }
-template<> EIGEN_STRONG_INLINE Packet4us pcast<Packet4f,Packet4us>(const Packet4f& a)
-{ return vmovn_u32(vreinterpretq_u32_s32(vcvtq_s32_f32(a))); }
-template<> EIGEN_STRONG_INLINE Packet4us pcast<Packet4c,Packet4us>(const Packet4c& a)
-{ return vget_low_u16(vreinterpretq_u16_s16(vmovl_s8(vreinterpret_s8_s32(vdup_n_s32(a))))); }
-template<> EIGEN_STRONG_INLINE Packet4us pcast<Packet4uc,Packet4us>(const Packet4uc& a)
-{ return vget_low_u16(vmovl_u8(vreinterpret_u8_u32(vdup_n_u32(a)))); }
-template<> EIGEN_STRONG_INLINE Packet4us pcast<Packet4s,Packet4us>(const Packet4s& a)
-{ return vreinterpret_u16_s16(a); }
-template<> EIGEN_STRONG_INLINE Packet4us pcast<Packet4i,Packet4us>(const Packet4i& a)
-{ return vmovn_u32(vreinterpretq_u32_s32(a)); }
-template<> EIGEN_STRONG_INLINE Packet4us pcast<Packet4ui,Packet4us>(const Packet4ui& a) { return vmovn_u32(a); }
-template<> EIGEN_STRONG_INLINE Packet8us pcast<Packet8c,Packet8us>(const Packet8c& a)
-{ return vreinterpretq_u16_s16(vmovl_s8(a)); }
-template<> EIGEN_STRONG_INLINE Packet8us pcast<Packet8uc,Packet8us>(const Packet8uc& a) { return vmovl_u8(a); }
-template<> EIGEN_STRONG_INLINE Packet8us pcast<Packet8s,Packet8us>(const Packet8s& a)
-{ return vreinterpretq_u16_s16(a); }
-template<> EIGEN_STRONG_INLINE Packet2i pcast<Packet2f,Packet2i>(const Packet2f& a) { return vcvt_s32_f32(a); }
-template<> EIGEN_STRONG_INLINE Packet2i pcast<Packet2ui,Packet2i>(const Packet2ui& a)
-{ return vreinterpret_s32_u32(a); }
-template<> EIGEN_STRONG_INLINE Packet2i pcast<Packet2l,Packet2i>(const Packet2l& a)
-{ return vmovn_s64(a); }
-template<> EIGEN_STRONG_INLINE Packet2i pcast<Packet2ul,Packet2i>(const Packet2ul& a)
-{ return vmovn_s64(vreinterpretq_s64_u64(a)); }
-template<> EIGEN_STRONG_INLINE Packet4i pcast<Packet4f,Packet4i>(const Packet4f& a) { return vcvtq_s32_f32(a); }
-template<> EIGEN_STRONG_INLINE Packet4i pcast<Packet4c,Packet4i>(const Packet4c& a)
-{ return vmovl_s16(vget_low_s16(vmovl_s8(vreinterpret_s8_s32(vdup_n_s32(a))))); }
-template<> EIGEN_STRONG_INLINE Packet4i pcast<Packet4uc,Packet4i>(const Packet4uc& a)
-{ return vreinterpretq_s32_u32(vmovl_u16(vget_low_u16(vmovl_u8(vreinterpret_u8_u32(vdup_n_u32(a)))))); }
-template<> EIGEN_STRONG_INLINE Packet4i pcast<Packet4s,Packet4i>(const Packet4s& a) { return vmovl_s16(a); }
-template<> EIGEN_STRONG_INLINE Packet4i pcast<Packet4us,Packet4i>(const Packet4us& a)
-{ return vreinterpretq_s32_u32(vmovl_u16(a)); }
-template<> EIGEN_STRONG_INLINE Packet4i pcast<Packet4ui,Packet4i>(const Packet4ui& a)
-{ return vreinterpretq_s32_u32(a); }
-template<> EIGEN_STRONG_INLINE Packet2ui pcast<Packet2f,Packet2ui>(const Packet2f& a) { return vcvt_u32_f32(a); }
-template<> EIGEN_STRONG_INLINE Packet2ui pcast<Packet2i,Packet2ui>(const Packet2i& a)
-{ return vreinterpret_u32_s32(a); }
-template<> EIGEN_STRONG_INLINE Packet2ui pcast<Packet2l,Packet2ui>(const Packet2l& a)
-{ return vmovn_u64(vreinterpretq_u64_s64(a)); }
-template<> EIGEN_STRONG_INLINE Packet2ui pcast<Packet2ul,Packet2ui>(const Packet2ul& a)
-{ return vmovn_u64(a); }
-template<> EIGEN_STRONG_INLINE Packet4ui pcast<Packet4f,Packet4ui>(const Packet4f& a) { return vcvtq_u32_f32(a); }
-template<> EIGEN_STRONG_INLINE Packet4ui pcast<Packet4c,Packet4ui>(const Packet4c& a)
-{ return vreinterpretq_u32_s32(vmovl_s16(vget_low_s16(vmovl_s8(vreinterpret_s8_s32(vdup_n_s32(a)))))); }
-template<> EIGEN_STRONG_INLINE Packet4ui pcast<Packet4uc,Packet4ui>(const Packet4uc& a)
-{ return vmovl_u16(vget_low_u16(vmovl_u8(vreinterpret_u8_u32(vdup_n_u32(a))))); }
-template<> EIGEN_STRONG_INLINE Packet4ui pcast<Packet4s,Packet4ui>(const Packet4s& a)
-{ return vreinterpretq_u32_s32(vmovl_s16(a)); }
-template<> EIGEN_STRONG_INLINE Packet4ui pcast<Packet4us,Packet4ui>(const Packet4us& a) { return vmovl_u16(a); }
-template<> EIGEN_STRONG_INLINE Packet4ui pcast<Packet4i,Packet4ui>(const Packet4i& a)
-{ return vreinterpretq_u32_s32(a); }
-template<> EIGEN_STRONG_INLINE Packet2l pcast<Packet2f,Packet2l>(const Packet2f& a)
-{ return vmovl_s32(vcvt_s32_f32(a)); }
-template<> EIGEN_STRONG_INLINE Packet2l pcast<Packet2i,Packet2l>(const Packet2i& a)
-{ return vmovl_s32(a); }
-template<> EIGEN_STRONG_INLINE Packet2l pcast<Packet2ui,Packet2l>(const Packet2ui& a)
-{ return vreinterpretq_s64_u64(vmovl_u32(a)); }
-template<> EIGEN_STRONG_INLINE Packet2l pcast<Packet2ul,Packet2l>(const Packet2ul& a)
-{ return vreinterpretq_s64_u64(a); }
-template<> EIGEN_STRONG_INLINE Packet2ul pcast<Packet2f,Packet2ul>(const Packet2f& a)
-{ return vmovl_u32(vcvt_u32_f32(a)); }
-template<> EIGEN_STRONG_INLINE Packet2ul pcast<Packet2i,Packet2ul>(const Packet2i& a)
-{ return vreinterpretq_u64_s64(vmovl_s32(a)); }
-template<> EIGEN_STRONG_INLINE Packet2ul pcast<Packet2ui,Packet2ul>(const Packet2ui& a)
-{ return vmovl_u32(a); }
-template<> EIGEN_STRONG_INLINE Packet2ul pcast<Packet2l,Packet2ul>(const Packet2l& a)
-{ return vreinterpretq_u64_s64(a); }
-
-template<> EIGEN_STRONG_INLINE Packet2f preinterpret<Packet2f,Packet2i>(const Packet2i& a)
-{ return vreinterpret_f32_s32(a); }
-template<> EIGEN_STRONG_INLINE Packet2f preinterpret<Packet2f,Packet2ui>(const Packet2ui& a)
-{ return vreinterpret_f32_u32(a); }
-template<> EIGEN_STRONG_INLINE Packet4f preinterpret<Packet4f,Packet4i>(const Packet4i& a)
-{ return vreinterpretq_f32_s32(a); }
-template<> EIGEN_STRONG_INLINE Packet4f preinterpret<Packet4f,Packet4ui>(const Packet4ui& a)
-{ return vreinterpretq_f32_u32(a); }
-template<> EIGEN_STRONG_INLINE Packet4c preinterpret<Packet4c,Packet4uc>(const Packet4uc& a)
-{ return static_cast<Packet4c>(a); }
-template<> EIGEN_STRONG_INLINE Packet8c preinterpret<Packet8c,Packet8uc>(const Packet8uc& a)
-{ return vreinterpret_s8_u8(a); }
-template<> EIGEN_STRONG_INLINE Packet16c preinterpret<Packet16c,Packet16uc>(const Packet16uc& a)
-{ return vreinterpretq_s8_u8(a); }
-template<> EIGEN_STRONG_INLINE Packet4uc preinterpret<Packet4uc,Packet4c>(const Packet4c& a)
-{ return static_cast<Packet4uc>(a); }
-template<> EIGEN_STRONG_INLINE Packet8uc preinterpret<Packet8uc,Packet8c>(const Packet8c& a)
-{ return vreinterpret_u8_s8(a); }
-template<> EIGEN_STRONG_INLINE Packet16uc preinterpret<Packet16uc,Packet16c>(const Packet16c& a)
-{ return vreinterpretq_u8_s8(a); }
-template<> EIGEN_STRONG_INLINE Packet4s preinterpret<Packet4s,Packet4us>(const Packet4us& a)
-{ return vreinterpret_s16_u16(a); }
-template<> EIGEN_STRONG_INLINE Packet8s preinterpret<Packet8s,Packet8us>(const Packet8us& a)
-{ return vreinterpretq_s16_u16(a); }
-template<> EIGEN_STRONG_INLINE Packet4us preinterpret<Packet4us,Packet4s>(const Packet4s& a)
-{ return vreinterpret_u16_s16(a); }
-template<> EIGEN_STRONG_INLINE Packet8us preinterpret<Packet8us,Packet8s>(const Packet8s& a)
-{ return vreinterpretq_u16_s16(a); }
-template<> EIGEN_STRONG_INLINE Packet2i preinterpret<Packet2i,Packet2f>(const Packet2f& a)
-{ return vreinterpret_s32_f32(a); }
-template<> EIGEN_STRONG_INLINE Packet2i preinterpret<Packet2i,Packet2ui>(const Packet2ui& a)
-{ return vreinterpret_s32_u32(a); }
-template<> EIGEN_STRONG_INLINE Packet4i preinterpret<Packet4i,Packet4f>(const Packet4f& a)
-{ return vreinterpretq_s32_f32(a); }
-template<> EIGEN_STRONG_INLINE Packet4i preinterpret<Packet4i,Packet4ui>(const Packet4ui& a)
-{ return vreinterpretq_s32_u32(a); }
-template<> EIGEN_STRONG_INLINE Packet2ui preinterpret<Packet2ui,Packet2f>(const Packet2f& a)
-{ return vreinterpret_u32_f32(a); }
-template<> EIGEN_STRONG_INLINE Packet2ui preinterpret<Packet2ui,Packet2i>(const Packet2i& a)
-{ return vreinterpret_u32_s32(a); }
-template<> EIGEN_STRONG_INLINE Packet4ui preinterpret<Packet4ui,Packet4f>(const Packet4f& a)
-{ return vreinterpretq_u32_f32(a); }
-template<> EIGEN_STRONG_INLINE Packet4ui preinterpret<Packet4ui,Packet4i>(const Packet4i& a)
-{ return vreinterpretq_u32_s32(a); }
-template<> EIGEN_STRONG_INLINE Packet2l preinterpret<Packet2l,Packet2ul>(const Packet2ul& a)
-{ return vreinterpretq_s64_u64(a); }
-template<> EIGEN_STRONG_INLINE Packet2ul preinterpret<Packet2ul,Packet2l>(const Packet2l& a)
-{ return vreinterpretq_u64_s64(a); }
+//==============================================================================
+// pcast, SrcType = float
+//==============================================================================
+template <>
+struct type_casting_traits<float, float> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4f pcast<Packet4f, Packet4f>(const Packet4f& a) {
+  return a;
+}
+template <>
+EIGEN_STRONG_INLINE Packet2f pcast<Packet2f, Packet2f>(const Packet2f& a) {
+  return a;
+}
 
+template <>
+struct type_casting_traits<float, numext::int64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+struct type_casting_traits<float, numext::uint64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+// If float64 exists, first convert to that to keep as much precision as possible.
 #if EIGEN_ARCH_ARM64
+template <>
+EIGEN_STRONG_INLINE Packet2l pcast<Packet4f, Packet2l>(const Packet4f& a) {
+  // Discard second half of input.
+  return vcvtq_s64_f64(vcvt_f64_f32(vget_low_f32(a)));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ul pcast<Packet4f, Packet2ul>(const Packet4f& a) {
+  // Discard second half of input.
+  return vcvtq_u64_f64(vcvt_f64_f32(vget_low_f32(a)));
+}
+#else
+template <>
+EIGEN_STRONG_INLINE Packet2l pcast<Packet4f, Packet2l>(const Packet4f& a) {
+  // Discard second half of input.
+  return vmovl_s32(vget_low_s32(vcvtq_s32_f32(a)));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ul pcast<Packet4f, Packet2ul>(const Packet4f& a) {
+  // Discard second half of input.
+  return vmovl_u32(vget_low_u32(vcvtq_u32_f32(a)));
+}
+#endif  // EIGEN_ARCH_ARM64
+
+template <>
+struct type_casting_traits<float, numext::int32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4i pcast<Packet4f, Packet4i>(const Packet4f& a) {
+  return vcvtq_s32_f32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2i pcast<Packet2f, Packet2i>(const Packet2f& a) {
+  return vcvt_s32_f32(a);
+}
+
+template <>
+struct type_casting_traits<float, numext::uint32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4ui pcast<Packet4f, Packet4ui>(const Packet4f& a) {
+  return vcvtq_u32_f32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ui pcast<Packet2f, Packet2ui>(const Packet2f& a) {
+  return vcvt_u32_f32(a);
+}
+
+template <>
+struct type_casting_traits<float, numext::int16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8s pcast<Packet4f, Packet8s>(const Packet4f& a, const Packet4f& b) {
+  return vcombine_s16(vmovn_s32(vcvtq_s32_f32(a)), vmovn_s32(vcvtq_s32_f32(b)));
+}
+template <>
+EIGEN_STRONG_INLINE Packet4s pcast<Packet2f, Packet4s>(const Packet2f& a, const Packet2f& b) {
+  return vmovn_s32(vcombine_s32(vcvt_s32_f32(a), vcvt_s32_f32(b)));
+}
+
+template <>
+struct type_casting_traits<float, numext::uint16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8us pcast<Packet4f, Packet8us>(const Packet4f& a, const Packet4f& b) {
+  return vcombine_u16(vmovn_u32(vcvtq_u32_f32(a)), vmovn_u32(vcvtq_u32_f32(b)));
+}
+template <>
+EIGEN_STRONG_INLINE Packet4us pcast<Packet2f, Packet4us>(const Packet2f& a, const Packet2f& b) {
+  return vmovn_u32(vcombine_u32(vcvt_u32_f32(a), vcvt_u32_f32(b)));
+}
+
+template <>
+struct type_casting_traits<float, numext::int8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16c pcast<Packet4f, Packet16c>(const Packet4f& a, const Packet4f& b, const Packet4f& c,
+                                                         const Packet4f& d) {
+  const int16x8_t ab_s16 = pcast<Packet4f, Packet8s>(a, b);
+  const int16x8_t cd_s16 = pcast<Packet4f, Packet8s>(c, d);
+  return vcombine_s8(vmovn_s16(ab_s16), vmovn_s16(cd_s16));
+}
+template <>
+EIGEN_STRONG_INLINE Packet8c pcast<Packet2f, Packet8c>(const Packet2f& a, const Packet2f& b, const Packet2f& c,
+                                                       const Packet2f& d) {
+  const int16x4_t ab_s16 = pcast<Packet2f, Packet4s>(a, b);
+  const int16x4_t cd_s16 = pcast<Packet2f, Packet4s>(c, d);
+  return vmovn_s16(vcombine_s16(ab_s16, cd_s16));
+}
+
+template <>
+struct type_casting_traits<float, numext::uint8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16uc pcast<Packet4f, Packet16uc>(const Packet4f& a, const Packet4f& b, const Packet4f& c,
+                                                           const Packet4f& d) {
+  const uint16x8_t ab_u16 = pcast<Packet4f, Packet8us>(a, b);
+  const uint16x8_t cd_u16 = pcast<Packet4f, Packet8us>(c, d);
+  return vcombine_u8(vmovn_u16(ab_u16), vmovn_u16(cd_u16));
+}
+template <>
+EIGEN_STRONG_INLINE Packet8uc pcast<Packet2f, Packet8uc>(const Packet2f& a, const Packet2f& b, const Packet2f& c,
+                                                         const Packet2f& d) {
+  const uint16x4_t ab_u16 = pcast<Packet2f, Packet4us>(a, b);
+  const uint16x4_t cd_u16 = pcast<Packet2f, Packet4us>(c, d);
+  return vmovn_u16(vcombine_u16(ab_u16, cd_u16));
+}
+
+//==============================================================================
+// pcast, SrcType = int8_t
+//==============================================================================
+template <>
+struct type_casting_traits<numext::int8_t, float> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4f pcast<Packet16c, Packet4f>(const Packet16c& a) {
+  // Discard all but first 4 bytes.
+  return vcvtq_f32_s32(vmovl_s16(vget_low_s16(vmovl_s8(vget_low_s8(a)))));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2f pcast<Packet8c, Packet2f>(const Packet8c& a) {
+  // Discard all but first 2 bytes.
+  return vcvt_f32_s32(vget_low_s32(vmovl_s16(vget_low_s16(vmovl_s8(a)))));
+}
+
+template <>
+struct type_casting_traits<numext::int8_t, numext::int64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 8 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2l pcast<Packet16c, Packet2l>(const Packet16c& a) {
+  // Discard all but first two bytes.
+  return vmovl_s32(vget_low_s32(vmovl_s16(vget_low_s16(vmovl_s8(vget_low_s8(a))))));
+}
+
+template <>
+struct type_casting_traits<numext::int8_t, numext::uint64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 8 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2ul pcast<Packet16c, Packet2ul>(const Packet16c& a) {
+  return vreinterpretq_u64_s64(pcast<Packet16c, Packet2l>(a));
+}
+
+template <>
+struct type_casting_traits<numext::int8_t, numext::int32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4i pcast<Packet16c, Packet4i>(const Packet16c& a) {
+  // Discard all but first 4 bytes.
+  return vmovl_s16(vget_low_s16(vmovl_s8(vget_low_s8(a))));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2i pcast<Packet8c, Packet2i>(const Packet8c& a) {
+  // Discard all but first 2 bytes.
+  return vget_low_s32(vmovl_s16(vget_low_s16(vmovl_s8(a))));
+}
+
+template <>
+struct type_casting_traits<numext::int8_t, numext::uint32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4ui pcast<Packet16c, Packet4ui>(const Packet16c& a) {
+  return vreinterpretq_u32_s32(pcast<Packet16c, Packet4i>(a));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ui pcast<Packet8c, Packet2ui>(const Packet8c& a) {
+  return vreinterpret_u32_s32(pcast<Packet8c, Packet2i>(a));
+}
+
+template <>
+struct type_casting_traits<numext::int8_t, numext::int16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8s pcast<Packet16c, Packet8s>(const Packet16c& a) {
+  // Discard second half of input.
+  return vmovl_s8(vget_low_s8(a));
+}
+template <>
+EIGEN_STRONG_INLINE Packet4s pcast<Packet8c, Packet4s>(const Packet8c& a) {
+  // Discard second half of input.
+  return vget_low_s16(vmovl_s8(a));
+}
+
+template <>
+struct type_casting_traits<numext::int8_t, numext::uint16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8us pcast<Packet16c, Packet8us>(const Packet16c& a) {
+  return vreinterpretq_u16_s16(pcast<Packet16c, Packet8s>(a));
+}
+template <>
+EIGEN_STRONG_INLINE Packet4us pcast<Packet8c, Packet4us>(const Packet8c& a) {
+  return vreinterpret_u16_s16(pcast<Packet8c, Packet4s>(a));
+}
+
+template <>
+struct type_casting_traits<numext::int8_t, numext::int8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16c pcast<Packet16c, Packet16c>(const Packet16c& a) {
+  return a;
+}
+template <>
+EIGEN_STRONG_INLINE Packet8c pcast<Packet8c, Packet8c>(const Packet8c& a) {
+  return a;
+}
+template <>
+EIGEN_STRONG_INLINE Packet4c pcast<Packet4c, Packet4c>(const Packet4c& a) {
+  return a;
+}
+
+template <>
+struct type_casting_traits<numext::int8_t, numext::uint8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16uc pcast<Packet16c, Packet16uc>(const Packet16c& a) {
+  return vreinterpretq_u8_s8(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet8uc pcast<Packet8c, Packet8uc>(const Packet8c& a) {
+  return vreinterpret_u8_s8(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet4uc pcast<Packet4c, Packet4uc>(const Packet4c& a) {
+  return static_cast<Packet4uc>(a);
+}
+
+//==============================================================================
+// pcast, SrcType = uint8_t
+//==============================================================================
+template <>
+struct type_casting_traits<numext::uint8_t, float> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4f pcast<Packet16uc, Packet4f>(const Packet16uc& a) {
+  // Discard all but first 4 bytes.
+  return vcvtq_f32_u32(vmovl_u16(vget_low_u16(vmovl_u8(vget_low_u8(a)))));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2f pcast<Packet8uc, Packet2f>(const Packet8uc& a) {
+  // Discard all but first 2 bytes.
+  return vcvt_f32_u32(vget_low_u32(vmovl_u16(vget_low_u16(vmovl_u8(a)))));
+}
+
+template <>
+struct type_casting_traits<numext::uint8_t, numext::uint64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 8 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2ul pcast<Packet16uc, Packet2ul>(const Packet16uc& a) {
+  // Discard all but first two bytes.
+  return vmovl_u32(vget_low_u32(vmovl_u16(vget_low_u16(vmovl_u8(vget_low_u8(a))))));
+}
+
+template <>
+struct type_casting_traits<numext::uint8_t, numext::int64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 8 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2l pcast<Packet16uc, Packet2l>(const Packet16uc& a) {
+  return vreinterpretq_s64_u64(pcast<Packet16uc, Packet2ul>(a));
+}
+
+template <>
+struct type_casting_traits<numext::uint8_t, numext::uint32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4ui pcast<Packet16uc, Packet4ui>(const Packet16uc& a) {
+  // Discard all but first 4 bytes.
+  return vmovl_u16(vget_low_u16(vmovl_u8(vget_low_u8(a))));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ui pcast<Packet8uc, Packet2ui>(const Packet8uc& a) {
+  // Discard all but first 2 bytes.
+  return vget_low_u32(vmovl_u16(vget_low_u16(vmovl_u8(a))));
+}
+
+template <>
+struct type_casting_traits<numext::uint8_t, numext::int32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4i pcast<Packet16uc, Packet4i>(const Packet16uc& a) {
+  return vreinterpretq_s32_u32(pcast<Packet16uc, Packet4ui>(a));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2i pcast<Packet8uc, Packet2i>(const Packet8uc& a) {
+  return vreinterpret_s32_u32(pcast<Packet8uc, Packet2ui>(a));
+}
+
+template <>
+struct type_casting_traits<numext::uint8_t, numext::uint16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8us pcast<Packet16uc, Packet8us>(const Packet16uc& a) {
+  // Discard second half of input.
+  return vmovl_u8(vget_low_u8(a));
+}
+template <>
+EIGEN_STRONG_INLINE Packet4us pcast<Packet8uc, Packet4us>(const Packet8uc& a) {
+  // Discard second half of input.
+  return vget_low_u16(vmovl_u8(a));
+}
+
+template <>
+struct type_casting_traits<numext::uint8_t, numext::int16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8s pcast<Packet16uc, Packet8s>(const Packet16uc& a) {
+  return vreinterpretq_s16_u16(pcast<Packet16uc, Packet8us>(a));
+}
+template <>
+EIGEN_STRONG_INLINE Packet4s pcast<Packet8uc, Packet4s>(const Packet8uc& a) {
+  return vreinterpret_s16_u16(pcast<Packet8uc, Packet4us>(a));
+}
+
+template <>
+struct type_casting_traits<numext::uint8_t, numext::uint8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16uc pcast<Packet16uc, Packet16uc>(const Packet16uc& a) {
+  return a;
+}
+template <>
+EIGEN_STRONG_INLINE Packet8uc pcast<Packet8uc, Packet8uc>(const Packet8uc& a) {
+  return a;
+}
+template <>
+EIGEN_STRONG_INLINE Packet4uc pcast<Packet4uc, Packet4uc>(const Packet4uc& a) {
+  return a;
+}
+
+template <>
+struct type_casting_traits<numext::uint8_t, numext::int8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16c pcast<Packet16uc, Packet16c>(const Packet16uc& a) {
+  return vreinterpretq_s8_u8(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet8c pcast<Packet8uc, Packet8c>(const Packet8uc& a) {
+  return vreinterpret_s8_u8(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet4c pcast<Packet4uc, Packet4c>(const Packet4uc& a) {
+  return static_cast<Packet4c>(a);
+}
+
+//==============================================================================
+// pcast, SrcType = int16_t
+//==============================================================================
+template <>
+struct type_casting_traits<numext::int16_t, float> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4f pcast<Packet8s, Packet4f>(const Packet8s& a) {
+  // Discard second half of input.
+  return vcvtq_f32_s32(vmovl_s16(vget_low_s16(a)));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2f pcast<Packet4s, Packet2f>(const Packet4s& a) {
+  // Discard second half of input.
+  return vcvt_f32_s32(vget_low_s32(vmovl_s16(a)));
+}
+
+template <>
+struct type_casting_traits<numext::int16_t, numext::int64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2l pcast<Packet8s, Packet2l>(const Packet8s& a) {
+  // Discard all but first two values.
+  return vmovl_s32(vget_low_s32(vmovl_s16(vget_low_s16(a))));
+}
+
+template <>
+struct type_casting_traits<numext::int16_t, numext::uint64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2ul pcast<Packet8s, Packet2ul>(const Packet8s& a) {
+  return vreinterpretq_u64_s64(pcast<Packet8s, Packet2l>(a));
+}
+
+template <>
+struct type_casting_traits<numext::int16_t, numext::int32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4i pcast<Packet8s, Packet4i>(const Packet8s& a) {
+  // Discard second half of input.
+  return vmovl_s16(vget_low_s16(a));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2i pcast<Packet4s, Packet2i>(const Packet4s& a) {
+  // Discard second half of input.
+  return vget_low_s32(vmovl_s16(a));
+}
+
+template <>
+struct type_casting_traits<numext::int16_t, numext::uint32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4ui pcast<Packet8s, Packet4ui>(const Packet8s& a) {
+  return vreinterpretq_u32_s32(pcast<Packet8s, Packet4i>(a));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ui pcast<Packet4s, Packet2ui>(const Packet4s& a) {
+  return vreinterpret_u32_s32(pcast<Packet4s, Packet2i>(a));
+}
+
+template <>
+struct type_casting_traits<numext::int16_t, numext::int16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8s pcast<Packet8s, Packet8s>(const Packet8s& a) {
+  return a;
+}
+template <>
+EIGEN_STRONG_INLINE Packet4s pcast<Packet4s, Packet4s>(const Packet4s& a) {
+  return a;
+}
+
+template <>
+struct type_casting_traits<numext::int16_t, numext::uint16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8us pcast<Packet8s, Packet8us>(const Packet8s& a) {
+  return vreinterpretq_u16_s16(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet4us pcast<Packet4s, Packet4us>(const Packet4s& a) {
+  return vreinterpret_u16_s16(a);
+}
+
+template <>
+struct type_casting_traits<numext::int16_t, numext::int8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16c pcast<Packet8s, Packet16c>(const Packet8s& a, const Packet8s& b) {
+  return vcombine_s8(vmovn_s16(a), vmovn_s16(b));
+}
+template <>
+EIGEN_STRONG_INLINE Packet8c pcast<Packet4s, Packet8c>(const Packet4s& a, const Packet4s& b) {
+  return vmovn_s16(vcombine_s16(a, b));
+}
+
+template <>
+struct type_casting_traits<numext::int16_t, numext::uint8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16uc pcast<Packet8s, Packet16uc>(const Packet8s& a, const Packet8s& b) {
+  return vcombine_u8(vmovn_u16(vreinterpretq_u16_s16(a)), vmovn_u16(vreinterpretq_u16_s16(b)));
+}
+template <>
+EIGEN_STRONG_INLINE Packet8uc pcast<Packet4s, Packet8uc>(const Packet4s& a, const Packet4s& b) {
+  return vmovn_u16(vcombine_u16(vreinterpret_u16_s16(a), vreinterpret_u16_s16(b)));
+}
+
+//==============================================================================
+// pcast, SrcType = uint16_t
+//==============================================================================
+template <>
+struct type_casting_traits<numext::uint16_t, float> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4f pcast<Packet8us, Packet4f>(const Packet8us& a) {
+  // Discard second half of input.
+  return vcvtq_f32_u32(vmovl_u16(vget_low_u16(a)));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2f pcast<Packet4us, Packet2f>(const Packet4us& a) {
+  // Discard second half of input.
+  return vcvt_f32_u32(vget_low_u32(vmovl_u16(a)));
+}
+
+template <>
+struct type_casting_traits<numext::uint16_t, numext::uint64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2ul pcast<Packet8us, Packet2ul>(const Packet8us& a) {
+  // Discard all but first two values.
+  return vmovl_u32(vget_low_u32(vmovl_u16(vget_low_u16(a))));
+}
+
+template <>
+struct type_casting_traits<numext::uint16_t, numext::int64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2l pcast<Packet8us, Packet2l>(const Packet8us& a) {
+  return vreinterpretq_s64_u64(pcast<Packet8us, Packet2ul>(a));
+}
+
+template <>
+struct type_casting_traits<numext::uint16_t, numext::uint32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4ui pcast<Packet8us, Packet4ui>(const Packet8us& a) {
+  // Discard second half of input.
+  return vmovl_u16(vget_low_u16(a));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ui pcast<Packet4us, Packet2ui>(const Packet4us& a) {
+  // Discard second half of input.
+  return vget_low_u32(vmovl_u16(a));
+}
+
+template <>
+struct type_casting_traits<numext::uint16_t, numext::int32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4i pcast<Packet8us, Packet4i>(const Packet8us& a) {
+  return vreinterpretq_s32_u32(pcast<Packet8us, Packet4ui>(a));
+}
+template <>
+EIGEN_STRONG_INLINE Packet2i pcast<Packet4us, Packet2i>(const Packet4us& a) {
+  return vreinterpret_s32_u32(pcast<Packet4us, Packet2ui>(a));
+}
+
+template <>
+struct type_casting_traits<numext::uint16_t, numext::uint16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8us pcast<Packet8us, Packet8us>(const Packet8us& a) {
+  return a;
+}
+template <>
+EIGEN_STRONG_INLINE Packet4us pcast<Packet4us, Packet4us>(const Packet4us& a) {
+  return a;
+}
+
+template <>
+struct type_casting_traits<numext::uint16_t, numext::int16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8s pcast<Packet8us, Packet8s>(const Packet8us& a) {
+  return vreinterpretq_s16_u16(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet4s pcast<Packet4us, Packet4s>(const Packet4us& a) {
+  return vreinterpret_s16_u16(a);
+}
+
+template <>
+struct type_casting_traits<numext::uint16_t, numext::uint8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16uc pcast<Packet8us, Packet16uc>(const Packet8us& a, const Packet8us& b) {
+  return vcombine_u8(vmovn_u16(a), vmovn_u16(b));
+}
+template <>
+EIGEN_STRONG_INLINE Packet8uc pcast<Packet4us, Packet8uc>(const Packet4us& a, const Packet4us& b) {
+  return vmovn_u16(vcombine_u16(a, b));
+}
+
+template <>
+struct type_casting_traits<numext::uint16_t, numext::int8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16c pcast<Packet8us, Packet16c>(const Packet8us& a, const Packet8us& b) {
+  return vreinterpretq_s8_u8(pcast<Packet8us, Packet16uc>(a, b));
+}
+template <>
+EIGEN_STRONG_INLINE Packet8c pcast<Packet4us, Packet8c>(const Packet4us& a, const Packet4us& b) {
+  return vreinterpret_s8_u8(pcast<Packet4us, Packet8uc>(a, b));
+}
+
+//==============================================================================
+// pcast, SrcType = int32_t
+//==============================================================================
+template <>
+struct type_casting_traits<numext::int32_t, float> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4f pcast<Packet4i, Packet4f>(const Packet4i& a) {
+  return vcvtq_f32_s32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2f pcast<Packet2i, Packet2f>(const Packet2i& a) {
+  return vcvt_f32_s32(a);
+}
+
+template <>
+struct type_casting_traits<numext::int32_t, numext::int64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2l pcast<Packet4i, Packet2l>(const Packet4i& a) {
+  // Discard second half of input.
+  return vmovl_s32(vget_low_s32(a));
+}
+
+template <>
+struct type_casting_traits<numext::int32_t, numext::uint64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2ul pcast<Packet4i, Packet2ul>(const Packet4i& a) {
+  return vreinterpretq_u64_s64(pcast<Packet4i, Packet2l>(a));
+}
+
+template <>
+struct type_casting_traits<numext::int32_t, numext::int32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4i pcast<Packet4i, Packet4i>(const Packet4i& a) {
+  return a;
+}
+template <>
+EIGEN_STRONG_INLINE Packet2i pcast<Packet2i, Packet2i>(const Packet2i& a) {
+  return a;
+}
+
+template <>
+struct type_casting_traits<numext::int32_t, numext::uint32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4ui pcast<Packet4i, Packet4ui>(const Packet4i& a) {
+  return vreinterpretq_u32_s32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ui pcast<Packet2i, Packet2ui>(const Packet2i& a) {
+  return vreinterpret_u32_s32(a);
+}
+
+template <>
+struct type_casting_traits<numext::int32_t, numext::int16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8s pcast<Packet4i, Packet8s>(const Packet4i& a, const Packet4i& b) {
+  return vcombine_s16(vmovn_s32(a), vmovn_s32(b));
+}
+template <>
+EIGEN_STRONG_INLINE Packet4s pcast<Packet2i, Packet4s>(const Packet2i& a, const Packet2i& b) {
+  return vmovn_s32(vcombine_s32(a, b));
+}
+
+template <>
+struct type_casting_traits<numext::int32_t, numext::uint16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8us pcast<Packet4i, Packet8us>(const Packet4i& a, const Packet4i& b) {
+  return vcombine_u16(vmovn_u32(vreinterpretq_u32_s32(a)), vmovn_u32(vreinterpretq_u32_s32(b)));
+}
+template <>
+EIGEN_STRONG_INLINE Packet4us pcast<Packet2i, Packet4us>(const Packet2i& a, const Packet2i& b) {
+  return vmovn_u32(vreinterpretq_u32_s32(vcombine_s32(a, b)));
+}
+
+template <>
+struct type_casting_traits<numext::int32_t, numext::int8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16c pcast<Packet4i, Packet16c>(const Packet4i& a, const Packet4i& b, const Packet4i& c,
+                                                         const Packet4i& d) {
+  const int16x8_t ab_s16 = pcast<Packet4i, Packet8s>(a, b);
+  const int16x8_t cd_s16 = pcast<Packet4i, Packet8s>(c, d);
+  return vcombine_s8(vmovn_s16(ab_s16), vmovn_s16(cd_s16));
+}
+template <>
+EIGEN_STRONG_INLINE Packet8c pcast<Packet2i, Packet8c>(const Packet2i& a, const Packet2i& b, const Packet2i& c,
+                                                       const Packet2i& d) {
+  const int16x4_t ab_s16 = vmovn_s32(vcombine_s32(a, b));
+  const int16x4_t cd_s16 = vmovn_s32(vcombine_s32(c, d));
+  return vmovn_s16(vcombine_s16(ab_s16, cd_s16));
+}
+
+template <>
+struct type_casting_traits<numext::int32_t, numext::uint8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16uc pcast<Packet4i, Packet16uc>(const Packet4i& a, const Packet4i& b, const Packet4i& c,
+                                                           const Packet4i& d) {
+  const uint16x8_t ab_u16 = pcast<Packet4i, Packet8us>(a, b);
+  const uint16x8_t cd_u16 = pcast<Packet4i, Packet8us>(c, d);
+  return vcombine_u8(vmovn_u16(ab_u16), vmovn_u16(cd_u16));
+}
+template <>
+EIGEN_STRONG_INLINE Packet8uc pcast<Packet2i, Packet8uc>(const Packet2i& a, const Packet2i& b, const Packet2i& c,
+                                                         const Packet2i& d) {
+  const uint16x4_t ab_u16 = pcast<Packet2i, Packet4us>(a, b);
+  const uint16x4_t cd_u16 = pcast<Packet2i, Packet4us>(c, d);
+  return vmovn_u16(vcombine_u16(ab_u16, cd_u16));
+}
+
+//==============================================================================
+// pcast, SrcType = uint32_t
+//==============================================================================
+template <>
+struct type_casting_traits<numext::uint32_t, float> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4f pcast<Packet4ui, Packet4f>(const Packet4ui& a) {
+  return vcvtq_f32_u32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2f pcast<Packet2ui, Packet2f>(const Packet2ui& a) {
+  return vcvt_f32_u32(a);
+}
+
+template <>
+struct type_casting_traits<numext::uint32_t, numext::uint64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2ul pcast<Packet4ui, Packet2ul>(const Packet4ui& a) {
+  // Discard second half of input.
+  return vmovl_u32(vget_low_u32(a));
+}
+
+template <>
+struct type_casting_traits<numext::uint32_t, numext::int64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2l pcast<Packet4ui, Packet2l>(const Packet4ui& a) {
+  return vreinterpretq_s64_u64(pcast<Packet4ui, Packet2ul>(a));
+}
+
+template <>
+struct type_casting_traits<numext::uint32_t, numext::uint32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4ui pcast<Packet4ui, Packet4ui>(const Packet4ui& a) {
+  return a;
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ui pcast<Packet2ui, Packet2ui>(const Packet2ui& a) {
+  return a;
+}
+
+template <>
+struct type_casting_traits<numext::uint32_t, numext::int32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4i pcast<Packet4ui, Packet4i>(const Packet4ui& a) {
+  return vreinterpretq_s32_u32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2i pcast<Packet2ui, Packet2i>(const Packet2ui& a) {
+  return vreinterpret_s32_u32(a);
+}
+
+template <>
+struct type_casting_traits<numext::uint32_t, numext::uint16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8us pcast<Packet4ui, Packet8us>(const Packet4ui& a, const Packet4ui& b) {
+  return vcombine_u16(vmovn_u32(a), vmovn_u32(b));
+}
+template <>
+EIGEN_STRONG_INLINE Packet4us pcast<Packet2ui, Packet4us>(const Packet2ui& a, const Packet2ui& b) {
+  return vmovn_u32(vcombine_u32(a, b));
+}
+
+template <>
+struct type_casting_traits<numext::uint32_t, numext::int16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8s pcast<Packet4ui, Packet8s>(const Packet4ui& a, const Packet4ui& b) {
+  return vreinterpretq_s16_u16(pcast<Packet4ui, Packet8us>(a, b));
+}
+template <>
+EIGEN_STRONG_INLINE Packet4s pcast<Packet2ui, Packet4s>(const Packet2ui& a, const Packet2ui& b) {
+  return vreinterpret_s16_u16(pcast<Packet2ui, Packet4us>(a, b));
+}
+
+template <>
+struct type_casting_traits<numext::uint32_t, numext::uint8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16uc pcast<Packet4ui, Packet16uc>(const Packet4ui& a, const Packet4ui& b, const Packet4ui& c,
+                                                            const Packet4ui& d) {
+  const uint16x8_t ab_u16 = vcombine_u16(vmovn_u32(a), vmovn_u32(b));
+  const uint16x8_t cd_u16 = vcombine_u16(vmovn_u32(c), vmovn_u32(d));
+  return vcombine_u8(vmovn_u16(ab_u16), vmovn_u16(cd_u16));
+}
+template <>
+EIGEN_STRONG_INLINE Packet8uc pcast<Packet2ui, Packet8uc>(const Packet2ui& a, const Packet2ui& b, const Packet2ui& c,
+                                                          const Packet2ui& d) {
+  const uint16x4_t ab_u16 = vmovn_u32(vcombine_u32(a, b));
+  const uint16x4_t cd_u16 = vmovn_u32(vcombine_u32(c, d));
+  return vmovn_u16(vcombine_u16(ab_u16, cd_u16));
+}
+
+template <>
+struct type_casting_traits<numext::uint32_t, numext::int8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16c pcast<Packet4ui, Packet16c>(const Packet4ui& a, const Packet4ui& b, const Packet4ui& c,
+                                                          const Packet4ui& d) {
+  return vreinterpretq_s8_u8(pcast<Packet4ui, Packet16uc>(a, b, c, d));
+}
+template <>
+EIGEN_STRONG_INLINE Packet8c pcast<Packet2ui, Packet8c>(const Packet2ui& a, const Packet2ui& b, const Packet2ui& c,
+                                                        const Packet2ui& d) {
+  return vreinterpret_s8_u8(pcast<Packet2ui, Packet8uc>(a, b, c, d));
+}
+
+//==============================================================================
+// pcast, SrcType = int64_t
+//==============================================================================
+template <>
+struct type_casting_traits<numext::int64_t, float> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4f pcast<Packet2l, Packet4f>(const Packet2l& a, const Packet2l& b) {
+  return vcvtq_f32_s32(vcombine_s32(vmovn_s64(a), vmovn_s64(b)));
+}
+
+template <>
+struct type_casting_traits<numext::int64_t, numext::int64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2l pcast<Packet2l, Packet2l>(const Packet2l& a) {
+  return a;
+}
+
+template <>
+struct type_casting_traits<numext::int64_t, numext::uint64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2ul pcast<Packet2l, Packet2ul>(const Packet2l& a) {
+  return vreinterpretq_u64_s64(a);
+}
+
+template <>
+struct type_casting_traits<numext::int64_t, numext::int32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4i pcast<Packet2l, Packet4i>(const Packet2l& a, const Packet2l& b) {
+  return vcombine_s32(vmovn_s64(a), vmovn_s64(b));
+}
+
+template <>
+struct type_casting_traits<numext::int64_t, numext::uint32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4ui pcast<Packet2l, Packet4ui>(const Packet2l& a, const Packet2l& b) {
+  return vcombine_u32(vmovn_u64(vreinterpretq_u64_s64(a)), vmovn_u64(vreinterpretq_u64_s64(b)));
+}
+
+template <>
+struct type_casting_traits<numext::int64_t, numext::int16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8s pcast<Packet2l, Packet8s>(const Packet2l& a, const Packet2l& b, const Packet2l& c,
+                                                       const Packet2l& d) {
+  const int32x4_t ab_s32 = pcast<Packet2l, Packet4i>(a, b);
+  const int32x4_t cd_s32 = pcast<Packet2l, Packet4i>(c, d);
+  return vcombine_s16(vmovn_s32(ab_s32), vmovn_s32(cd_s32));
+}
+
+template <>
+struct type_casting_traits<numext::int64_t, numext::uint16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8us pcast<Packet2l, Packet8us>(const Packet2l& a, const Packet2l& b, const Packet2l& c,
+                                                         const Packet2l& d) {
+  const uint32x4_t ab_u32 = pcast<Packet2l, Packet4ui>(a, b);
+  const uint32x4_t cd_u32 = pcast<Packet2l, Packet4ui>(c, d);
+  return vcombine_u16(vmovn_u32(ab_u32), vmovn_u32(cd_u32));
+}
+
+template <>
+struct type_casting_traits<numext::int64_t, numext::int8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 8, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16c pcast<Packet2l, Packet16c>(const Packet2l& a, const Packet2l& b, const Packet2l& c,
+                                                         const Packet2l& d, const Packet2l& e, const Packet2l& f,
+                                                         const Packet2l& g, const Packet2l& h) {
+  const int16x8_t abcd_s16 = pcast<Packet2l, Packet8s>(a, b, c, d);
+  const int16x8_t efgh_s16 = pcast<Packet2l, Packet8s>(e, f, g, h);
+  return vcombine_s8(vmovn_s16(abcd_s16), vmovn_s16(efgh_s16));
+}
+
+template <>
+struct type_casting_traits<numext::int64_t, numext::uint8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 8, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16uc pcast<Packet2l, Packet16uc>(const Packet2l& a, const Packet2l& b, const Packet2l& c,
+                                                           const Packet2l& d, const Packet2l& e, const Packet2l& f,
+                                                           const Packet2l& g, const Packet2l& h) {
+  const uint16x8_t abcd_u16 = pcast<Packet2l, Packet8us>(a, b, c, d);
+  const uint16x8_t efgh_u16 = pcast<Packet2l, Packet8us>(e, f, g, h);
+  return vcombine_u8(vmovn_u16(abcd_u16), vmovn_u16(efgh_u16));
+}
+
+//==============================================================================
+// pcast, SrcType = uint64_t
+//==============================================================================
+template <>
+struct type_casting_traits<numext::uint64_t, float> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4f pcast<Packet2ul, Packet4f>(const Packet2ul& a, const Packet2ul& b) {
+  return vcvtq_f32_u32(vcombine_u32(vmovn_u64(a), vmovn_u64(b)));
+}
+
+template <>
+struct type_casting_traits<numext::uint64_t, numext::uint64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2ul pcast<Packet2ul, Packet2ul>(const Packet2ul& a) {
+  return a;
+}
+
+template <>
+struct type_casting_traits<numext::uint64_t, numext::int64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2l pcast<Packet2ul, Packet2l>(const Packet2ul& a) {
+  return vreinterpretq_s64_u64(a);
+}
+
+template <>
+struct type_casting_traits<numext::uint64_t, numext::uint32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4ui pcast<Packet2ul, Packet4ui>(const Packet2ul& a, const Packet2ul& b) {
+  return vcombine_u32(vmovn_u64(a), vmovn_u64(b));
+}
+
+template <>
+struct type_casting_traits<numext::uint64_t, numext::int32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4i pcast<Packet2ul, Packet4i>(const Packet2ul& a, const Packet2ul& b) {
+  return vreinterpretq_s32_u32(pcast<Packet2ul, Packet4ui>(a, b));
+}
+
+template <>
+struct type_casting_traits<numext::uint64_t, numext::uint16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8us pcast<Packet2ul, Packet8us>(const Packet2ul& a, const Packet2ul& b, const Packet2ul& c,
+                                                          const Packet2ul& d) {
+  const uint16x4_t ab_u16 = vmovn_u32(vcombine_u32(vmovn_u64(a), vmovn_u64(b)));
+  const uint16x4_t cd_u16 = vmovn_u32(vcombine_u32(vmovn_u64(c), vmovn_u64(d)));
+  return vcombine_u16(ab_u16, cd_u16);
+}
+
+template <>
+struct type_casting_traits<numext::uint64_t, numext::int16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8s pcast<Packet2ul, Packet8s>(const Packet2ul& a, const Packet2ul& b, const Packet2ul& c,
+                                                        const Packet2ul& d) {
+  return vreinterpretq_s16_u16(pcast<Packet2ul, Packet8us>(a, b, c, d));
+}
+
+template <>
+struct type_casting_traits<numext::uint64_t, numext::uint8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 8, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16uc pcast<Packet2ul, Packet16uc>(const Packet2ul& a, const Packet2ul& b, const Packet2ul& c,
+                                                            const Packet2ul& d, const Packet2ul& e, const Packet2ul& f,
+                                                            const Packet2ul& g, const Packet2ul& h) {
+  const uint16x8_t abcd_u16 = pcast<Packet2ul, Packet8us>(a, b, c, d);
+  const uint16x8_t efgh_u16 = pcast<Packet2ul, Packet8us>(e, f, g, h);
+  return vcombine_u8(vmovn_u16(abcd_u16), vmovn_u16(efgh_u16));
+}
+
+template <>
+struct type_casting_traits<numext::uint64_t, numext::int8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 8, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16c pcast<Packet2ul, Packet16c>(const Packet2ul& a, const Packet2ul& b, const Packet2ul& c,
+                                                          const Packet2ul& d, const Packet2ul& e, const Packet2ul& f,
+                                                          const Packet2ul& g, const Packet2ul& h) {
+  return vreinterpretq_s8_u8(pcast<Packet2ul, Packet16uc>(a, b, c, d, e, f, g, h));
+}
+
+//==============================================================================
+// preinterpret
+//==============================================================================
+template <>
+EIGEN_STRONG_INLINE Packet2f preinterpret<Packet2f, Packet2i>(const Packet2i& a) {
+  return vreinterpret_f32_s32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2f preinterpret<Packet2f, Packet2ui>(const Packet2ui& a) {
+  return vreinterpret_f32_u32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet4f preinterpret<Packet4f, Packet4i>(const Packet4i& a) {
+  return vreinterpretq_f32_s32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet4f preinterpret<Packet4f, Packet4ui>(const Packet4ui& a) {
+  return vreinterpretq_f32_u32(a);
+}
+
+template <>
+EIGEN_STRONG_INLINE Packet4c preinterpret<Packet4c, Packet4uc>(const Packet4uc& a) {
+  return static_cast<Packet4c>(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet8c preinterpret<Packet8c, Packet8uc>(const Packet8uc& a) {
+  return vreinterpret_s8_u8(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet16c preinterpret<Packet16c, Packet16uc>(const Packet16uc& a) {
+  return vreinterpretq_s8_u8(a);
+}
+
+template <>
+EIGEN_STRONG_INLINE Packet4uc preinterpret<Packet4uc, Packet4c>(const Packet4c& a) {
+  return static_cast<Packet4uc>(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet8uc preinterpret<Packet8uc, Packet8c>(const Packet8c& a) {
+  return vreinterpret_u8_s8(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet16uc preinterpret<Packet16uc, Packet16c>(const Packet16c& a) {
+  return vreinterpretq_u8_s8(a);
+}
+
+template <>
+EIGEN_STRONG_INLINE Packet4s preinterpret<Packet4s, Packet4us>(const Packet4us& a) {
+  return vreinterpret_s16_u16(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet8s preinterpret<Packet8s, Packet8us>(const Packet8us& a) {
+  return vreinterpretq_s16_u16(a);
+}
+
+template <>
+EIGEN_STRONG_INLINE Packet4us preinterpret<Packet4us, Packet4s>(const Packet4s& a) {
+  return vreinterpret_u16_s16(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet8us preinterpret<Packet8us, Packet8s>(const Packet8s& a) {
+  return vreinterpretq_u16_s16(a);
+}
+
+template <>
+EIGEN_STRONG_INLINE Packet2i preinterpret<Packet2i, Packet2f>(const Packet2f& a) {
+  return vreinterpret_s32_f32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2i preinterpret<Packet2i, Packet2ui>(const Packet2ui& a) {
+  return vreinterpret_s32_u32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet4i preinterpret<Packet4i, Packet4f>(const Packet4f& a) {
+  return vreinterpretq_s32_f32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet4i preinterpret<Packet4i, Packet4ui>(const Packet4ui& a) {
+  return vreinterpretq_s32_u32(a);
+}
+
+template <>
+EIGEN_STRONG_INLINE Packet2ui preinterpret<Packet2ui, Packet2f>(const Packet2f& a) {
+  return vreinterpret_u32_f32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ui preinterpret<Packet2ui, Packet2i>(const Packet2i& a) {
+  return vreinterpret_u32_s32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet4ui preinterpret<Packet4ui, Packet4f>(const Packet4f& a) {
+  return vreinterpretq_u32_f32(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet4ui preinterpret<Packet4ui, Packet4i>(const Packet4i& a) {
+  return vreinterpretq_u32_s32(a);
+}
+
+template <>
+EIGEN_STRONG_INLINE Packet2l preinterpret<Packet2l, Packet2ul>(const Packet2ul& a) {
+  return vreinterpretq_s64_u64(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ul preinterpret<Packet2ul, Packet2l>(const Packet2l& a) {
+  return vreinterpretq_u64_s64(a);
+}
+
+#if EIGEN_ARCH_ARM64
+
+//==============================================================================
+// pcast/preinterpret, Double
+//==============================================================================
+
+template <>
+struct type_casting_traits<double, double> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2d pcast<Packet2d, Packet2d>(const Packet2d& a) {
+  return a;
+}
+
+template <>
+struct type_casting_traits<double, float> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4f pcast<Packet2d, Packet4f>(const Packet2d& a, const Packet2d& b) {
+  return vcombine_f32(vcvt_f32_f64(a), vcvt_f32_f64(b));
+}
+
+template <>
+struct type_casting_traits<double, numext::int64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2l pcast<Packet2d, Packet2l>(const Packet2d& a) {
+  return vcvtq_s64_f64(a);
+}
+
+template <>
+struct type_casting_traits<double, numext::uint64_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2ul pcast<Packet2d, Packet2ul>(const Packet2d& a) {
+  return vcvtq_u64_f64(a);
+}
+
+template <>
+struct type_casting_traits<double, numext::int32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4i pcast<Packet2d, Packet4i>(const Packet2d& a, const Packet2d& b) {
+  return vcombine_s32(vmovn_s64(vcvtq_s64_f64(a)), vmovn_s64(vcvtq_s64_f64(b)));
+}
+
+template <>
+struct type_casting_traits<double, numext::uint32_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 2, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet4ui pcast<Packet2d, Packet4ui>(const Packet2d& a, const Packet2d& b) {
+  return vcombine_u32(vmovn_u64(vcvtq_u64_f64(a)), vmovn_u64(vcvtq_u64_f64(b)));
+}
+
+template <>
+struct type_casting_traits<double, numext::int16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8s pcast<Packet2d, Packet8s>(const Packet2d& a, const Packet2d& b, const Packet2d& c,
+                                                       const Packet2d& d) {
+  const int32x4_t ab_s32 = pcast<Packet2d, Packet4i>(a, b);
+  const int32x4_t cd_s32 = pcast<Packet2d, Packet4i>(c, d);
+  return vcombine_s16(vmovn_s32(ab_s32), vmovn_s32(cd_s32));
+}
+
+template <>
+struct type_casting_traits<double, numext::uint16_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 4, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet8us pcast<Packet2d, Packet8us>(const Packet2d& a, const Packet2d& b, const Packet2d& c,
+                                                         const Packet2d& d) {
+  const uint32x4_t ab_u32 = pcast<Packet2d, Packet4ui>(a, b);
+  const uint32x4_t cd_u32 = pcast<Packet2d, Packet4ui>(c, d);
+  return vcombine_u16(vmovn_u32(ab_u32), vmovn_u32(cd_u32));
+}
+
+template <>
+struct type_casting_traits<double, numext::int8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 8, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16c pcast<Packet2d, Packet16c>(const Packet2d& a, const Packet2d& b, const Packet2d& c,
+                                                         const Packet2d& d, const Packet2d& e, const Packet2d& f,
+                                                         const Packet2d& g, const Packet2d& h) {
+  const int16x8_t abcd_s16 = pcast<Packet2d, Packet8s>(a, b, c, d);
+  const int16x8_t efgh_s16 = pcast<Packet2d, Packet8s>(e, f, g, h);
+  return vcombine_s8(vmovn_s16(abcd_s16), vmovn_s16(efgh_s16));
+}
+
+template <>
+struct type_casting_traits<double, numext::uint8_t> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 8, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet16uc pcast<Packet2d, Packet16uc>(const Packet2d& a, const Packet2d& b, const Packet2d& c,
+                                                           const Packet2d& d, const Packet2d& e, const Packet2d& f,
+                                                           const Packet2d& g, const Packet2d& h) {
+  const uint16x8_t abcd_u16 = pcast<Packet2d, Packet8us>(a, b, c, d);
+  const uint16x8_t efgh_u16 = pcast<Packet2d, Packet8us>(e, f, g, h);
+  return vcombine_u8(vmovn_u16(abcd_u16), vmovn_u16(efgh_u16));
+}
+
+template <>
+struct type_casting_traits<float, double> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2d pcast<Packet4f, Packet2d>(const Packet4f& a) {
+  // Discard second-half of input.
+  return vcvt_f64_f32(vget_low_f32(a));
+}
+
+template <>
+struct type_casting_traits<numext::int8_t, double> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 8 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2d pcast<Packet16c, Packet2d>(const Packet16c& a) {
+  // Discard all but first two values.
+  return vcvt_f64_f32(pcast<Packet8c, Packet2f>(vget_low_s8(a)));
+}
+
+template <>
+struct type_casting_traits<numext::uint8_t, double> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 8 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2d pcast<Packet16uc, Packet2d>(const Packet16uc& a) {
+  // Discard all but first two values.
+  return vcvt_f64_f32(pcast<Packet8uc, Packet2f>(vget_low_u8(a)));
+}
+
+template <>
+struct type_casting_traits<numext::int16_t, double> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2d pcast<Packet8s, Packet2d>(const Packet8s& a) {
+  // Discard all but first two values.
+  return vcvt_f64_f32(pcast<Packet4s, Packet2f>(vget_low_s16(a)));
+}
+
+template <>
+struct type_casting_traits<numext::uint16_t, double> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 4 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2d pcast<Packet8us, Packet2d>(const Packet8us& a) {
+  // Discard all but first two values.
+  return vcvt_f64_f32(pcast<Packet4us, Packet2f>(vget_low_u16(a)));
+}
+
+template <>
+struct type_casting_traits<numext::int32_t, double> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2d pcast<Packet4i, Packet2d>(const Packet4i& a) {
+  // Discard second half of input.
+  return vcvtq_f64_s64(vmovl_s32(vget_low_s32(a)));
+}
+
+template <>
+struct type_casting_traits<numext::uint32_t, double> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 2 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2d pcast<Packet4ui, Packet2d>(const Packet4ui& a) {
+  // Discard second half of input.
+  return vcvtq_f64_u64(vmovl_u32(vget_low_u32(a)));
+}
+
+template <>
+struct type_casting_traits<numext::int64_t, double> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2d pcast<Packet2l, Packet2d>(const Packet2l& a) {
+  return vcvtq_f64_s64(a);
+}
+
+template <>
+struct type_casting_traits<numext::uint64_t, double> {
+  enum { VectorizedCast = 1, SrcCoeffRatio = 1, TgtCoeffRatio = 1 };
+};
+template <>
+EIGEN_STRONG_INLINE Packet2d pcast<Packet2ul, Packet2d>(const Packet2ul& a) {
+  return vcvtq_f64_u64(a);
+}
+
+template <>
+EIGEN_STRONG_INLINE Packet2d preinterpret<Packet2d, Packet2l>(const Packet2l& a) {
+  return vreinterpretq_f64_s64(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2d preinterpret<Packet2d, Packet2ul>(const Packet2ul& a) {
+  return vreinterpretq_f64_u64(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2l preinterpret<Packet2l, Packet2d>(const Packet2d& a) {
+  return vreinterpretq_s64_f64(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet2ul preinterpret<Packet2ul, Packet2d>(const Packet2d& a) {
+  return vreinterpretq_u64_f64(a);
+}
+
+#endif  // EIGEN_ARCH_ARM64
+
+}  // end namespace internal
+
+}  // end namespace Eigen
 
-template<> EIGEN_STRONG_INLINE Packet2f pcast<Packet2d,Packet2f>(const Packet2d& a) { return vcvt_f32_f64(a); }
-template<> EIGEN_STRONG_INLINE Packet2d pcast<Packet2f,Packet2d>(const Packet2f& a) { return vcvt_f64_f32(a); }
-template<> EIGEN_STRONG_INLINE Packet2d pcast<Packet2i,Packet2d>(const Packet2i& a) { return vcvtq_f64_s64(vmovl_s32(a)); }
-template<> EIGEN_STRONG_INLINE Packet2d pcast<Packet2ui,Packet2d>(const Packet2ui& a) { return vcvtq_f64_u64(vmovl_u32(a)); }
-template<> EIGEN_STRONG_INLINE Packet2d pcast<Packet2l,Packet2d>(const Packet2l& a) { return vcvtq_f64_s64(a); }
-template<> EIGEN_STRONG_INLINE Packet2d pcast<Packet2ul,Packet2d>(const Packet2ul& a) { return vcvtq_f64_u64(a); }
-template<> EIGEN_STRONG_INLINE Packet2i pcast<Packet2d,Packet2i>(const Packet2d& a) { return vcvt_s32_f32(vcvt_f32_f64(a)); }
-template<> EIGEN_STRONG_INLINE Packet2ui pcast<Packet2d,Packet2ui>(const Packet2d& a) { return vcvt_u32_f32(vcvt_f32_f64(a)); }
-template<> EIGEN_STRONG_INLINE Packet2l pcast<Packet2d,Packet2l>(const Packet2d& a) { return vcvtq_s64_f64(a); }
-template<> EIGEN_STRONG_INLINE Packet2ul pcast<Packet2d,Packet2ul>(const Packet2d& a) { return vcvtq_u64_f64(a); }
-
-template<> EIGEN_STRONG_INLINE Packet2d preinterpret<Packet2d,Packet2l>(const Packet2l& a)
-{ return vreinterpretq_f64_s64(a); }
-template<> EIGEN_STRONG_INLINE Packet2d preinterpret<Packet2d,Packet2ul>(const Packet2ul& a)
-{ return vreinterpretq_f64_u64(a); }
-template<> EIGEN_STRONG_INLINE Packet2l preinterpret<Packet2l,Packet2d>(const Packet2d& a)
-{ return vreinterpretq_s64_f64(a); }
-template<> EIGEN_STRONG_INLINE Packet2ul preinterpret<Packet2ul,Packet2d>(const Packet2d& a)
-{ return vreinterpretq_u64_f64(a); }
-
-#endif // EIGEN_ARCH_ARM64
-
-} // end namespace internal
-
-} // end namespace Eigen
-
-#endif // EIGEN_TYPE_CASTING_NEON_H
+#endif  // EIGEN_TYPE_CASTING_NEON_H
diff --git a/test/packetmath.cpp b/test/packetmath.cpp
index a82b2b87a..a57a3b6a5 100644
--- a/test/packetmath.cpp
+++ b/test/packetmath.cpp
@@ -11,95 +11,264 @@
 #include "packetmath_test_shared.h"
 
 template <typename T>
-inline T REF_ADD(const T& a, const T& b) { return a + b;}
+inline T REF_ADD(const T& a, const T& b) {
+  return a + b;
+}
 template <typename T>
-inline T REF_SUB(const T& a, const T& b) { return a - b;}
+inline T REF_SUB(const T& a, const T& b) {
+  return a - b;
+}
 template <typename T>
-inline T REF_MUL(const T& a, const T& b) { return a * b;}
+inline T REF_MUL(const T& a, const T& b) {
+  return a * b;
+}
 template <typename T>
-inline T REF_DIV(const T& a, const T& b) { return a / b;}
+inline T REF_DIV(const T& a, const T& b) {
+  return a / b;
+}
 template <typename T>
-inline T REF_ABS_DIFF(const T& a, const T& b) { return a>b ? a - b : b-a;}
+inline T REF_ABS_DIFF(const T& a, const T& b) {
+  return a > b ? a - b : b - a;
+}
 
 // Specializations for bool.
 template <>
-inline bool REF_ADD(const bool& a, const bool& b) { return a || b;}
+inline bool REF_ADD(const bool& a, const bool& b) {
+  return a || b;
+}
 template <>
-inline bool REF_SUB(const bool& a, const bool& b) { return a ^ b;}
+inline bool REF_SUB(const bool& a, const bool& b) {
+  return a ^ b;
+}
 template <>
-inline bool REF_MUL(const bool& a, const bool& b) { return a && b;}
+inline bool REF_MUL(const bool& a, const bool& b) {
+  return a && b;
+}
 
 // Uses pcast to cast from one array to another.
-template<typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio>
-struct pcast_array {
-  static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size, typename internal::unpacket_traits<TgtPacket>::type* dst) {
+template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio>
+struct pcast_array;
+
+template <typename SrcPacket, typename TgtPacket, int TgtCoeffRatio>
+struct pcast_array<SrcPacket, TgtPacket, 1, TgtCoeffRatio> {
+  typedef typename internal::unpacket_traits<SrcPacket>::type SrcScalar;
+  typedef typename internal::unpacket_traits<TgtPacket>::type TgtScalar;
+  static void cast(const SrcScalar* src, size_t size, TgtScalar* dst) {
+    static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
     static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
-    for (size_t i=0; i<size; i+=TgtPacketSize) {
-      internal::pstoreu(dst+i, internal::pcast<SrcPacket,TgtPacket>(internal::ploadu<SrcPacket>(src+i)));
+    size_t i;
+    for (i = 0; i < size && i + SrcPacketSize <= size; i += TgtPacketSize) {
+      internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(internal::ploadu<SrcPacket>(src + i)));
+    }
+    // Leftovers that cannot be loaded into a packet.
+    for (; i < size; ++i) {
+      dst[i] = static_cast<TgtScalar>(src[i]);
     }
   }
 };
 
-template<typename SrcPacket,typename TgtPacket>
-struct pcast_array<SrcPacket, TgtPacket, 2, 1>{
-  static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size, typename internal::unpacket_traits<TgtPacket>::type* dst) {
+template <typename SrcPacket, typename TgtPacket>
+struct pcast_array<SrcPacket, TgtPacket, 2, 1> {
+  static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size,
+                   typename internal::unpacket_traits<TgtPacket>::type* dst) {
     static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
     static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
-    for (size_t i=0; i<size; i+=TgtPacketSize) {
-      SrcPacket a = internal::ploadu<SrcPacket>(src+i);
-      SrcPacket b = internal::ploadu<SrcPacket>(src+i+SrcPacketSize);
-      internal::pstoreu(dst+i, internal::pcast<SrcPacket,TgtPacket>(a, b));
+    for (size_t i = 0; i < size; i += TgtPacketSize) {
+      SrcPacket a = internal::ploadu<SrcPacket>(src + i);
+      SrcPacket b = internal::ploadu<SrcPacket>(src + i + SrcPacketSize);
+      internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(a, b));
     }
   }
 };
 
-template<typename SrcPacket, typename TgtPacket>
-struct pcast_array<SrcPacket, TgtPacket, 4, 1>{
-  static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size, typename internal::unpacket_traits<TgtPacket>::type* dst) {
+template <typename SrcPacket, typename TgtPacket>
+struct pcast_array<SrcPacket, TgtPacket, 4, 1> {
+  static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size,
+                   typename internal::unpacket_traits<TgtPacket>::type* dst) {
     static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
     static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
-    for (size_t i=0; i<size; i+=TgtPacketSize) {
-      SrcPacket a = internal::ploadu<SrcPacket>(src+i);
-      SrcPacket b = internal::ploadu<SrcPacket>(src+i+SrcPacketSize);
-      SrcPacket c = internal::ploadu<SrcPacket>(src+i+2*SrcPacketSize);
-      SrcPacket d = internal::ploadu<SrcPacket>(src+i+3*SrcPacketSize);
-      internal::pstoreu(dst+i, internal::pcast<SrcPacket,TgtPacket>(a, b));
+    for (size_t i = 0; i < size; i += TgtPacketSize) {
+      SrcPacket a = internal::ploadu<SrcPacket>(src + i);
+      SrcPacket b = internal::ploadu<SrcPacket>(src + i + SrcPacketSize);
+      SrcPacket c = internal::ploadu<SrcPacket>(src + i + 2 * SrcPacketSize);
+      SrcPacket d = internal::ploadu<SrcPacket>(src + i + 3 * SrcPacketSize);
+      internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(a, b, c, d));
     }
   }
 };
 
-template<typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio, bool CanCast = false>
+template <typename SrcPacket, typename TgtPacket>
+struct pcast_array<SrcPacket, TgtPacket, 8, 1> {
+  static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size,
+                   typename internal::unpacket_traits<TgtPacket>::type* dst) {
+    static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
+    static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
+    for (size_t i = 0; i < size; i += TgtPacketSize) {
+      SrcPacket a = internal::ploadu<SrcPacket>(src + i);
+      SrcPacket b = internal::ploadu<SrcPacket>(src + i + SrcPacketSize);
+      SrcPacket c = internal::ploadu<SrcPacket>(src + i + 2 * SrcPacketSize);
+      SrcPacket d = internal::ploadu<SrcPacket>(src + i + 3 * SrcPacketSize);
+      SrcPacket e = internal::ploadu<SrcPacket>(src + i + 4 * SrcPacketSize);
+      SrcPacket f = internal::ploadu<SrcPacket>(src + i + 5 * SrcPacketSize);
+      SrcPacket g = internal::ploadu<SrcPacket>(src + i + 6 * SrcPacketSize);
+      SrcPacket h = internal::ploadu<SrcPacket>(src + i + 7 * SrcPacketSize);
+      internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(a, b, c, d, e, f, g, h));
+    }
+  }
+};
+
+template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio, bool CanCast = false>
 struct test_cast_helper;
 
-template<typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio>
+template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio>
 struct test_cast_helper<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio, false> {
   static void run() {}
 };
 
-template<typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio>
+// Generates random values that fit in both SrcScalar and TgtScalar without
+// overflowing when cast.
+template <typename SrcScalar, typename TgtScalar, typename EnableIf = void>
+struct random_without_cast_overflow {
+  static SrcScalar value() { return internal::random<SrcScalar>(); }
+};
+
+// Widening integer cast signed to unsigned.
+template <typename SrcScalar, typename TgtScalar>
+struct random_without_cast_overflow<
+    SrcScalar, TgtScalar,
+    typename internal::enable_if<NumTraits<SrcScalar>::IsInteger && NumTraits<TgtScalar>::IsInteger &&
+                                 !NumTraits<TgtScalar>::IsSigned &&
+                                 (std::numeric_limits<SrcScalar>::digits < std::numeric_limits<TgtScalar>::digits ||
+                                  (std::numeric_limits<SrcScalar>::digits == std::numeric_limits<TgtScalar>::digits &&
+                                   NumTraits<SrcScalar>::IsSigned))>::type> {
+  static SrcScalar value() {
+    SrcScalar a = internal::random<SrcScalar>();
+    return a < SrcScalar(0) ? -(a + 1) : a;
+  }
+};
+
+// Narrowing integer cast to unsigned.
+template <typename SrcScalar, typename TgtScalar>
+struct random_without_cast_overflow<
+    SrcScalar, TgtScalar,
+    typename internal::enable_if<
+        NumTraits<SrcScalar>::IsInteger && NumTraits<TgtScalar>::IsInteger && !NumTraits<SrcScalar>::IsSigned &&
+        (std::numeric_limits<SrcScalar>::digits > std::numeric_limits<TgtScalar>::digits)>::type> {
+  static SrcScalar value() {
+    TgtScalar b = internal::random<TgtScalar>();
+    return static_cast<SrcScalar>(b < TgtScalar(0) ? -(b + 1) : b);
+  }
+};
+
+// Narrowing integer cast to signed.
+template <typename SrcScalar, typename TgtScalar>
+struct random_without_cast_overflow<
+    SrcScalar, TgtScalar,
+    typename internal::enable_if<
+        NumTraits<SrcScalar>::IsInteger && NumTraits<TgtScalar>::IsInteger && NumTraits<SrcScalar>::IsSigned &&
+        (std::numeric_limits<SrcScalar>::digits > std::numeric_limits<TgtScalar>::digits)>::type> {
+  static SrcScalar value() {
+    TgtScalar b = internal::random<TgtScalar>();
+    return static_cast<SrcScalar>(b);
+  }
+};
+
+// Unsigned to signed narrowing cast.
+template <typename SrcScalar, typename TgtScalar>
+struct random_without_cast_overflow<
+    SrcScalar, TgtScalar,
+    typename internal::enable_if<NumTraits<SrcScalar>::IsInteger && NumTraits<TgtScalar>::IsInteger &&
+                                 !NumTraits<SrcScalar>::IsSigned && NumTraits<TgtScalar>::IsSigned &&
+                                 (std::numeric_limits<SrcScalar>::digits ==
+                                  std::numeric_limits<TgtScalar>::digits)>::type> {
+  static SrcScalar value() { return internal::random<SrcScalar>() / 2; }
+};
+
+template <typename Scalar>
+struct is_floating_point {
+  enum { value = 0 };
+};
+template <>
+struct is_floating_point<float> {
+  enum { value = 1 };
+};
+template <>
+struct is_floating_point<double> {
+  enum { value = 1 };
+};
+template <>
+struct is_floating_point<half> {
+  enum { value = 1 };
+};
+template <>
+struct is_floating_point<bfloat16> {
+  enum { value = 1 };
+};
+
+// Floating-point to integer, full precision.
+template <typename SrcScalar, typename TgtScalar>
+struct random_without_cast_overflow<
+    SrcScalar, TgtScalar,
+    typename internal::enable_if<is_floating_point<SrcScalar>::value && NumTraits<TgtScalar>::IsInteger &&
+                                 (std::numeric_limits<TgtScalar>::digits <=
+                                  std::numeric_limits<SrcScalar>::digits)>::type> {
+  static SrcScalar value() { return static_cast<SrcScalar>(internal::random<TgtScalar>()); }
+};
+
+// Floating-point to integer, narrowing precision.
+template <typename SrcScalar, typename TgtScalar>
+struct random_without_cast_overflow<
+    SrcScalar, TgtScalar,
+    typename internal::enable_if<is_floating_point<SrcScalar>::value && NumTraits<TgtScalar>::IsInteger &&
+                                 (std::numeric_limits<TgtScalar>::digits >
+                                  std::numeric_limits<SrcScalar>::digits)>::type> {
+  static SrcScalar value() {
+    static const int BitShift = std::numeric_limits<TgtScalar>::digits - std::numeric_limits<SrcScalar>::digits;
+    return static_cast<SrcScalar>(internal::random<TgtScalar>() >> BitShift);
+  }
+};
+
+// Floating-point target from integer, re-use above logic.
+template <typename SrcScalar, typename TgtScalar>
+struct random_without_cast_overflow<
+    SrcScalar, TgtScalar,
+    typename internal::enable_if<NumTraits<SrcScalar>::IsInteger && is_floating_point<TgtScalar>::value>::type> {
+  static SrcScalar value() {
+    return static_cast<SrcScalar>(random_without_cast_overflow<TgtScalar, SrcScalar>::value());
+  }
+};
+
+// Floating-point narrowing conversion.
+template <typename SrcScalar, typename TgtScalar>
+struct random_without_cast_overflow<
+    SrcScalar, TgtScalar,
+    typename internal::enable_if<is_floating_point<SrcScalar>::value && is_floating_point<TgtScalar>::value &&
+                                 (std::numeric_limits<SrcScalar>::digits >
+                                  std::numeric_limits<TgtScalar>::digits)>::type> {
+  static SrcScalar value() { return static_cast<SrcScalar>(internal::random<TgtScalar>()); }
+};
+
+template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio>
 struct test_cast_helper<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio, true> {
   static void run() {
     typedef typename internal::unpacket_traits<SrcPacket>::type SrcScalar;
     typedef typename internal::unpacket_traits<TgtPacket>::type TgtScalar;
     static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
     static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
-    static const int DataSize = SrcPacketSize*SrcCoeffRatio;
-    VERIFY(DataSize == TgtPacketSize*TgtCoeffRatio && "Packet sizes and cast ratios are mismatched.");
+    static const int BlockSize = SrcPacketSize * SrcCoeffRatio;
+    eigen_assert(BlockSize == TgtPacketSize * TgtCoeffRatio && "Packet sizes and cast ratios are mismatched.");
 
+    static const int DataSize = 10 * BlockSize;
     EIGEN_ALIGN_MAX SrcScalar data1[DataSize];
     EIGEN_ALIGN_MAX TgtScalar data2[DataSize];
     EIGEN_ALIGN_MAX TgtScalar ref[DataSize];
 
     // Construct a packet of scalars that will not overflow when casting
-    for (int i=0; i<DataSize; ++i) {
-      const SrcScalar a = Array<SrcScalar,1,1>::Random().value();
-      const TgtScalar b = Array<TgtScalar,1,1>::Random().value();
-      const SrcScalar c = sizeof(TgtScalar) > sizeof(SrcScalar) ? static_cast<SrcScalar>(b) : a;
-      data1[i] = (NumTraits<SrcScalar>::IsSigned && !NumTraits<TgtScalar>::IsSigned) ? numext::abs(c) : c;
+    for (int i = 0; i < DataSize; ++i) {
+      data1[i] = random_without_cast_overflow<SrcScalar, TgtScalar>::value();
     }
 
-    for (int i=0; i<DataSize; ++i)
-      ref[i] = static_cast<const TgtScalar>(data1[i]);
+    for (int i = 0; i < DataSize; ++i) ref[i] = static_cast<const TgtScalar>(data1[i]);
 
     pcast_array<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio>::cast(data1, DataSize, data2);
 
@@ -107,48 +276,72 @@ struct test_cast_helper<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio, true
   }
 };
 
-template<typename SrcScalar, typename TgtScalar>
-void test_cast() {
-  typedef typename internal::packet_traits<SrcScalar> SrcPacketTraits;
-  typedef typename internal::packet_traits<TgtScalar> TgtPacketTraits;
-  typedef typename internal::type_casting_traits<SrcScalar, TgtScalar> TypeCastingTraits;
-  static const int SrcCoeffRatio = TypeCastingTraits::SrcCoeffRatio;
-  static const int TgtCoeffRatio = TypeCastingTraits::TgtCoeffRatio;
+template <typename SrcPacket, typename TgtPacket>
+struct test_cast {
+  static void run() {
+    typedef typename internal::unpacket_traits<SrcPacket>::type SrcScalar;
+    typedef typename internal::unpacket_traits<TgtPacket>::type TgtScalar;
+    typedef typename internal::type_casting_traits<SrcScalar, TgtScalar> TypeCastingTraits;
+    static const int SrcCoeffRatio = TypeCastingTraits::SrcCoeffRatio;
+    static const int TgtCoeffRatio = TypeCastingTraits::TgtCoeffRatio;
+    static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
+    static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
+    static const bool HasCast =
+        internal::unpacket_traits<SrcPacket>::vectorizable && internal::unpacket_traits<TgtPacket>::vectorizable &&
+        TypeCastingTraits::VectorizedCast && (SrcPacketSize * SrcCoeffRatio == TgtPacketSize * TgtCoeffRatio);
+    test_cast_helper<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio, HasCast>::run();
+  }
+};
 
-  static const bool HasFullCast = TypeCastingTraits::VectorizedCast;
-  static const bool HasHalfCast = HasFullCast && internal::packet_traits<SrcScalar>::HasHalfPacket && internal::packet_traits<TgtScalar>::HasHalfPacket;
+template <typename SrcPacket, typename TgtScalar,
+          typename TgtPacket = typename internal::packet_traits<TgtScalar>::type,
+          bool Vectorized = internal::packet_traits<TgtScalar>::Vectorizable,
+          bool HasHalf = !internal::is_same<typename internal::unpacket_traits<TgtPacket>::half, TgtPacket>::value>
+struct test_cast_runner;
 
-  test_cast_helper<typename SrcPacketTraits::type, typename TgtPacketTraits::type, SrcCoeffRatio, TgtCoeffRatio, HasFullCast>::run();
-  test_cast_helper<typename SrcPacketTraits::half, typename TgtPacketTraits::half, SrcCoeffRatio, TgtCoeffRatio, HasHalfCast>::run();
-}
+template <typename SrcPacket, typename TgtScalar, typename TgtPacket>
+struct test_cast_runner<SrcPacket, TgtScalar, TgtPacket, true, false> {
+  static void run() { test_cast<SrcPacket, TgtPacket>::run(); }
+};
 
-template<typename Scalar, typename Packet> void packetmath_pcast_ops() {
-  const static bool IsFullPacket = internal::is_same<typename internal::packet_traits<Scalar>::type,Packet>::value;
-  if (IsFullPacket) {
-    test_cast<Scalar, float>();
-    test_cast<Scalar, double>();
-    test_cast<Scalar, int8_t>();
-    test_cast<Scalar, uint8_t>();
-    test_cast<Scalar, int16_t>();
-    test_cast<Scalar, uint16_t>();
-    test_cast<Scalar, int32_t>();
-    test_cast<Scalar, uint32_t>();
-    test_cast<Scalar, int64_t>();
-    test_cast<Scalar, uint64_t>();
+template <typename SrcPacket, typename TgtScalar, typename TgtPacket>
+struct test_cast_runner<SrcPacket, TgtScalar, TgtPacket, true, true> {
+  static void run() {
+    test_cast<SrcPacket, TgtPacket>::run();
+    test_cast_runner<SrcPacket, TgtScalar, typename internal::unpacket_traits<TgtPacket>::half>::run();
   }
+};
+
+template <typename SrcPacket, typename TgtScalar, typename TgtPacket>
+struct test_cast_runner<SrcPacket, TgtScalar, TgtPacket, false, false> {
+  static void run() {}
+};
+
+template <typename Scalar, typename Packet>
+void packetmath_pcast_ops() {
+  test_cast_runner<Packet, float>::run();
+  test_cast_runner<Packet, double>::run();
+  test_cast_runner<Packet, int8_t>::run();
+  test_cast_runner<Packet, uint8_t>::run();
+  test_cast_runner<Packet, int16_t>::run();
+  test_cast_runner<Packet, uint16_t>::run();
+  test_cast_runner<Packet, int32_t>::run();
+  test_cast_runner<Packet, uint32_t>::run();
+  test_cast_runner<Packet, int64_t>::run();
+  test_cast_runner<Packet, uint64_t>::run();
+  test_cast_runner<Packet, bool>::run();
+  test_cast_runner<Packet, half>::run();
 }
 
-template<typename Scalar,typename Packet>
-void packetmath_boolean_mask_ops()
-{
+template <typename Scalar, typename Packet>
+void packetmath_boolean_mask_ops() {
   const int PacketSize = internal::unpacket_traits<Packet>::size;
-  const int size = 2*PacketSize;
+  const int size = 2 * PacketSize;
   EIGEN_ALIGN_MAX Scalar data1[size];
   EIGEN_ALIGN_MAX Scalar data2[size];
   EIGEN_ALIGN_MAX Scalar ref[size];
 
-  for (int i=0; i<size; ++i)
-  {
+  for (int i = 0; i < size; ++i) {
     data1[i] = internal::random<Scalar>();
   }
   CHECK_CWISE1(internal::ptrue, internal::ptrue);
@@ -163,82 +356,68 @@ void packetmath_boolean_mask_ops()
 // Packet16b representing bool does not support ptrue, pandnot or pcmp_eq, since the scalar path
 // (for some compilers) compute the bitwise and with 0x1 of the results to keep the value in [0,1].
 #ifdef EIGEN_PACKET_MATH_SSE_H
-template<>
-void packetmath_boolean_mask_ops<bool, internal::Packet16b>()
-{
-}
+template <>
+void packetmath_boolean_mask_ops<bool, internal::Packet16b>() {}
 #endif
 
-template<typename Scalar,typename Packet> void packetmath()
-{
+template <typename Scalar, typename Packet>
+void packetmath() {
   typedef internal::packet_traits<Scalar> PacketTraits;
   const int PacketSize = internal::unpacket_traits<Packet>::size;
   typedef typename NumTraits<Scalar>::Real RealScalar;
 
   if (g_first_pass)
-    std::cerr << "=== Testing packet of type '" << typeid(Packet).name()
-              << "' and scalar type '" << typeid(Scalar).name()
-              << "' and size '" << PacketSize << "' ===\n" ;
+    std::cerr << "=== Testing packet of type '" << typeid(Packet).name() << "' and scalar type '"
+              << typeid(Scalar).name() << "' and size '" << PacketSize << "' ===\n";
 
   const int max_size = PacketSize > 4 ? PacketSize : 4;
-  const int size = PacketSize*max_size;
+  const int size = PacketSize * max_size;
   EIGEN_ALIGN_MAX Scalar data1[size];
   EIGEN_ALIGN_MAX Scalar data2[size];
   EIGEN_ALIGN_MAX Scalar data3[size];
   EIGEN_ALIGN_MAX Scalar ref[size];
   RealScalar refvalue = RealScalar(0);
-  for (int i=0; i<size; ++i)
-  {
-    data1[i] = internal::random<Scalar>()/RealScalar(PacketSize);
-    data2[i] = internal::random<Scalar>()/RealScalar(PacketSize);
+  for (int i = 0; i < size; ++i) {
+    data1[i] = internal::random<Scalar>() / RealScalar(PacketSize);
+    data2[i] = internal::random<Scalar>() / RealScalar(PacketSize);
     refvalue = (std::max)(refvalue, numext::abs(data1[i]));
   }
 
   internal::pstore(data2, internal::pload<Packet>(data1));
   VERIFY(test::areApprox(data1, data2, PacketSize) && "aligned load/store");
 
-  for (int offset=0; offset<PacketSize; ++offset)
-  {
-    internal::pstore(data2, internal::ploadu<Packet>(data1+offset));
-    VERIFY(test::areApprox(data1+offset, data2, PacketSize) && "internal::ploadu");
+  for (int offset = 0; offset < PacketSize; ++offset) {
+    internal::pstore(data2, internal::ploadu<Packet>(data1 + offset));
+    VERIFY(test::areApprox(data1 + offset, data2, PacketSize) && "internal::ploadu");
   }
 
-  for (int offset=0; offset<PacketSize; ++offset)
-  {
-    internal::pstoreu(data2+offset, internal::pload<Packet>(data1));
-    VERIFY(test::areApprox(data1, data2+offset, PacketSize) && "internal::pstoreu");
+  for (int offset = 0; offset < PacketSize; ++offset) {
+    internal::pstoreu(data2 + offset, internal::pload<Packet>(data1));
+    VERIFY(test::areApprox(data1, data2 + offset, PacketSize) && "internal::pstoreu");
   }
 
-  if (internal::unpacket_traits<Packet>::masked_load_available)
-  {
+  if (internal::unpacket_traits<Packet>::masked_load_available) {
     test::packet_helper<internal::unpacket_traits<Packet>::masked_load_available, Packet> h;
     unsigned long long max_umask = (0x1ull << PacketSize);
 
-    for (int offset=0; offset<PacketSize; ++offset)
-    {
-      for (unsigned long long umask=0; umask<max_umask; ++umask)
-      {
-        h.store(data2, h.load(data1+offset, umask));
-        for (int k=0; k<PacketSize; ++k)
-          data3[k] = ((umask & ( 0x1ull << k )) >> k) ? data1[k+offset] : Scalar(0);
+    for (int offset = 0; offset < PacketSize; ++offset) {
+      for (unsigned long long umask = 0; umask < max_umask; ++umask) {
+        h.store(data2, h.load(data1 + offset, umask));
+        for (int k = 0; k < PacketSize; ++k) data3[k] = ((umask & (0x1ull << k)) >> k) ? data1[k + offset] : Scalar(0);
         VERIFY(test::areApprox(data3, data2, PacketSize) && "internal::ploadu masked");
       }
     }
   }
 
-  if (internal::unpacket_traits<Packet>::masked_store_available)
-  {
+  if (internal::unpacket_traits<Packet>::masked_store_available) {
     test::packet_helper<internal::unpacket_traits<Packet>::masked_store_available, Packet> h;
     unsigned long long max_umask = (0x1ull << PacketSize);
 
-    for (int offset=0; offset<PacketSize; ++offset)
-    {
-      for (unsigned long long umask=0; umask<max_umask; ++umask)
-      {
+    for (int offset = 0; offset < PacketSize; ++offset) {
+      for (unsigned long long umask = 0; umask < max_umask; ++umask) {
         internal::pstore(data2, internal::pset1<Packet>(Scalar(0)));
-        h.store(data2, h.loadu(data1+offset), umask);
-        for (int k=0; k<PacketSize; ++k)
-          data3[k] = ((umask & ( 0x1ull << k )) >> k) ? data1[k+offset] : Scalar(0);
+        h.store(data2, h.loadu(data1 + offset), umask);
+        for (int k = 0; k < PacketSize; ++k) data3[k] = ((umask & (0x1ull << k)) >> k) ? data1[k + offset] : Scalar(0);
         VERIFY(test::areApprox(data3, data2, PacketSize) && "internal::pstoreu masked");
       }
     }
@@ -248,110 +427,95 @@ template<typename Scalar,typename Packet> void packetmath()
   VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasSub);
   VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasMul);
 
-  CHECK_CWISE2_IF(PacketTraits::HasAdd, REF_ADD,  internal::padd);
-  CHECK_CWISE2_IF(PacketTraits::HasSub, REF_SUB,  internal::psub);
-  CHECK_CWISE2_IF(PacketTraits::HasMul, REF_MUL,  internal::pmul);
+  CHECK_CWISE2_IF(PacketTraits::HasAdd, REF_ADD, internal::padd);
+  CHECK_CWISE2_IF(PacketTraits::HasSub, REF_SUB, internal::psub);
+  CHECK_CWISE2_IF(PacketTraits::HasMul, REF_MUL, internal::pmul);
   CHECK_CWISE2_IF(PacketTraits::HasDiv, REF_DIV, internal::pdiv);
 
-  if (PacketTraits::HasNegate)
-    CHECK_CWISE1(internal::negate, internal::pnegate);
+  if (PacketTraits::HasNegate) CHECK_CWISE1(internal::negate, internal::pnegate);
   CHECK_CWISE1(numext::conj, internal::pconj);
 
-  for(int offset=0;offset<3;++offset)
-  {
-    for (int i=0; i<PacketSize; ++i)
-      ref[i] = data1[offset];
+  for (int offset = 0; offset < 3; ++offset) {
+    for (int i = 0; i < PacketSize; ++i) ref[i] = data1[offset];
     internal::pstore(data2, internal::pset1<Packet>(data1[offset]));
     VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::pset1");
   }
 
   {
-    for (int i=0; i<PacketSize*4; ++i)
-      ref[i] = data1[i/PacketSize];
+    for (int i = 0; i < PacketSize * 4; ++i) ref[i] = data1[i / PacketSize];
     Packet A0, A1, A2, A3;
     internal::pbroadcast4<Packet>(data1, A0, A1, A2, A3);
-    internal::pstore(data2+0*PacketSize, A0);
-    internal::pstore(data2+1*PacketSize, A1);
-    internal::pstore(data2+2*PacketSize, A2);
-    internal::pstore(data2+3*PacketSize, A3);
-    VERIFY(test::areApprox(ref, data2, 4*PacketSize) && "internal::pbroadcast4");
+    internal::pstore(data2 + 0 * PacketSize, A0);
+    internal::pstore(data2 + 1 * PacketSize, A1);
+    internal::pstore(data2 + 2 * PacketSize, A2);
+    internal::pstore(data2 + 3 * PacketSize, A3);
+    VERIFY(test::areApprox(ref, data2, 4 * PacketSize) && "internal::pbroadcast4");
   }
 
   {
-    for (int i=0; i<PacketSize*2; ++i)
-      ref[i] = data1[i/PacketSize];
+    for (int i = 0; i < PacketSize * 2; ++i) ref[i] = data1[i / PacketSize];
     Packet A0, A1;
     internal::pbroadcast2<Packet>(data1, A0, A1);
-    internal::pstore(data2+0*PacketSize, A0);
-    internal::pstore(data2+1*PacketSize, A1);
-    VERIFY(test::areApprox(ref, data2, 2*PacketSize) && "internal::pbroadcast2");
+    internal::pstore(data2 + 0 * PacketSize, A0);
+    internal::pstore(data2 + 1 * PacketSize, A1);
+    VERIFY(test::areApprox(ref, data2, 2 * PacketSize) && "internal::pbroadcast2");
   }
 
   VERIFY(internal::isApprox(data1[0], internal::pfirst(internal::pload<Packet>(data1))) && "internal::pfirst");
 
-  if(PacketSize>1)
-  {
+  if (PacketSize > 1) {
     // apply different offsets to check that ploaddup is robust to unaligned inputs
-    for(int offset=0;offset<4;++offset)
-    {
-      for(int i=0;i<PacketSize/2;++i)
-        ref[2*i+0] = ref[2*i+1] = data1[offset+i];
-      internal::pstore(data2,internal::ploaddup<Packet>(data1+offset));
+    for (int offset = 0; offset < 4; ++offset) {
+      for (int i = 0; i < PacketSize / 2; ++i) ref[2 * i + 0] = ref[2 * i + 1] = data1[offset + i];
+      internal::pstore(data2, internal::ploaddup<Packet>(data1 + offset));
       VERIFY(test::areApprox(ref, data2, PacketSize) && "ploaddup");
     }
   }
 
-  if(PacketSize>2)
-  {
+  if (PacketSize > 2) {
     // apply different offsets to check that ploadquad is robust to unaligned inputs
-    for(int offset=0;offset<4;++offset)
-    {
-      for(int i=0;i<PacketSize/4;++i)
-        ref[4*i+0] = ref[4*i+1] = ref[4*i+2] = ref[4*i+3] = data1[offset+i];
-      internal::pstore(data2,internal::ploadquad<Packet>(data1+offset));
+    for (int offset = 0; offset < 4; ++offset) {
+      for (int i = 0; i < PacketSize / 4; ++i)
+        ref[4 * i + 0] = ref[4 * i + 1] = ref[4 * i + 2] = ref[4 * i + 3] = data1[offset + i];
+      internal::pstore(data2, internal::ploadquad<Packet>(data1 + offset));
       VERIFY(test::areApprox(ref, data2, PacketSize) && "ploadquad");
     }
   }
 
   ref[0] = Scalar(0);
-  for (int i=0; i<PacketSize; ++i)
-    ref[0] += data1[i];
+  for (int i = 0; i < PacketSize; ++i) ref[0] += data1[i];
   VERIFY(test::isApproxAbs(ref[0], internal::predux(internal::pload<Packet>(data1)), refvalue) && "internal::predux");
 
-  if(PacketSize==8 && internal::unpacket_traits<typename internal::unpacket_traits<Packet>::half>::size ==4) // so far, predux_half_downto4 is only required in such a case
+  if (PacketSize == 8 && internal::unpacket_traits<typename internal::unpacket_traits<Packet>::half>::size ==
+                             4)  // so far, predux_half_downto4 is only required in such a case
   {
-    int HalfPacketSize = PacketSize>4 ? PacketSize/2 : PacketSize;
-    for (int i=0; i<HalfPacketSize; ++i)
-      ref[i] = Scalar(0);
-    for (int i=0; i<PacketSize; ++i)
-      ref[i%HalfPacketSize] += data1[i];
+    int HalfPacketSize = PacketSize > 4 ? PacketSize / 2 : PacketSize;
+    for (int i = 0; i < HalfPacketSize; ++i) ref[i] = Scalar(0);
+    for (int i = 0; i < PacketSize; ++i) ref[i % HalfPacketSize] += data1[i];
     internal::pstore(data2, internal::predux_half_dowto4(internal::pload<Packet>(data1)));
     VERIFY(test::areApprox(ref, data2, HalfPacketSize) && "internal::predux_half_dowto4");
   }
 
   ref[0] = Scalar(1);
-  for (int i=0; i<PacketSize; ++i)
-    ref[0] = REF_MUL(ref[0], data1[i]);
+  for (int i = 0; i < PacketSize; ++i) ref[0] = REF_MUL(ref[0], data1[i]);
   VERIFY(internal::isApprox(ref[0], internal::predux_mul(internal::pload<Packet>(data1))) && "internal::predux_mul");
 
-  for (int i=0; i<PacketSize; ++i)
-    ref[i] = data1[PacketSize-i-1];
+  for (int i = 0; i < PacketSize; ++i) ref[i] = data1[PacketSize - i - 1];
   internal::pstore(data2, internal::preverse(internal::pload<Packet>(data1)));
   VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::preverse");
 
   internal::PacketBlock<Packet> kernel;
-  for (int i=0; i<PacketSize; ++i) {
-    kernel.packet[i] = internal::pload<Packet>(data1+i*PacketSize);
+  for (int i = 0; i < PacketSize; ++i) {
+    kernel.packet[i] = internal::pload<Packet>(data1 + i * PacketSize);
   }
   ptranspose(kernel);
-  for (int i=0; i<PacketSize; ++i) {
+  for (int i = 0; i < PacketSize; ++i) {
     internal::pstore(data2, kernel.packet[i]);
     for (int j = 0; j < PacketSize; ++j) {
-      VERIFY(test::isApproxAbs(data2[j], data1[i+j*PacketSize], refvalue) && "ptranspose");
+      VERIFY(test::isApproxAbs(data2[j], data1[i + j * PacketSize], refvalue) && "ptranspose");
     }
   }
 
-
   if (PacketTraits::HasBlend) {
     Packet thenPacket = internal::pload<Packet>(data1);
     Packet elsePacket = internal::pload<Packet>(data2);
@@ -372,22 +536,21 @@ template<typename Scalar,typename Packet> void packetmath()
     for (int i = 0; i < PacketSize; ++i) {
       // "if" mask
       unsigned char v = internal::random<bool>() ? 0xff : 0;
-      char* bytes = (char*)(data1+i);
-      for(int k=0; k<int(sizeof(Scalar)); ++k) {
+      char* bytes = (char*)(data1 + i);
+      for (int k = 0; k < int(sizeof(Scalar)); ++k) {
         bytes[k] = v;
       }
       // "then" packet
-      data1[i+PacketSize] = internal::random<Scalar>();
+      data1[i + PacketSize] = internal::random<Scalar>();
       // "else" packet
-      data1[i+2*PacketSize] = internal::random<Scalar>();
+      data1[i + 2 * PacketSize] = internal::random<Scalar>();
     }
     CHECK_CWISE3_IF(true, internal::pselect, internal::pselect);
   }
 
   CHECK_CWISE1_IF(PacketTraits::HasSqrt, numext::sqrt, internal::psqrt);
 
-  for (int i=0; i<size; ++i)
-  {
+  for (int i = 0; i < size; ++i) {
     data1[i] = internal::random<Scalar>();
   }
   CHECK_CWISE1(internal::pzero, internal::pzero);
@@ -399,34 +562,29 @@ template<typename Scalar,typename Packet> void packetmath()
   packetmath_pcast_ops<Scalar, Packet>();
 }
 
-
-
-template<typename Scalar,typename Packet> void packetmath_real()
-{
+template <typename Scalar, typename Packet>
+void packetmath_real() {
   typedef internal::packet_traits<Scalar> PacketTraits;
   const int PacketSize = internal::unpacket_traits<Packet>::size;
 
-  const int size = PacketSize*4;
-  EIGEN_ALIGN_MAX Scalar data1[PacketSize*4];
-  EIGEN_ALIGN_MAX Scalar data2[PacketSize*4];
-  EIGEN_ALIGN_MAX Scalar ref[PacketSize*4];
+  const int size = PacketSize * 4;
+  EIGEN_ALIGN_MAX Scalar data1[PacketSize * 4];
+  EIGEN_ALIGN_MAX Scalar data2[PacketSize * 4];
+  EIGEN_ALIGN_MAX Scalar ref[PacketSize * 4];
 
-  for (int i=0; i<size; ++i)
-  {
-    data1[i] = internal::random<Scalar>(0,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6));
-    data2[i] = internal::random<Scalar>(0,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6));
+  for (int i = 0; i < size; ++i) {
+    data1[i] = internal::random<Scalar>(0, 1) * std::pow(Scalar(10), internal::random<Scalar>(-6, 6));
+    data2[i] = internal::random<Scalar>(0, 1) * std::pow(Scalar(10), internal::random<Scalar>(-6, 6));
   }
 
-  if(internal::random<float>(0,1)<0.1f)
-     data1[internal::random<int>(0, PacketSize)] = 0;
+  if (internal::random<float>(0, 1) < 0.1f) data1[internal::random<int>(0, PacketSize)] = 0;
 
   CHECK_CWISE1_IF(PacketTraits::HasLog, std::log, internal::plog);
-  CHECK_CWISE1_IF(PacketTraits::HasRsqrt, Scalar(1)/std::sqrt, internal::prsqrt);
+  CHECK_CWISE1_IF(PacketTraits::HasRsqrt, Scalar(1) / std::sqrt, internal::prsqrt);
 
-  for (int i=0; i<size; ++i)
-  {
-    data1[i] = internal::random<Scalar>(-1,1) * std::pow(Scalar(10), internal::random<Scalar>(-3,3));
-    data2[i] = internal::random<Scalar>(-1,1) * std::pow(Scalar(10), internal::random<Scalar>(-3,3));
+  for (int i = 0; i < size; ++i) {
+    data1[i] = internal::random<Scalar>(-1, 1) * std::pow(Scalar(10), internal::random<Scalar>(-3, 3));
+    data2[i] = internal::random<Scalar>(-1, 1) * std::pow(Scalar(10), internal::random<Scalar>(-3, 3));
   }
   CHECK_CWISE1_IF(PacketTraits::HasSin, std::sin, internal::psin);
   CHECK_CWISE1_IF(PacketTraits::HasCos, std::cos, internal::pcos);
@@ -438,40 +596,35 @@ template<typename Scalar,typename Packet> void packetmath_real()
   CHECK_CWISE1_IF(PacketTraits::HasRint, numext::rint, internal::print);
 
   // See bug 1785.
-  for (int i=0; i<size; ++i)
-   {
-     data1[i] = -1.5 + i;
-     data2[i] = -1.5 + i;
-   }
+  for (int i = 0; i < size; ++i) {
+    data1[i] = -1.5 + i;
+    data2[i] = -1.5 + i;
+  }
   CHECK_CWISE1_IF(PacketTraits::HasRound, numext::round, internal::pround);
   CHECK_CWISE1_IF(PacketTraits::HasRint, numext::rint, internal::print);
 
-  for (int i=0; i<size; ++i)
-  {
-    data1[i] = internal::random<Scalar>(-1,1);
-    data2[i] = internal::random<Scalar>(-1,1);
+  for (int i = 0; i < size; ++i) {
+    data1[i] = internal::random<Scalar>(-1, 1);
+    data2[i] = internal::random<Scalar>(-1, 1);
   }
   CHECK_CWISE1_IF(PacketTraits::HasASin, std::asin, internal::pasin);
   CHECK_CWISE1_IF(PacketTraits::HasACos, std::acos, internal::pacos);
 
-  for (int i=0; i<size; ++i)
-  {
-    data1[i] = internal::random<Scalar>(-87,88);
-    data2[i] = internal::random<Scalar>(-87,88);
+  for (int i = 0; i < size; ++i) {
+    data1[i] = internal::random<Scalar>(-87, 88);
+    data2[i] = internal::random<Scalar>(-87, 88);
   }
   CHECK_CWISE1_IF(PacketTraits::HasExp, std::exp, internal::pexp);
-  for (int i=0; i<size; ++i)
-  {
-    data1[i] = internal::random<Scalar>(-1,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6));
-    data2[i] = internal::random<Scalar>(-1,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6));
+  for (int i = 0; i < size; ++i) {
+    data1[i] = internal::random<Scalar>(-1, 1) * std::pow(Scalar(10), internal::random<Scalar>(-6, 6));
+    data2[i] = internal::random<Scalar>(-1, 1) * std::pow(Scalar(10), internal::random<Scalar>(-6, 6));
   }
   data1[0] = 1e-20;
   CHECK_CWISE1_IF(PacketTraits::HasTanh, std::tanh, internal::ptanh);
-  if(PacketTraits::HasExp && PacketSize>=2)
-  {
+  if (PacketTraits::HasExp && PacketSize >= 2) {
     data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
     data1[1] = std::numeric_limits<Scalar>::epsilon();
-    test::packet_helper<PacketTraits::HasExp,Packet> h;
+    test::packet_helper<PacketTraits::HasExp, Packet> h;
     h.store(data2, internal::pexp(h.load(data1)));
     VERIFY((numext::isnan)(data2[0]));
     VERIFY_IS_EQUAL(std::exp(std::numeric_limits<Scalar>::epsilon()), data2[1]);
@@ -498,24 +651,23 @@ template<typename Scalar,typename Packet> void packetmath_real()
   if (PacketTraits::HasTanh) {
     // NOTE this test migh fail with GCC prior to 6.3, see MathFunctionsImpl.h for details.
     data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
-    test::packet_helper<internal::packet_traits<Scalar>::HasTanh,Packet> h;
+    test::packet_helper<internal::packet_traits<Scalar>::HasTanh, Packet> h;
     h.store(data2, internal::ptanh(h.load(data1)));
     VERIFY((numext::isnan)(data2[0]));
   }
 
   if (PacketTraits::HasExp) {
     internal::scalar_logistic_op<Scalar> logistic;
-    for (int i=0; i<size; ++i)
-    {
-      data1[i] = internal::random<Scalar>(-20,20);
+    for (int i = 0; i < size; ++i) {
+      data1[i] = internal::random<Scalar>(-20, 20);
     }
 
-    test::packet_helper<PacketTraits::HasExp,Packet> h;
+    test::packet_helper<PacketTraits::HasExp, Packet> h;
     h.store(data2, logistic.packetOp(h.load(data1)));
-    for (int i=0; i<PacketSize; ++i) {
-      VERIFY_IS_APPROX(data2[i],logistic(data1[i]));
-#ifdef EIGEN_VECTORIZE // don't check for exactness when using the i387 FPU
-      VERIFY_IS_EQUAL(data2[i],logistic(data1[i]));
+    for (int i = 0; i < PacketSize; ++i) {
+      VERIFY_IS_APPROX(data2[i], logistic(data1[i]));
+#ifdef EIGEN_VECTORIZE  // don't check for exactness when using the i387 FPU
+      VERIFY_IS_EQUAL(data2[i], logistic(data1[i]));
 #endif
     }
   }
@@ -529,13 +681,11 @@ template<typename Scalar,typename Packet> void packetmath_real()
   CHECK_CWISE1_IF(PacketTraits::HasExpm1, std::expm1, internal::pexpm1);
 #endif
 
-  if(PacketSize>=2)
-  {
+  if (PacketSize >= 2) {
     data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
     data1[1] = std::numeric_limits<Scalar>::epsilon();
-    if(PacketTraits::HasLog)
-    {
-      test::packet_helper<PacketTraits::HasLog,Packet> h;
+    if (PacketTraits::HasLog) {
+      test::packet_helper<PacketTraits::HasLog, Packet> h;
       h.store(data2, internal::plog(h.load(data1)));
       VERIFY((numext::isnan)(data2[0]));
       VERIFY_IS_EQUAL(std::log(std::numeric_limits<Scalar>::epsilon()), data2[1]);
@@ -566,45 +716,41 @@ template<typename Scalar,typename Packet> void packetmath_real()
       h.store(data2, internal::plog(h.load(data1)));
       VERIFY((numext::isinf)(data2[0]));
     }
-    if(PacketTraits::HasLog1p) {
-      test::packet_helper<PacketTraits::HasLog1p,Packet> h;
+    if (PacketTraits::HasLog1p) {
+      test::packet_helper<PacketTraits::HasLog1p, Packet> h;
       data1[0] = Scalar(-2);
       data1[1] = -std::numeric_limits<Scalar>::infinity();
       h.store(data2, internal::plog1p(h.load(data1)));
       VERIFY((numext::isnan)(data2[0]));
       VERIFY((numext::isnan)(data2[1]));
     }
-    if(PacketTraits::HasSqrt)
-    {
-      test::packet_helper<PacketTraits::HasSqrt,Packet> h;
+    if (PacketTraits::HasSqrt) {
+      test::packet_helper<PacketTraits::HasSqrt, Packet> h;
       data1[0] = Scalar(-1.0f);
       data1[1] = -std::numeric_limits<Scalar>::denorm_min();
       h.store(data2, internal::psqrt(h.load(data1)));
       VERIFY((numext::isnan)(data2[0]));
       VERIFY((numext::isnan)(data2[1]));
     }
-    if(PacketTraits::HasCos)
-    {
-      test::packet_helper<PacketTraits::HasCos,Packet> h;
-      for(Scalar k = 1; k<Scalar(10000)/std::numeric_limits<Scalar>::epsilon(); k*=2)
-      {
-        for(int k1=0;k1<=1; ++k1)
-        {
-          data1[0] = (2*k+k1  )*Scalar(EIGEN_PI)/2 * internal::random<Scalar>(0.8,1.2);
-          data1[1] = (2*k+2+k1)*Scalar(EIGEN_PI)/2 * internal::random<Scalar>(0.8,1.2);
-          h.store(data2,            internal::pcos(h.load(data1)));
-          h.store(data2+PacketSize, internal::psin(h.load(data1)));
-          VERIFY(data2[0]<=Scalar(1.) && data2[0]>=Scalar(-1.));
-          VERIFY(data2[1]<=Scalar(1.) && data2[1]>=Scalar(-1.));
-          VERIFY(data2[PacketSize+0]<=Scalar(1.) && data2[PacketSize+0]>=Scalar(-1.));
-          VERIFY(data2[PacketSize+1]<=Scalar(1.) && data2[PacketSize+1]>=Scalar(-1.));
-
-          VERIFY_IS_APPROX(numext::abs2(data2[0])+numext::abs2(data2[PacketSize+0]), Scalar(1));
-          VERIFY_IS_APPROX(numext::abs2(data2[1])+numext::abs2(data2[PacketSize+1]), Scalar(1));
+    if (PacketTraits::HasCos) {
+      test::packet_helper<PacketTraits::HasCos, Packet> h;
+      for (Scalar k = 1; k < Scalar(10000) / std::numeric_limits<Scalar>::epsilon(); k *= 2) {
+        for (int k1 = 0; k1 <= 1; ++k1) {
+          data1[0] = (2 * k + k1) * Scalar(EIGEN_PI) / 2 * internal::random<Scalar>(0.8, 1.2);
+          data1[1] = (2 * k + 2 + k1) * Scalar(EIGEN_PI) / 2 * internal::random<Scalar>(0.8, 1.2);
+          h.store(data2, internal::pcos(h.load(data1)));
+          h.store(data2 + PacketSize, internal::psin(h.load(data1)));
+          VERIFY(data2[0] <= Scalar(1.) && data2[0] >= Scalar(-1.));
+          VERIFY(data2[1] <= Scalar(1.) && data2[1] >= Scalar(-1.));
+          VERIFY(data2[PacketSize + 0] <= Scalar(1.) && data2[PacketSize + 0] >= Scalar(-1.));
+          VERIFY(data2[PacketSize + 1] <= Scalar(1.) && data2[PacketSize + 1] >= Scalar(-1.));
+
+          VERIFY_IS_APPROX(numext::abs2(data2[0]) + numext::abs2(data2[PacketSize + 0]), Scalar(1));
+          VERIFY_IS_APPROX(numext::abs2(data2[1]) + numext::abs2(data2[PacketSize + 1]), Scalar(1));
         }
       }
 
-      data1[0] =  std::numeric_limits<Scalar>::infinity();
+      data1[0] = std::numeric_limits<Scalar>::infinity();
       data1[1] = -std::numeric_limits<Scalar>::infinity();
       h.store(data2, internal::psin(h.load(data1)));
       VERIFY((numext::isnan)(data2[0]));
@@ -614,7 +760,7 @@ template<typename Scalar,typename Packet> void packetmath_real()
       VERIFY((numext::isnan)(data2[0]));
       VERIFY((numext::isnan)(data2[1]));
 
-      data1[0] =  std::numeric_limits<Scalar>::quiet_NaN();
+      data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
       h.store(data2, internal::psin(h.load(data1)));
       VERIFY((numext::isnan)(data2[0]));
       h.store(data2, internal::pcos(h.load(data1)));
@@ -622,27 +768,26 @@ template<typename Scalar,typename Packet> void packetmath_real()
 
       data1[0] = -Scalar(0.);
       h.store(data2, internal::psin(h.load(data1)));
-      VERIFY( internal::biteq(data2[0], data1[0]) );
+      VERIFY(internal::biteq(data2[0], data1[0]));
       h.store(data2, internal::pcos(h.load(data1)));
       VERIFY_IS_EQUAL(data2[0], Scalar(1));
     }
   }
 }
 
-template<typename Scalar,typename Packet> void packetmath_notcomplex()
-{
+template <typename Scalar, typename Packet>
+void packetmath_notcomplex() {
   typedef internal::packet_traits<Scalar> PacketTraits;
   const int PacketSize = internal::unpacket_traits<Packet>::size;
 
-  EIGEN_ALIGN_MAX Scalar data1[PacketSize*4];
-  EIGEN_ALIGN_MAX Scalar data2[PacketSize*4];
-  EIGEN_ALIGN_MAX Scalar ref[PacketSize*4];
+  EIGEN_ALIGN_MAX Scalar data1[PacketSize * 4];
+  EIGEN_ALIGN_MAX Scalar data2[PacketSize * 4];
+  EIGEN_ALIGN_MAX Scalar ref[PacketSize * 4];
 
-  Array<Scalar,Dynamic,1>::Map(data1, PacketSize*4).setRandom();
+  Array<Scalar, Dynamic, 1>::Map(data1, PacketSize * 4).setRandom();
 
   ref[0] = data1[0];
-  for (int i=0; i<PacketSize; ++i)
-    ref[0] = (std::min)(ref[0],data1[i]);
+  for (int i = 0; i < PacketSize; ++i) ref[0] = (std::min)(ref[0], data1[i]);
   VERIFY(internal::isApprox(ref[0], internal::predux_min(internal::pload<Packet>(data1))) && "internal::predux_min");
 
   VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasMin);
@@ -654,12 +799,10 @@ template<typename Scalar,typename Packet> void packetmath_notcomplex()
   CHECK_CWISE2_IF(PacketTraits::HasAbsDiff, REF_ABS_DIFF, internal::pabsdiff);
 
   ref[0] = data1[0];
-  for (int i=0; i<PacketSize; ++i)
-    ref[0] = (std::max)(ref[0],data1[i]);
+  for (int i = 0; i < PacketSize; ++i) ref[0] = (std::max)(ref[0], data1[i]);
   VERIFY(internal::isApprox(ref[0], internal::predux_max(internal::pload<Packet>(data1))) && "internal::predux_max");
 
-  for (int i=0; i<PacketSize; ++i)
-    ref[i] = data1[0]+Scalar(i);
+  for (int i = 0; i < PacketSize; ++i) ref[i] = data1[0] + Scalar(i);
   internal::pstore(data2, internal::plset<Packet>(data1[0]));
   VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::plset");
 
@@ -676,167 +819,158 @@ template<typename Scalar,typename Packet> void packetmath_notcomplex()
     // }
 
     // predux_any
-    for (unsigned int i=0; i<PacketSize*sizeof(Scalar); ++i) data1_bits[i] = 0x0;
-    VERIFY( (!internal::predux_any(internal::pload<Packet>(data1))) && "internal::predux_any(0000)");
-    for(int k=0; k<PacketSize; ++k)
-    {
-      for (unsigned int i=0; i<sizeof(Scalar); ++i) data1_bits[k*sizeof(Scalar)+i] = 0xff;
-      VERIFY( internal::predux_any(internal::pload<Packet>(data1)) && "internal::predux_any(0101)");
-      for (unsigned int i=0; i<sizeof(Scalar); ++i) data1_bits[k*sizeof(Scalar)+i] = 0x00;
+    for (unsigned int i = 0; i < PacketSize * sizeof(Scalar); ++i) data1_bits[i] = 0x0;
+    VERIFY((!internal::predux_any(internal::pload<Packet>(data1))) && "internal::predux_any(0000)");
+    for (int k = 0; k < PacketSize; ++k) {
+      for (unsigned int i = 0; i < sizeof(Scalar); ++i) data1_bits[k * sizeof(Scalar) + i] = 0xff;
+      VERIFY(internal::predux_any(internal::pload<Packet>(data1)) && "internal::predux_any(0101)");
+      for (unsigned int i = 0; i < sizeof(Scalar); ++i) data1_bits[k * sizeof(Scalar) + i] = 0x00;
     }
   }
 }
 
-template<typename Scalar,typename Packet,bool ConjLhs,bool ConjRhs> void test_conj_helper(Scalar* data1, Scalar* data2, Scalar* ref, Scalar* pval)
-{
+template <typename Scalar, typename Packet, bool ConjLhs, bool ConjRhs>
+void test_conj_helper(Scalar* data1, Scalar* data2, Scalar* ref, Scalar* pval) {
   const int PacketSize = internal::unpacket_traits<Packet>::size;
 
   internal::conj_if<ConjLhs> cj0;
   internal::conj_if<ConjRhs> cj1;
-  internal::conj_helper<Scalar,Scalar,ConjLhs,ConjRhs> cj;
-  internal::conj_helper<Packet,Packet,ConjLhs,ConjRhs> pcj;
+  internal::conj_helper<Scalar, Scalar, ConjLhs, ConjRhs> cj;
+  internal::conj_helper<Packet, Packet, ConjLhs, ConjRhs> pcj;
 
-  for(int i=0;i<PacketSize;++i)
-  {
+  for (int i = 0; i < PacketSize; ++i) {
     ref[i] = cj0(data1[i]) * cj1(data2[i]);
-    VERIFY(internal::isApprox(ref[i], cj.pmul(data1[i],data2[i])) && "conj_helper pmul");
+    VERIFY(internal::isApprox(ref[i], cj.pmul(data1[i], data2[i])) && "conj_helper pmul");
   }
-  internal::pstore(pval,pcj.pmul(internal::pload<Packet>(data1),internal::pload<Packet>(data2)));
+  internal::pstore(pval, pcj.pmul(internal::pload<Packet>(data1), internal::pload<Packet>(data2)));
   VERIFY(test::areApprox(ref, pval, PacketSize) && "conj_helper pmul");
 
-  for(int i=0;i<PacketSize;++i)
-  {
+  for (int i = 0; i < PacketSize; ++i) {
     Scalar tmp = ref[i];
     ref[i] += cj0(data1[i]) * cj1(data2[i]);
-    VERIFY(internal::isApprox(ref[i], cj.pmadd(data1[i],data2[i],tmp)) && "conj_helper pmadd");
+    VERIFY(internal::isApprox(ref[i], cj.pmadd(data1[i], data2[i], tmp)) && "conj_helper pmadd");
   }
-  internal::pstore(pval,pcj.pmadd(internal::pload<Packet>(data1),internal::pload<Packet>(data2),internal::pload<Packet>(pval)));
+  internal::pstore(
+      pval, pcj.pmadd(internal::pload<Packet>(data1), internal::pload<Packet>(data2), internal::pload<Packet>(pval)));
   VERIFY(test::areApprox(ref, pval, PacketSize) && "conj_helper pmadd");
 }
 
-template<typename Scalar,typename Packet> void packetmath_complex()
-{
+template <typename Scalar, typename Packet>
+void packetmath_complex() {
   const int PacketSize = internal::unpacket_traits<Packet>::size;
 
-  const int size = PacketSize*4;
-  EIGEN_ALIGN_MAX Scalar data1[PacketSize*4];
-  EIGEN_ALIGN_MAX Scalar data2[PacketSize*4];
-  EIGEN_ALIGN_MAX Scalar ref[PacketSize*4];
-  EIGEN_ALIGN_MAX Scalar pval[PacketSize*4];
+  const int size = PacketSize * 4;
+  EIGEN_ALIGN_MAX Scalar data1[PacketSize * 4];
+  EIGEN_ALIGN_MAX Scalar data2[PacketSize * 4];
+  EIGEN_ALIGN_MAX Scalar ref[PacketSize * 4];
+  EIGEN_ALIGN_MAX Scalar pval[PacketSize * 4];
 
-  for (int i=0; i<size; ++i)
-  {
+  for (int i = 0; i < size; ++i) {
     data1[i] = internal::random<Scalar>() * Scalar(1e2);
     data2[i] = internal::random<Scalar>() * Scalar(1e2);
   }
 
-  test_conj_helper<Scalar,Packet,false,false> (data1,data2,ref,pval);
-  test_conj_helper<Scalar,Packet,false,true>  (data1,data2,ref,pval);
-  test_conj_helper<Scalar,Packet,true,false>  (data1,data2,ref,pval);
-  test_conj_helper<Scalar,Packet,true,true>   (data1,data2,ref,pval);
+  test_conj_helper<Scalar, Packet, false, false>(data1, data2, ref, pval);
+  test_conj_helper<Scalar, Packet, false, true>(data1, data2, ref, pval);
+  test_conj_helper<Scalar, Packet, true, false>(data1, data2, ref, pval);
+  test_conj_helper<Scalar, Packet, true, true>(data1, data2, ref, pval);
 
   {
-    for(int i=0;i<PacketSize;++i)
-      ref[i] = Scalar(std::imag(data1[i]),std::real(data1[i]));
-    internal::pstore(pval,internal::pcplxflip(internal::pload<Packet>(data1)));
+    for (int i = 0; i < PacketSize; ++i) ref[i] = Scalar(std::imag(data1[i]), std::real(data1[i]));
+    internal::pstore(pval, internal::pcplxflip(internal::pload<Packet>(data1)));
     VERIFY(test::areApprox(ref, pval, PacketSize) && "pcplxflip");
   }
 }
 
-template<typename Scalar,typename Packet> void packetmath_scatter_gather()
-{
+template <typename Scalar, typename Packet>
+void packetmath_scatter_gather() {
   typedef typename NumTraits<Scalar>::Real RealScalar;
   const int PacketSize = internal::unpacket_traits<Packet>::size;
   EIGEN_ALIGN_MAX Scalar data1[PacketSize];
   RealScalar refvalue = 0;
-  for (int i=0; i<PacketSize; ++i) {
-    data1[i] = internal::random<Scalar>()/RealScalar(PacketSize);
+  for (int i = 0; i < PacketSize; ++i) {
+    data1[i] = internal::random<Scalar>() / RealScalar(PacketSize);
   }
 
-  int stride = internal::random<int>(1,20);
+  int stride = internal::random<int>(1, 20);
 
-  EIGEN_ALIGN_MAX Scalar buffer[PacketSize*20];
-  memset(buffer, 0, 20*PacketSize*sizeof(Scalar));
+  EIGEN_ALIGN_MAX Scalar buffer[PacketSize * 20];
+  memset(buffer, 0, 20 * PacketSize * sizeof(Scalar));
   Packet packet = internal::pload<Packet>(data1);
   internal::pscatter<Scalar, Packet>(buffer, packet, stride);
 
-  for (int i = 0; i < PacketSize*20; ++i) {
-    if ((i%stride) == 0 && i<stride*PacketSize) {
-      VERIFY(
-          test::isApproxAbs(buffer[i], data1[i/stride], refvalue) && "pscatter");
+  for (int i = 0; i < PacketSize * 20; ++i) {
+    if ((i % stride) == 0 && i < stride * PacketSize) {
+      VERIFY(test::isApproxAbs(buffer[i], data1[i / stride], refvalue) && "pscatter");
     } else {
-      VERIFY(
-          test::isApproxAbs(buffer[i], Scalar(0), refvalue) && "pscatter");
+      VERIFY(test::isApproxAbs(buffer[i], Scalar(0), refvalue) && "pscatter");
     }
   }
 
-  for (int i=0; i<PacketSize*7; ++i) {
-    buffer[i] = internal::random<Scalar>()/RealScalar(PacketSize);
+  for (int i = 0; i < PacketSize * 7; ++i) {
+    buffer[i] = internal::random<Scalar>() / RealScalar(PacketSize);
   }
   packet = internal::pgather<Scalar, Packet>(buffer, 7);
   internal::pstore(data1, packet);
   for (int i = 0; i < PacketSize; ++i) {
-    VERIFY(test::isApproxAbs(data1[i], buffer[i*7], refvalue) && "pgather");
+    VERIFY(test::isApproxAbs(data1[i], buffer[i * 7], refvalue) && "pgather");
   }
 }
 
 namespace Eigen {
 namespace test {
 
-template<typename Scalar,typename PacketType>
-struct runall<Scalar,PacketType,false,false> { // i.e. float or double
+template <typename Scalar, typename PacketType>
+struct runall<Scalar, PacketType, false, false> {  // i.e. float or double
   static void run() {
-    packetmath<Scalar,PacketType>();
-    packetmath_scatter_gather<Scalar,PacketType>();
-    packetmath_notcomplex<Scalar,PacketType>();
-    packetmath_real<Scalar,PacketType>();
+    packetmath<Scalar, PacketType>();
+    packetmath_scatter_gather<Scalar, PacketType>();
+    packetmath_notcomplex<Scalar, PacketType>();
+    packetmath_real<Scalar, PacketType>();
   }
 };
 
-template<typename Scalar,typename PacketType>
-struct runall<Scalar,PacketType,false,true> { // i.e. int
+template <typename Scalar, typename PacketType>
+struct runall<Scalar, PacketType, false, true> {  // i.e. int
   static void run() {
-    packetmath<Scalar,PacketType>();
-    packetmath_scatter_gather<Scalar,PacketType>();
-    packetmath_notcomplex<Scalar,PacketType>();
+    packetmath<Scalar, PacketType>();
+    packetmath_scatter_gather<Scalar, PacketType>();
+    packetmath_notcomplex<Scalar, PacketType>();
   }
 };
 
-template<typename Scalar,typename PacketType>
-struct runall<Scalar,PacketType,true,false> { // i.e. complex
+template <typename Scalar, typename PacketType>
+struct runall<Scalar, PacketType, true, false> {  // i.e. complex
   static void run() {
-    packetmath<Scalar,PacketType>();
-    packetmath_scatter_gather<Scalar,PacketType>();
-    packetmath_complex<Scalar,PacketType>();
+    packetmath<Scalar, PacketType>();
+    packetmath_scatter_gather<Scalar, PacketType>();
+    packetmath_complex<Scalar, PacketType>();
   }
 };
 
-}
-}
-
+}  // namespace test
+}  // namespace Eigen
 
-EIGEN_DECLARE_TEST(packetmath)
-{
+EIGEN_DECLARE_TEST(packetmath) {
   g_first_pass = true;
-  for(int i = 0; i < g_repeat; i++) {
-
-    CALL_SUBTEST_1( test::runner<float>::run() );
-    CALL_SUBTEST_2( test::runner<double>::run() );
-    CALL_SUBTEST_3( test::runner<int8_t>::run() );
-    CALL_SUBTEST_4( test::runner<uint8_t>::run() );
-    CALL_SUBTEST_5( test::runner<int16_t>::run() );
-    CALL_SUBTEST_6( test::runner<uint16_t>::run() );
-    CALL_SUBTEST_7( test::runner<int32_t>::run() );
-    CALL_SUBTEST_8( test::runner<uint32_t>::run() );
-    CALL_SUBTEST_9( test::runner<int64_t>::run() );
-    CALL_SUBTEST_10( test::runner<uint64_t>::run() );
-    CALL_SUBTEST_11( test::runner<std::complex<float> >::run() );
-    CALL_SUBTEST_12( test::runner<std::complex<double> >::run() );
-    CALL_SUBTEST_13(( packetmath<half,internal::packet_traits<half>::type>() ));
+  for (int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST_1(test::runner<float>::run());
+    CALL_SUBTEST_2(test::runner<double>::run());
+    CALL_SUBTEST_3(test::runner<int8_t>::run());
+    CALL_SUBTEST_4(test::runner<uint8_t>::run());
+    CALL_SUBTEST_5(test::runner<int16_t>::run());
+    CALL_SUBTEST_6(test::runner<uint16_t>::run());
+    CALL_SUBTEST_7(test::runner<int32_t>::run());
+    CALL_SUBTEST_8(test::runner<uint32_t>::run());
+    CALL_SUBTEST_9(test::runner<int64_t>::run());
+    CALL_SUBTEST_10(test::runner<uint64_t>::run());
+    CALL_SUBTEST_11(test::runner<std::complex<float> >::run());
+    CALL_SUBTEST_12(test::runner<std::complex<double> >::run());
+    CALL_SUBTEST_13((packetmath<half, internal::packet_traits<half>::type>()));
 #ifdef EIGEN_PACKET_MATH_SSE_H
-    CALL_SUBTEST_14(( packetmath<bool,internal::packet_traits<bool>::type>() ));
+    CALL_SUBTEST_14((packetmath<bool, internal::packet_traits<bool>::type>()));
 #endif
-    CALL_SUBTEST_15(( packetmath<bfloat16,internal::packet_traits<bfloat16>::type>() ));
+    CALL_SUBTEST_15((packetmath<bfloat16, internal::packet_traits<bfloat16>::type>()));
     g_first_pass = false;
   }
 }