From: Szabolcs Nagy <nsz@port70.net>
To: musl@lists.openwall.com
Subject: [musl] [PATCH] math: remove the x86 specific expm1l with accuracy problem
Date: Mon, 10 May 2021 22:54:04 +0200 [thread overview]
Message-ID: <20210510205404.GP2799122@port70.net> (raw)
the x86 expm1l used exp2l to compute e^x - 1 = 2^(x*log2(e)) - 1, but
the rounding error of the mul in the exponent is amplified by exp2l:
2^(x*log2(e)*(1+r)) = e^x e^(x*r) ~= e^x (1 + x*r)
so r relative error in the input causes x*r relative error in the
output, which can be large for large x.
see commits a8f73bb1a685dd7d67669c6f6ceb255cfa967790 and
58bba42d1bd14e1ab01f3249ffc98afdbf841a6a for similar issue in expl.
fixing the expm1l asm is not obvious (large x can use expl, but mid
range still needs precise mul and reconstructing the result from
multiple parts), maintaining x87 asm is error prone so the code is
just removed.
---
src/math/i386/exp_ld.s | 36 +-----------------------------------
src/math/i386/expm1l.s | 1 -
src/math/x86_64/exp2l.s | 31 +------------------------------
src/math/x86_64/expm1l.s | 1 -
4 files changed, 2 insertions(+), 67 deletions(-)
delete mode 100644 src/math/i386/expm1l.s
delete mode 100644 src/math/x86_64/expm1l.s
diff --git a/src/math/i386/exp_ld.s b/src/math/i386/exp_ld.s
index 99cba01f..9c782ad6 100644
--- a/src/math/i386/exp_ld.s
+++ b/src/math/i386/exp_ld.s
@@ -1,37 +1,3 @@
-.global expm1l
-.type expm1l,@function
-expm1l:
- fldt 4(%esp)
- fldl2e
- fmulp
- mov $0xc2820000,%eax
- push %eax
- flds (%esp)
- pop %eax
- fucomp %st(1)
- fnstsw %ax
- sahf
- fld1
- jb 1f
- # x*log2e < -65, return -1 without underflow
- fstp %st(1)
- fchs
- ret
-1: fld %st(1)
- fabs
- fucom %st(1)
- fnstsw %ax
- fstp %st(0)
- fstp %st(0)
- sahf
- ja 1f
- f2xm1
- ret
-1: call 1f
- fld1
- fsubrp
- ret
-
.global exp2l
.global __exp2l
.hidden __exp2l
@@ -39,7 +5,7 @@ expm1l:
exp2l:
__exp2l:
fldt 4(%esp)
-1: sub $12,%esp
+ sub $12,%esp
fld %st(0)
fstpt (%esp)
mov 8(%esp),%ax
diff --git a/src/math/i386/expm1l.s b/src/math/i386/expm1l.s
deleted file mode 100644
index 8125761d..00000000
--- a/src/math/i386/expm1l.s
+++ /dev/null
@@ -1 +0,0 @@
-# see exp_ld.s
diff --git a/src/math/x86_64/exp2l.s b/src/math/x86_64/exp2l.s
index effab2bd..a48fb0f9 100644
--- a/src/math/x86_64/exp2l.s
+++ b/src/math/x86_64/exp2l.s
@@ -1,37 +1,8 @@
-.global expm1l
-.type expm1l,@function
-expm1l:
- fldt 8(%rsp)
- fldl2e
- fmulp
- movl $0xc2820000,-4(%rsp)
- flds -4(%rsp)
- fucomip %st(1),%st
- fld1
- jb 1f
- # x*log2e <= -65, return -1 without underflow
- fstp %st(1)
- fchs
- ret
-1: fld %st(1)
- fabs
- fucomip %st(1),%st
- fstp %st(0)
- ja 1f
- f2xm1
- ret
-1: push %rax
- call 1f
- pop %rax
- fld1
- fsubrp
- ret
-
.global exp2l
.type exp2l,@function
exp2l:
fldt 8(%rsp)
-1: fld %st(0)
+ fld %st(0)
sub $16,%rsp
fstpt (%rsp)
mov 8(%rsp),%ax
diff --git a/src/math/x86_64/expm1l.s b/src/math/x86_64/expm1l.s
deleted file mode 100644
index e773f080..00000000
--- a/src/math/x86_64/expm1l.s
+++ /dev/null
@@ -1 +0,0 @@
-# see exp2l.s
--
2.29.2
reply other threads:[~2021-05-10 20:54 UTC|newest]
Thread overview: [no followups] expand[flat|nested] mbox.gz Atom feed
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=20210510205404.GP2799122@port70.net \
--to=nsz@port70.net \
--cc=musl@lists.openwall.com \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
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).