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