Here's the code of atanh2_impl in MathFunctions.h (starts at line 362) /**************************************************************************** * 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) #if (__cplusplus >= 201103L) && !defined(__CYGWIN__) using std::log1p; return log1p(2 * x / (r - x)) / 2; #else using std::abs; 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; #endif } }; Here's a compilation error: workdir/Eigen/src/Core/MathFunctions.h>: In static member function 'static Scalar Eigen::internal::atanh2_impl<Scalar>::run(const Scalar&, const Scalar&)': workdir/Eigen/src/Core/MathFunctions.h>:373:18: error: 'std::log1p' has not been declared using std::log1p; I guess it's not enough to just check the C++ version, maybe not-android check should be done as well.
Yes, we should do that. Patch welcome! (Or just tell us the flag for which we can check). However, maybe it would be safer to check for systems which do support std::log1p instead for systems which don't. OTOH, that would silently generate inefficient code for such systems. OTOOH, how bad is our fall-back implementation, compared to the log1p implementation? In this case, we might simply use it as default.
yes, let's be conservative, we could enable it for "C++11 & (gcc|clang|msvc|icc) & (x86|x86_64)" only.
Actually, it seems to me that such a atanh2 function is rather non-standard. If I'm correct, this function has been introduced for the MathFunctions module, and it is its unique use case. In this module, given two numbers a and b, atanh2 is called like this: 2*atanh2(a-b,a+b) Then atanh2(x,r) recovers a and b by computing 2b=r-x and 2a=r+x and finally computes log(1+(a-b)/b)/2. All these intermediate operations are a waste of time and accuracy. So what about directly calling: log1p((a-b)/b) ? and thus replace atanh2 by a portable log1p function.
(In reply to Gael Guennebaud from comment #3) Indeed. 2*atanh2(a-b,a+b) == log(a/b) == log(1+(a-b)/b) I don't even understand how log1p((a-b)/b) gains any accuracy over directly calling log(a/b)? I.e. in what way would log1p((a-b)/b) be better than log1p(a/b-1.0)? And for the latter call manually subtracting 1 and calling log1p is evidently pointless. However, what could make sense here is having a Taylor approximation for log(a/b)/(a-b) = 1/b - (a-b)/(2*b*b) + O{(a-b)^2} if a and b are close. http://www.wolframalpha.com/input/?i=series+log%28a%2Fb%29%2F%28a-b%29+at+a%3Db This should also seamlessly integrate the A(0,0)==A(1,1) branch. Admittedly, I did not entirely check the complex cases of matrix_log_compute_2x2 yet.
Using log1p remains more accurate when a approx b because it is as accurate as the first order Taylor approx of log(a/b), and log(a/b) is also needed in MatrixPower.h (without the division by a-b). Here is a selfcontained example comparing different implementations of log(a/b) and log(a/b)/(a-b): #include <cmath> #include <iostream> template<typename Scalar> Scalar atanh2(const Scalar& x, const Scalar& r) { using std::abs; using std::log; using std::sqrt; Scalar z = x / r; if (r == 0 || abs(z) > 1e-8) return log((r + x) / (r - x)) / 2; else return z + z*z*z / 3; } void foo(double a, double b) { double r1 = std::log(a/b); double r2 = std::log1p((a-b)/b); double r3 = 2.0*atanh2(a-b,a+b); double r4 = (a-b)/b; // Taylor of log(a/b) for a=b double r5 = 1.0/b - (a-b)/(2*b*b); // Taylor of log(a/b)/(a-b) for a=b std::cout.precision(20); std::cout << "a=" << a << " b=" << b << "\n"; std::cout << "log(a/b): " << r1 << "\t " << r1/(a-b) << "\n" << "log1p((a-b)/b) " << r2 << "\t " << r2/(a-b) << "\n" << "atanh2(a-b,a+b) " << r3 << "\t " << r3/(a-b) << "\n" << "series log(a/b) " << r4 << "\t " << r4/(a-b) << "\n" << "series log(a/b)/(a-b) " << r5*(a-b) << "\t " << r5 << "\n\n"; } int main() { foo(1,1+1e-15); foo(1,1+2e-16); foo(1,1+1e-16); } and the results are: a=1 b=1.0000000000000011102 log(a/b): -1.1102230246251571321e-15 1.0000000000000004441 log1p((a-b)/b) -1.1102230246251559488e-15 0.99999999999999944489 atanh2(a-b,a+b) -1.110223024625156146e-15 0.99999999999999966693 series log(a/b) -1.1102230246251553571e-15 0.99999999999999888978 series log(a/b)/(a-b) -1.1102230246251559488e-15 0.99999999999999944489 a=1 b=1.000000000000000222 log(a/b): -2.2204460492503130808e-16 1 log1p((a-b)/b) -2.2204460492503128343e-16 0.99999999999999988898 atanh2(a-b,a+b) -2.2204460492503130808e-16 1 series log(a/b) -2.2204460492503125878e-16 0.99999999999999977796 series log(a/b)/(a-b) -2.2204460492503128343e-16 0.99999999999999988898 a=1 b=1 log(a/b): 0 nan log1p((a-b)/b) 0 nan atanh2(a-b,a+b) 0 nan series log(a/b) 0 nan series log(a/b)/(a-b) 0 1 This also shows that the current atanh2 is clearly not very good.
(In reply to Gael Guennebaud from comment #5) > Using log1p remains more accurate when a approx b because it is as accurate You're right! I only considered that the relative errors of (a-b)/b and a/b-1.0 are the same magnitude. However, if we divide by (a-b) afterwards, we will have exactly the same relative error as in (a-b)/b, whereas a/b-1.0 can round off in a different direction. [Actually, with x87 math there could still be cases where both (a-b) expressions are rounded differently, but I guess that's not too severe] Even more impressive case against log(a/b) is foo((1-1e-15)*std::sqrt(2.0),std::sqrt(2.0)); a=1.4142135623730938132 b=1.4142135623730951455 log(a/b): -8.8817841970012562677e-16 0.66666666666666696273 log1p((a-b)/b) -9.4205547521026534655e-16 0.70710678118654779478 atanh2(a-b,a+b) -9.4205547521026534655e-16 0.70710678118654779478 series log(a/b) -9.4205547521026495212e-16 0.70710678118654746172 series log(a/b)/(a-b) -9.4205547521026534655e-16 0.70710678118654779478 --> about 5% error for log(a/b)/(a-b). Actually, MatrixPower also divides by (a-b), however it calculates sinh(p*w)/( curr - prev ); // with w=atanh2(curr-prev, curr+prev); (and sinh(x) is close to x, for small x) http://www.wolframalpha.com/input/?i=series+sinh%28log%28a%2Fb%29*p%29%2F%28a-b%29+at+a%3Db
Even if we have a (new) portable log1p, we still need atanh2 for MatrixLogarithm, which is the origin of atanh2. \sa internal::matrix_log_compute_2x2()
(In reply to Chen-Pang He from comment #7) > Even if we have a (new) portable log1p, we still need atanh2 for > MatrixLogarithm, which is the origin of atanh2. > > \sa internal::matrix_log_compute_2x2() We have (neglecting the complex unwinding for brevity): Scalar y = A(1,1) - A(0,0), x = A(1,1) + A(0,0); result(0,1) = A(0,1) * (Scalar(2) * numext::atanh2(y,x) ) / y; That is res = A01*2*atanh2(A11-A00, A11+A00)/(A11-A00) = A01*log(A11/A00)/(A11-A00) = A01*log(1+(A11-A00)/A00))/(A11-A00) --> Most likely less rounding errors with log1p!
It seems that libstdc++ included with the Android NDK does just not pull log1p() into the std:: namespace; the function itself is supported, though. This problem occurs only when using libstdc++ (e.g. when specifying APP_STL := gnustl_[static|shared] in Application.mk); there are no issues with libc++ (e.g. when specifying APP_STL := c++_[static|shared] in Application.mk). A valid fix therefore is adding the following check for 1) Android, and 2) libstdc++. (Whether that nesting of #ifs is desirable is a different question, of course.) template<typename Scalar> struct atanh2_impl { static inline Scalar run(const Scalar& x, const Scalar& r) { EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) #if (__cplusplus >= 201103L) && !defined(__CYGWIN__) #if !(defined(__ANDROID__) && defined(__GLIBCXX__)) using std::log1p; #endif return log1p(2 * x / (r - x)) / 2; #else using std::abs; 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; #endif } };
This issue will likely be fixed in the future, so we need to check for the version of GLIBC++
For the case of NDK r10c, this would make the check #if !(defined(__ANDROID__) && defined(__GLIBCXX__) && (__GLIBCXX__ <= 20140827)) I'm not so sure when this will be fixed, however, so I'm not sure how safe a bet the addition of the version check is. A similar bug was already filed in August in the Android bug tracker, and acknowledged though not fixed yet: https://code.google.com/p/android/issues/detail?id=74835 Seems that a bunch of other functions are also affected, such as copysign().
I implemented a custom version of log1p so that we can be very conservative to enable the use of the C++11 standard one. I also removed atanh2. Feel free to re-open if the tests is not conservative enough. https://bitbucket.org/eigen/eigen/commits/56cb933e3b96/ Changeset: 56cb933e3b96 User: ggael Date: 2014-12-08 15:26:53+00:00 Summary: Bug 876: implement a portable log1p function Affected #: 1 file https://bitbucket.org/eigen/eigen/commits/1f5ee4502698/ Changeset: 1f5ee4502698 User: ggael Date: 2014-12-08 15:28:06+00:00 Summary: Bug 876, matrix_log_compute_2x2: directly use logp1 instead of atanh2 Affected #: 1 file https://bitbucket.org/eigen/eigen/commits/668f38f99419/ Changeset: 668f38f99419 User: ggael Date: 2014-12-08 15:44:05+00:00 Summary: Bug 876: remove usage of atanh2 in matrix power Affected #: 1 file https://bitbucket.org/eigen/eigen/commits/96c5f6a3e844/ Changeset: 96c5f6a3e844 User: ggael Date: 2014-12-08 15:44:34+00:00 Summary: Remove useless and non standard numext::atanh2 function.
This issue seems to have resurfaced when compiling using Android NDK/GCC. (NDK/Clang seems unaffected.) Using r10e/GCC 4.9, I get error messages such as: /eigen/Eigen/src/Core/MathFunctions.h: In static member function 'static Scalar Eigen::internal::round_impl<Scalar>::run(const Scalar&)': /eigen/Eigen/src/Core/MathFunctions.h:388:18: error: 'std::round' has not been declared using std::round; /eigen/Eigen/src/Core/MathFunctions.h: In static member function 'static Scalar Eigen::internal::log1p_impl<Scalar, false>::run(const Scalar&)': /eigen/Eigen/src/Core/MathFunctions.h:479:16: error: 'std::log1p' has not been declared using std::log1p; A working fix is to change: diff -r 98fcfcb99dee Eigen/src/Core/util/Macros.h --- a/Eigen/src/Core/util/Macros.h Wed Jan 13 14:24:37 2016 -0800 +++ b/Eigen/src/Core/util/Macros.h Thu Jan 14 10:50:26 2016 +0100 @@ -372,8 +372,9 @@ // Does the compiler support C++11 math? // Let's be conservative and enable the default C++11 implementation only if we are sure it exists #ifndef EIGEN_HAS_CXX11_MATH - #if (__cplusplus > 201103L) || (__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) + #if (__cplusplus > 201103L && !(EIGEN_COMP_GNUC_STRICT && EIGEN_OS_ANDROID)) \ + || (__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) #define EIGEN_HAS_CXX11_MATH 1 #else #define EIGEN_HAS_CXX11_MATH 0
That's a shame that the compiler claims C++14 without even providing a C++11 standard library. Maybe there is an analogue to __cplusplus but for the STL?
hm, maybe you have to specify the correct STL to use using APP_STL and LOCAL_C_INCLUDES, e.g.: APP_STL:= gnustl_static LOCAL_C_INCLUDES += ${ANDROID_NDK}/sources/cxx-stl/gnu-libstdc++/4.9/include
Adding APP_STL:= gnustl_static works fine for me with clang3.5 or gcc4.9.
-- GitLab Migration Automatic Message -- This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity. You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/876.