Factories code between numext::hypot and scalar_hyot_op functor.

(grafted from 4213b63f5ce33d3f904674ee7b0cabd6934dda6b
)
This commit is contained in:
Gael Guennebaud 2018-04-04 15:12:43 +02:00
parent b18e2d422b
commit 90cd199d4b
2 changed files with 21 additions and 24 deletions

View File

@ -71,6 +71,17 @@ T generic_fast_tanh_float(const T& a_x)
return pdiv(p, q); return pdiv(p, q);
} }
template<typename RealScalar>
EIGEN_STRONG_INLINE
RealScalar positive_real_hypot(const RealScalar& x, const RealScalar& y)
{
EIGEN_USING_STD_MATH(sqrt);
RealScalar p, qp;
p = numext::maxi(x,y);
if(p==RealScalar(0)) return RealScalar(0);
qp = numext::mini(y,x) / p;
return p * sqrt(RealScalar(1) + qp*qp);
}
template<typename Scalar> template<typename Scalar>
struct hypot_impl struct hypot_impl
@ -79,14 +90,7 @@ struct hypot_impl
static inline RealScalar run(const Scalar& x, const Scalar& y) static inline RealScalar run(const Scalar& x, const Scalar& y)
{ {
EIGEN_USING_STD_MATH(abs); EIGEN_USING_STD_MATH(abs);
EIGEN_USING_STD_MATH(sqrt); return positive_real_hypot(abs(x), abs(y));
RealScalar _x = abs(x);
RealScalar _y = abs(y);
RealScalar p, qp;
p = numext::maxi(_x,_y);
if(p==RealScalar(0)) return RealScalar(0);
qp = numext::mini(_y,_x) / p;
return p * sqrt(RealScalar(1) + qp*qp);
} }
}; };

View File

@ -255,7 +255,7 @@ struct scalar_cmp_op<LhsScalar,RhsScalar, cmp_NEQ> : binary_op_base<LhsScalar,Rh
/** \internal /** \internal
* \brief Template functor to compute the hypot of two scalars * \brief Template functor to compute the hypot of two \b positive \b and \b real scalars
* *
* \sa MatrixBase::stableNorm(), class Redux * \sa MatrixBase::stableNorm(), class Redux
*/ */
@ -263,22 +263,15 @@ template<typename Scalar>
struct scalar_hypot_op<Scalar,Scalar> : binary_op_base<Scalar,Scalar> struct scalar_hypot_op<Scalar,Scalar> : binary_op_base<Scalar,Scalar>
{ {
EIGEN_EMPTY_STRUCT_CTOR(scalar_hypot_op) EIGEN_EMPTY_STRUCT_CTOR(scalar_hypot_op)
// typedef typename NumTraits<Scalar>::Real result_type;
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& _x, const Scalar& _y) const EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar &x, const Scalar &y) const
{ {
EIGEN_USING_STD_MATH(sqrt) // This functor is used by hypotNorm only for which it is faster to first apply abs
Scalar p, qp; // on all coefficients prior to reduction through hypot.
if(_x>_y) // This way we avoid calling abs on positive and real entries, and this also permits
{ // to seamlessly handle complexes. Otherwise we would have to handle both real and complexes
p = _x; // through the same functor...
qp = _y / p; return internal::positive_real_hypot(x,y);
}
else
{
p = _y;
qp = _x / p;
}
return p * sqrt(Scalar(1) + qp*qp);
} }
}; };
template<typename Scalar> template<typename Scalar>