From 94b19dc5f24f8cc649be7deb849c5ad38e4ce60e Mon Sep 17 00:00:00 2001 From: Charles Schlosser Date: Wed, 15 Feb 2023 21:33:06 +0000 Subject: [PATCH] Add CArg --- Eigen/src/Core/GenericPacketMath.h | 44 +++++++++++++++++-------- Eigen/src/Core/GlobalFunctions.h | 1 + Eigen/src/Core/functors/UnaryFunctors.h | 22 +++++++++++++ Eigen/src/plugins/ArrayCwiseUnaryOps.h | 5 +++ Eigen/src/plugins/MatrixCwiseUnaryOps.h | 5 +++ test/array_cwise.cpp | 3 ++ 6 files changed, 67 insertions(+), 13 deletions(-) diff --git a/Eigen/src/Core/GenericPacketMath.h b/Eigen/src/Core/GenericPacketMath.h index d65c63a10..499a94b22 100644 --- a/Eigen/src/Core/GenericPacketMath.h +++ b/Eigen/src/Core/GenericPacketMath.h @@ -1192,27 +1192,27 @@ Packet prsqrt(const Packet& a) { } template ::value, - bool IsInteger = NumTraits::type>::IsInteger> - struct psignbit_impl; + bool IsInteger = NumTraits::type>::IsInteger> +struct psignbit_impl; template struct psignbit_impl { - EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE static constexpr Packet run(const Packet& a) { return numext::signbit(a); } + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE static constexpr Packet run(const Packet& a) { return numext::signbit(a); } }; template struct psignbit_impl { - // generic implementation if not specialized in PacketMath.h - // slower than arithmetic shift - typedef typename unpacket_traits::type Scalar; - EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE static Packet run(const Packet& a) { - const Packet cst_pos_one = pset1(Scalar(1)); - const Packet cst_neg_one = pset1(Scalar(-1)); - return pcmp_eq(por(pand(a, cst_neg_one), cst_pos_one), cst_neg_one); - } + // generic implementation if not specialized in PacketMath.h + // slower than arithmetic shift + typedef typename unpacket_traits::type Scalar; + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE static Packet run(const Packet& a) { + const Packet cst_pos_one = pset1(Scalar(1)); + const Packet cst_neg_one = pset1(Scalar(-1)); + return pcmp_eq(por(pand(a, cst_neg_one), cst_pos_one), cst_neg_one); + } }; template struct psignbit_impl { - // generic implementation for integer packets - EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE static constexpr Packet run(const Packet& a) { return pcmp_lt(a, pzero(a)); } + // generic implementation for integer packets + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE static constexpr Packet run(const Packet& a) { return pcmp_lt(a, pzero(a)); } }; /** \internal \returns the sign bit of \a a as a bitmask*/ template @@ -1256,6 +1256,24 @@ EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet patan2(const Packet& y, const Packe return result; } +/** \internal \returns the argument of \a a as a complex number */ +template ::value, int> = 0> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet pcarg(const Packet& a) { + return Packet(numext::arg(a)); +} + +/** \internal \returns the argument of \a a as a complex number */ +template ::value, int> = 0> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet pcarg(const Packet& a) { + using Scalar = typename unpacket_traits::type; + EIGEN_STATIC_ASSERT(NumTraits::IsComplex, THIS METHOD IS FOR COMPLEX TYPES ONLY) + using RealPacket = typename unpacket_traits::as_real; + // a // r i r i ... + RealPacket aflip = pcplxflip(a).v; // i r i r ... + RealPacket result = patan2(aflip, a.v); // atan2 crap atan2 crap ... + return (Packet)pand(result, peven_mask(result)); // atan2 0 atan2 0 ... +} + } // end namespace internal } // end namespace Eigen diff --git a/Eigen/src/Core/GlobalFunctions.h b/Eigen/src/Core/GlobalFunctions.h index 18792cb12..30f2a38af 100644 --- a/Eigen/src/Core/GlobalFunctions.h +++ b/Eigen/src/Core/GlobalFunctions.h @@ -86,6 +86,7 @@ namespace Eigen EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(abs,scalar_abs_op,absolute value,\sa ArrayBase::abs DOXCOMMA MatrixBase::cwiseAbs) EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(abs2,scalar_abs2_op,squared absolute value,\sa ArrayBase::abs2 DOXCOMMA MatrixBase::cwiseAbs2) EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(arg,scalar_arg_op,complex argument,\sa ArrayBase::arg DOXCOMMA MatrixBase::cwiseArg) + EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(carg, scalar_carg_op, complex argument, \sa ArrayBase::carg DOXCOMMA MatrixBase::cwiseCArg) EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(sqrt,scalar_sqrt_op,square root,\sa ArrayBase::sqrt DOXCOMMA MatrixBase::cwiseSqrt) EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(rsqrt,scalar_rsqrt_op,reciprocal square root,\sa ArrayBase::rsqrt) EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(square,scalar_square_op,square (power 2),\sa Eigen::abs2 DOXCOMMA Eigen::pow DOXCOMMA ArrayBase::square) diff --git a/Eigen/src/Core/functors/UnaryFunctors.h b/Eigen/src/Core/functors/UnaryFunctors.h index 348536968..a76d8d264 100644 --- a/Eigen/src/Core/functors/UnaryFunctors.h +++ b/Eigen/src/Core/functors/UnaryFunctors.h @@ -146,6 +146,28 @@ struct functor_traits > PacketAccess = packet_traits::HasArg }; }; + +/** \internal + * \brief Template functor to compute the complex argument, returned as a complex type + * + * \sa class CwiseUnaryOp, Cwise::carg + */ +template +struct scalar_carg_op { + using result_type = Scalar; + EIGEN_DEVICE_FUNC + EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& a) const { return Scalar(numext::arg(a)); } + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const { + return pcarg(a); + } +}; +template +struct functor_traits> { + using RealScalar = typename NumTraits::Real; + enum { Cost = functor_traits>::Cost, PacketAccess = packet_traits::HasATan }; +}; + /** \internal * \brief Template functor to cast a scalar to another type * diff --git a/Eigen/src/plugins/ArrayCwiseUnaryOps.h b/Eigen/src/plugins/ArrayCwiseUnaryOps.h index d8c1a849b..807a31325 100644 --- a/Eigen/src/plugins/ArrayCwiseUnaryOps.h +++ b/Eigen/src/plugins/ArrayCwiseUnaryOps.h @@ -2,6 +2,7 @@ typedef CwiseUnaryOp, const Derived> AbsReturnType; typedef CwiseUnaryOp, const Derived> ArgReturnType; +typedef CwiseUnaryOp, const Derived> CArgReturnType; typedef CwiseUnaryOp, const Derived> Abs2ReturnType; typedef CwiseUnaryOp, const Derived> SqrtReturnType; typedef CwiseUnaryOp, const Derived> RsqrtReturnType; @@ -66,6 +67,10 @@ arg() const return ArgReturnType(derived()); } +EIGEN_DEVICE_FUNC +EIGEN_STRONG_INLINE const CArgReturnType +carg() const { return CArgReturnType(derived()); } + /** \returns an expression of the coefficient-wise squared absolute value of \c *this * * Example: \include Cwise_abs2.cpp diff --git a/Eigen/src/plugins/MatrixCwiseUnaryOps.h b/Eigen/src/plugins/MatrixCwiseUnaryOps.h index 98d925dd2..0222137ae 100644 --- a/Eigen/src/plugins/MatrixCwiseUnaryOps.h +++ b/Eigen/src/plugins/MatrixCwiseUnaryOps.h @@ -15,6 +15,7 @@ typedef CwiseUnaryOp, const Derived> CwiseAbsReturnType; typedef CwiseUnaryOp, const Derived> CwiseAbs2ReturnType; typedef CwiseUnaryOp, const Derived> CwiseArgReturnType; +typedef CwiseUnaryOp, const Derived> CwiseCArgReturnType; typedef CwiseUnaryOp, const Derived> CwiseSqrtReturnType; typedef CwiseUnaryOp, const Derived> CwiseSignReturnType; typedef CwiseUnaryOp, const Derived> CwiseInverseReturnType; @@ -94,6 +95,10 @@ EIGEN_DEVICE_FUNC inline const CwiseArgReturnType cwiseArg() const { return CwiseArgReturnType(derived()); } +EIGEN_DEVICE_FUNC +EIGEN_STRONG_INLINE const CwiseCArgReturnType +cwiseCArg() const { return CwiseCArgReturnType(derived()); } + template using CwisePowReturnType = std::enable_if_t::Real>::value, diff --git a/test/array_cwise.cpp b/test/array_cwise.cpp index 451b1ccb6..6cd6469b4 100644 --- a/test/array_cwise.cpp +++ b/test/array_cwise.cpp @@ -773,6 +773,8 @@ template void array_complex(const ArrayType& m) VERIFY_IS_APPROX(m1.tanh(), tanh(m1)); VERIFY_IS_APPROX(m1.logistic(), logistic(m1)); VERIFY_IS_APPROX(m1.arg(), arg(m1)); + VERIFY_IS_APPROX(m1.carg(), carg(m1)); + VERIFY_IS_APPROX(arg(m1), carg(m1)); VERIFY((m1.isNaN() == (Eigen::isnan)(m1)).all()); VERIFY((m1.isInf() == (Eigen::isinf)(m1)).all()); VERIFY((m1.isFinite() == (Eigen::isfinite)(m1)).all()); @@ -806,6 +808,7 @@ template void array_complex(const ArrayType& m) for (Index j = 0; j < m.cols(); ++j) m3(i,j) = std::atan2(m1(i,j).imag(), m1(i,j).real()); VERIFY_IS_APPROX(arg(m1), m3); + VERIFY_IS_APPROX(carg(m1), m3); std::complex zero(0.0,0.0); VERIFY((Eigen::isnan)(m1*zero/zero).all());