GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 11:49:02 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #876: implement a portable log1p function

Browse Source
This commit is contained in:
Gael Guennebaud 2014-12-08 16:26:53 +01:00
parent 7f7a712062
commit bea36925db
1 changed files with 52 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

85
Eigen/src/Core/MathFunctions.h
View File

@ -14,7 +14,7 @@ namespace Eigen {
// On WINCE, std::abs is defined for int only, so let's defined our own overloads:
// This issue has been confirmed with MSVC 2008 only, but the issue might exist for more recent versions too.
#if defined(_WIN32_WCE) && defined(_MSC_VER) && _MSC_VER<=1500
#if EIGEN_OS_WINCE && EIGEN_COMP_MSVC && EIGEN_COMP_MSVC<=1500
long abs(long x) { return (labs(x)); }
double abs(double x) { return (fabs(x)); }
float abs(float x) { return (fabsf(x)); }
@ -360,50 +360,31 @@ inline NewType cast(const OldType& x)
}
/****************************************************************************
* Implementation of atanh2 *
* Implementation of logp1 *
****************************************************************************/
template<typename Scalar>
struct atanh2_impl
struct log1p_impl
{
static inline Scalar run(const Scalar& x, const Scalar& r)
static inline Scalar run(const Scalar& x)
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
#if (__cplusplus >= 201103L) && !defined(__CYGWIN__)
// Let's be conservative and enable the default C++11 implementation only if we are sure it exists
#if (__cplusplus >= 201103L) && (EIGEN_COMP_GNUC_STRICT || EIGEN_COMP_CLANG || EIGEN_COMP_MSVC || EIGEN_COMP_ICC) \
&& (EIGEN_ARCH_i386_OR_x86_64) && (EIGEN_OS_GNULINUX || EIGEN_OS_WIN_STRICT || EIGEN_OS_MAC)
using std::log1p;
return log1p(2 * x / (r - x)) / 2;
return log1p(x);
#else
using std::abs;
typedef typename NumTraits<Scalar>::Real RealScalar;
using std::log;
using std::sqrt;
Scalar z = x / r;
if (r == 0 || abs(z) > sqrt(NumTraits<Scalar>::epsilon()))
return log((r + x) / (r - x)) / 2;
else
return z + z*z*z / 3;
Scalar x1p = RealScalar(1) + x;
return ( x1p == Scalar(1) ) ? x : x * ( log(x1p) / (x1p - RealScalar(1)) );
#endif
}
};
template<typename RealScalar>
struct atanh2_impl<std::complex<RealScalar> >
{
typedef std::complex<RealScalar> Scalar;
static inline Scalar run(const Scalar& x, const Scalar& r)
{
using std::log;
using std::norm;
using std::sqrt;
Scalar z = x / r;
if (r == Scalar(0) || norm(z) > NumTraits<RealScalar>::epsilon())
return RealScalar(0.5) * log((r + x) / (r - x));
else
return z + z*z*z / RealScalar(3);
}
};
template<typename Scalar>
struct atanh2_retval
struct log1p_retval
{
typedef Scalar type;
};
@ -680,9 +661,9 @@ inline EIGEN_MATHFUNC_RETVAL(hypot, Scalar) hypot(const Scalar& x, const Scalar&
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(atanh2, Scalar) atanh2(const Scalar& x, const Scalar& y)
inline EIGEN_MATHFUNC_RETVAL(log1p, Scalar) log1p(const Scalar& x)
{
return EIGEN_MATHFUNC_IMPL(atanh2, Scalar)::run(x, y);
return EIGEN_MATHFUNC_IMPL(log1p, Scalar)::run(x);
}
template<typename Scalar>
@ -727,6 +708,44 @@ inline int log2(int x)
} // end namespace numext
namespace internal {
/****************************************************************************
* Implementation of atanh2 *
****************************************************************************/
template<typename Scalar>
struct atanh2_impl
{
static inline Scalar run(const Scalar& x, const Scalar& r)
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
typedef typename NumTraits<Scalar>::Real RealScalar;
return numext::log1p(RealScalar(2) * x / (r - x)) / RealScalar(2);
}
};
template<typename Scalar>
struct atanh2_retval
{
typedef Scalar type;
};
} // end namespace internal
namespace numext {
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(atanh2, Scalar) atanh2(const Scalar& x, const Scalar& y)
{
return EIGEN_MATHFUNC_IMPL(atanh2, Scalar)::run(x, y);
}
} // end namespace numext
namespace internal {
/****************************************************************************