Power-Programmierung/CD2.mdf/c/library/dos/math/cephes/qfloat/divn.asm

Assembly Source File | 1992-11-17 | 5.1 KB | 429 lines

;    Static Name Aliases
;
    TITLE   divm

include qhead.asm
extrn    _mulm:proc
extrn    mdnorm:proc

_DATA    SEGMENT
;extrn    SC:word
B    DW    0
count    DW    0
quot    DW    NQ+2 DUP (?)
_square    DW    NQ+2 DUP (?)
_ans    DW    NQ+2 DUP (?)
_DATA    ENDS

_TEXT      SEGMENT


; Macro to accumulate twice the partial product
z2mul    MACRO    a,b
ao = 6 + 2*a
bo = 6 + 2*b
co = 8 + 2*a + 2*b
do = 6 + 2*a + 2*b
eo = 4 + 2*a + 2*b
    mov    ax, word ptr [di]+ao
    mul    word ptr [si]+bo
    xor    cx,cx
    sal    ax,1
    rcl    dx,1
    rcl    cx,1
    add    word ptr [bx]+co,ax
    adc    word ptr [bx]+do,dx
    adc    word ptr [bx]+eo,cx
    ENDM




extrn    _shdn1:PROC
extrn    _shup1:PROC
extrn    _normlz:PROC
;
; Floating point divide
; @stack into @stack+1
;
;
; The program makes use of 80286 16-bit unsigned
; MUL and DIV instructions.
; -Copyright 1988 by Stephen L. Moshier


; zdiv( source, dest )
; dest = numerator
    PUBLIC _divm
_divm    PROC    NEAR
    push    bp
    mov    bp,sp
    push    si
    push    di
    push    bx
    push    cx
    push    es

    push    ds
    pop    es



    mov    si, word ptr [bp]+4    ;source, denominator
    mov    di,si
    add    si,8

; Test if the low order denominator is zero.
    lodsw
    or    ax,ax
    je    tiszer
    jmp    nzero

tiszer:
    mov    cx,OMG-2
$iszer:
    lodsw
    or    ax,ax
    jne    nnzero
    loop    $iszer

    jmp    short yzer


nnzero:
    jmp    nzero


yzer:
; Denominator only has one word of significant bits.
; Do a single precision divide.

    xor    bx,bx
    mov    dx, word ptr [di]+6    ; denominator value
; get numerator and shift down 1
    mov    si, word ptr [bp]+6    ;numerator
    add    si,4
    lea    di, quot+4
    mov    cx, OMG
    clc
$gnum:
    lodsw
    rcr    ax,1
    stosw
    loop    $gnum

    mov    ax,bx
    rcr    ax,1
    stosw

    xor    ax,ax
    stosw

    lea    si, quot+6        ; numerator
    lea    di,_ans+6
    mov    cx,dx            ;denominator

    lodsw
    mov    dx,ax
    lodsw
    div    cx
    stosw

REPT OMG-1
    lodsw
    div    cx
    stosw
ENDM

    mov    di, word ptr [bp]+6    ; answer
    sub    word ptr [di]+2,1
; This is the exact quotient since the low order denominator is zero.
    jmp    divfin


nzero:


    mov    si,di    ;get denominator pointer

; clear temp area
    lea    di,quot
    mov    cx,NQ+2
    xor    ax,ax
    rep    stosw

    mov    cx, word ptr [si]+6    ; denominator value

; divide numerator = 1 by most significant word of denominator
; compute quotient 1/B to 2 word precision
    lea    di,quot+6    ; significand of quotient
    mov    dx,04000H    ; 1.0

REPT 2
    xor    ax,ax
    div    cx
    stosw
ENDM

; Calculate double precision quotient using Taylor series

; clear out temporary storage
    xor    ax,ax
    lea    di,_square
    mov    cx,2*NQ+4
    rep    stosw

; square the quotient
    lea    di,quot
    mov    si,di
    lea    bx,_square
    z2mul    0,1
    zmul    1,1
    zmul    0,0

; multiply squared quotient by low order denominator
    mov    si, word ptr [bp]+4
    lea    di,_square
    lea    bx,_ans
;        di si
    zmul    2,1
    zmul    1,1
    zmul    0,1
; shift up 2
    sal    word ptr [bx]+10,1
    rcl    word ptr [bx]+8,1
    rcl    word ptr [bx]+6,1
    sal    word ptr [bx]+10,1
    rcl    word ptr [bx]+8,1
    rcl    word ptr [bx]+6,1

; subtract from quotient
    lea    di,quot
    lea    si,_ans
    mov    ax, word ptr [si]+8
    sub    word ptr [di]+8,ax
    mov    ax, word ptr [si]+6
    sbb    word ptr [di]+6,ax

; do first Newton-Raphson iteration to four word precision
prec = 4
;    mov    count,3
;newtlup:

; Loop is done by in-line code with increasing
; arithmetic precision each iteration.

; Start of Newton-Raphson iterations.
rept NQDIV

; limit precision to maximum available
if prec gt OMG-1
prec = OMG-1
endif

; clear out accumulators
    lea    di,_square
    mov    cx,2*NQ+4
    xor    ax,ax
    rep    stosw

; Form the square of the approximate quotient

    lea    si,quot
    mov    di,si
    lea    bx,_square


hh = prec

rept prec
g = 0
h = hh
rept (hh+1)/2
if g LE OMG-1
if h LE OMG-1
    z2mul    g,h
endif
endif
g = g+1
h = h-1
endm

if g EQ h
    xor    cx,cx
    zmul    g,g
endif

hh = hh-1

endm
    xor    cx,cx
    zmul    0,0


; multiply the square of the quotient by the denominator

    mov    si, word ptr [bp]+4    ;source, denominator
    lea    di,_square
    lea    bx,_ans

hh = prec

rept prec
g = 0
h = hh
rept (hh+1)/2
if g LE OMG-1
if h LE OMG-2
    zmul    g,h
endif
endif
if h LE OMG-1
if g LE OMG-2
    zmul    h,g
endif
endif
g = g+1
h = h-1
endm

if g EQ h
    zmul    g,g
endif

hh = hh-1

endm
    zmul    0,0


; shift the product up 1 bit

    lea    bx,_ans

h = 2*prec+8
    sal    word ptr [bx]+h,1

rept prec+2
h = h-2
    rcl    word ptr [bx]+h,1
endm

; subtract quot - ans

i = 2*prec+8
    std
    lea    si,_ans+i
    lea    di,quot

    lodsw
    sub    word ptr [di]+i,ax

rept prec+2
i=i-2
    lodsw
    sbb    word ptr [di]+i,ax
endm
    cld

; shift result up 1 bit

    lea    bx,quot

h = 2*prec+8
    sal    word ptr [bx]+h,1

rept prec+2
h = h-2
    rcl    word ptr [bx]+h,1
endm

; increase the arithmetic precision for next loop
prec = prec+prec

; end of Newton-Raphson iteration loop
endm

;    sub    count,1
;    je    divdon
;    jmp    newtlup

divdon:
    cld

; multiply 1/x by the numerator

; temp area for the product
    lea    bx, _ans
; clear the temp area
    xor    ax,ax
    mov    di,bx
    mov    cx,NQ+2
    rep    stosw

    lea    si,quot
;    mov    si,WORD PTR [bp+4]    ;source
    mov    di,WORD PTR [bp+6]    ;dest

    xor    cx,cx    ; for adc 0

hh = OMG-1

rept OMG-1

g = 0
h = hh
rept (hh+1)/2
if g LE OMG-2
if h LE OMG-1
    zmul    g,h
endif
endif
if h LE OMG-2
if g LE OMG-1
    zmul    h,g
endif
endif
g = g+1
h = h-1
endm

if g EQ h
    zmul    g,g
endif

hh = hh-1

endm

    zmul    0,0



divfin:


; normalize
    mov    di, word ptr [bp]+6
    lea    bx, _ans
    call    mdnorm

; move out the answer
    lea    si, _ans+4
    add    di,4
    mov    cx,OMG
    rep    movsw

    pop    es
    pop    cx
    pop    bx
    pop    di
    pop    si
    pop    bp
    ret
_divm    ENDP


_TEXT    ENDS
END