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&#39;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&#39;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&#39;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&#39;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==--