# 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@PLT
                # 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