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 &lt;<a href=3D"mailt=
o:[email protected]">[email protected]</a>&gt; 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&#39=
;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&#39;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&#39;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&#39;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==--