Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : DIVN.ASM

 
Output of file : DIVN.ASM contained in archive : CEPHES22.ZIP
; 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


  3 Responses to “Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : DIVN.ASM

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/