Re: Unexpected overflow in acosh and asinh
Raymond Toy <[email protected]> Mon, 10 Jan 2022 09:05:18 -0800
| Newsgroups | gmane.lisp.clisp.devel |
|---|---|
| Message-ID | <CAG14z1FbG1Qz8gcejrJRYyJoep2djbF8j+B9HCZ6_L1todxvew@mail.gmail.com> |
--===============2951907146292039776== Content-Type: multipart/alternative; boundary="00000000000004037705d53d58af" --00000000000004037705d53d58af Content-Type: text/plain; charset="UTF-8" One more thing I just noticed. (acosh 1d50) gives a division by 0. But asinh gives the expected answer. On Mon, Jan 10, 2022 at 8:59 AM Raymond Toy <[email protected]> wrote: > While working on some maxima stuff, I found out that clisp can't compute > acosh(most-positive-double-float). I get an overflow. Same happens with > asinh. The actual answer should be about 710.475. > > I'm guessing this happens because clisp is using the log definition for > these functions. For example, asinh(x) = log(x+sqrt(x^2+1)). > > So for any value above sqrt(most-positive), there's an overflow. > > In this case, for large enough x, sqrt(x^2+1) = x so the final answer is > log(2*x). However, that would overflow if x is greater than > most-positive/2. In that case, we can use log(x)+log(x). > > This will allow clisp to compute acosh and asinh for the entire range of > floats as expected. > > Oh, just noticed that clisp also fails to compute sinh(710.475d0). This > probably happens because sinh(x) is probably computed using > (exp(x)-exp(-x))/2. exp(710.475) overflows. One possible way to avoid > this is to compute exp(x/2)/sqrt(2). This won't overflow. Then square the > result to get the final answer. > > -- > Ray > -- Ray --00000000000004037705d53d58af Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"ltr"><div class=3D"gmail_default" style=3D"font-family:arial,he= lvetica,sans-serif;font-size:small">One more thing I just noticed.=C2=A0 (a= cosh 1d50) gives a division by 0.=C2=A0 But asinh gives the expected answer= .<br></div></div><br><div class=3D"gmail_quote"><div dir=3D"ltr" class=3D"g= mail_attr">On Mon, Jan 10, 2022 at 8:59 AM Raymond Toy <<a href=3D"mailt= o:[email protected]">[email protected]</a>> wrote:<br></div><blo= ckquote class=3D"gmail_quote" style=3D"margin:0px 0px 0px 0.8ex;border-left= :1px solid rgb(204,204,204);padding-left:1ex"><div dir=3D"ltr"><div class= =3D"gmail_default" style=3D"font-family:arial,helvetica,sans-serif;font-siz= e:small">While working on some maxima stuff, I found out that clisp can'= ;t compute acosh(most-positive-double-float).=C2=A0 I get an overflow.=C2= =A0 Same happens with asinh.=C2=A0 The actual answer should be about 710.47= 5.<br></div><div class=3D"gmail_default" style=3D"font-family:arial,helveti= ca,sans-serif;font-size:small"><br></div><div class=3D"gmail_default" style= =3D"font-family:arial,helvetica,sans-serif;font-size:small">I'm guessin= g this happens because clisp is using the log definition for these function= s.=C2=A0 For example, asinh(x) =3D log(x+sqrt(x^2+1)).</div><div class=3D"g= mail_default" style=3D"font-family:arial,helvetica,sans-serif;font-size:sma= ll"><br></div><div class=3D"gmail_default" style=3D"font-family:arial,helve= tica,sans-serif;font-size:small">So for any value above sqrt(most-positive)= , there's an overflow.</div><div class=3D"gmail_default" style=3D"font-= family:arial,helvetica,sans-serif;font-size:small"><br></div><div class=3D"= gmail_default" style=3D"font-family:arial,helvetica,sans-serif;font-size:sm= all">In this case, for large enough x, sqrt(x^2+1) =3D x so the final answe= r is log(2*x).=C2=A0 However, that would overflow if x is greater than most= -positive/2.=C2=A0 In that case, we can use log(x)+log(x).</div><div class= =3D"gmail_default" style=3D"font-family:arial,helvetica,sans-serif;font-siz= e:small"><br></div><div class=3D"gmail_default" style=3D"font-family:arial,= helvetica,sans-serif;font-size:small">This will allow clisp to compute acos= h and asinh for the entire range of floats as expected.</div><div class=3D"= gmail_default" style=3D"font-family:arial,helvetica,sans-serif;font-size:sm= all"><br></div><div class=3D"gmail_default" style=3D"font-family:arial,helv= etica,sans-serif;font-size:small">Oh, just noticed that clisp also fails to= compute sinh(710.475d0).=C2=A0 This probably happens because sinh(x) is pr= obably computed using (exp(x)-exp(-x))/2.=C2=A0 exp(710.475) overflows.=C2= =A0 One possible way to avoid this is to compute exp(x/2)/sqrt(2).=C2=A0 Th= is won't overflow.=C2=A0 Then square the result to get the final answer= .<br></div><br>-- <br><div dir=3D"ltr"><div dir=3D"ltr"><div></div>Ray<br><= /div></div></div> </blockquote></div><br clear=3D"all"><br>-- <br><div dir=3D"ltr" class=3D"g= mail_signature"><div dir=3D"ltr"><div></div>Ray<br></div></div> --00000000000004037705d53d58af-- --===============2951907146292039776== Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline --===============2951907146292039776== Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline _______________________________________________ clisp-devel mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/clisp-devel --===============2951907146292039776==--