 ; pow.asm
comment ^

        This sample program shows one method of calculating X raised to
        the Y power where Y is a relatively small positive integer.
        This code requires a numeric coprocessor.

        The algorithm:
        Starting from the most significant bit (bit number i) in our
        exponent, we calculate Z = X raised to the Y power as in the
        following psuedocode:

        Z = X
        while (i > 0)
        {
                Z = Z * Z
                if (bit i in exponent is 1)
                {
                        Z = Z * X
                }
                i = i - 1
        }


        released to the public domain by the author
^
        IDEAL
        MODEL small
        p486

        STACK   400h

        DATASEG
X           dt  23.125          ; an arbitrary X
Y           db  12              ; a positive integral exponent
result      dt  ?               ; the answer, as everyone knows,
                              ; is 23387343638035088 and change

        CODESEG
proc  main
        STARTUPCODE                     ;
        finit                           ; initialize the NPU
        mov     cx,8                    ; get set for eight bits max
        mov     al,[Y]                  ; load our exponent
        or      al,al                   ; Q: is our exponent zero?
        jnz     LoadX                   ;   N: we continue
        fld1                            ;   Y: we know the answer
        jmp     Exit                    ; and exit
LoadX:
        fld     [X]                     ; ST(0) = X
        fld     [X]                     ; ST(0) = X, ST(1) = X
Slider:
        dec     cx                      ; decrement our count
        shl     al,1                    ; shift exponent left until..
        jnc     Slider                  ;  ..we find our first set bit
CalcExp:
        fmul    st(0), st(0)            ; square our result
        shl     al,1                    ; test high bit
        jnc     SkipMul                 ; if it's clear, skip mul
        fmul    st(0), st(1)            ; multiply by original X
SkipMul:
        loop    CalcExp                 ; continue with all remaining bits
Exit:
        fstp    [result]                ; store X**Y in result
        ffree   st(0)                   ; clean up
        EXITCODE                        ;
endp  main
        END
