;	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