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