math: fix i386/expl.s with more precise x*log2e

with naive exp2l(x*log2e) the last 12bits of the result was incorrect
for x with large absolute value

with hi + lo = x*log2e is caluclated to 128 bits precision and then
  expl(x) = exp2l(hi) + exp2l(hi) * f2xm1(lo)
this gives <1.5ulp measured error everywhere in nearest rounding mode

2 files changed, 107 insertions, 7 deletions

diff --git a/src/math/i386/exp.s b/src/math/i386/exp.s
index c7f5ad0f..e3b42af5 100644
--- a/src/math/i386/exp.s
+++ b/src/math/i386/exp.s
@@ -50,12 +50,6 @@ expf:
 	flds 4(%esp)
 	jmp 2f
 
-.global expl
-.type expl,@function
-expl:
-	fldt 4(%esp)
-	jmp 2f
-
 .global exp
 .type exp,@function
 exp:
diff --git a/src/math/i386/expl.s b/src/math/i386/expl.s
index f335a3e5..3f2f707d 100644
--- a/src/math/i386/expl.s
+++ b/src/math/i386/expl.s
@@ -1 +1,107 @@
-# see exp.s
+# exp(x) = 2^hi + 2^hi (2^lo - 1)
+# where hi+lo = log2e*x with 128bit precision
+# exact log2e*x calculation depends on nearest rounding mode
+
+.global expl
+.type expl,@function
+expl:
+	fldt 4(%esp)
+
+		# special cases: 2*x is +-inf, nan or |x| < 0x1p-32
+		# check (exponent|0x8000)+2 < 0xbfff+2-32
+	movw 12(%esp), %ax
+	movw %ax, %dx
+	orw $0x8000, %dx
+	addw $2, %dx
+	cmpw $0xbfff-30, %dx
+	jnb 3f
+	cmpw $1, %dx
+	jbe 1f
+		# if |x|<0x1p-32 return 1+x
+	fld1
+	jmp 2f
+1:	testw %ax, %ax
+	jns 1f
+		# if 2*x == -inf,-nan return -0/x
+	fldz
+	fchs
+	fdivp
+	ret
+		# if 2*x == inf,nan return 2*x
+1:	fld %st(0)
+2:	faddp
+	ret
+
+		# should be 0x1.71547652b82fe178p0 == 0x3fff b8aa3b29 5c17f0bc
+		# it will be wrong on non-nearest rounding mode
+3:	fldl2e
+#	subl $32, %esp
+	subl $44, %esp
+		# hi = log2e_hi*x
+		# 2^hi = exp2l(hi)
+	fmul %st(1),%st
+	fld %st(0)
+	fstpt (%esp)
+	fstpt 16(%esp)
+	fstpt 32(%esp)
+	call exp2l
+		# if 2^hi == inf return 2^hi
+	fld %st(0)
+	fstpt (%esp)
+	cmpw $0x7fff, 8(%esp)
+	je 1f
+	fldt 32(%esp)
+	fldt 16(%esp)
+		# fpu stack: 2^hi x hi
+		# exact mult: x*log2e
+	fld %st(1) # x
+		# c = 0x1p32+1
+	pushl $0x41f00000
+	pushl $0x00100000
+	fldl (%esp)
+		# xh = x - c*x + c*x
+		# xl = x - xh
+	fmulp
+	fld %st(2)
+	fsub %st(1), %st
+	faddp
+	fld %st(2)
+	fsub %st(1), %st
+		# yh = log2e_hi - c*log2e_hi + c*log2e_hi
+	pushl $0x3ff71547
+	pushl $0x65200000
+	fldl (%esp)
+		# fpu stack: 2^hi x hi xh xl yh
+		# lo = hi - xh*yh + xl*yh
+	fld %st(2)
+	fmul %st(1), %st
+	fsubp %st, %st(4)
+	fmul %st(1), %st
+	faddp %st, %st(3)
+		# yl = log2e_hi - yh
+	pushl $0x3de705fc
+	pushl $0x2f000000
+	fldl (%esp)
+		# fpu stack: 2^hi x lo xh xl yl
+		# lo += xh*yl + xl*yl
+	fmul %st, %st(2)
+	fmulp %st, %st(1)
+	fxch %st(2)
+	faddp
+	faddp
+		# log2e_lo
+	pushl $0xbfbe
+	pushl $0x82f0025f
+	pushl $0x2dc582ee
+	fldt (%esp)
+	addl $36,%esp
+		# fpu stack: 2^hi x lo log2e_lo
+		# lo += log2e_lo*x
+		# return 2^hi + 2^hi (2^lo - 1)
+	fmulp %st, %st(2)
+	faddp
+	f2xm1
+	fmul %st(1), %st
+	faddp
+1:	addl $44, %esp
+	ret