From 3ea33411d7cbf32c20940cc72ca07df6830eeed7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= Date: Wed, 28 Mar 2018 15:51:16 +0200 Subject: [PATCH] Stable area hyperbolic sine for `Double` and `Float`. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- 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 c534bafa07..d60c660bd0 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 9f2f52a0fb..37fff44bde 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, ['']) -- GitLab