[<prev] [next>] [day] [month] [year] [list]
Date: Mon, 10 May 2021 22:54:04 +0200
From: Szabolcs Nagy <nsz@...t70.net>
To: musl@...ts.openwall.com
Subject: [PATCH] math: remove the x86 specific expm1l with accuracy problem
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
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
