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 |