[GHC] #14927: Hyperbolic area sine is unstable for (even moderately) big negative arguments.
#14927: Hyperbolic area sine is unstable for (even moderately) big negative arguments. -------------------------------------+------------------------------------- Reporter: | Owner: (none) leftaroundabout | Type: bug | Status: new Priority: normal | Milestone: Component: | Version: 8.2.1 libraries/base | Keywords: | Operating System: Unknown/Multiple Architecture: | Type of failure: Incorrect result Unknown/Multiple | at runtime Test Case: | Blocked By: Blocking: | Related Tickets: Differential Rev(s): | Wiki Page: -------------------------------------+------------------------------------- `asinh` is supposed to be the inverse of `sinh`, and this works pretty reliable for positive arguments. However, for negative arguments, the [http://hackage.haskell.org/package/base-4.10.1.0/docs/src/GHC.Float.html#lin... 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)) }}} -- Ticket URL: http://ghc.haskell.org/trac/ghc/ticket/14927 GHC http://www.haskell.org/ghc/ The Glasgow Haskell Compiler
#14927: Hyperbolic area sine is unstable for (even moderately) big negative arguments. -------------------------------------+------------------------------------- Reporter: leftaroundabout | Owner: (none) Type: bug | Status: new Priority: normal | Milestone: Component: libraries/base | Version: 8.2.1 Resolution: | Keywords: Operating System: Unknown/Multiple | Architecture: Type of failure: Incorrect result | Unknown/Multiple at runtime | Test Case: Blocked By: | Blocking: Related Tickets: | Differential Rev(s): Wiki Page: | -------------------------------------+------------------------------------- Comment (by leftaroundabout): Alternatives to the fix suggestion above: - {{{asinh x = signum x * log (abs x + sqrt (1 + x*x))}}} More concise form of the same thing; quite possibly faster performance because a branch can be avoided. - Actually, the Python version has also a higher ''upper'' limit: {{{ Prelude> asinh 1e200 In [3]: asinh(1e200) Infinity Out[3]: 461.2101657793691 }}} The reason is that `x*x` overflows even though the function as a whole could still be well-defined. I consider this not nearly as problematic as the much earlier problems I showed for the negative side, but it's still not great. This issue can be fixed by noting that for `x > 1/ε²`, we have anyways `1 + x*x == x*x` and thus `asinh x` gives the same result as `log (2*x) = log 2 + log x`. So, a ''really'' stable area hyperbolic sine would be {{{ asinh x | x > huge = log 2 + log x | x < 0 = -asinh (-x) | otherwise = log (1 + sqrt (1 + x*x)) }}} with {{{ huge = 1 / last (takeWhile ((>1).(+1)) $ recip . (2^) <$> [0..]) {- == 1 / Numeric.IEEE.epsilon -} }}} Mind, the exact value is completely uncritical, we could as well hard- code `huge = 1e20`. That implementation accurately left-inverts `sinh` (for `Double`) for all the range from `-710` to `710` (in other words, it right-inverts `sinh` for all the range from `-1.17e308` to `1.17e308`), which is basically the same as in Python. The current implementation only reaches up to `sinh 355`. -- Ticket URL: http://ghc.haskell.org/trac/ghc/ticket/14927#comment:1 GHC http://www.haskell.org/ghc/ The Glasgow Haskell Compiler
#14927: Hyperbolic area sine is unstable for (even moderately) big negative arguments. -------------------------------------+------------------------------------- Reporter: leftaroundabout | Owner: (none) Type: bug | Status: new Priority: normal | Milestone: Component: libraries/base | Version: 8.2.1 Resolution: | Keywords: Floating | IEEE754 trigonometric Operating System: Unknown/Multiple | Architecture: Type of failure: Incorrect result | Unknown/Multiple at runtime | Test Case: Blocked By: | Blocking: Related Tickets: | Differential Rev(s): Wiki Page: | -------------------------------------+------------------------------------- Changes (by leftaroundabout): * keywords: => Floating IEEE754 trigonometric -- Ticket URL: http://ghc.haskell.org/trac/ghc/ticket/14927#comment:2 GHC http://www.haskell.org/ghc/ The Glasgow Haskell Compiler
#14927: Hyperbolic area sine is unstable for (even moderately) big negative arguments. -------------------------------------+------------------------------------- Reporter: leftaroundabout | Owner: (none) Type: bug | Status: new Priority: normal | Milestone: Component: libraries/base | Version: 8.2.1 Resolution: | Keywords: Floating | IEEE754 trigonometric Operating System: Unknown/Multiple | Architecture: Type of failure: Incorrect result | Unknown/Multiple at runtime | Test Case: Blocked By: | Blocking: Related Tickets: | Differential Rev(s): Wiki Page: | -------------------------------------+------------------------------------- Comment (by bgamari): Either of these options sound reasonable to me. Would you be willing to offer a patch? -- Ticket URL: http://ghc.haskell.org/trac/ghc/ticket/14927#comment:3 GHC http://www.haskell.org/ghc/ The Glasgow Haskell Compiler
#14927: Hyperbolic area sine is unstable for (even moderately) big negative arguments. -------------------------------------+------------------------------------- Reporter: leftaroundabout | Owner: (none) Type: bug | Status: new Priority: normal | Milestone: Component: libraries/base | Version: 8.2.1 Resolution: | Keywords: Floating | IEEE754 trigonometric Operating System: Unknown/Multiple | Architecture: Type of failure: Incorrect result | Unknown/Multiple at runtime | Test Case: Blocked By: | Blocking: Related Tickets: | Differential Rev(s): Wiki Page: | -------------------------------------+------------------------------------- Comment (by leftaroundabout): Fix available at [https://github.com/leftaroundabout/ghc/tree/testing /float-inverses]. -- Ticket URL: http://ghc.haskell.org/trac/ghc/ticket/14927#comment:4 GHC http://www.haskell.org/ghc/ The Glasgow Haskell Compiler
#14927: Hyperbolic area sine is unstable for (even moderately) big negative arguments. -------------------------------------+------------------------------------- Reporter: leftaroundabout | Owner: (none) Type: bug | Status: new Priority: normal | Milestone: Component: libraries/base | Version: 8.2.1 Resolution: | Keywords: Floating | IEEE754 trigonometric Operating System: Unknown/Multiple | Architecture: Type of failure: Incorrect result | Unknown/Multiple at runtime | Test Case: Blocked By: | Blocking: Related Tickets: | Differential Rev(s): Wiki Page: | -------------------------------------+------------------------------------- Changes (by dfeuer): * cc: dfeuer (added) Comment: I don't really understand how you get your calculation of `huge`. Shouldn't that be `1/sqrt epsilon` rather than `1/epsilon`? If so, the cut-offs come out quite different. -- Ticket URL: http://ghc.haskell.org/trac/ghc/ticket/14927#comment:5 GHC http://www.haskell.org/ghc/ The Glasgow Haskell Compiler
#14927: Hyperbolic area sine is unstable for (even moderately) big negative arguments. -------------------------------------+------------------------------------- Reporter: leftaroundabout | Owner: (none) Type: bug | Status: new Priority: normal | Milestone: Component: libraries/base | Version: 8.2.1 Resolution: | Keywords: Floating | IEEE754 trigonometric Operating System: Unknown/Multiple | Architecture: Type of failure: Incorrect result | Unknown/Multiple at runtime | Test Case: Blocked By: | Blocking: Related Tickets: | Differential Rev(s): Wiki Page: | -------------------------------------+------------------------------------- Comment (by leftaroundabout): Yes, the cut-offs come out different, but it doesn't really matter: `sqrt (1 + x*x)` is still stable and accurate for a long way even when `x` is so big that `1 + x*x == x*x`. So, whereas it's important that `huge` should not be ''smaller'' than `1/sqrt epsilon ≈ 6.7e7`, it ''can'' actually be much bigger: `1e100` would also work in case of `Double`. Therefore, `1/epsilon ≈ 4.5e15` is a safe compromise that's presumably less likely to give new surprising corner-cases. -- Ticket URL: http://ghc.haskell.org/trac/ghc/ticket/14927#comment:6 GHC http://www.haskell.org/ghc/ The Glasgow Haskell Compiler
#14927: Hyperbolic area sine is unstable for (even moderately) big negative arguments. -------------------------------------+------------------------------------- Reporter: leftaroundabout | Owner: (none) Type: bug | Status: patch Priority: normal | Milestone: Component: libraries/base | Version: 8.2.1 Resolution: | Keywords: Floating | IEEE754 trigonometric Operating System: Unknown/Multiple | Architecture: Type of failure: Incorrect result | Unknown/Multiple at runtime | Test Case: Blocked By: | Blocking: Related Tickets: | Differential Rev(s): Phab:D4672, Wiki Page: | Phab:D4671, Phab:D4670 -------------------------------------+------------------------------------- Changes (by bgamari): * status: new => patch * differential: => Phab:D4672, Phab:D4671, Phab:D4670 Comment: Thanks leftaroundabout! I've pushed the patches on this branch to Phab:D4672, Phab:D4671, and Phab:D4670. -- Ticket URL: http://ghc.haskell.org/trac/ghc/ticket/14927#comment:7 GHC http://www.haskell.org/ghc/ The Glasgow Haskell Compiler
#14927: Hyperbolic area sine is unstable for (even moderately) big negative arguments. -------------------------------------+------------------------------------- Reporter: leftaroundabout | Owner: (none) Type: bug | Status: closed Priority: normal | Milestone: Component: libraries/base | Version: 8.2.1 Resolution: fixed | Keywords: Floating | IEEE754 trigonometric Operating System: Unknown/Multiple | Architecture: Type of failure: Incorrect result | Unknown/Multiple at runtime | Test Case: Blocked By: | Blocking: Related Tickets: | Differential Rev(s): Phab:D4672, Wiki Page: | Phab:D4671, Phab:D4670 -------------------------------------+------------------------------------- Changes (by bgamari): * status: patch => closed * resolution: => fixed Comment: Pushed as 3ea33411d7cbf32c20940cc72ca07df6830eeed7. -- Ticket URL: http://ghc.haskell.org/trac/ghc/ticket/14927#comment:8 GHC http://www.haskell.org/ghc/ The Glasgow Haskell Compiler
#14927: Hyperbolic area sine is unstable for (even moderately) big negative
arguments.
-------------------------------------+-------------------------------------
Reporter: leftaroundabout | Owner: (none)
Type: bug | Status: closed
Priority: normal | Milestone:
Component: libraries/base | Version: 8.2.1
Resolution: fixed | Keywords: Floating
| IEEE754 trigonometric
Operating System: Unknown/Multiple | Architecture:
Type of failure: Incorrect result | Unknown/Multiple
at runtime | Test Case:
Blocked By: | Blocking:
Related Tickets: | Differential Rev(s): Phab:D4672,
Wiki Page: | Phab:D4671, Phab:D4670
-------------------------------------+-------------------------------------
Comment (by Ben Gamari
participants (1)
-
GHC