1 # exp(x) = 2^hi + 2^hi (2^lo - 1) 2 # where hi+lo = log2e*x with 128bit precision 3 # exact log2e*x calculation depends on nearest rounding mode 4 # using the exact multiplication method of Dekker and Veltkamp 5 6 .global expl 7 .type expl,@function 8 expl: 9 fldt 8(%esp) 10 11 # interesting case: 0x1p-32 <= |x| < 16384 12 # check if (exponent|0x8000) is in [0xbfff-32, 0xbfff+13] 13 mov 16(%esp), %ax 14 or $0x8000, %ax 15 sub $0xbfdf, %ax 16 cmp $45, %ax 17 jbe 2f 18 test %ax, %ax 19 fld1 20 js 1f 21 # if |x|>=0x1p14 or nan return 2^trunc(x) 22 fscale 23 fstp %st(1) 24 ret 25 # if |x|<0x1p-32 return 1+x 26 1: faddp 27 ret 28 29 # should be 0x1.71547652b82fe178p0L == 0x3fff b8aa3b29 5c17f0bc 30 # it will be wrong on non-nearest rounding mode 31 2: fldl2e 32 sub $48, %esp 33 # hi = log2e_hi*x 34 # 2^hi = exp2l(hi) 35 fmul %st(1),%st 36 fld %st(0) 37 fstpt (%esp) 38 fstpt 16(%esp) 39 fstpt 32(%esp) 40 call exp2l@PLT 41 # if 2^hi == inf return 2^hi 42 fld %st(0) 43 fstpt (%esp) 44 cmpw $0x7fff, 8(%esp) 45 je 1f 46 fldt 32(%esp) 47 fldt 16(%esp) 48 # fpu stack: 2^hi x hi 49 # exact mult: x*log2e 50 fld %st(1) 51 # c = 0x1p32+1 52 movq $0x41f0000000100000,%rax 53 pushq %rax 54 fldl (%esp) 55 # xh = x - c*x + c*x 56 # xl = x - xh 57 fmulp 58 fld %st(2) 59 fsub %st(1), %st 60 faddp 61 fld %st(2) 62 fsub %st(1), %st 63 # yh = log2e_hi - c*log2e_hi + c*log2e_hi 64 movq $0x3ff7154765200000,%rax 65 pushq %rax 66 fldl (%esp) 67 # fpu stack: 2^hi x hi xh xl yh 68 # lo = hi - xh*yh + xl*yh 69 fld %st(2) 70 fmul %st(1), %st 71 fsubp %st, %st(4) 72 fmul %st(1), %st 73 faddp %st, %st(3) 74 # yl = log2e_hi - yh 75 movq $0x3de705fc2f000000,%rax 76 pushq %rax 77 fldl (%esp) 78 # fpu stack: 2^hi x lo xh xl yl 79 # lo += xh*yl + xl*yl 80 fmul %st, %st(2) 81 fmulp %st, %st(1) 82 fxch %st(2) 83 faddp 84 faddp 85 # log2e_lo 86 movq $0xbfbe,%rax 87 pushq %rax 88 movq $0x82f0025f2dc582ee,%rax 89 pushq %rax 90 fldt (%esp) 91 add $40,%esp 92 # fpu stack: 2^hi x lo log2e_lo 93 # lo += log2e_lo*x 94 # return 2^hi + 2^hi (2^lo - 1) 95 fmulp %st, %st(2) 96 faddp 97 f2xm1 98 fmul %st(1), %st 99 faddp 100 1: add $48, %esp 101 ret 102