*
*   $VER: pow.s 33.1 (19.1.97)
*
*   calculates the result of a to the power of b.
*
*   Version history:
*
*   33.1    19.1.97 laukkanen
*
*           - created (using the fact that pow(a,b) == exp(b*log(a))
*
*   33.2    20.1.97 laukkanen
*
*           - added pow2 and pow10
*
*   33.3    25.1.97 laukkanen
*
*           - hopefully fixed all the special cases in pow()
*           - further optimization is possible if the cases
*             where y is integer are handled separately.
*
*

    machine 68040
    fpu     1

    XDEF    _pow
    XDEF    @pow
    XDEF    _pow2
    XDEF    @pow2
    XDEF    _pow10
    XDEF    @pow10

    XREF    @exp
    XREF    @log

        cnop            0,4
_pow
        fmove.d         (4,sp),fp0
        fmove.d         (12,sp),fp1
@pow
        fmove.s         fp0,-(sp)               ; Let's check if
        move.w          (sp),d0                 ; x is negative or zero
        bmi.s           .neg                    ; and branch accordingly
        addq.l          #4,sp
        and.w           #$7f80,d0
        beq.s           .zero
.pos
        fmove.x         fp2,-(sp)               ; x**y is calculated
        fmove.x         fp1,fp2                 ; with exp(y*ln(x))
        jsr             @log
        fmul.x          fp2,fp0
        jsr             @exp
        fmove.x         (sp)+,fp2
        rts
.zero                                           ; x is zero
        fmove.s         fp1,-(sp)
        move.w          (sp),d0                 ; if y is zero
        addq.l          #4,sp                   ;   return 1;
        and.w           #$7f80,d0               ; else return 0;
        beq             .zero1                  ; My calculator returns
        rts                                     ; error for 0**0 but
.zero1                                          ; "cephes" math package
        fmove.s         #1,fp0                  ; did so. OTOH SAS/C
        rts                                     ; scm881.lib just returns
                                                ; undefined number.
.neg
        fabs.x          fp0
        addq.l          #4,sp
        fmove.x         fp2,-(sp)
        moveq           #0,d1
        fmove.x         fp1,fp2
        fabs.x          fp1
        fmove.l         fp1,-(sp)
        move.l          (sp)+,d0
        and.b           #$01,d0
        beq             .even
        moveq           #1,d1
.even
        move.l          d1,-(sp)
        jsr             @log
        fmul.x          fp2,fp0
        jsr             @exp
        move.l          (sp)+,d1
        beq             .even2
        fneg.x          fp0
.even2
        fmove.x         (sp)+,fp2
        rts

        cnop            0,4
LOGTABLE



_pow2
        fmove.d         (4,sp),fp0
@pow2
        fmul.d          (.ln2,pc),fp0
        jsr             @exp
        rts

        cnop            0,4

.ln2    dc.d            0.69314718055995

_pow10
        fmove.d         (4,sp),fp0
@pow10
        fmul.d          (.ln10,pc),fp0
        jsr             @exp
        rts

        cnop            0,4

.ln10   dc.d            2.302585092994