/ pow.s (emx+gcc) -- Copyright (c) 1991-1993 by Kolja Elsaesser
/                    Modified 1992-1993 by Eberhard Mattes

#include <libm.h>

        .globl  FUN(pow)

        .text

        .align  2, 0x90

/ double pow (double x, double y)

#define cw1           -8(%ebp)
#define cw2           -6(%ebp)
#define int_exponent  -4(%ebp)
/define saved_ebp      0(%ebp)
/define ret_addr       4(%ebp)
#define x              8(%ebp)
#if defined (LONG_DOUBLE)
#define y             20(%ebp)
#else
#define y             16(%ebp)
#endif

DEF(pow)
        pushl   %ebp
        movl    %esp, %ebp
        subl    $8, %esp
        pushf                           / Save flags
        fstcww  cw1
        movw    cw1, %ax
        andw    $0xf3ff, %ax            / round nearest zero
        movw    %ax, cw2
        FLD     y
        FLD     x
        fcoml   __const_ZERO
        fnstsww %ax
        sahf
        jz      zero
        jc      negative

/ x > 0
positive:
        fyl2x                           / y * ln(x)
        fld     %st
        fldcww  cw2
        frndint
        fldcww  cw1
        fsubr   %st, %st(1)
        fxch    %st(1)
        f2xm1
        faddl   __const_ONE
        fscale
        fxch    %st(1)
        ffree   %st
        fincstp
        jmp     return

/ x = 0
        .align  2, 0x90
zero:
        fld     %st(1)
        fcompl  __const_ZERO
        fnstsww %ax
        sahf
        jz      positive        / 0^0           -> NAN
        jc      positive        / 0^negative    -> NAN
        ffree   %st
        fincstp
        ffree   %st
        fincstp
        fldz
        jmp     return

/ x < 0
        .align  2, 0x90
negative:
        fld     %st(1)
        fldcww  cw2
        frndint
/ Copy integral part of exponent to memory
        fld     %st
        fistpl  int_exponent
        fldcww  cw1
        fsub    %st(2), %st
        fcompl  __const_ZERO
        fstsww  %ax
        sahf
/ Non-integral exponent -> error
        jnz     positive
/ Integral exponent -> special case
/ Compute x^y where x is positive
        fchs
        fyl2x                           / y * ln(x)
        fld     %st
        fldcww  cw2
        frndint
        fldcww  cw1
        fsubr   %st, %st(1)
        fxch    %st(1)
        f2xm1
        faddl   __const_ONE
        fscale
        fxch    %st(1)
        ffree   %st
        fincstp
/ Correction for odd exponent and negative base
        testl   $1, int_exponent
        jz      even_exp
        fchs
        .align  4, 0x90
even_exp:
return: CONV(x)                         / convert to double
        popf
        leave
        _xam
        j_inf   1f
        j_nan   2f
        ret

        .align  2, 0x90
1:      SETERRNO($ERANGE)
        ret

        .align  2, 0x90
2:      SETERRNO($EDOM)
        ret
