accuracy of asinh and atanh

Hi devs, I tried to use asinh :: Double -> Double and discovered that it's inaccurate compared to my system library (GNU libm), even returning -Infinity in place of finite values in the neighborhood of -22 for large negative arguments. `atanh` is also inaccurate compared to the system library. I wrote up a more detailed description of the problem including plots in the README file at https://github.com/peddie/ghc-inverse-hyperbolic -- this repository is package that can help you examine the error for yourself or generate the plots, and it also contains accurate pure-Haskell translations of the system library's implementation for these functions. What's the next step to fixing this in GHC? Cheers Matt Peddie

Note: From skimming your readme it is worth noting that log1p _is_ in base now (alongside expm1, log1pexp, and log1mexp). We added them all a couple of years back as a result of the very thread linked in your README. You need to `import Numeric` to see them, though. Switching to more accurate functions for doubles and floats for asinh, atanh, etc. to exploit this sort of functionality at least seems to make a lot of sense. That can be done locally without any user API impact as the current definitions aren't supplied as defaults, merely as pointwise implementations instance by instance. Things will just become more accurate. In that same spirit, we can probably crib a better version for complex numbers from somewhere as well, as it follows the same general simplistic formula right now, even if it can't be plugged directly into the equations you've given. For that matter, the log1p definition we're using for complex numbers was the best I could come up with, but there may well be a more accurate version you can find down in the mines of libm or another math library written by real analysts. log1p https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Float.html#log... x https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#...@(a https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#... :+ https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#... b https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#...) | abs https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Num.html#abs a https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#... < 0.5 && abs https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Num.html#abs b https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#... < 0.5 , u https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#... <- 2* https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Num.html#%2Aa https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#... + https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Num.html#%2B a https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#...* https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Num.html#%2Aa https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#... + https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Num.html#%2B b https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#...* https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Num.html#%2Ab https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#... = log1p https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Float.html#log... (u https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#.../ https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Real.html#%2F(1 + https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Num.html#%2B sqrt https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Float.html#sqr...(u https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#...+ https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Num.html#%2B1))) :+ https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#... atan2 https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Float.html#ata... (1 + https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Num.html#%2B a https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#...) b https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#... | otherwise https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Base.html#othe... = log https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Float.html#log (1 + https://hackage.haskell.org/package/base-4.9.1.0/docs/src/GHC.Num.html#%2B x https://hackage.haskell.org/package/base-4.9.1.0/docs/src/Data.Complex.html#...) So, here's a +1 from the libraries committee side towards improving the situation.
From there, it's a small matter of implementation.
Here's where I'd usually get Ben involved. Hi Ben!
-Edward
On Sat, Jun 2, 2018 at 1:23 AM, Matt Peddie
Hi devs,
I tried to use asinh :: Double -> Double and discovered that it's inaccurate compared to my system library (GNU libm), even returning -Infinity in place of finite values in the neighborhood of -22 for large negative arguments. `atanh` is also inaccurate compared to the system library. I wrote up a more detailed description of the problem including plots in the README file at https://github.com/peddie/ghc-inverse-hyperbolic -- this repository is package that can help you examine the error for yourself or generate the plots, and it also contains accurate pure-Haskell translations of the system library's implementation for these functions. What's the next step to fixing this in GHC?
Cheers
Matt Peddie _______________________________________________ ghc-devs mailing list ghc-devs@haskell.org http://mail.haskell.org/cgi-bin/mailman/listinfo/ghc-devs

Sending again, this time including ghc-devs.
Edward Kmett
Note: From skimming your readme it is worth noting that log1p _is_ in base now (alongside expm1, log1pexp, and log1mexp). We added them all a couple of years back as a result of the very thread linked in your README.
You need to `import Numeric` to see them, though.
Switching to more accurate functions for doubles and floats for asinh, atanh, etc. to exploit this sort of functionality at least seems to make a lot of sense.
That can be done locally without any user API impact as the current definitions aren't supplied as defaults, merely as pointwise implementations instance by instance. Things will just become more accurate.
In that same spirit, we can probably crib a better version for complex numbers from somewhere as well, as it follows the same general simplistic formula right now, even if it can't be plugged directly into the equations you've given. For that matter, the log1p definition we're using for complex numbers was the best I could come up with, but there may well be a more accurate version you can find down in the mines of libm or another math library written by real analysts.
[snip]
So, here's a +1 from the libraries committee side towards improving the situation.
From there, it's a small matter of implementation.
Here's where I'd usually get Ben involved. Hi Ben!
Hi Edward and Matt! Indeed the sinh sinh situation should already be improved in 8.6. See #14927 and GHC commit 3ea33411d7cbf32c20940cc72ca07df6830eeed7. I suspect this should fix Matt's issue. Regarding log1p, I'd happily accept patches improving on the status quo. Cheers, - Ben

Hello Mat
Just curious, why the preferred solution isn't to call the system math
library? As it says in the README you reference below,
- One good solution would be to always call the system math library for
these functions.
Hope this is isn't a stupid question.
Thanks
George
On Sat, Jun 2, 2018 at 2:23 AM Matt Peddie
Hi devs,
I tried to use asinh :: Double -> Double and discovered that it's inaccurate compared to my system library (GNU libm), even returning -Infinity in place of finite values in the neighborhood of -22 for large negative arguments. `atanh` is also inaccurate compared to the system library. I wrote up a more detailed description of the problem including plots in the README file at https://github.com/peddie/ghc-inverse-hyperbolic -- this repository is package that can help you examine the error for yourself or generate the plots, and it also contains accurate pure-Haskell translations of the system library's implementation for these functions. What's the next step to fixing this in GHC?
Cheers
Matt Peddie _______________________________________________ ghc-devs mailing list ghc-devs@haskell.org http://mail.haskell.org/cgi-bin/mailman/listinfo/ghc-devs

Hi George,
Not a stupid question. I don't have a single source at hand, but I
think I read in a few places on the wiki that calling out to the
system math library is not an option due to the variety of system math
libraries on the platforms GHC supports. It'd be great if I got the
wrong impression and this could just be a call to C. Can anyone set
me straight on this point?
Matt
On Thu, Aug 2, 2018 at 5:13 AM, George Colpitts
Hello Mat
Just curious, why the preferred solution isn't to call the system math library? As it says in the README you reference below,
One good solution would be to always call the system math library for these functions.
Hope this is isn't a stupid question.
Thanks
George
On Sat, Jun 2, 2018 at 2:23 AM Matt Peddie
wrote: Hi devs,
I tried to use asinh :: Double -> Double and discovered that it's inaccurate compared to my system library (GNU libm), even returning -Infinity in place of finite values in the neighborhood of -22 for large negative arguments. `atanh` is also inaccurate compared to the system library. I wrote up a more detailed description of the problem including plots in the README file at https://github.com/peddie/ghc-inverse-hyperbolic -- this repository is package that can help you examine the error for yourself or generate the plots, and it also contains accurate pure-Haskell translations of the system library's implementation for these functions. What's the next step to fixing this in GHC?
Cheers
Matt Peddie _______________________________________________ ghc-devs mailing list ghc-devs@haskell.org http://mail.haskell.org/cgi-bin/mailman/listinfo/ghc-devs

Matt Peddie
Hi George,
Not a stupid question. I don't have a single source at hand, but I think I read in a few places on the wiki that calling out to the system math library is not an option due to the variety of system math libraries on the platforms GHC supports. It'd be great if I got the wrong impression and this could just be a call to C. Can anyone set me straight on this point?
Indeed it's not a stupid question at all. Indeed this is precisely what we do for the simpler transcendentals (e.g. sin, asin, log). We very well could move in this direction in the case of asinh/atanh as well. I believe the reason we don't currently is that atanh was only standardized in C99, which we only started requiring a few releases ago. Perhaps this is ultimately the right direction. Cheers, - Ben

Thanks, Ben, for chiming in. I think calling out to C for these
functions is the way to go if it's now feasible. (Calling out to libm
is the workaround I'm using in the application that led me to discover
the inaccuracy.)
Matt
On Thu, Aug 2, 2018 at 11:33 AM, Ben Gamari
Matt Peddie
writes: Hi George,
Not a stupid question. I don't have a single source at hand, but I think I read in a few places on the wiki that calling out to the system math library is not an option due to the variety of system math libraries on the platforms GHC supports. It'd be great if I got the wrong impression and this could just be a call to C. Can anyone set me straight on this point?
Indeed it's not a stupid question at all. Indeed this is precisely what we do for the simpler transcendentals (e.g. sin, asin, log). We very well could move in this direction in the case of asinh/atanh as well. I believe the reason we don't currently is that atanh was only standardized in C99, which we only started requiring a few releases ago. Perhaps this is ultimately the right direction.
Cheers,
- Ben

I'd be willing to do this.
--
Best wishes,
Artem
On Thu, 2 Aug 2018, 04:38 Matt Peddie,
Thanks, Ben, for chiming in. I think calling out to C for these functions is the way to go if it's now feasible. (Calling out to libm is the workaround I'm using in the application that led me to discover the inaccuracy.)
Matt
On Thu, Aug 2, 2018 at 11:33 AM, Ben Gamari
wrote: Matt Peddie
writes: Hi George,
Not a stupid question. I don't have a single source at hand, but I think I read in a few places on the wiki that calling out to the system math library is not an option due to the variety of system math libraries on the platforms GHC supports. It'd be great if I got the wrong impression and this could just be a call to C. Can anyone set me straight on this point?
Indeed it's not a stupid question at all. Indeed this is precisely what we do for the simpler transcendentals (e.g. sin, asin, log). We very well could move in this direction in the case of asinh/atanh as well. I believe the reason we don't currently is that atanh was only standardized in C99, which we only started requiring a few releases ago. Perhaps this is ultimately the right direction.
Cheers,
- Ben
ghc-devs mailing list ghc-devs@haskell.org http://mail.haskell.org/cgi-bin/mailman/listinfo/ghc-devs

Here is the patch: https://phabricator.haskell.org/D5034
--
Best, Artem
On Thu, 2 Aug 2018 at 06:26 Artem Pelenitsyn
I'd be willing to do this.
-- Best wishes, Artem
On Thu, 2 Aug 2018, 04:38 Matt Peddie,
wrote: Thanks, Ben, for chiming in. I think calling out to C for these functions is the way to go if it's now feasible. (Calling out to libm is the workaround I'm using in the application that led me to discover the inaccuracy.)
Matt
On Thu, Aug 2, 2018 at 11:33 AM, Ben Gamari
wrote: Matt Peddie
writes: Hi George,
Not a stupid question. I don't have a single source at hand, but I think I read in a few places on the wiki that calling out to the system math library is not an option due to the variety of system math libraries on the platforms GHC supports. It'd be great if I got the wrong impression and this could just be a call to C. Can anyone set me straight on this point?
Indeed it's not a stupid question at all. Indeed this is precisely what we do for the simpler transcendentals (e.g. sin, asin, log). We very well could move in this direction in the case of asinh/atanh as well. I believe the reason we don't currently is that atanh was only standardized in C99, which we only started requiring a few releases ago. Perhaps this is ultimately the right direction.
Cheers,
- Ben
ghc-devs mailing list ghc-devs@haskell.org http://mail.haskell.org/cgi-bin/mailman/listinfo/ghc-devs

Wasn't there a very recent commit to improve these functions, by leftaroundabout? On Thursday, August 2, 2018 8:16:10 AM EDT Artem Pelenitsyn wrote:
Here is the patch: https://phabricator.haskell.org/D5034
-- Best, Artem
On Thu, 2 Aug 2018 at 06:26 Artem Pelenitsyn
wrote: I'd be willing to do this.
-- Best wishes, Artem
On Thu, 2 Aug 2018, 04:38 Matt Peddie,
wrote: Thanks, Ben, for chiming in. I think calling out to C for these functions is the way to go if it's now feasible. (Calling out to libm is the workaround I'm using in the application that led me to discover the inaccuracy.)
Matt
On Thu, Aug 2, 2018 at 11:33 AM, Ben Gamari
wrote: Matt Peddie
writes: Hi George,
Not a stupid question. I don't have a single source at hand, but I think I read in a few places on the wiki that calling out to the system math library is not an option due to the variety of system math libraries on the platforms GHC supports. It'd be great if I got the wrong impression and this could just be a call to C. Can anyone set me straight on this point?
Indeed it's not a stupid question at all. Indeed this is precisely what we do for the simpler transcendentals (e.g. sin, asin, log). We very well could move in this direction in the case of asinh/atanh as well. I believe the reason we don't currently is that atanh was only standardized in C99, which we only started requiring a few releases ago. Perhaps this is ultimately the right direction.
Cheers,
- Ben
ghc-devs mailing list ghc-devs@haskell.org http://mail.haskell.org/cgi-bin/mailman/listinfo/ghc-devs
-- David Feuer Well-Typed

Thanks David! Indeed, here is the commit and ticket:
https://github.com/ghc/ghc/commit/3ea33411d7cbf32c20940cc72ca07df6830eeed7
https://ghc.haskell.org/trac/ghc/ticket/14927
This concerns only `asinh` though. The implementation is closer to what
Matt proposes in his package but simpler. Nevertheless, the original issue
about `Infinity` on large negative numbers seems to be fixed with this.
So, I guess, feel free to kill the patch.
--
Best, Artem
On Thu, 2 Aug 2018 at 23:19 David Feuer
Wasn't there a very recent commit to improve these functions, by leftaroundabout?
Here is the patch: https://phabricator.haskell.org/D5034
-- Best, Artem
On Thu, 2 Aug 2018 at 06:26 Artem Pelenitsyn
wrote: I'd be willing to do this.
-- Best wishes, Artem
On Thu, 2 Aug 2018, 04:38 Matt Peddie,
wrote: Thanks, Ben, for chiming in. I think calling out to C for these functions is the way to go if it's now feasible. (Calling out to libm is the workaround I'm using in the application that led me to discover the inaccuracy.)
Matt
On Thu, Aug 2, 2018 at 11:33 AM, Ben Gamari
wrote: Matt Peddie
writes: Hi George,
Not a stupid question. I don't have a single source at hand, but I think I read in a few places on the wiki that calling out to the system math library is not an option due to the variety of system math libraries on the platforms GHC supports. It'd be great if I got
On Thursday, August 2, 2018 8:16:10 AM EDT Artem Pelenitsyn wrote: the
wrong impression and this could just be a call to C. Can anyone set me straight on this point?
Indeed it's not a stupid question at all. Indeed this is precisely what we do for the simpler transcendentals (e.g. sin, asin, log). We very well could move in this direction in the case of asinh/atanh as well. I believe the reason we don't currently is that atanh was only standardized in C99, which we only started requiring a few releases ago. Perhaps this is ultimately the right direction.
Cheers,
- Ben
ghc-devs mailing list ghc-devs@haskell.org http://mail.haskell.org/cgi-bin/mailman/listinfo/ghc-devs
-- David Feuer Well-Typed

I believe there was but IMHO calling the libm library function is the
better solution, as discussed in the previous emails below. Less
maintenance, less testing, and possibly better performance, e.g. as
explained here
https://software.intel.com/en-us/articles/using-avx-without-writing-avx-code:
the C compiler used to produce libm may generate appropriate 128 and
256-bit Intel AVX VEX-encoded instructions, generating multiple,
processor-specific, auto-dispatched code paths when there is a performance
benefit. The most appropriate code would be executed at run time.
Cheers
George
On Thu, Aug 2, 2018 at 5:19 PM David Feuer
Wasn't there a very recent commit to improve these functions, by leftaroundabout?
Here is the patch: https://phabricator.haskell.org/D5034
-- Best, Artem
On Thu, 2 Aug 2018 at 06:26 Artem Pelenitsyn
wrote: I'd be willing to do this.
-- Best wishes, Artem
On Thu, 2 Aug 2018, 04:38 Matt Peddie,
wrote: Thanks, Ben, for chiming in. I think calling out to C for these functions is the way to go if it's now feasible. (Calling out to libm is the workaround I'm using in the application that led me to discover the inaccuracy.)
Matt
On Thu, Aug 2, 2018 at 11:33 AM, Ben Gamari
wrote: Matt Peddie
writes: Hi George,
Not a stupid question. I don't have a single source at hand, but I think I read in a few places on the wiki that calling out to the system math library is not an option due to the variety of system math libraries on the platforms GHC supports. It'd be great if I got
On Thursday, August 2, 2018 8:16:10 AM EDT Artem Pelenitsyn wrote: the
wrong impression and this could just be a call to C. Can anyone set me straight on this point?
Indeed it's not a stupid question at all. Indeed this is precisely what we do for the simpler transcendentals (e.g. sin, asin, log). We very well could move in this direction in the case of asinh/atanh as well. I believe the reason we don't currently is that atanh was only standardized in C99, which we only started requiring a few releases ago. Perhaps this is ultimately the right direction.
Cheers,
- Ben
ghc-devs mailing list ghc-devs@haskell.org http://mail.haskell.org/cgi-bin/mailman/listinfo/ghc-devs
-- David Feuer Well-Typed _______________________________________________ ghc-devs mailing list ghc-devs@haskell.org http://mail.haskell.org/cgi-bin/mailman/listinfo/ghc-devs
participants (7)
-
Artem Pelenitsyn
-
Ben Gamari
-
Ben Gamari
-
David Feuer
-
Edward Kmett
-
George Colpitts
-
Matt Peddie