3ec7c072409b07a0c1cbfa379e50c47e0a023775
[reactos.git] / reactos / lib / rtl / i386 / math_asm.S
1 /*
2 * COPYRIGHT: See COPYING in the top level directory
3 * PROJECT: ReactOS kernel
4 * PURPOSE: Run-Time Library
5 * FILE: lib/rtl/i386/math.S
6 * PROGRAMER: Alex Ionescu (alex@relsoft.net)
7 * Eric Kohl (ekohl@rz-online.de)
8 *
9 * Copyright (C) 2002 Michael Ringgaard.
10 * All rights reserved.
11 *
12 * Redistribution and use in source and binary forms, with or without
13 * modification, are permitted provided that the following conditions
14 * are met:
15 *
16 * 1. Redistributions of source code must retain the above copyright
17 * notice, this list of conditions and the following disclaimer.
18 * 2. Redistributions in binary form must reproduce the above copyright
19 * notice, this list of conditions and the following disclaimer in the
20 * documentation and/or other materials provided with the distribution.
21 * 3. Neither the name of the project nor the names of its contributors
22 * may be used to endorse or promote products derived from this software
23 * without specific prior written permission.
24
25 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
26 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
28 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
29 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
30 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
31 * OR SERVICES// LOSS OF USE, DATA, OR PROFITS// OR BUSINESS INTERRUPTION)
32 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
33 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
34 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
35 * SUCH DAMAGE.
36 */
37
38 /* GLOBALS ****************************************************************/
39
40 .globl __ftol
41 .globl __aullshr
42 .globl __allrem
43 .globl __aulldiv
44 .globl __allshr
45 .globl __allshl
46 .globl __aullrem
47 .globl __allmul
48 .globl __alldiv
49 .globl __aulldvrm
50 .globl __alldvrm
51 .globl _atan
52 .globl _ceil
53 .globl _cos
54 .globl _fabs
55 .globl _floor
56 .globl _log
57 .globl _pow
58 .globl _sin
59 .globl _sqrt
60 .globl _tan
61 .globl __fltused
62
63 /* DATA ********************************************************************/
64
65 fzero:
66 .long 0 // Floating point zero
67 .long 0 // Floating point zero
68
69 __fltused:
70 .long 0x9875
71
72 /* FUNCTIONS ***************************************************************/
73
74 /*
75 * long long
76 * __alldiv(long long Dividend, long long Divisor)//
77 *
78 * Parameters:
79 * [ESP+04h] - long long Dividend
80 * [ESP+0Ch] - long long Divisor
81 * Registers:
82 * Unknown
83 * Returns:
84 * EDX:EAX - long long quotient (Dividend/Divisor)
85 * Notes:
86 * Routine removes the arguments from the stack.
87 */
88 __alldiv:
89 call ___divdi3
90 ret $0x10
91
92 /*
93 * long long
94 * __allmul(long long Multiplier, long long Multiplicand)//
95 *
96 * Parameters:
97 * [ESP+04h] - long long Multiplier
98 * [ESP+0Ch] - long long Multiplicand
99 * Registers:
100 * Unknown
101 * Returns:
102 * EDX:EAX - long long product (Multiplier*Multiplicand)
103 * Notes:
104 * Routine removes the arguments from the stack.
105 */
106 __allmul:
107 pushl %ebp
108 movl %esp, %ebp
109 pushl %edi
110 pushl %esi
111 pushl %ebx
112 subl $12, %esp
113 movl 16(%ebp), %ebx
114 movl 8(%ebp), %eax
115 mull %ebx
116 movl 20(%ebp), %ecx
117 movl %eax, -24(%ebp)
118 movl 8(%ebp), %eax
119 movl %edx, %esi
120 imull %ecx, %eax
121 addl %eax, %esi
122 movl 12(%ebp), %eax
123 imull %eax, %ebx
124 leal (%ebx,%esi), %eax
125 movl %eax, -20(%ebp)
126 movl -24(%ebp), %eax
127 movl -20(%ebp), %edx
128 addl $12, %esp
129 popl %ebx
130 popl %esi
131 popl %edi
132 popl %ebp
133 ret $0x10
134
135 /*
136 * unsigned long long
137 * __aullrem(unsigned long long Dividend, unsigned long long Divisor)//
138 *
139 * Parameters:
140 * [ESP+04h] - unsigned long long Dividend
141 * [ESP+0Ch] - unsigned long long Divisor
142 * Registers:
143 * Unknown
144 * Returns:
145 * EDX:EAX - unsigned long long remainder (Dividend%Divisor)
146 * Notes:
147 * Routine removes the arguments from the stack.
148 */
149 __aullrem:
150 call ___umoddi3
151 ret $16
152
153 /*
154 * long long
155 * __allshl(long long Value, unsigned char Shift)//
156 *
157 * Parameters:
158 * EDX:EAX - signed long long value to be shifted left
159 * CL - number of bits to shift by
160 * Registers:
161 * Destroys CL
162 * Returns:
163 * EDX:EAX - shifted value
164 */
165 __allshl:
166 shldl %cl, %eax, %edx
167 sall %cl, %eax
168 andl $32, %ecx
169 je 1f
170 movl %eax, %edx
171 xorl %eax, %eax
172 1:
173 ret
174
175 /*
176 * long long
177 * __allshr(long long Value, unsigned char Shift)//
178 *
179 * Parameters:
180 * EDX:EAX - signed long long value to be shifted right
181 * CL - number of bits to shift by
182 * Registers:
183 * Destroys CL
184 * Returns:
185 * EDX:EAX - shifted value
186 */
187 __allshr:
188 shrdl %cl, %edx, %eax
189 sarl %cl, %edx
190 andl $32, %ecx
191 je 1f
192 movl %edx, %eax
193 sarl $31, %edx
194 1:
195 ret
196
197 /*
198 * unsigned long long
199 * __aulldiv(unsigned long long Dividend, unsigned long long Divisor)//
200 *
201 * Parameters:
202 * [ESP+04h] - unsigned long long Dividend
203 * [ESP+0Ch] - unsigned long long Divisor
204 * Registers:
205 * Unknown
206 * Returns:
207 * EDX:EAX - unsigned long long quotient (Dividend/Divisor)
208 * Notes:
209 * Routine removes the arguments from the stack.
210 */
211 __aulldiv:
212 call ___udivdi3
213 ret $16
214
215 /*
216 * unsigned long long
217 * __aullshr(unsigned long long Value, unsigned char Shift)//
218 *
219 * Parameters:
220 * EDX:EAX - unsigned long long value to be shifted right
221 * CL - number of bits to shift by
222 * Registers:
223 * Destroys CL
224 * Returns:
225 * EDX:EAX - shifted value
226 */
227 __aullshr:
228 shrdl %cl, %edx, %eax
229 shrl %cl, %edx
230 andl $32, %ecx
231 je 1f
232 movl %edx, %eax
233 1:
234 ret
235
236 /*
237 * long long
238 * __allrem(long long Dividend, long long Divisor)//
239 *
240 * Parameters:
241 * [ESP+04h] - long long Dividend
242 * [ESP+0Ch] - long long Divisor
243 * Registers:
244 * Unknown
245 * Returns:
246 * EDX:EAX - long long remainder (Dividend/Divisor)
247 * Notes:
248 * Routine removes the arguments from the stack.
249 */
250 __allrem:
251 call ___moddi3
252 ret $16
253
254 .intel_syntax noprefix
255
256 /*
257 * This routine is called by MSVC-generated code to convert from floating point
258 * to integer representation. The floating point number to be converted is
259 * on the top of the floating point stack.
260 */
261 __ftol:
262 /* Set up stack frame */
263 push ebp
264 mov ebp, esp
265
266 /* Set "round towards zero" mode */
267 fstcw [ebp-2]
268 wait
269 mov ax, [ebp-2]
270 or ah, 0xC
271 mov [ebp-4], ax
272 fldcw [ebp-4]
273
274 /* Do the conversion */
275 fistp qword ptr [ebp-12]
276
277 /* Restore rounding mode */
278 fldcw [ebp-2]
279
280 /* Return value */
281 mov eax, [ebp-12]
282 mov edx, [ebp-8]
283
284 /* Remove stack frame and return*/
285 leave
286 ret
287
288 __alldvrm:
289 push edi
290 push esi
291 push ebp
292
293 // Set up the local stack and save the index registers. When this is done
294 // the stack frame will look as follows (assuming that the expression a/b will
295 // generate a call to alldvrm(a, b)):
296 //
297 // -----------------
298 // | |
299 // |---------------|
300 // | |
301 // |--divisor (b)--|
302 // | |
303 // |---------------|
304 // | |
305 // |--dividend (a)-|
306 // | |
307 // |---------------|
308 // | return addr** |
309 // |---------------|
310 // | EDI |
311 // |---------------|
312 // | ESI |
313 // |---------------|
314 // ESP---->| EBP |
315 // -----------------
316 //
317
318 #define DVNDLO [esp + 16] // stack address of dividend (a)
319 #define DVNDHI [esp + 20] // stack address of dividend (a)
320 #define DVSRLO [esp + 24] // stack address of divisor (b)
321 #define DVSRHI [esp + 28] // stack address of divisor (b)
322
323 // Determine sign of the quotient (edi = 0 if result is positive, non-zero
324 // otherwise) and make operands positive.
325 // Sign of the remainder is kept in ebp.
326
327 xor edi,edi // result sign assumed positive
328 xor ebp,ebp // result sign assumed positive
329
330 mov eax,DVNDHI // hi word of a
331 or eax,eax // test to see if signed
332 jge short L1 // skip rest if a is already positive
333 inc edi // complement result sign flag
334 inc ebp // complement result sign flag
335 mov edx,DVNDLO // lo word of a
336 neg eax // make a positive
337 neg edx
338 sbb eax,0
339 mov DVNDHI,eax // save positive value
340 mov DVNDLO,edx
341 L1:
342 mov eax,DVSRHI // hi word of b
343 or eax,eax // test to see if signed
344 jge short L2 // skip rest if b is already positive
345 inc edi // complement the result sign flag
346 mov edx,DVSRLO // lo word of a
347 neg eax // make b positive
348 neg edx
349 sbb eax,0
350 mov DVSRHI,eax // save positive value
351 mov DVSRLO,edx
352 L2:
353
354 //
355 // Now do the divide. First look to see if the divisor is less than 4194304K.
356 // If so, then we can use a simple algorithm with word divides, otherwise
357 // things get a little more complex.
358 //
359 // NOTE - eax currently contains the high order word of DVSR
360 //
361
362 or eax,eax // check to see if divisor < 4194304K
363 jnz short L3 // nope, gotta do this the hard way
364 mov ecx,DVSRLO // load divisor
365 mov eax,DVNDHI // load high word of dividend
366 xor edx,edx
367 div ecx // eax <- high order bits of quotient
368 mov ebx,eax // save high bits of quotient
369 mov eax,DVNDLO // edx:eax <- remainder:lo word of dividend
370 div ecx // eax <- low order bits of quotient
371 mov esi,eax // ebx:esi <- quotient
372 //
373 // Now we need to do a multiply so that we can compute the remainder.
374 //
375 mov eax,ebx // set up high word of quotient
376 mul dword ptr DVSRLO // HIWORD(QUOT) * DVSR
377 mov ecx,eax // save the result in ecx
378 mov eax,esi // set up low word of quotient
379 mul dword ptr DVSRLO // LOWORD(QUOT) * DVSR
380 add edx,ecx // EDX:EAX = QUOT * DVSR
381 jmp short L4 // complete remainder calculation
382
383 //
384 // Here we do it the hard way. Remember, eax contains the high word of DVSR
385 //
386
387 L3:
388 mov ebx,eax // ebx:ecx <- divisor
389 mov ecx,DVSRLO
390 mov edx,DVNDHI // edx:eax <- dividend
391 mov eax,DVNDLO
392 L5:
393 shr ebx,1 // shift divisor right one bit
394 rcr ecx,1
395 shr edx,1 // shift dividend right one bit
396 rcr eax,1
397 or ebx,ebx
398 jnz short L5 // loop until divisor < 4194304K
399 div ecx // now divide, ignore remainder
400 mov esi,eax // save quotient
401
402 //
403 // We may be off by one, so to check, we will multiply the quotient
404 // by the divisor and check the result against the orignal dividend
405 // Note that we must also check for overflow, which can occur if the
406 // dividend is close to 2**64 and the quotient is off by 1.
407 //
408
409 mul dword ptr DVSRHI // QUOT * DVSRHI
410 mov ecx,eax
411 mov eax,DVSRLO
412 mul esi // QUOT * DVSRLO
413 add edx,ecx // EDX:EAX = QUOT * DVSR
414 jc short L6 // carry means Quotient is off by 1
415
416 //
417 // do long compare here between original dividend and the result of the
418 // multiply in edx:eax. If original is larger or equal, we are ok, otherwise
419 // subtract one (1) from the quotient.
420 //
421
422 cmp edx,DVNDHI // compare hi words of result and original
423 ja short L6 // if result > original, do subtract
424 jb short L7 // if result < original, we are ok
425 cmp eax,DVNDLO // hi words are equal, compare lo words
426 jbe short L7 // if less or equal we are ok, else subtract
427 L6:
428 dec esi // subtract 1 from quotient
429 sub eax,DVSRLO // subtract divisor from result
430 sbb edx,DVSRHI
431 L7:
432 xor ebx,ebx // ebx:esi <- quotient
433
434 L4:
435 //
436 // Calculate remainder by subtracting the result from the original dividend.
437 // Since the result is already in a register, we will do the subtract in the
438 // opposite direction and negate the result if necessary.
439 //
440
441 sub eax,DVNDLO // subtract dividend from result
442 sbb edx,DVNDHI
443
444 //
445 // Now check the result sign flag to see if the result is supposed to be positive
446 // or negative. It is currently negated (because we subtracted in the 'wrong'
447 // direction), so if the sign flag is set we are done, otherwise we must negate
448 // the result to make it positive again.
449 //
450
451 dec ebp // check result sign flag
452 jns short L9 // result is ok, set up the quotient
453 neg edx // otherwise, negate the result
454 neg eax
455 sbb edx,0
456
457 //
458 // Now we need to get the quotient into edx:eax and the remainder into ebx:ecx.
459 //
460 L9:
461 mov ecx,edx
462 mov edx,ebx
463 mov ebx,ecx
464 mov ecx,eax
465 mov eax,esi
466
467 //
468 // Just the cleanup left to do. edx:eax contains the quotient. Set the sign
469 // according to the save value, cleanup the stack, and return.
470 //
471
472 dec edi // check to see if result is negative
473 jnz short L8 // if EDI == 0, result should be negative
474 neg edx // otherwise, negate the result
475 neg eax
476 sbb edx,0
477
478 //
479 // Restore the saved registers and return.
480 //
481
482 L8:
483 pop ebp
484 pop esi
485 pop edi
486
487 ret 16
488
489 __aulldvrm:
490
491 // ulldvrm - unsigned long divide and remainder
492 //
493 // Purpose:
494 // Does a unsigned long divide and remainder of the arguments. Arguments
495 // are not changed.
496 //
497 // Entry:
498 // Arguments are passed on the stack:
499 // 1st pushed: divisor (QWORD)
500 // 2nd pushed: dividend (QWORD)
501 //
502 // Exit:
503 // EDX:EAX contains the quotient (dividend/divisor)
504 // EBX:ECX contains the remainder (divided % divisor)
505 // NOTE: this routine removes the parameters from the stack.
506 //
507 // Uses:
508 // ECX
509 //
510 push esi
511
512 // Set up the local stack and save the index registers. When this is done
513 // the stack frame will look as follows (assuming that the expression a/b will
514 // generate a call to aulldvrm(a, b)):
515 //
516 // -----------------
517 // | |
518 // |---------------|
519 // | |
520 // |--divisor (b)--|
521 // | |
522 // |---------------|
523 // | |
524 // |--dividend (a)-|
525 // | |
526 // |---------------|
527 // | return addr** |
528 // |---------------|
529 // ESP---->| ESI |
530 // -----------------
531 //
532
533 #undef DVNDLO
534 #undef DVNDHI
535 #undef DVSRLO
536 #undef DVSRHI
537 #define DVNDLO [esp + 8] // stack address of dividend (a)
538 #define DVNDHI [esp + 8] // stack address of dividend (a)
539 #define DVSRLO [esp + 16] // stack address of divisor (b)
540 #define DVSRHI [esp + 20] // stack address of divisor (b)
541
542 //
543 // Now do the divide. First look to see if the divisor is less than 4194304K.
544 // If so, then we can use a simple algorithm with word divides, otherwise
545 // things get a little more complex.
546 //
547
548 mov eax,DVSRHI // check to see if divisor < 4194304K
549 or eax,eax
550 jnz short .L1 // nope, gotta do this the hard way
551 mov ecx,DVSRLO // load divisor
552 mov eax,DVNDHI // load high word of dividend
553 xor edx,edx
554 div ecx // get high order bits of quotient
555 mov ebx,eax // save high bits of quotient
556 mov eax,DVNDLO // edx:eax <- remainder:lo word of dividend
557 div ecx // get low order bits of quotient
558 mov esi,eax // ebx:esi <- quotient
559
560 //
561 // Now we need to do a multiply so that we can compute the remainder.
562 //
563 mov eax,ebx // set up high word of quotient
564 mul dword ptr DVSRLO // HIWORD(QUOT) * DVSR
565 mov ecx,eax // save the result in ecx
566 mov eax,esi // set up low word of quotient
567 mul dword ptr DVSRLO // LOWORD(QUOT) * DVSR
568 add edx,ecx // EDX:EAX = QUOT * DVSR
569 jmp short .L2 // complete remainder calculation
570
571 //
572 // Here we do it the hard way. Remember, eax contains DVSRHI
573 //
574
575 .L1:
576 mov ecx,eax // ecx:ebx <- divisor
577 mov ebx,DVSRLO
578 mov edx,DVNDHI // edx:eax <- dividend
579 mov eax,DVNDLO
580 .L3:
581 shr ecx,1 // shift divisor right one bit// hi bit <- 0
582 rcr ebx,1
583 shr edx,1 // shift dividend right one bit// hi bit <- 0
584 rcr eax,1
585 or ecx,ecx
586 jnz short .L3 // loop until divisor < 4194304K
587 div ebx // now divide, ignore remainder
588 mov esi,eax // save quotient
589
590 //
591 // We may be off by one, so to check, we will multiply the quotient
592 // by the divisor and check the result against the orignal dividend
593 // Note that we must also check for overflow, which can occur if the
594 // dividend is close to 2**64 and the quotient is off by 1.
595 //
596
597 mul dword ptr DVSRHI // QUOT * DVSRHI
598 mov ecx,eax
599 mov eax,DVSRLO
600 mul esi // QUOT * DVSRLO
601 add edx,ecx // EDX:EAX = QUOT * DVSR
602 jc short .L4 // carry means Quotient is off by 1
603
604 //
605 // do long compare here between original dividend and the result of the
606 // multiply in edx:eax. If original is larger or equal, we are ok, otherwise
607 // subtract one (1) from the quotient.
608 //
609
610 cmp edx,DVNDHI // compare hi words of result and original
611 ja short .L4 // if result > original, do subtract
612 jb short .L5 // if result < original, we are ok
613 cmp eax,DVNDLO // hi words are equal, compare lo words
614 jbe short .L5 // if less or equal we are ok, else subtract
615 .L4:
616 dec esi // subtract 1 from quotient
617 sub eax,DVSRLO // subtract divisor from result
618 sbb edx,DVSRHI
619 .L5:
620 xor ebx,ebx // ebx:esi <- quotient
621
622 .L2:
623 //
624 // Calculate remainder by subtracting the result from the original dividend.
625 // Since the result is already in a register, we will do the subtract in the
626 // opposite direction and negate the result.
627 //
628
629 sub eax,DVNDLO // subtract dividend from result
630 sbb edx,DVNDHI
631 neg edx // otherwise, negate the result
632 neg eax
633 sbb edx,0
634
635 //
636 // Now we need to get the quotient into edx:eax and the remainder into ebx:ecx.
637 //
638 mov ecx,edx
639 mov edx,ebx
640 mov ebx,ecx
641 mov ecx,eax
642 mov eax,esi
643 //
644 // Just the cleanup left to do. edx:eax contains the quotient.
645 // Restore the saved registers and return.
646 //
647
648 pop esi
649
650 ret 16
651
652 _atan:
653 push ebp
654 mov ebp,esp
655 fld qword ptr [ebp+8] // Load real from stack
656 fld1 // Load constant 1
657 fpatan // Take the arctangent
658 pop ebp
659 ret
660
661 _ceil:
662 push ebp
663 mov ebp,esp
664 sub esp,4 // Allocate temporary space
665 fld qword ptr [ebp+8] // Load real from stack
666 fstcw [ebp-2] // Save control word
667 fclex // Clear exceptions
668 mov word ptr [ebp-4],0xb63 // Rounding control word
669 fldcw [ebp-4] // Set new rounding control
670 frndint // Round to integer
671 fclex // Clear exceptions
672 fldcw [ebp-2] // Restore control word
673 mov esp,ebp // Deallocate temporary space
674 pop ebp
675 ret
676
677 _cos:
678 push ebp
679 mov ebp,esp // Point to the stack frame
680 fld qword ptr [ebp+8] // Load real from stack
681 fcos // Take the cosine
682 pop ebp
683 ret
684
685 _fabs:
686 push ebp
687 mov ebp,esp
688 fld qword ptr [ebp+8] // Load real from stack
689 fabs // Take the absolute value
690 pop ebp
691 ret
692
693 _floor:
694 push ebp
695 mov ebp,esp
696 sub esp,4 // Allocate temporary space
697 fld qword ptr [ebp+8] // Load real from stack
698 fstcw [ebp-2] // Save control word
699 fclex // Clear exceptions
700 mov word ptr [ebp-4],0x763 // Rounding control word
701 fldcw [ebp-4] // Set new rounding control
702 frndint // Round to integer
703 fclex // Clear exceptions
704 fldcw [ebp-2] // Restore control word
705 mov esp,ebp
706 pop ebp
707 ret
708
709 _log:
710 push ebp
711 mov ebp,esp
712 fld qword ptr [ebp+8] // Load real from stack
713 fldln2 // Load log base e of 2
714 fxch st(1) // Exchange st, st(1)
715 fyl2x // Compute the natural log(x)
716 pop ebp
717 ret
718
719 _pow:
720 push ebp
721 mov ebp,esp
722 sub esp,12 // Allocate temporary space
723 push edi // Save register edi
724 push eax // Save register eax
725 mov dword ptr [ebp-12],0 // Set negation flag to zero
726 fld qword ptr [ebp+16] // Load real from stack
727 fld qword ptr [ebp+8] // Load real from stack
728 mov edi,offset flat:fzero // Point to real zero
729 fcom qword ptr [edi] // Compare x with zero
730 fstsw ax // Get the FPU status word
731 mov al,ah // Move condition flags to AL
732 lahf // Load Flags into AH
733 and al, 0b01000101 // Isolate C0, C2 and C3
734 and ah,not 0b01000101 // Turn off CF, PF and ZF
735 or ah,al // Set new CF, PF and ZF
736 sahf // Store AH into Flags
737 jb __fpow1 // Re-direct if x < 0
738 ja __fpow3 // Re-direct if x > 0
739 fxch // Swap st, st(1)
740 fcom qword ptr [edi] // Compare y with zero
741 fxch // Restore x as top of stack
742 fstsw ax // Get the FPU status word
743 mov al,ah // Move condition flags to AL
744 lahf // Load Flags into AH
745 and al, 0b01000101 // Isolate C0, C2 and C3
746 and ah,not 0b01000101 // Turn off CF, PF and ZF
747 or ah,al // Set new CF, PF and ZF
748 sahf // Store AH into Flags
749 ja __fpow3 // Re-direct if y > 0
750 fstp st(1) // Set new stack top and pop
751 mov eax,1 // Set domain error (EDOM)
752 jmp __fpow5 // End of case
753 __fpow1: fxch // Put y on top of stack
754 fld st // Duplicate y as st(1)
755 frndint // Round to integer
756 fxch // Put y on top of stack
757 fcomp // y = int(y) ?
758 fstsw ax // Get the FPU status word
759 mov al,ah // Move condition flags to AL
760 lahf // Load Flags into AH
761 and al, 0b01000101 // Isolate C0, C2 and C3
762 and ah,not 0b01000101 // Turn off CF, PF and ZF
763 or ah,al // Set new CF, PF and ZF
764 sahf // Store AH into Flags
765 je __fpow2 // Proceed if y = int(y)
766 fstp st(1) // Set new stack top and pop
767 fldz // Set result to zero
768 fstp st(1) // Set new stack top and pop
769 mov eax,1 // Set domain error (EDOM)
770 jmp __fpow5 // End of case
771 __fpow2: fist dword ptr [ebp-12] // Store y as integer
772 and dword ptr [ebp-12],1 // Set bit if y is odd
773 fxch // Put x on top of stack
774 fabs // x = |x|
775 __fpow3: fldln2 // Load log base e of 2
776 fxch st(1) // Exchange st, st(1)
777 fyl2x // Compute the natural log(x)
778 fmulp // Compute y * ln(x)
779 fldl2e // Load log base 2(e)
780 fmulp st(1),st // Multiply x * log base 2(e)
781 fst st(1) // Push result
782 frndint // Round to integer
783 fsub st(1),st // Subtract
784 fxch // Exchange st, st(1)
785 f2xm1 // Compute 2 to the (x - 1)
786 fld1 // Load real number 1
787 faddp // 2 to the x
788 fscale // Scale by power of 2
789 fstp st(1) // Set new stack top and pop
790 test dword ptr [ebp-12],1 // Negation required ?
791 jz __fpow4 // No, re-direct
792 fchs // Negate the result
793 __fpow4: fstp qword ptr [ebp-8] // Save (double)pow(x, y)
794 fld qword ptr [ebp-8] // Load (double)pow(x, y)
795 fxam // Examine st
796 fstsw ax // Get the FPU status word
797 cmp ah,5 // Infinity ?
798 jne __fpow6 // No, end of case
799 mov eax,2 // Set range error (ERANGE)
800 // Get errno pointer offset
801 __fpow5: int 3
802 mov edi,0 // TODO: offset flat:__crt_errno
803 mov edi,[edi] // Get C errno variable pointer
804 mov dword ptr [edi],eax // Set errno
805 __fpow6: pop eax // Restore register eax
806 pop edi // Restore register edi
807 mov esp,ebp // Deallocate temporary space
808 pop ebp
809 ret
810
811 _sin:
812 push ebp // Save register bp
813 mov ebp,esp // Point to the stack frame
814 fld qword ptr [ebp+8] // Load real from stack
815 fsin // Take the sine
816 pop ebp // Restore register bp
817 ret
818
819 _sqrt:
820 push ebp
821 mov ebp,esp
822 fld qword ptr [ebp+8] // Load real from stack
823 fsqrt // Take the square root
824 pop ebp
825 ret
826
827 _tan:
828 push ebp
829 mov ebp,esp
830 sub esp,4 // Allocate temporary space
831 fld qword ptr [ebp+8] // Load real from stack
832 fptan // Take the tangent
833 fstp dword ptr [ebp-4] // Throw away the constant 1
834 mov esp,ebp // Deallocate temporary space
835 pop ebp
836 ret