1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
|
# 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
.global expl
.type expl,@function
expl:
fldt 4(%esp)
# special cases: 2*x is +-inf, nan or |x| < 0x1p-32
# check (exponent|0x8000)+2 < 0xbfff+2-32
movw 12(%esp), %ax
movw %ax, %dx
orw $0x8000, %dx
addw $2, %dx
cmpw $0xbfff-30, %dx
jnb 3f
cmpw $1, %dx
jbe 1f
# if |x|<0x1p-32 return 1+x
fld1
jmp 2f
1: testw %ax, %ax
jns 1f
# if 2*x == -inf,-nan return -0/x
fldz
fchs
fdivp
ret
# if 2*x == inf,nan return 2*x
1: fld %st(0)
2: faddp
ret
# should be 0x1.71547652b82fe178p0 == 0x3fff b8aa3b29 5c17f0bc
# it will be wrong on non-nearest rounding mode
3: fldl2e
# subl $32, %esp
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)
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) # x
# 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
|
