Commit 3ea33411 authored by Justus Sagemüller's avatar Justus Sagemüller Committed by Ben Gamari

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.
parent d814dd38
...@@ -367,7 +367,11 @@ instance Floating Float where ...@@ -367,7 +367,11 @@ instance Floating Float where
(**) x y = powerFloat x y (**) x y = powerFloat x y
logBase x y = log y / log x 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))) 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)) atanh x = 0.5 * log ((1.0+x) / (1.0-x))
...@@ -492,7 +496,11 @@ instance Floating Double where ...@@ -492,7 +496,11 @@ instance Floating Double where
(**) x y = powerDouble x y (**) x y = powerDouble x y
logBase x y = log y / log x 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))) 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)) atanh x = 0.5 * log ((1.0+x) / (1.0-x))
......
...@@ -41,7 +41,7 @@ test('arith018', normal, compile_and_run, ['']) ...@@ -41,7 +41,7 @@ test('arith018', normal, compile_and_run, [''])
test('arith019', normal, compile_and_run, ['']) test('arith019', normal, compile_and_run, [''])
test('expfloat', 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('T1603', skip, compile_and_run, [''])
test('T3676', expect_broken(3676), compile_and_run, ['']) test('T3676', expect_broken(3676), compile_and_run, [''])
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment