mailing list of musl libc
 help / color / mirror / code / Atom feed
* [musl] [PATCH] math: avoid internal overflow in expm1l.c
@ 2021-05-10 20:53 Szabolcs Nagy
  0 siblings, 0 replies; only message in thread
From: Szabolcs Nagy @ 2021-05-10 20:53 UTC (permalink / raw)
  To: musl

the ld80 expm1l implementation in c tries to compute

  2^k q + 2^k - 1

but 2^k can overflow (when k==16384) while q is small and negative.
then the result should be finite, but the code gives nan. to avoid the
overflow use

  2 (2^(k-1) q + 2^(k-1) - 0.5).

note: this implementation is not used on x86 currently.
---
 src/math/expm1l.c | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/math/expm1l.c b/src/math/expm1l.c
index d1715078..79795dc0 100644
--- a/src/math/expm1l.c
+++ b/src/math/expm1l.c
@@ -109,9 +109,9 @@ long double expm1l(long double x)
 
 	/* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
 	 We have qx = exp(remainder ln 2) - 1, so
-	 exp(x) - 1  =  2^k (qx + 1) - 1  =  2^k qx + 2^k - 1.  */
-	px = scalbnl(1.0, k);
-	x = px * qx + (px - 1.0);
+	 exp(x) - 1  =  2^k (qx + 1) - 1  =  2 (2^(k-1) qx + 2^(k-1) - 0.5).  */
+	px = scalbnl(1.0, k-1);
+	x = 2*(px * qx + (px - 0.5));
 	return x;
 }
 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
-- 
2.29.2


^ permalink raw reply	[flat|nested] only message in thread

only message in thread, other threads:[~2021-05-10 20:53 UTC | newest]

Thread overview: (only message) (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2021-05-10 20:53 [musl] [PATCH] math: avoid internal overflow in expm1l.c Szabolcs Nagy

Code repositories for project(s) associated with this public inbox

	https://git.vuxu.org/mirror/musl/

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for NNTP newsgroup(s).