[musl] [PATCH] math: remove the x86 specific expm1l with accuracy problem

[musl] [PATCH] math: remove the x86 specific expm1l with accuracy problem

[Szabolcs Nagy]

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