home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Stars of Shareware: Programmierung
/
SOURCE.mdf
/
programm
/
msdos
/
c
/
cephes22
/
qfloat
/
divn.asm
< prev
next >
Wrap
Assembly Source File
|
1992-11-17
|
5KB
|
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