;*DDK*************************************************************************/
;
; COPYRIGHT    Copyright (C) 1995 IBM Corporation
;
;    The following IBM OS/2 WARP source code is provided to you solely for
;    the purpose of assisting you in your development of OS/2 WARP device
;    drivers. You may use this code in accordance with the IBM License
;    Agreement provided in the IBM Device Driver Source Kit for OS/2. This
;    Copyright statement may not be removed.;
;*****************************************************************************/
        page    ,132
;/*****************************************************************************
;*
;* SOURCE FILE NAME = MATH.ASM
;*
;*
;* VERSION      V2.0
;*
;* DATE         11/20/90
;*
;* DESCRIPTION  Contains math routines for 64- by 32-bit division.
;*
;* NOTES        NONE
;*
;* STRUCTURES   NONE
;*
;* EXTERNAL REFERENCES
;*
;*              NONE
;*
;* EXTERNAL FUNCTIONS
;*
;*              NONE
;*
;* CHANGE ACTIVITY =
;*   DATE      FLAG        APAR   CHANGE DESCRIPTION
;*   --------  ----------  -----  --------------------------------------
;*   mm/dd/yy  @Vr.mpppxx  xxxxx  xxxxxxx
;*
;*****************************************************************************/

        .386

        .XLIST

        include pmgre.inc
        include driver.inc
        include assert.mac

        .MODEL FLAT

        ASSUME  CS:FLAT,SS:FLAT,DS:FLAT,ES:FLAT

        .CODE

page

;---------------------------Public-Routine------------------------------;
; get_vector_length
;
; finds the length of a vector in an array of vectors
;
; Entry:
;   esi --> point
;
; Returns:
;   eax --> the length of the point
; Error:
;   OF set
; History:
;  Mon 28-Mar-1988 17:17:13  -by-  Paul Klingler [paulk]
;-----------------------------------------------------------------------;

get_vector_length PROC SYSCALL uses ECX

        mov     eax,[esi].POINTL.ptl_x
        imul    eax                 ; edx:eax = x^2
        mov     ebx,eax
        mov     ecx,edx
        mov     eax,[esi].POINTL.ptl_y
        imul    eax                 ; edx:eax = y^2
        add     eax,ebx
        adc     edx,ecx
        call    big_fix_square_root     ; EAX = sqrt(EDX:EAX)
        ret
get_vector_length ENDP

page

;---------------------------Public-Routine------------------------------;
; ulNormalize
;
; Normalizes a ULONG so that the highest order bit is 1.  Returns the
; number of shifts done.  Also returns ZF=1 if the ULONG was zero.
;
; Entry:
;       DX:AX = ULONG
; Returns:
;       DX:AX = normalized ULONG
;       CX    = shift count
;       ZF  = 1 if the ULONG is zero, 0 otherwise
; Registers Destroyed:
;       none
; History:
;  Mon 25-Jan-1988 22:07:03  -by-  Charles Whitmer [chuckwh]
; Wrote it.
;-----------------------------------------------------------------------;

ulNormalize PROC SYSCALL

; shift by words

        xor     cx,cx
        or      dx,dx
        js      ulNormalize_exit
        jnz     top_word_ok
        xchg    ax,dx
        or      dx,dx
        jz      ulNormalize_exit        ; the zero exit
        mov     cl,16
        js      ulNormalize_exit
top_word_ok:

; shift by bytes

        or      dh,dh
        jnz     top_byte_ok
        mov     dh,dl
        mov     dl,ah
        mov     ah,al
        xor     al,al
        add     cl,8
        or      dh,dh
        js      ulNormalize_exit
top_byte_ok:

; do the rest by bits

        inc     cx
        add     ax,ax
        adc     dx,dx
        js      ulNormalize_exit
        inc     cx
        add     ax,ax
        adc     dx,dx
        js      ulNormalize_exit
        inc     cx
        add     ax,ax
        adc     dx,dx
        js      ulNormalize_exit
        inc     cx
        add     ax,ax
        adc     dx,dx
        js      ulNormalize_exit
        inc     cx
        add     ax,ax
        adc     dx,dx
        js      ulNormalize_exit
        inc     cx
        add     ax,ax
        adc     dx,dx
        js      ulNormalize_exit
        inc     cx
        add     ax,ax
        adc     dx,dx
ulNormalize_exit:
        ret
ulNormalize ENDP

page

;---------------------------Public-Routine------------------------------;
; uqNormalize
;
; Normalizes a UQUAD (unsigned quadword) so that the highest order bit
; is 1.  Returns only the top 32 bits of the normalized UQUAD.  Returns
; the number of shifts done.  Retuns ZF=1 if the UQUAD was zero.
;
; Entry:
;       EDX:EAX = UQUAD
; Returns:
;       EAX   = normalized top 32 bits
;       ECX   = shift count
;       ZF    = 1 if the UQUAD is zero, 0 otherwise
; Registers Destroyed:
;       EBX
; History:
;  Mon 25-Jan-1988 22:07:03  -by-  Charles Whitmer [chuckwh]
; Wrote it.
;-----------------------------------------------------------------------;

uqNormalize PROC SYSCALL uses ESI EDI

; tmp until move without any trouble

        mov     ebx,eax
        shr     ebx,16
        mov     cx,dx
        shr     edx,16

; use SI to count the shifts

        xor     esi,esi

; get the data into the top three words to free CX

        xchg    ax,cx                   ; CX = lowest order
        or      dx,dx
        js      uqNormalize_exit
        jnz     top_word_good
        mov     dx,ax
        mov     ax,bx
        mov     bx,cx
        mov     si,16
        or      dx,dx
        js      uqNormalize_exit
        jnz     top_word_good

; in this case, all the data fits in DX:AX, let ulNormalize do it

        mov     dx,ax
        mov     ax,bx
        INVOKE  ulNormalize
        jz      uqNormalize_final_exit
        add     cx,32                   ; how much we'd shifted ourselves
        jmp     uqNormalize_final_exit

; all data is in DX:AX:BX, let ulNormalize help

top_word_good:
        INVOKE  ulNormalize
        add     si,cx

; DX:AX is normalized, now we should move the top bits of BX into AX

        mov     di,-1
        rol     bx,cl           ; abcdexxxxxxxxxxx => xxxxxxxxxxxabcde
        shl     di,cl           ; 1111111111111111 => 1111111111100000
        not     di              ; 1111111111100000 => 0000000000011111
        and     bx,di           ; xxxxxxxxxxxabcde => 00000000000abcde
        or      ax,bx
        or      dx,dx           ; ZF = 0
uqNormalize_exit:
        mov     ecx,esi
uqNormalize_final_exit:
        and     eax,0ffffh      ; tmp until move
        shl     edx,16
        or      eax,edx
        ret
uqNormalize ENDP

page

;---------------------------Public-Routine------------------------------;
; square_root
;
; Computes the Fixed square root of a Fixed using Newton's method.
;
; This routine will occasionally give a different result from
; the old_square_root routine.  This is caused by a very small lack of
; precission.  For instance, the square root of 1.0001h is approximately
; 1.00007FFh.  This routine calculates it to be 1.0000800.  Hence when
; we round it off and put it in a fixed point we get 1.0001 instead
; of 1.0000.
;
; When firewalls are set this routine is called by another routine
; which borrows this routine's name.
;
; Entry:
;       DX:AX = number to square root
; Returns:
;       DX:AX = square root of input
;       Zero flag set if the output (and the input) are zero.
; Error Returns:
;       OF = 1, DX:AX = 0       if input is negative
; Registers Destroyed:
;       BX,CX
; Registers Preserved:
;       DS,ES,SI,DI
; History:
;  Tue 02-Aug-1988 12:50:30 -by-  Lee A. Newberg [leen]
; Wrote it.
;-----------------------------------------------------------------------;

square_root PROC SYSCALL uses ESI EDI

        INVOKE  ulNormalize
        jnz     input_non_zero          ; Input wasn't zero
        jmp     lee_root_exit
negative_input:
        xor     ax,ax
        shr     dx,1                    ; Sets Overflow
        cwd
;        assert  O
        jmp     lee_root_exit

input_non_zero:
        jcxz    negative_input          ; Input was negative
        shr     cx,1
        jc      even_shifts
        dec     cx
        shr     dx,1
        rcr     ax,1
even_shifts:
        shr     dx,1
        rcr     ax,1

; Let SI:DI = the normalized square

        push    cx              ; Save the number of shifts
        mov     di,ax
        mov     si,dx           ; SI:DI = the input

; First find the square root of the *integer* DX:AX.
; We will use Newton's method.

; Zeroth guess is 8000

        mov     bx,dx
        add     bh,40h

; SI:DI = the normalized square
; BX = first guess = ((square / zeroth_guess) + zeroth_guess) / 2

; second_guess = ((square / first_guess) + first_guess) / 2

        div     bx
        add     bx,ax
        rcr     bx,1

; SI:DI = the normalized square
; BX = second guess
; third_guess = ((square / second_guess) + second_guess) / 2

        mov     ax,di
        mov     dx,si
        div     bx
        add     bx,ax
        rcr     bx,1

; SI:DI = the normalized square
; BX = third guess.  Now we'll go for some more precision.
; That is, we'll remember that we're working with fixed point
; numbers rather than integers.

        mov     ax,di
        mov     dx,si
        div     bx
        mov     cx,ax           ; CX = high part of quotient
        xor     ax,ax
        div     bx              ; AX = low part of quotient
        mov     dx,cx           ; DX = high part of quotient

; Now average BX:0000 and DX:AX

        add     dx,bx
        rcr     dx,1
        rcr     ax,1

unnormalize:
        pop     cx              ; Restore the number of shifts.

; 7 = half the number of bits in the fraction of a fixed - 1

        add     cl,7
        test    cl,10h
        jz      @F
        mov     ax,dx
        xor     dx,dx
@@:
        test    cl,8
        jz      @F
        mov     al,ah
        mov     ah,dl
        mov     dl,dh
        xor     dh,dh
@@:
        and     cl,7
        jz      @F
        mov     bx,-1
        shr     bx,cl
        not     bx
        shr     ax,cl
        ror     dx,cl
        and     bx,dx
        xor     dx,bx
        or      ax,bx
@@:

; DX:AX = double our answer.  Half it with round off

        rcr     dx,1
        rcr     ax,1
        adc     ax,0
        adc     dx,0
;       assert  NO,msg_Overflow
lee_root_exit:
        ret
square_root ENDP

page

;-------------------------------Public-Routine---------------------------;
; big_fix_square_root
;
;   Approximates the square root of an unsigned BIGFIXED
;
; Entry:
;     EDX:EAX = unsigned BIGFIXED
; Returns:
;     EAX = square root
; Error Returns:
;     none
; Registers Destroyed:
;     EBX,ECX
; Calls:
;     none
; History:
;  Fri 29-Jul-1988 09:18:55 -by-  Lee A. Newberg [leen]
; Rounds answer to the nearest 0.0001h.  Fractions of 0.00008h or greater
; are rounded up.
;
;  Sun 06-Mar-1988            -by-   Martin Picha        [martinpi]
; Wrote it.  (pulled out of solve_quadratic)
;------------------------------------------------------------------------;

big_fix_square_root PROC SYSCALL

; normalize the BIGFIXED for a sqrt
        INVOKE  uqNormalize

; tmp until move      - top -

        mov     edx,eax
        shr     edx,16

; tmp until move      - bot -

        shr     cx,1
        jc      it_will_be_even
        shr     dx,1
        rcr     ax,1
        dec     cx
it_will_be_even:
        shr     dx,1
        rcr     ax,1

; call the square root with almost 32 bits of precision

        push    cx
        invoke  square_root
        pop     cx
;        assert  NO

; remove the normalization
;                                                                16
        mov     dh,dl                           ; we lied about 2
        mov     dl,ah
        mov     ah,al
        xor     al,al

        jcxz    have_sqrt
        dec     cl                              ; leave answer doubled
        cmp     cl,16
        jb      remove_final_16
        mov     ax,dx                           ; shift right by 16
        xor     dx,dx
        sub     cl,16
        jz      have_double
remove_final_16:
        mov     bx,-1
        shr     bx,cl
        not     bx
        shr     ax,cl
        ror     dx,cl
        and     bx,dx
        or      ax,bx
        xor     dx,bx
have_double:

; Half this number but round it off.

        add     ax,1
        adc     dx,0
        shr     dx,1
        rcr     ax,1
        or      ax,ax                   ; clears overflow
have_sqrt:
;       assert  NO,msg_Overflow

; tmp until move

        and     eax,0ffffh      ; tmp until move
        shl     edx,16
        or      eax,edx
        ret
big_fix_square_root ENDP

        end
