Unexpected overflow in acosh and asinh
Raymond Toy <[email protected]> Mon, 10 Jan 2022 08:59:49 -0800
| Newsgroups | gmane.lisp.clisp.devel |
|---|---|
| Message-ID | <CAG14z1G9Rb4j1tYcB_KOnnsJR8GraLMY7_uNLeEJX3QLSejQaQ@mail.gmail.com> |
--===============4437602794182144883== Content-Type: multipart/alternative; boundary="0000000000006426ad05d53d449a" --0000000000006426ad05d53d449a Content-Type: text/plain; charset="UTF-8" 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 --0000000000006426ad05d53d449a 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">While working on some maxima stuff, I f= ound 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 answ= er should be about 710.475.<br></div><div class=3D"gmail_default" style=3D"= font-family:arial,helvetica,sans-serif;font-size:small"><br></div><div clas= s=3D"gmail_default" style=3D"font-family:arial,helvetica,sans-serif;font-si= ze:small">I'm guessing this happens because clisp is using the log defi= nition for these functions.=C2=A0 For example, asinh(x) =3D log(x+sqrt(x^2+= 1)).</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:small">So for any valu= e above sqrt(most-positive), there's an overflow.</div><div class=3D"gm= ail_default" style=3D"font-family:arial,helvetica,sans-serif;font-size:smal= l"><br></div><div class=3D"gmail_default" style=3D"font-family:arial,helvet= ica,sans-serif;font-size:small">In this case, for large enough x, sqrt(x^2+= 1) =3D x so the final answer is log(2*x).=C2=A0 However, that would overflo= w 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,hel= vetica,sans-serif;font-size:small"><br></div><div class=3D"gmail_default" s= tyle=3D"font-family:arial,helvetica,sans-serif;font-size:small">This will a= llow clisp to compute acosh and asinh for the entire range of floats as exp= ected.</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">Oh, just notice= d that clisp also fails to compute sinh(710.475d0).=C2=A0 This probably hap= pens because sinh(x) is probably computed using (exp(x)-exp(-x))/2.=C2=A0 e= xp(710.475) overflows.=C2=A0 One possible way to avoid this is to compute e= xp(x/2)/sqrt(2).=C2=A0 This won't overflow.=C2=A0 Then square the resul= t to get the final answer.<br></div><br>-- <br><div dir=3D"ltr" class=3D"gm= ail_signature" data-smartmail=3D"gmail_signature"><div dir=3D"ltr"><div></d= iv>Ray<br></div></div></div> --0000000000006426ad05d53d449a-- --===============4437602794182144883== Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline --===============4437602794182144883== 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 --===============4437602794182144883==--