| Message ID | CAEFO=4C+BLJYGHNwuCbjdeY2Xj5KJkBeRRCKmSxqTUqwVokkag@mail.gmail.com |
|---|---|
| State | New |
| Headers | show |
| Series | Add sinh(tanh(x)) and cosh(tanh(x)) rules | expand |
> On Aug 7, 2018, at 4:00 PM, Giuliano Augusto Faulin Belinassi <giuliano.belinassi@usp.br> wrote: > > Related with bug 86829, but for hyperbolic trigonometric functions. > This patch adds substitution rules to both sinh(tanh(x)) -> x / sqrt(1 > - x*x) and cosh(tanh(x)) -> 1 / sqrt(1 - x*x). Notice that the both > formulas has division by 0, but it causes no harm because 1/(+0) -> > +infinity, thus the math is still safe. What about non-IEEE targets that don't have "infinite" in their float representation? paul
That is a good question because I didn't know that such targets exists. Any suggestion? On Tue, Aug 7, 2018 at 5:29 PM, Paul Koning <paulkoning@comcast.net> wrote: > > >> On Aug 7, 2018, at 4:00 PM, Giuliano Augusto Faulin Belinassi <giuliano.belinassi@usp.br> wrote: >> >> Related with bug 86829, but for hyperbolic trigonometric functions. >> This patch adds substitution rules to both sinh(tanh(x)) -> x / sqrt(1 >> - x*x) and cosh(tanh(x)) -> 1 / sqrt(1 - x*x). Notice that the both >> formulas has division by 0, but it causes no harm because 1/(+0) -> >> +infinity, thus the math is still safe. > > What about non-IEEE targets that don't have "infinite" in their float representation? > > paul > >
Now I'm puzzled. I don't see how an infinite would show up in the original expression. I don't know hyperbolic functions, so I just constructed a small test program, and the original vs. the substitution you mention are not at all similar. paul > On Aug 7, 2018, at 4:42 PM, Giuliano Augusto Faulin Belinassi <giuliano.belinassi@usp.br> wrote: > > That is a good question because I didn't know that such targets > exists. Any suggestion? > > > On Tue, Aug 7, 2018 at 5:29 PM, Paul Koning <paulkoning@comcast.net> wrote: >> >> >>> On Aug 7, 2018, at 4:00 PM, Giuliano Augusto Faulin Belinassi <giuliano.belinassi@usp.br> wrote: >>> >>> Related with bug 86829, but for hyperbolic trigonometric functions. >>> This patch adds substitution rules to both sinh(tanh(x)) -> x / sqrt(1 >>> - x*x) and cosh(tanh(x)) -> 1 / sqrt(1 - x*x). Notice that the both >>> formulas has division by 0, but it causes no harm because 1/(+0) -> >>> +infinity, thus the math is still safe. >> >> What about non-IEEE targets that don't have "infinite" in their float representation? >> >> paul >> >>
Sorry about that. In the e-mail text field I wrote sinh(tanh(x)) and cosh(tanh(x)) where it was supposed to be sinh(atanh(x)) and cosh(atanh(x)), thus I am talking about the inverse hyperbolic tangent function. The patch code and comments are still correct. On Wed, Aug 8, 2018 at 10:58 AM, Paul Koning <paulkoning@comcast.net> wrote: > Now I'm puzzled. > > I don't see how an infinite would show up in the original expression. I don't know hyperbolic functions, so I just constructed a small test program, and the original vs. the substitution you mention are not at all similar. > > paul > > >> On Aug 7, 2018, at 4:42 PM, Giuliano Augusto Faulin Belinassi <giuliano.belinassi@usp.br> wrote: >> >> That is a good question because I didn't know that such targets >> exists. Any suggestion? >> >> >> On Tue, Aug 7, 2018 at 5:29 PM, Paul Koning <paulkoning@comcast.net> wrote: >>> >>> >>>> On Aug 7, 2018, at 4:00 PM, Giuliano Augusto Faulin Belinassi <giuliano.belinassi@usp.br> wrote: >>>> >>>> Related with bug 86829, but for hyperbolic trigonometric functions. >>>> This patch adds substitution rules to both sinh(tanh(x)) -> x / sqrt(1 >>>> - x*x) and cosh(tanh(x)) -> 1 / sqrt(1 - x*x). Notice that the both >>>> formulas has division by 0, but it causes no harm because 1/(+0) -> >>>> +infinity, thus the math is still safe. >>> >>> What about non-IEEE targets that don't have "infinite" in their float representation? >>> >>> paul >>> >>> >
Thanks. Ok, so the expressions you gave are undefined for x==1, which says that substituting something that is also undefined for x==1 is permitted. You can argue from "undefined" rather than relying on IEEE features like NaN or infinite. paul > On Aug 8, 2018, at 2:57 PM, Giuliano Augusto Faulin Belinassi <giuliano.belinassi@usp.br> wrote: > > Sorry about that. In the e-mail text field I wrote sinh(tanh(x)) and > cosh(tanh(x)) where it was supposed to be sinh(atanh(x)) and > cosh(atanh(x)), thus I am talking about the inverse hyperbolic tangent > function. The patch code and comments are still correct. > > On Wed, Aug 8, 2018 at 10:58 AM, Paul Koning <paulkoning@comcast.net> wrote: >> Now I'm puzzled. >> >> I don't see how an infinite would show up in the original expression. I don't know hyperbolic functions, so I just constructed a small test program, and the original vs. the substitution you mention are not at all similar. >> >> paul >> >> >>> On Aug 7, 2018, at 4:42 PM, Giuliano Augusto Faulin Belinassi <giuliano.belinassi@usp.br> wrote: >>> >>> That is a good question because I didn't know that such targets >>> exists. Any suggestion? >>> >>> >>> On Tue, Aug 7, 2018 at 5:29 PM, Paul Koning <paulkoning@comcast.net> wrote: >>>> >>>> >>>>> On Aug 7, 2018, at 4:00 PM, Giuliano Augusto Faulin Belinassi <giuliano.belinassi@usp.br> wrote: >>>>> >>>>> Related with bug 86829, but for hyperbolic trigonometric functions. >>>>> This patch adds substitution rules to both sinh(tanh(x)) -> x / sqrt(1 >>>>> - x*x) and cosh(tanh(x)) -> 1 / sqrt(1 - x*x). Notice that the both >>>>> formulas has division by 0, but it causes no harm because 1/(+0) -> >>>>> +infinity, thus the math is still safe. >>>> >>>> What about non-IEEE targets that don't have "infinite" in their float representation? >>>> >>>> paul >>>> >>>> >>
On 8/7/18 2:00 PM, Giuliano Augusto Faulin Belinassi wrote: > Related with bug 86829, but for hyperbolic trigonometric functions. > This patch adds substitution rules to both sinh(tanh(x)) -> x / sqrt(1 > - x*x) and cosh(tanh(x)) -> 1 / sqrt(1 - x*x). Notice that the both > formulas has division by 0, but it causes no harm because 1/(+0) -> > +infinity, thus the math is still safe. > > Changelog: > 2018-08-07 Giuliano Belinassi <giuliano.belinassi@usp.br> > > * match.pd: add simplification rules to sinh(atanh(x)) and cosh(atanh(x)). > > All tests added by this patch runs without errors in trunk, however, > there are tests unrelated with this patch that fails in my x86_64 > Ubuntu 18.04. I think these are going to need similar handling because the x*x can overflow. Are the domains constrained in a way that is helpful? jeff
Hello.
I don't think there is a need for overflow handling here because
the argument is bound by the argument of the sqrt function :-)
Since we have to compute sqrt (1 - x*x), the input is only valid
if 1 - x*x >= 0, implying that -1 <= x <= 1. For any x outside of this
set, the sqrt will return a invalid value, as imaginary numbers are
required to represent the answer.
One can also argue a problem regarding division by 0, however in
the extremes x -> -1 by the right and x -> 1 by the left we have:
sinh(atanh(-1)) = -1 / sqrt (0) = -inf
sinh(atanh( 1)) = 1 / sqrt (0) = +inf
cos(atanh(-1)) = 1 / sqrt (0) = +inf
cos(atanh( 1)) = 1 / sqrt (0) = +inf
Therefore it seems that the target has to support infinity anyway.
Well, I think I can take a look about how glibc handles such cases
on targets where infinity is not supported to try to keep
compatibility, but I think this is safe :-).
On Fri, Oct 12, 2018 at 1:09 AM Jeff Law <law@redhat.com> wrote:
>
> On 8/7/18 2:00 PM, Giuliano Augusto Faulin Belinassi wrote:
> > Related with bug 86829, but for hyperbolic trigonometric functions.
> > This patch adds substitution rules to both sinh(tanh(x)) -> x / sqrt(1
> > - x*x) and cosh(tanh(x)) -> 1 / sqrt(1 - x*x). Notice that the both
> > formulas has division by 0, but it causes no harm because 1/(+0) ->
> > +infinity, thus the math is still safe.
> >
> > Changelog:
> > 2018-08-07 Giuliano Belinassi <giuliano.belinassi@usp.br>
> >
> > * match.pd: add simplification rules to sinh(atanh(x)) and cosh(atanh(x)).
> >
> > All tests added by this patch runs without errors in trunk, however,
> > there are tests unrelated with this patch that fails in my x86_64
> > Ubuntu 18.04.
> I think these are going to need similar handling because the x*x can
> overflow. Are the domains constrained in a way that is helpful?
>
> jeff
On 10/12/18 8:36 AM, Giuliano Augusto Faulin Belinassi wrote: > Hello. > I don't think there is a need for overflow handling here because > the argument is bound by the argument of the sqrt function :-) Yea, I guess you're right. The domain of arctanh is -1 to 1, so I guess we're safe there. Except for the case where the input is -1 or 1 in which case I think you just set the output to +- INF as appropriate. Hmm, do we have problems as we get close to -1 or 1 where the outputs of the two forms might diverge? Jeff
> Hmm, do we have problems as we get close to -1 or 1 where the outputs of > the two forms might diverge? Well, I did some minor testing with that with input x around nextafter(1, -1); There are a minor imprecision when comparing directly with sinh(atanh(x)) and cosh(atanh(x)). * On 32-bits floats, for such x the error is about 10^-4 * On 64-bits floats, for such x the error is about 10^-7 * On 80-bits floats, for such x the error is about 10^-9 here are the code that I used for the test: https://pastebin.com/JzYZyigQ I can create a testcase based on this if needed :-)
On 10/17/18 3:25 PM, Giuliano Augusto Faulin Belinassi wrote: >> Hmm, do we have problems as we get close to -1 or 1 where the outputs of >> the two forms might diverge? > > Well, I did some minor testing with that with input x around nextafter(1, -1); > There are a minor imprecision when comparing directly with > sinh(atanh(x)) and cosh(atanh(x)). > * On 32-bits floats, for such x the error is about 10^-4 > * On 64-bits floats, for such x the error is about 10^-7 > * On 80-bits floats, for such x the error is about 10^-9 > > here are the code that I used for the test: https://pastebin.com/JzYZyigQ > > I can create a testcase based on this if needed :-) My gut instinct is those errors are too significant in practice. It also just occurred to me that we may have problems as X approaches X from either direction. Clearly when x^2 is indistinguishable from 0 or 1, then the result has to be +-0 or +-1. But I'm not sure if figuring out where those points are is sufficient to avoid the imprecisions noted above. This is *well* outside my areas of expertise. jeff
Oh, please note that the error that I'm talking about is the comparison with the result obtained before and after the simplification. It is possible that the result obtained after the simplification be more precise when compared to an arbitrary precise value (example, a 30 digits precise approximation). Well, I will try check that. But yes, with regard to compatibility this may be a problem. On Wed, Oct 17, 2018 at 6:42 PM Jeff Law <law@redhat.com> wrote: > > On 10/17/18 3:25 PM, Giuliano Augusto Faulin Belinassi wrote: > >> Hmm, do we have problems as we get close to -1 or 1 where the outputs of > >> the two forms might diverge? > > > > Well, I did some minor testing with that with input x around nextafter(1, -1); > > There are a minor imprecision when comparing directly with > > sinh(atanh(x)) and cosh(atanh(x)). > > * On 32-bits floats, for such x the error is about 10^-4 > > * On 64-bits floats, for such x the error is about 10^-7 > > * On 80-bits floats, for such x the error is about 10^-9 > > > > here are the code that I used for the test: https://pastebin.com/JzYZyigQ > > > > I can create a testcase based on this if needed :-) > My gut instinct is those errors are too significant in practice. > > It also just occurred to me that we may have problems as X approaches X > from either direction. > > Clearly when x^2 is indistinguishable from 0 or 1, then the result has > to be +-0 or +-1. But I'm not sure if figuring out where those points > are is sufficient to avoid the imprecisions noted above. This is *well* > outside my areas of expertise. > > jeff
On 10/17/18 4:21 PM, Giuliano Augusto Faulin Belinassi wrote: > Oh, please note that the error that I'm talking about is the > comparison with the result obtained before and after the > simplification. It is possible that the result obtained after the > simplification be more precise when compared to an arbitrary precise > value (example, a 30 digits precise approximation). Well, I will try > check that. That would be helpful. Obviously if we're getting more precise, then that's a good thing :-) jeff
On 10/18, Jeff Law wrote: > On 10/17/18 4:21 PM, Giuliano Augusto Faulin Belinassi wrote: > > Oh, please note that the error that I'm talking about is the > > comparison with the result obtained before and after the > > simplification. It is possible that the result obtained after the > > simplification be more precise when compared to an arbitrary precise > > value (example, a 30 digits precise approximation). Well, I will try > > check that. > That would be helpful. Obviously if we're getting more precise, then > that's a good thing :-) > > jeff Well, I compared the results before and after the simplifications with a 512-bit precise mpfr value. Unfortunately, I found that sometimes the error is very noticeable :-( . For example, using floats and comparing with a 512 precision mpfr calculation with input : = 9.99966979026794433593750000000000000000000000000000e-01 cosh: before : = 1.23053413391113281250000000000000000000000000000000e+02 cosh: after : = 1.23052398681640625000000000000000000000000000000000e+02 cosh: mpfr512: = 1.23053409952258504358633865742873246642102963529577e+02 error before : = 3.43885477689136613425712675335789703647042270993727e-06 error after : = 1.01127061787935863386574287324664210296352957729006e-03 There are also some significant loss of precision with long doubles: with input : = 9.99999999999996799706237365912286918501195032149553e-01 cosh: before : = 1.24994262843556815705596818588674068450927734375000e+07 cosh: after : = 1.24994262843556715697559411637485027313232421875000e+07 cosh: mpfr512: = 1.24994262843556815704069193408098058772318248178348e+07 error before : = 1.52762518057600967860948619665184971612393688891101e-13 error after : = 1.00006509781770613031459085826303348150283876063111e-08 So yes, precision may be a problem here.
Index: gcc/match.pd =================================================================== --- gcc/match.pd (revisão 263359) +++ gcc/match.pd (cópia de trabalho) @@ -4219,6 +4219,25 @@ (mult:c (TAN:s @0) (COS:s @0)) (SIN @0)) + /* Simplify sinh(atanh(x)) -> x / sqrt(1 - x*x). */ + (for sins (SINH) + atans (ATANH) + sqrts (SQRT) + (simplify + (sins (atans:s @0)) + (rdiv @0 (sqrts (minus {build_one_cst (type);} + (mult @0 @0)))))) + + /* Simplify cosh(atanh(x)) -> 1 / sqrt(1 - x*x). */ + (for coss (COSH) + atans (ATANH) + sqrts (SQRT) + (simplify + (coss (atans:s @0)) + (rdiv {build_one_cst (type);} + (sqrts (minus {build_one_cst (type);} + (mult @0 @0)))))) + /* Simplify x * pow(x,c) -> pow(x,c+1). */ (simplify (mult:c @0 (POW:s @0 REAL_CST@1)) Index: gcc/testsuite/gcc.dg/sinhtanh-1.c =================================================================== --- gcc/testsuite/gcc.dg/sinhtanh-1.c (nonexistent) +++ gcc/testsuite/gcc.dg/sinhtanh-1.c (cópia de trabalho) @@ -0,0 +1,15 @@ +/* { dg-do compile } */ +/* { dg-options "-O3 -ffast-math -fdump-tree-optimized" } */ + +extern double sinh(double x); +extern double atanh(double x); + +double __attribute__ ((noinline)) +sinhatanh_(double x) +{ + return sinh(atanh(x)); +} + +/* There should be no calls to sinh nor atanh */ +/* { dg-final { scan-tree-dump-not "sinh " "optimized" } } */ +/* { dg-final { scan-tree-dump-not "atanh " "optimized" } } */ Index: gcc/testsuite/gcc.dg/sinhtanh-2.c =================================================================== --- gcc/testsuite/gcc.dg/sinhtanh-2.c (nonexistent) +++ gcc/testsuite/gcc.dg/sinhtanh-2.c (cópia de trabalho) @@ -0,0 +1,15 @@ +/* { dg-do compile } */ +/* { dg-options "-O3 -ffast-math -fdump-tree-optimized" } */ + +extern double cosh(double x); +extern double atanh(double x); + +double __attribute__ ((noinline)) +coshatanh_(double x) +{ + return cosh(atanh(x)); +} + +/* There should be no calls to cosh nor atanh */ +/* { dg-final { scan-tree-dump-not "cosh " "optimized" } } */ +/* { dg-final { scan-tree-dump-not "atanh " "optimized" } } */ Index: gcc/testsuite/gcc.dg/sinhtanh-3.c =================================================================== --- gcc/testsuite/gcc.dg/sinhtanh-3.c (nonexistent) +++ gcc/testsuite/gcc.dg/sinhtanh-3.c (cópia de trabalho) @@ -0,0 +1,16 @@ +/* { dg-do compile } */ +/* { dg-options "-O3 -ffast-math -fdump-tree-optimized" } */ + +extern double sinh(double x); +extern double atanh(double x); + +double __attribute__ ((noinline)) +sinhatanh_(double x) +{ + double atgh = atanh(x); + return sinh(atgh) + atgh; +} + +/* There should be calls to both sinh and atanh */ +/* { dg-final { scan-tree-dump "sinh " "optimized" } } */ +/* { dg-final { scan-tree-dump "atanh " "optimized" } } */ Index: gcc/testsuite/gcc.dg/sinhtanh-4.c =================================================================== --- gcc/testsuite/gcc.dg/sinhtanh-4.c (nonexistent) +++ gcc/testsuite/gcc.dg/sinhtanh-4.c (cópia de trabalho) @@ -0,0 +1,16 @@ +/* { dg-do compile } */ +/* { dg-options "-O3 -ffast-math -fdump-tree-optimized" } */ + +extern double cosh(double x); +extern double atanh(double x); + +double __attribute__ ((noinline)) +coshatanh_(double x) +{ + double atgh = atanh(x); + return cosh(atgh) + atgh; +} + +/* There should be calls to both cosh and atanh */ +/* { dg-final { scan-tree-dump "cosh " "optimized" } } */ +/* { dg-final { scan-tree-dump "atanh " "optimized" } } */
Related with bug 86829, but for hyperbolic trigonometric functions. This patch adds substitution rules to both sinh(tanh(x)) -> x / sqrt(1 - x*x) and cosh(tanh(x)) -> 1 / sqrt(1 - x*x). Notice that the both formulas has division by 0, but it causes no harm because 1/(+0) -> +infinity, thus the math is still safe. Changelog: 2018-08-07 Giuliano Belinassi <giuliano.belinassi@usp.br> * match.pd: add simplification rules to sinh(atanh(x)) and cosh(atanh(x)). All tests added by this patch runs without errors in trunk, however, there are tests unrelated with this patch that fails in my x86_64 Ubuntu 18.04.