Hi,

tested x86_64-linux, committed to mainline.

Paolo.

//////////////////
2011-11-02  Richard B. Kreckel  <krec...@ginac.de>
            Paolo Carlini  <paolo.carl...@oracle.com>

        PR libstdc++/50880
        * include/std/complex (__complex_acosh): Fix in a better way,
        use Kahan's formula.
        * include/tr1/complex (__complex_acosh): Likewise.
Index: include/std/complex
===================================================================
--- include/std/complex (revision 180785)
+++ include/std/complex (working copy)
@@ -1686,14 +1686,9 @@
     std::complex<_Tp>
     __complex_acosh(const std::complex<_Tp>& __z)
     {
-      std::complex<_Tp> __t((__z.real() - __z.imag())
-                           * (__z.real() + __z.imag()) - _Tp(1.0),
-                           _Tp(2.0) * __z.real() * __z.imag());
-      __t = std::sqrt(__t);
-      if (__z.real() < _Tp())
-       __t = -__t;
-
-      return std::log(__t + __z);
+      // Kahan's formula.
+      return _Tp(2.0) * std::log(std::sqrt(_Tp(0.5) * (__z + _Tp(1.0)))
+                                + std::sqrt(_Tp(0.5) * (__z - _Tp(1.0))));
     }
 
 #if _GLIBCXX_USE_C99_COMPLEX_TR1
Index: include/tr1/complex
===================================================================
--- include/tr1/complex (revision 180785)
+++ include/tr1/complex (working copy)
@@ -185,14 +185,9 @@
     std::complex<_Tp>
     __complex_acosh(const std::complex<_Tp>& __z)
     {
-      std::complex<_Tp> __t((__z.real() - __z.imag())
-                           * (__z.real() + __z.imag()) - _Tp(1.0),
-                           _Tp(2.0) * __z.real() * __z.imag());
-      __t = std::sqrt(__t);
-      if (__z.real() < _Tp())
-       __t = -__t;
-
-      return std::log(__t + __z);
+      // Kahan's formula.
+      return _Tp(2.0) * std::log(std::sqrt(_Tp(0.5) * (__z + _Tp(1.0)))
+                                + std::sqrt(_Tp(0.5) * (__z - _Tp(1.0))));
     }
 
 #if _GLIBCXX_USE_C99_COMPLEX_TR1

Reply via email to