diff options
Diffstat (limited to 'Eigen/src/Core/arch/Default')
-rw-r--r-- | Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h | 20 | ||||
-rw-r--r-- | Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h | 6 |
2 files changed, 26 insertions, 0 deletions
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h index c9fbaf68b..f1e10c898 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h @@ -757,6 +757,26 @@ Packet pcos_float(const Packet& x) return psincos_float<false>(x); } +template<typename Packet> +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS +EIGEN_UNUSED Packet pdiv_complex(const Packet& x, const Packet& y) { + typedef typename unpacket_traits<Packet>::as_real RealPacket; + // In the following we annotate the code for the case where the inputs + // are a pair length-2 SIMD vectors representing a single pair of complex + // numbers x = a + i*b, y = c + i*d. + const RealPacket y_abs = pabs(y.v); // |c|, |d| + const RealPacket y_abs_flip = pcplxflip(Packet(y_abs)).v; // |d|, |c| + const RealPacket y_max = pmax(y_abs, y_abs_flip); // max(|c|, |d|), max(|c|, |d|) + const RealPacket y_scaled = pdiv(y.v, y_max); // c / max(|c|, |d|), d / max(|c|, |d|) + // Compute scaled denominator. + const RealPacket y_scaled_sq = pmul(y_scaled, y_scaled); // c'**2, d'**2 + const RealPacket denom = y_scaled_sq + pcplxflip(Packet(y_scaled_sq)).v; + Packet result_scaled = pmul(x, pconj(Packet(y_scaled))); // a * c' + b * d', -a * d + b * c + // Divide elementwise by denom. + result_scaled = Packet(pdiv(result_scaled.v, denom)); + // Rescale result + return Packet(pdiv(result_scaled.v, y_max)); +} template<typename Packet> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h index 177a04e93..730cc7395 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h @@ -101,6 +101,12 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet psqrt_complex(const Packet& a); +/** \internal \returns x / y for complex types */ +template<typename Packet> +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS +EIGEN_UNUSED +Packet pdiv_complex(const Packet& x, const Packet& y); + template <typename Packet, int N> struct ppolevl; |