This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

Bug 876

Summary: Android NDK C++11 build fails to compile on std::log1p in MathFunctions.h (dev branch)
Product: Eigen Reporter: Tim Chernov (SE) <chernov.tim.y>
Component: Core - generalAssignee: Nobody <eigen.nobody>
Status: RESOLVED WORKSFORME    
Severity: Compilation Problem CC: chtz, gael.guennebaud, jacob.benoit.1, jdh8, jitseniesen, kmhofmann
Priority: Normal Keywords: Compatibility
Version: unspecified   
Hardware: All   
OS: Android   
Whiteboard:
Bug Depends on:    
Bug Blocks: 558    

Description Tim Chernov (SE) 2014-09-12 15:35:09 UTC
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.
Comment 1 Christoph Hertzberg 2014-11-03 14:55:58 UTC
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.
Comment 2 Gael Guennebaud 2014-11-04 15:36:46 UTC
yes, let's be conservative, we could enable it for "C++11 & (gcc|clang|msvc|icc) & (x86|x86_64)" only.
Comment 3 Gael Guennebaud 2014-11-05 15:07:47 UTC
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.
Comment 4 Christoph Hertzberg 2014-11-05 15:53:52 UTC
(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.
Comment 5 Gael Guennebaud 2014-11-05 17:34:10 UTC
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.
Comment 6 Christoph Hertzberg 2014-11-05 18:30:23 UTC
(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
Comment 7 Chen-Pang He 2014-11-06 17:51:17 UTC
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()
Comment 8 Christoph Hertzberg 2014-11-06 19:01:59 UTC
(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!
Comment 9 Michael Hofmann 2014-11-18 18:27:15 UTC
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
  }
};
Comment 10 Gael Guennebaud 2014-11-18 20:39:12 UTC
This issue will likely be fixed in the future, so we need to check for the version of GLIBC++
Comment 11 Michael Hofmann 2014-11-19 10:51:57 UTC
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().
Comment 12 Gael Guennebaud 2014-12-08 16:48:42 UTC
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.
Comment 13 Michael Hofmann 2016-01-14 09:51:52 UTC
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
Comment 14 Gael Guennebaud 2016-01-25 10:12:03 UTC
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?
Comment 15 Gael Guennebaud 2016-01-25 10:20:21 UTC
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
Comment 16 Gael Guennebaud 2016-02-12 21:53:50 UTC
Adding APP_STL:= gnustl_static works fine for me with clang3.5 or gcc4.9.
Comment 17 Nobody 2019-12-04 13:43:20 UTC
-- 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.