;*DDK*************************************************************************/
;
; COPYRIGHT (C) Microsoft Corporation, 1989
; 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.;
;*****************************************************************************/
;/*****************************************************************************
;*
;* SOURCE FILE NAME = PRDMATH.ASM
;*
;* DESCRIPTIVE NAME = FIXED Point math routines
;*
;*
;* DATE         08/30/87
;*
;* DESCRIPTION
;*
;* FUNCTIONS    frmul
;*              frdiv
;*              frsqrt
;*
;* NOTES        NONE
;*
;* STRUCTURES   UQUAD, LONG
;*
;* EXTERNAL REFERENCES
;*
;*              NONE
;*
;* EXTERNAL FUNCTIONS
;*
;*              NONE
;*
;* CHANGE ACTIVITY =
;*   DATE      FLAG        APAR   CHANGE DESCRIPTION
;*   --------  ----------  -----  --------------------------------------
;*   mm/dd/yy  @Vr.mpppxx  xxxxx  xxxxxxx
;*   08/30/87                     Created
;*                                Charles Whitmer [chuckwh]
;*
;*   03/22/88                     Wrote frdiv
;*                                Tues Mar 22, 1988           -by- Larry French
;*                                [larryf]
;*
;*   03/23/88                     Converted the code from the engine's
;*                                MATH.ASM so that it takes parameters
;*                                from the stack.  This used to be a
;*                                macro.
;*                                Wed Mar 23, 1988           -by- Larry French
;*                                [larryf]
;*                                frmul
;*
;*   09/13/88                     Wrote frsqrt
;*                                Tue 13-Sep-1988 11:42:00 -by-  Roger Schrag
;*                                [rogersc]
;*
;*   09/13/88                     Replaced explicit code with invocation
;*                                of multiply routine from MATH.ASM
;*                                Tue Sep 13, 1988           -by- Roger Schrag
;*                                [rogersc]
;*                                frmul
;*
;*   08/21/89                     Eliminated extraneous code
;*                                Mon 21 Aug 1989 -by- Mark Owen [c-marko]
;*   05/20/93                     Modified for FLAT MODEL
;*                                Kran
;*****************************************************************************/

        .386p
        .xlist
        include 32cmacro.inc
        include profile.inc
        include  os2def.inc
        .list
?PLM = CC_SYSCALL

INCL_NOCOMMON  EQU 1

UQUAD   struc
uq0     dw      ?
uq1     dw      ?
uq2     dw      ?
uq3     dw      ?
UQUAD   ends

assert  macro
        endm

FIREWALLS     equ     0
MININT        equ     8000h
wp            equ     word ptr

LONG    struc
lo      dw      ?
hi      dw      ?
LONG    ends

multiply_fixed  macro
        local   multiply_fixed_done
        local   no_ovflow

        mov      edx, ebx
        imul     edx
        mov      ebx, edx
        and      ebx, 0FFFF8000h
        jz       no_ovflow       ;no overflow
        cmp      ebx, 0FFFF8000h
        jz       no_ovflow       ;no overflow
        jmp      multiply_fixed_done   ;overflow
no_ovflow:
        mov      dword ptr qtmp.uq0, eax    ; move result to app reg
        mov      dword ptr qtmp.uq2, edx
        mov      eax, dword ptr qtmp.uq1

multiply_fixed_done:
        endm

; PAGE_TUNE
; Modify segment to SEG35 to fit in with new tuned code.
; Also, remove 'SEG35 ENDS' immediately following next line.
SEG35     SEGMENT dword use32 public 'CODE'
          ASSUME  CS:FLAT, DS:FLAT, ES:FLAT, SS:FLAT



;/***************************************************************************
;*
;* FUNCTION NAME = ulNormalize
;*
;* DESCRIPTION   = 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.
;*
;* INPUT         = EAX = ULONG
;*
;* OUTPUT        = EAX = normalized ULONG
;*                 ECX    = shift count
;*                 ZF    = 1 if the ULONG is zero, 0 otherwise
;*
;* RETURN-NORMAL = NONE
;* RETURN-ERROR  = NONE
;*
;**************************************************************************/

cProc   ulNormalize,<PUBLIC,NEAR>
cBegin
        xor       ecx,ecx
        or        eax,eax
        jnz       No_Zero
        jmp       ulNormalize_exit
No_Zero:
        bsr       ecx,eax       ; scan right for the first set bit
        neg       ecx
        add       ecx,31
        shl       eax,cl
ulNormalize_exit:
cEnd


;/***************************************************************************
;*
;* FUNCTION NAME = frmul
;*
;* DESCRIPTION   = Computes the FIXED product of two FIXED values.
;*
;* INPUT         = x  = FIXED value
;*                 y  = FIXED value
;*
;* OUTPUT        = EAX = FIXED x*y
;*
;* RETURN-NORMAL = OF = 0
;* RETURN-ERROR  = OF = 1
;*
;**************************************************************************/

cProc   frmul,<PUBLIC,NODATA,NONWIN>,<esi,edi,ebx>
        parmD   x
        parmD   y
        localV  qtmp,%(size UQUAD)
cBegin
        mov     eax,x               ; eax = x
        mov     ebx,y               ; ebx = y
        multiply_fixed              ; get result in eax
cEnd


;/***************************************************************************
;*
;* FUNCTION NAME = frdiv
;*
;* DESCRIPTION   = Computes the FIXED quotient of two longs.  Works by
;*                 noting that all we want is 16 bits of fraction
;*                 tacked onto the quotient.  To get this, we multiply
;*                 the numerator by 10000h, and let iqdiv do the rest.
;*
;* INPUT         = param1 = FIXED numerator
;*                 param2 = FIXED denominator
;*
;* OUTPUT        = EAX = FIXED quotient
;*
;* RETURN-NORMAL = OF = 0
;* RETURN-ERROR  = OF = 1
;*
;**************************************************************************/

cProc   frdiv,<PUBLIC,NODATA,NONWIN>,<edi,esi,ebx>
        parmD   numerator
        parmD   denominator
cBegin
        mov     eax,numerator
        cdq                             ; edx:eax = numerator * 10000h
        shld    edx,eax,16
        shl     eax, 16
        mov     edi, denominator        ; edi = denominator
        call    iqdiv
cEnd

;/***************************************************************************
;*
;* FUNCTION NAME =  iq_div
;*
;* DESCRIPTION   =
;*                 This routine just keeps track of the signs and calls
;*                 idiv to do the real work.
;*
;*                 Warning:  This routine uses BP as a local variable.
;*
;* INPUT         = EDX:EAX = QUAD Numerator
;*                 EDI     = LONG Denominator
;*
;* OUTPUT        = EAX = quotient
;*                 EDX = remainder
;*
;* RETURN-NORMAL = NONE
;* RETURN-ERROR  = NONE
;*
;**************************************************************************/


IQDIV_RESULT_SIGN       equ     1
IQDIV_REM_SIGN          equ     2


cProc   iqdiv,<PUBLIC>
cBegin
;/*
;** ecx will store bits indicating if the quotient and remainder need to be
;** negated.
;*/

        xor     ecx,ecx

;/*
;** take the absolute value of the denominator
;*/

        or      edi,edi
        jns     denominator_is_cool
        xor     cl,IQDIV_RESULT_SIGN
        neg     edi

;/*
;** EDI is now unsigned so possible overflow is unimportant.
;*/

denominator_is_cool:

;/*
;** take the absolute value of the numerator
;*/

        or      edx,edx
        jns     numerator_is_cool
        xor     cl,IQDIV_RESULT_SIGN + IQDIV_REM_SIGN
        not     edx
        neg     eax
        sbb     edx, -1

;/*
;** EDX:EAX is now unsigned so possible overflow is unimportant.
;*/

numerator_is_cool:

;/*
;** do the unsigned division
;*/
        ; check for overflow

        cmp     edi,edx
        jbe     short iqdiv_overflow

        div     edi

        ;check for overflow
        or      edx,edx
        jns     have_a_bit_to_spare

iqdiv_overflow:

        mov     ah,80h
        dec     ah
        jmp     short iqdiv_exit

have_a_bit_to_spare:

;/*
;** negate the result, if required
;*/

        test    cl,IQDIV_RESULT_SIGN
        jz      result_is_done
        neg     eax
result_is_done:

;/*
;** negate the remainder, if required
;*/

        test    cl,IQDIV_REM_SIGN
        jz      remainder_is_done
        neg     edx

;/*
;** Can't overflow because this number is closer to zero than the
;** denominator.
;*/

        assert  NO,msg_Overflow
remainder_is_done:
iqdiv_exit:

cEnd


;/***************************************************************************
;*
;* FUNCTION NAME = frsqrt
;*
;* DESCRIPTION   = Computes the Fixed square root of a Fixed.  This
;*                 routine simply pulls the parameter off the stack,
;*                 puts it in EAX, and calls the square_root
;*                 function.
;*
;* INPUT         = fxValue = number to square root
;*
;* OUTPUT        = EAX = square root of input
;*
;* RETURN-NORMAL = OF = 0
;* RETURN-ERROR  = OF = 1
;*
;**************************************************************************/
cProc   frsqrt,<PUBLIC,NODATA,NONWIN>,<ESI,EDI,EBX>
        parmD   fxValue
cBegin
        mov     eax,fxValue           ; Load the parameter into EAX
        cCall   square_root             ; Get the square root in EAX
cEnd


;/***************************************************************************
;*
;* FUNCTION NAME = square_root
;*
;* DESCRIPTION   = Computes the Fixed square root of a Fixed using
;*                 Newton's method.
;*
;* INPUT         = EAX = number to square root
;*
;* OUTPUT        = EAX = square root of input
;*                 Zero flag set if the output (and the input) are zero.
;*
;* RETURN-NORMAL = OF = 0
;* RETURN-ERROR  = OF = 1, EAX = 0       if input is negative
;*
;**************************************************************************/


cProc   square_root,<PUBLIC,FAR,NODATA,NONWIN>,<SI,DI>
         localD dwnumber
cBegin
        mov    dwnumber,eax
;;;
;;; Computes the square root of a fixed point integer using a
;;; variation of newtons method.
;;;
;;; History:
;;;     Totally accurate for
;;;     16.16 fixed integers.
        ;;
        ;; verify that dwNumber exists in Z+
        ;;
        or      eax, eax
        jns     @f
        xor     eax, eax                ; if dwNumber < 0, return 0
        jmp     square_root_exit
@@:
        cmp     eax, 0
        jne     @f
        jmp     square_root_exit        ; if dwNumber = 0, return 0
@@:
        cmp     eax, 1
        jne     @f
        mov     eax, 0100h              ; if dwNumber = 1, return 0100h
        jmp     square_root_exit
@@:
        ;;
        ;; at this point the value in EAX exists in Z+
        ;;
        ;; The following variation on newtons algorithm
        ;; will be used:
        ;;
        ;;      F (n+1) := ((square/F(n) + F(n)) / 2
        ;;
        ;; This calculation will be carried out iteratively
        ;; 8 times.  This should ensure a fairly high level
        ;; of accuracy.
        ;;
        shr     eax, 1          ;initial guess = dwNumber / 2
        mov     ebx, eax        ;EBX = F(n)

        mov     ecx, 8                  ;; set loop counter

next_guess:
        xor     edx, edx
        mov     eax, dwNumber   ;convert dwNumber to 32:32 in EDX:EAX
        shld    edx, eax, 16
        shl     eax, 16

        div     ebx             ; EAX = dwNumber(32:32) / Fn(16:16)
        add     eax, ebx        ; EAX = Sq/Fn + Fn
        shr     eax, 1          ; EAX = (Sq/Fn + Fn) / 2
        cmp     ebx, eax        ; are they the same
        je      square_root_exit
        mov     ebx, eax        ; EBX = F(n+1)

        dec     ecx             ; iterate
        jnz     next_guess

square_root_exit:
                                ;*; result in EAX
cEnd

SEG35   ENDS
        END
