Index: fusl/src/math/i386/expl.s |
diff --git a/fusl/src/math/i386/expl.s b/fusl/src/math/i386/expl.s |
deleted file mode 100644 |
index b5124e8f1f5f4565de61d6221ff9652d9f60260b..0000000000000000000000000000000000000000 |
--- a/fusl/src/math/i386/expl.s |
+++ /dev/null |
@@ -1,101 +0,0 @@ |
-# 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 |
-# using the exact multiplication method of Dekker and Veltkamp |
- |
-.global expl |
-.type expl,@function |
-expl: |
- fldt 4(%esp) |
- |
- # interesting case: 0x1p-32 <= |x| < 16384 |
- # check if (exponent|0x8000) is in [0xbfff-32, 0xbfff+13] |
- mov 12(%esp), %ax |
- or $0x8000, %ax |
- sub $0xbfdf, %ax |
- cmp $45, %ax |
- jbe 2f |
- test %ax, %ax |
- fld1 |
- js 1f |
- # if |x|>=0x1p14 or nan return 2^trunc(x) |
- fscale |
- fstp %st(1) |
- ret |
- # if |x|<0x1p-32 return 1+x |
-1: faddp |
- ret |
- |
- # should be 0x1.71547652b82fe178p0L == 0x3fff b8aa3b29 5c17f0bc |
- # it will be wrong on non-nearest rounding mode |
-2: fldl2e |
- 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) |
-.hidden __exp2l |
- 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) |
- # 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 |