home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / CLIPPER / MISC / EMXLIB8F.ZIP / EMX / LIB / MATH / POW.S < prev    next >
Encoding:
Text File  |  1993-01-02  |  2.8 KB  |  125 lines
  1. / pow.s (emx+gcc) -- Copyright (c) 1991-1993 by Kolja Elsaesser
  2. /                    Modified 1992-1993 by Eberhard Mattes
  3.  
  4. #include <libm.h>
  5.  
  6.         .globl  _pow
  7.  
  8.         .text
  9.  
  10.         .align  2, 0x90
  11.  
  12. / double pow (double x, double y)
  13.  
  14. #define cw1           -8(%ebp)
  15. #define cw2           -6(%ebp)
  16. #define int_exponent  -4(%ebp)
  17. /define saved_ebp      0(%ebp)
  18. /define ret_addr       4(%ebp)
  19. #define x              8(%ebp)
  20. #define y             16(%ebp)
  21.  
  22. _pow:   pushl   %ebp
  23.         movl    %esp, %ebp
  24.         subl    $8, %esp
  25.         pushf                           / Save flags
  26.         fstcww  cw1
  27.         movw    cw1, %ax
  28.         andw    $0xf3ff, %ax            / round nearest zero
  29.         movw    %ax, cw2
  30.         fldl    y
  31.         fldl    x
  32.         fcoml   __const_ZERO
  33.         fnstsww %ax
  34.         sahf
  35.         jz      zero
  36.         jc      negative
  37.  
  38. / x > 0
  39. positive:
  40.         fyl2x                           / y * ln(x)
  41.         fld     %st
  42.         fldcww  cw2
  43.         frndint
  44.         fldcww  cw1
  45.         fsubr   %st, %st(1)
  46.         fxch    %st(1)
  47.         f2xm1
  48.         faddl   __const_ONE
  49.         fscale
  50.         fxch    %st(1)
  51.         ffree   %st
  52.         fincstp
  53.         jmp     return
  54.  
  55. / x = 0
  56.         .align  2, 0x90
  57. zero:
  58.         fldl    %st(1)
  59.         fcompl  __const_ZERO
  60.         fnstsww %ax
  61.         sahf
  62.         jz      positive        / 0^0           -> NAN
  63.         jc      positive        / 0^negative    -> NAN
  64.         ffree   %st
  65.         fincstp
  66.         ffree   %st
  67.         fincstp
  68.         fldz
  69.         jmp     return
  70.  
  71. / x < 0
  72.         .align  2, 0x90
  73. negative:
  74.         fld     %st(1)
  75.         fldcww  cw2
  76.         frndint
  77. / Copy integral part of exponent to memory
  78.         fld     %st
  79.         fistpll int_exponent
  80.         fldcww  cw1
  81.         fsub    %st(2), %st
  82.         fcompl  __const_ZERO
  83.         fstsww  %ax
  84.         sahf
  85. / Non-integral exponent -> error
  86.         jnz     positive
  87. / Integral exponent -> special case
  88. / Compute x^y where x is positive
  89.         fchs
  90.         fyl2x                           / y * ln(x)
  91.         fld     %st
  92.         fldcww  cw2
  93.         frndint
  94.         fldcww  cw1
  95.         fsubr   %st, %st(1)
  96.         fxch    %st(1)
  97.         f2xm1
  98.         faddl   __const_ONE
  99.         fscale
  100.         fxch    %st(1)
  101.         ffree   %st
  102.         fincstp
  103. / Correction for odd exponent and negative base
  104.         testl   $1, int_exponent
  105.         jz      even_exp
  106.         fchs
  107.         .align  4, 0x90
  108. even_exp:
  109. return: fstpl   x                       / convert to double
  110.         fldl    x
  111.         popf
  112.         leave
  113.         _xam
  114.         j_inf   1f
  115.         j_nan   2f
  116.         ret
  117.  
  118.         .align  2, 0x90
  119. 1:      SETERRNO($ERANGE)
  120.         ret
  121.  
  122.         .align  2, 0x90
  123. 2:      SETERRNO($EDOM)
  124.         ret
  125.