[commit: ghc] master: Stable area hyperbolic sine for `Double` and `Float`. (3ea3341)

[git at git.haskell.org]

Repository : ssh://git@git.haskell.org/ghc

On branch  : master
Link       : http://ghc.haskell.org/trac/ghc/changeset/3ea33411d7cbf32c20940cc72ca07df6830eeed7/ghc

>---------------------------------------------------------------

commit 3ea33411d7cbf32c20940cc72ca07df6830eeed7
Author: Justus Sagemüller <sagemueller at geo.uni-koeln.de>
Date:   Wed Mar 28 15:51:16 2018 +0200

    Stable area hyperbolic sine for `Double` and `Float`.
    
    This function was unstable, in particular for negative arguments.
    
    https://ghc.haskell.org/trac/ghc/ticket/14927
    
    The reason is that the formula `log (x + sqrt (1 + x*x))` is dominated
    by the numerical error of the `sqrt` function when x is strongly negative
    (and thus the summands in the `log` mostly cancel). However, the area
    hyperbolic sine is an odd function, thus the negative side can as well
    be calculated by flipping over the positive side, which avoids this instability.
    
    Furthermore, for _very_ big arguments, the `x*x` subexpression overflows. However,
    long before that happens, the square root is anyways completely dominated
    by that term, so we can neglect the `1 +` and get
    
        sqrt (1 + x*x) ≈ sqrt (x*x) = x
    
    and therefore
    
        asinh x ≈ log (x + x) = log (2*x) = log 2 + log x
    
    which does not overflow for any normal-finite positive argument, but
    perfectly matches the exact formula within the floating-point accuracy.


>---------------------------------------------------------------

3ea33411d7cbf32c20940cc72ca07df6830eeed7
 libraries/base/GHC/Float.hs              | 12 ++++++++++--
 testsuite/tests/numeric/should_run/all.T |  2 +-
 2 files changed, 11 insertions(+), 3 deletions(-)

diff --git a/libraries/base/GHC/Float.hs b/libraries/base/GHC/Float.hs
index c534baf..d60c660 100644
--- a/libraries/base/GHC/Float.hs
+++ b/libraries/base/GHC/Float.hs
@@ -367,7 +367,11 @@ instance  Floating Float  where
     (**) x y            =  powerFloat x y
     logBase x y         =  log y / log x
 
-    asinh x = log (x + sqrt (1.0+x*x))
+    asinh x
+      | x > huge   = log 2 + log x
+      | x < 0      = -asinh (-x)
+      | otherwise  = log (x + sqrt (1 + x*x))
+     where huge = 1e10
     acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
     atanh x = 0.5 * log ((1.0+x) / (1.0-x))
 
@@ -492,7 +496,11 @@ instance  Floating Double  where
     (**) x y            =  powerDouble x y
     logBase x y         =  log y / log x
 
-    asinh x = log (x + sqrt (1.0+x*x))
+    asinh x
+      | x > huge   = log 2 + log x
+      | x < 0      = -asinh (-x)
+      | otherwise  = log (x + sqrt (1 + x*x))
+     where huge = 1e20
     acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
     atanh x = 0.5 * log ((1.0+x) / (1.0-x))
 
diff --git a/testsuite/tests/numeric/should_run/all.T b/testsuite/tests/numeric/should_run/all.T
index 9f2f52a..37fff44 100644
--- a/testsuite/tests/numeric/should_run/all.T
+++ b/testsuite/tests/numeric/should_run/all.T
@@ -41,7 +41,7 @@ test('arith018', normal, compile_and_run, [''])
 test('arith019', normal, compile_and_run, [''])
 test('expfloat', normal, compile_and_run, [''])
 
-test('FloatFnInverses', expect_broken(14927), compile_and_run, [''])
+test('FloatFnInverses', normal, compile_and_run, [''])
 
 test('T1603', skip, compile_and_run, [''])
 test('T3676', expect_broken(3676), compile_and_run, [''])