[CRT]
[reactos.git] / reactos / lib / sdk / crt / math / i386 / pow_asm.s
1 /* ix87 specific implementation of pow function.
2 Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004, 2005, 2007
3 Free Software Foundation, Inc.
4 This file is part of the GNU C Library.
5 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
6
7 The GNU C Library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public
9 License as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11
12 The GNU C Library is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 Lesser General Public License for more details.
16
17 You should have received a copy of the GNU Lesser General Public
18 License along with the GNU C Library; if not, write to the Free
19 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
20 02111-1307 USA. */
21
22 /* Reacros modifications */
23 #define ALIGNARG(log2) log2
24 #define ASM_TYPE_DIRECTIVE(name,typearg)
25 #define ASM_SIZE_DIRECTIVE(name)
26 #define cfi_adjust_cfa_offset(x)
27 #define ENTRY(x)
28 #define END(x)
29 .global _pow
30
31 .text
32
33 .align ALIGNARG(4)
34 ASM_TYPE_DIRECTIVE(infinity,@object)
35 inf_zero:
36 infinity:
37 .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
38 ASM_SIZE_DIRECTIVE(infinity)
39 ASM_TYPE_DIRECTIVE(zero,@object)
40 zero: .double 0.0
41 ASM_SIZE_DIRECTIVE(zero)
42 ASM_TYPE_DIRECTIVE(minf_mzero,@object)
43 minf_mzero:
44 minfinity:
45 .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
46 mzero:
47 .byte 0, 0, 0, 0, 0, 0, 0, 0x80
48 ASM_SIZE_DIRECTIVE(minf_mzero)
49 ASM_TYPE_DIRECTIVE(one,@object)
50 one: .double 1.0
51 ASM_SIZE_DIRECTIVE(one)
52 ASM_TYPE_DIRECTIVE(limit,@object)
53 limit: .double 0.29
54 ASM_SIZE_DIRECTIVE(limit)
55 ASM_TYPE_DIRECTIVE(p63,@object)
56 p63: .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
57 ASM_SIZE_DIRECTIVE(p63)
58
59 #ifdef PIC
60 #define MO(op) op##@GOTOFF(%ecx)
61 #define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
62 #else
63 #define MO(op) op
64 #define MOX(op,x,f) op(,x,f)
65 #endif
66
67 .text
68 _pow:
69 ENTRY(__ieee754_pow)
70 fldl 12(%esp) // y
71 fxam
72
73 #ifdef PIC
74 LOAD_PIC_REG (cx)
75 #endif
76
77 fnstsw
78 movb %ah, %dl
79 andb $0x45, %ah
80 cmpb $0x40, %ah // is y == 0 ?
81 je 11f
82
83 cmpb $0x05, %ah // is y == ±inf ?
84 je 12f
85
86 cmpb $0x01, %ah // is y == NaN ?
87 je 30f
88
89 fldl 4(%esp) // x : y
90
91 subl $8,%esp
92 cfi_adjust_cfa_offset (8)
93
94 fxam
95 fnstsw
96 movb %ah, %dh
97 andb $0x45, %ah
98 cmpb $0x40, %ah
99 je 20f // x is ±0
100
101 cmpb $0x05, %ah
102 je 15f // x is ±inf
103
104 fxch // y : x
105
106 /* fistpll raises invalid exception for |y| >= 1L<<63. */
107 fld %st // y : y : x
108 fabs // |y| : y : x
109 fcompl MO(p63) // y : x
110 fnstsw
111 sahf
112 jnc 2f
113
114 /* First see whether `y' is a natural number. In this case we
115 can use a more precise algorithm. */
116 fld %st // y : y : x
117 fistpll (%esp) // y : x
118 fildll (%esp) // int(y) : y : x
119 fucomp %st(1) // y : x
120 fnstsw
121 sahf
122 jne 2f
123
124 /* OK, we have an integer value for y. */
125 popl %eax
126 cfi_adjust_cfa_offset (-4)
127 popl %edx
128 cfi_adjust_cfa_offset (-4)
129 orl $0, %edx
130 fstp %st(0) // x
131 jns 4f // y >= 0, jump
132 fdivrl MO(one) // 1/x (now referred to as x)
133 negl %eax
134 adcl $0, %edx
135 negl %edx
136 4: fldl MO(one) // 1 : x
137 fxch
138
139 6: shrdl $1, %edx, %eax
140 jnc 5f
141 fxch
142 fmul %st(1) // x : ST*x
143 fxch
144 5: fmul %st(0), %st // x*x : ST*x
145 shrl $1, %edx
146 movl %eax, %ecx
147 orl %edx, %ecx
148 jnz 6b
149 fstp %st(0) // ST*x
150 ret
151
152 /* y is ±NAN */
153 30: fldl 4(%esp) // x : y
154 fldl MO(one) // 1.0 : x : y
155 fucomp %st(1) // x : y
156 fnstsw
157 sahf
158 je 31f
159 fxch // y : x
160 31: fstp %st(1)
161 ret
162
163 cfi_adjust_cfa_offset (8)
164 .align ALIGNARG(4)
165 2: /* y is a real number. */
166 fxch // x : y
167 fldl MO(one) // 1.0 : x : y
168 fldl MO(limit) // 0.29 : 1.0 : x : y
169 fld %st(2) // x : 0.29 : 1.0 : x : y
170 fsub %st(2) // x-1 : 0.29 : 1.0 : x : y
171 fabs // |x-1| : 0.29 : 1.0 : x : y
172 fucompp // 1.0 : x : y
173 fnstsw
174 fxch // x : 1.0 : y
175 sahf
176 ja 7f
177 fsub %st(1) // x-1 : 1.0 : y
178 fyl2xp1 // log2(x) : y
179 jmp 8f
180
181 7: fyl2x // log2(x) : y
182 8: fmul %st(1) // y*log2(x) : y
183 fst %st(1) // y*log2(x) : y*log2(x)
184 frndint // int(y*log2(x)) : y*log2(x)
185 fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x))
186 fxch // fract(y*log2(x)) : int(y*log2(x))
187 f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x))
188 faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x))
189 fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
190 addl $8, %esp
191 cfi_adjust_cfa_offset (-8)
192 fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x))
193 ret
194
195
196 // pow(x,±0) = 1
197 .align ALIGNARG(4)
198 11: fstp %st(0) // pop y
199 fldl MO(one)
200 ret
201
202 // y == ±inf
203 .align ALIGNARG(4)
204 12: fstp %st(0) // pop y
205 fldl MO(one) // 1
206 fldl 4(%esp) // x : 1
207 fabs // abs(x) : 1
208 fucompp // < 1, == 1, or > 1
209 fnstsw
210 andb $0x45, %ah
211 cmpb $0x45, %ah
212 je 13f // jump if x is NaN
213
214 cmpb $0x40, %ah
215 je 14f // jump if |x| == 1
216
217 shlb $1, %ah
218 xorb %ah, %dl
219 andl $2, %edx
220 fldl MOX(inf_zero, %edx, 4)
221 ret
222
223 .align ALIGNARG(4)
224 14: fldl MO(one)
225 ret
226
227 .align ALIGNARG(4)
228 13: fldl 4(%esp) // load x == NaN
229 ret
230
231 cfi_adjust_cfa_offset (8)
232 .align ALIGNARG(4)
233 // x is ±inf
234 15: fstp %st(0) // y
235 testb $2, %dh
236 jz 16f // jump if x == +inf
237
238 // We must find out whether y is an odd integer.
239 fld %st // y : y
240 fistpll (%esp) // y
241 fildll (%esp) // int(y) : y
242 fucompp // <empty>
243 fnstsw
244 sahf
245 jne 17f
246
247 // OK, the value is an integer, but is the number of bits small
248 // enough so that all are coming from the mantissa?
249 popl %eax
250 cfi_adjust_cfa_offset (-4)
251 popl %edx
252 cfi_adjust_cfa_offset (-4)
253 andb $1, %al
254 jz 18f // jump if not odd
255 movl %edx, %eax
256 orl %edx, %edx
257 jns 155f
258 negl %eax
259 155: cmpl $0x00200000, %eax
260 ja 18f // does not fit in mantissa bits
261 // It's an odd integer.
262 shrl $31, %edx
263 fldl MOX(minf_mzero, %edx, 8)
264 ret
265
266 cfi_adjust_cfa_offset (8)
267 .align ALIGNARG(4)
268 16: fcompl MO(zero)
269 addl $8, %esp
270 cfi_adjust_cfa_offset (-8)
271 fnstsw
272 shrl $5, %eax
273 andl $8, %eax
274 fldl MOX(inf_zero, %eax, 1)
275 ret
276
277 cfi_adjust_cfa_offset (8)
278 .align ALIGNARG(4)
279 17: shll $30, %edx // sign bit for y in right position
280 addl $8, %esp
281 cfi_adjust_cfa_offset (-8)
282 18: shrl $31, %edx
283 fldl MOX(inf_zero, %edx, 8)
284 ret
285
286 cfi_adjust_cfa_offset (8)
287 .align ALIGNARG(4)
288 // x is ±0
289 20: fstp %st(0) // y
290 testb $2, %dl
291 jz 21f // y > 0
292
293 // x is ±0 and y is < 0. We must find out whether y is an odd integer.
294 testb $2, %dh
295 jz 25f
296
297 fld %st // y : y
298 fistpll (%esp) // y
299 fildll (%esp) // int(y) : y
300 fucompp // <empty>
301 fnstsw
302 sahf
303 jne 26f
304
305 // OK, the value is an integer, but is the number of bits small
306 // enough so that all are coming from the mantissa?
307 popl %eax
308 cfi_adjust_cfa_offset (-4)
309 popl %edx
310 cfi_adjust_cfa_offset (-4)
311 andb $1, %al
312 jz 27f // jump if not odd
313 cmpl $0xffe00000, %edx
314 jbe 27f // does not fit in mantissa bits
315 // It's an odd integer.
316 // Raise divide-by-zero exception and get minus infinity value.
317 fldl MO(one)
318 fdivl MO(zero)
319 fchs
320 ret
321
322 cfi_adjust_cfa_offset (8)
323 25: fstp %st(0)
324 26: addl $8, %esp
325 cfi_adjust_cfa_offset (-8)
326 27: // Raise divide-by-zero exception and get infinity value.
327 fldl MO(one)
328 fdivl MO(zero)
329 ret
330
331 cfi_adjust_cfa_offset (8)
332 .align ALIGNARG(4)
333 // x is ±0 and y is > 0. We must find out whether y is an odd integer.
334 21: testb $2, %dh
335 jz 22f
336
337 fld %st // y : y
338 fistpll (%esp) // y
339 fildll (%esp) // int(y) : y
340 fucompp // <empty>
341 fnstsw
342 sahf
343 jne 23f
344
345 // OK, the value is an integer, but is the number of bits small
346 // enough so that all are coming from the mantissa?
347 popl %eax
348 cfi_adjust_cfa_offset (-4)
349 popl %edx
350 cfi_adjust_cfa_offset (-4)
351 andb $1, %al
352 jz 24f // jump if not odd
353 cmpl $0xffe00000, %edx
354 jae 24f // does not fit in mantissa bits
355 // It's an odd integer.
356 fldl MO(mzero)
357 ret
358
359 cfi_adjust_cfa_offset (8)
360 22: fstp %st(0)
361 23: addl $8, %esp // Don't use 2 x pop
362 cfi_adjust_cfa_offset (-8)
363 24: fldl MO(zero)
364 ret
365
366 END(__ieee754_pow)
367
368