summaryrefslogtreecommitdiff
path: root/src/math/x32/expl.s
diff options
context:
space:
mode:
Diffstat (limited to 'src/math/x32/expl.s')
-rw-r--r--src/math/x32/expl.s101
1 files changed, 101 insertions, 0 deletions
diff --git a/src/math/x32/expl.s b/src/math/x32/expl.s
new file mode 100644
index 00000000..3add810d
--- /dev/null
+++ b/src/math/x32/expl.s
@@ -0,0 +1,101 @@
+# 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 8(%rsp)
+
+ # interesting case: 0x1p-32 <= |x| < 16384
+ # check if (exponent|0x8000) is in [0xbfff-32, 0xbfff+13]
+ mov 16(%rsp), %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
+ subq $48, %rsp
+ # hi = log2e_hi*x
+ # 2^hi = exp2l(hi)
+ fmul %st(1),%st
+ fld %st(0)
+ fstpt (%rsp)
+ fstpt 16(%rsp)
+ fstpt 32(%rsp)
+ call exp2l
+ # if 2^hi == inf return 2^hi
+ fld %st(0)
+ fstpt (%rsp)
+ cmpw $0x7fff, 8(%rsp)
+ je 1f
+ fldt 32(%rsp)
+ fldt 16(%rsp)
+ # fpu stack: 2^hi x hi
+ # exact mult: x*log2e
+ fld %st(1)
+ # c = 0x1p32+1
+ movq $0x41f0000000100000,%rax
+ pushq %rax
+ fldl (%rsp)
+ # 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
+ movq $0x3ff7154765200000,%rax
+ pushq %rax
+ fldl (%rsp)
+ # 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
+ movq $0x3de705fc2f000000,%rax
+ pushq %rax
+ fldl (%rsp)
+ # 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
+ movq $0xbfbe,%rax
+ pushq %rax
+ movq $0x82f0025f2dc582ee,%rax
+ pushq %rax
+ fldt (%rsp)
+ addq $40,%rsp
+ # 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: addq $48, %rsp
+ ret