summaryrefslogtreecommitdiff
path: root/src/math/i386/exp.s
blob: abb9036957ad04bf9c4579ec5aa3251e69165e3f (plain) (blame)
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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
.global expm1f
.type expm1f,@function
expm1f:
	flds 4(%esp)
	mov 4(%esp),%eax
	add %eax,%eax
	cmp $0x01000000,%eax
	jae 1f
		# subnormal x, return x with underflow
	fnstsw %ax
	and $16,%ax
	jnz 2f
	fld %st(0)
	fmul %st(1)
	fstps 4(%esp)
2:	ret

.global expm1l
.type expm1l,@function
expm1l:
	fldt 4(%esp)
	jmp 1f

.global expm1
.type expm1,@function
expm1:
	fldl 4(%esp)
	mov 8(%esp),%eax
	add %eax,%eax
	cmp $0x00200000,%eax
	jae 1f
		# subnormal x, return x with underflow
	fnstsw %ax
	and $16,%ax
	jnz 2f
	fsts 4(%esp)
2:	ret
1:	fldl2e
	fmulp
	mov $0xc2820000,%eax
	push %eax
	flds (%esp)
	pop %eax
	fucomp %st(1)
	fnstsw %ax
	sahf
	fld1
	jb 1f
		# x*log2e < -65, return -1 without underflow
	fstp %st(1)
	fchs
	ret
1:	fld %st(1)
	fabs
	fucom %st(1)
	fnstsw %ax
	fstp %st(0)
	fstp %st(0)
	sahf
	ja 1f
	f2xm1
	ret
1:	call 1f
	fld1
	fsubrp
	ret

.global exp2f
.type exp2f,@function
exp2f:
	flds 4(%esp)
	jmp 1f

.global exp2l
.type exp2l,@function
exp2l:
	fldt 4(%esp)
	jmp 1f

.global expf
.type expf,@function
expf:
	flds 4(%esp)
	jmp 2f

.global exp
.type exp,@function
exp:
	fldl 4(%esp)
2:	fldl2e
	fmulp
	jmp 1f

.global exp2
.type exp2,@function
exp2:
	fldl 4(%esp)
1:	sub $12,%esp
	fld %st(0)
	fstpt (%esp)
	mov 8(%esp),%ax
	and $0x7fff,%ax
	cmp $0x3fff+13,%ax
	jb 4f             # |x| < 8192
	cmp $0x3fff+15,%ax
	jae 3f            # |x| >= 32768
	fsts (%esp)
	cmpl $0xc67ff800,(%esp)
	jb 2f             # x > -16382
	movl $0x5f000000,(%esp)
	flds (%esp)       # 0x1p63
	fld %st(1)
	fsub %st(1)
	faddp
	fucomp %st(1)
	fnstsw
	sahf
	je 2f             # x - 0x1p63 + 0x1p63 == x
	movl $1,(%esp)
	flds (%esp)       # 0x1p-149
	fdiv %st(1)
	fstps (%esp)      # raise underflow
2:	fld1
	fld %st(1)
	frndint
	fxch %st(2)
	fsub %st(2)       # st(0)=x-rint(x), st(1)=1, st(2)=rint(x)
	f2xm1
	faddp             # 2^(x-rint(x))
1:	fscale
	fstp %st(1)
	add $12,%esp
	ret
3:	xor %eax,%eax
4:	cmp $0x3fff-64,%ax
	fld1
	jb 1b             # |x| < 0x1p-64
	fstpt (%esp)
	fistl 8(%esp)
	fildl 8(%esp)
	fsubrp %st(1)
	addl $0x3fff,8(%esp)
	f2xm1
	fld1
	faddp             # 2^(x-rint(x))
	fldt (%esp)       # 2^rint(x)
	fmulp
	add $12,%esp
	ret