author Justus Sagemüller Wed, 28 Mar 2018 13:51:16 +0000 (15:51 +0200) committer Ben Gamari Sun, 6 May 2018 00:28:18 +0000 (20:28 -0400)
This function was unstable, in particular for negative arguments.

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.

index c534baf..d60c660 100644 (file)
@@ -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))

index 9f2f52a..37fff44 100644 (file)
@@ -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, [''])