Hyperbolic area sine is unstable for (even moderately) big negative arguments.

asinh is supposed to be the inverse of sinh, and this works pretty reliable for positive arguments. However, for negative arguments, the currently used formula

asinh x = log (x + sqrt (1.0+x*x))

gets unstable much earlier than necessary, namely when the summands in the logarithm cancel almost to zero, dominated by the numerical error of the square root.

This is particularly troubling because mathematically **a)** asinh is a very “inert” function (i.e. you can carelessly put in huge numbers and – as long as they're not outright Infinity – always get a somewhat sane result), pseudo-sigmoidal as it were **b)** it is an odd function, fulfilling asinh (-x) = -asinh x.

Both is reflected in other implementations, e.g. Python, but not in GHC Haskell:

GHCi, version 8.2.1         Python 3.5.2 (default, Nov 23 2017, 16:37:01) 

                            In [1]: from math import *

Prelude> asinh 1e6          In [2]: asinh(1e6)
14.50865773852447           Out[2]: 14.50865773852447

Prelude> asinh (-1e6)       In [3]: asinh(-1e6)
-14.50865012405984          Out[3]: -14.50865773852447

Prelude> asinh 1e9          In [4]: asinh(1e9)
21.416413017506358          Out[4]: 21.416413017506354

Prelude> asinh (-1e9)       In [5]: asinh(-1e9)
-Infinity                   Out[5]: -21.416413017506354

Prelude> asinh 1e76         In [6]: asinh(1e76)
175.6896142481074           Out[6]: 175.68961424810743

Prelude> asinh (-1e76)      In [7]: asinh(-1e76)
-Infinity                   Out[7]: -175.68961424810743

Demo of non-inverse property:

Prelude> [(x, asinh $ sinh x) | x <- [-25..25]]
[(-25.0,-Infinity)
,(-24.0,-Infinity)
,(-23.0,-Infinity)
,(-22.0,-Infinity)
,(-21.0,-Infinity)
,(-20.0,-Infinity)
,(-19.0,-18.021826694558577)
,(-18.0,-18.021826694558577)
,(-17.0,-17.0102257828801)
,(-16.0,-15.998624871201619)
,(-15.0,-14.999878578873695)
,(-14.0,-13.999968823323222)
,(-13.0,-12.999991335176079)
,(-12.0,-12.000000137072186)
,(-11.0,-10.999999903206444)
,(-10.0,-10.000000013503529)
,(-9.0,-9.000000000551713)
,(-8.0,-8.00000000017109)
,(-7.0,-7.000000000036329)
,(-6.0,-5.999999999998066)
,(-5.0,-5.000000000000641)
,(-4.0,-4.000000000000046)
,(-3.0,-2.999999999999989)
,(-2.0,-1.9999999999999991)
,(-1.0,-1.0)
,(0.0,0.0)
,(1.0,1.0)
,(2.0,2.0)
,(3.0,3.0)
,(4.0,4.0)
,(5.0,5.0)
,(6.0,6.0)
,(7.0,7.0)
,(8.0,8.0)
,(9.0,9.0)
,(10.0,10.0)
,(11.0,11.0)
,(12.0,12.0)
,(13.0,13.0)
,(14.0,14.0)
,(15.0,15.0)
,(16.0,16.0)
,(17.0,17.0)
,(18.0,18.0)
,(19.0,19.0)
,(20.0,20.0)
,(21.0,21.0)
,(22.0,22.0)
,(23.0,23.0)
,(24.0,24.0)
,(25.0,25.0)]

Those results are less than satisfying, even for inputs that aren't astronomically big at all.

A simple fix would be to “copy” the sane positive-number behaviour to the negative side:

    asinh x
      | x < 0      = -asinh (-x)
      | otherwise  = log (x + sqrt (1.0+x*x))

Trac metadata

Trac field	Value
Version	8.2.1
Type	Bug
TypeOfFailure	OtherFailure
Priority	normal
Resolution	Unresolved
Component	libraries/base
Test case
Differential revisions
BlockedBy
Related
Blocking
CC
Operating system
Architecture