subttl emfdiv.asm - Division
page
;*******************************************************************************
; Copyright (c) Microsoft Corporation 1991
; All Rights Reserved
;
;emfdiv.asm - long double divide
; by Tim Paterson
;
;Purpose:
; Long double division.
;Inputs:
; ebx:esi = op1 mantissa
; ecx = op1 sign in bit 15, exponent in high half
; edi = pointer to op2 and result location
; [Result] = edi
;
; Exponents are unbiased. Denormals have been normalized using
; this expanded exponent range. Neither operand is allowed to be zero.
;Outputs:
; Jumps to [RoundMode] to round and store result.
;
;Revision History:
;
; [] 09/05/91 TP Initial 32-bit version.
;
;*******************************************************************************
;Dispatch tables for division
;
;One operand has been loaded into ecx:ebx:esi ("source"), the other is
;pointed to by edi ("dest"). edi points to dividend for fdiv,
;to divisor for fdivr.
;
;Tag of source is shifted. Tag values are as follows:
;
.erre TAG_SNGL eq 0 ;SINGLE: low 32 bits are zero
.erre TAG_VALID eq 1
.erre TAG_ZERO eq 2
.erre TAG_SPCL eq 3 ;NAN, Infinity, Denormal, Empty
;dest = dest / source
tFdivDisp label dword ;Source (reg) Dest (*[di])
dd DivSingle ;single single
dd DivSingle ;single double
dd XorDestSign ;single zero
dd DivSpclDest ;single special
dd DivDouble ;double single
dd DivDouble ;double double
dd XorDestSign ;double zero
dd DivSpclDest ;double special
dd DivideByZero ;zero single
dd DivideByZero ;zero double
dd ReturnIndefinite ;zero zero
dd DivSpclDest ;zero special
dd DivSpclSource ;special single
dd DivSpclSource ;special double
dd DivSpclSource ;special zero
dd TwoOpBothSpcl ;special special
dd ReturnIndefinite ;Two infinities
;dest = source / dest
tFdivrDisp label dword ;Source (reg) Dest (*[di])
dd DivrSingle ;single single
dd DivrDouble ;single double
dd DivideByZero ;single zero
dd DivrSpclDest ;single special
dd DivrSingle ;double single
dd DivrDouble ;double double
dd DivideByZero ;double zero
dd DivrSpclDest ;double special
dd XorSourceSign ;zero single
dd XorSourceSign ;zero double
dd ReturnIndefinite ;zero zero
dd DivrSpclDest ;zero special
dd DivrSpclSource ;special single
dd DivrSpclSource ;special double
dd DivrSpclSource ;special zero
dd TwoOpBothSpcl ;special special
dd ReturnIndefinite ;Two infinities
EM_ENTRY eFIDIV16
eFIDIV16:
push offset DivSetResult
jmp Load16Int ;Returns to DivSetResult
EM_ENTRY eFIDIVR16
eFIDIVR16:
push offset DivrSetResult
jmp Load16Int
EM_ENTRY eFIDIV32
eFIDIV32:
push offset DivSetResult
jmp Load32Int
EM_ENTRY eFIDIVR32
eFIDIVR32:
push offset DivrSetResult
jmp Load32Int
EM_ENTRY eFDIV32
eFDIV32:
push offset DivSetResult
jmp Load32Real ;Returns to DivSetResult
EM_ENTRY eFDIVR32
eFDIVR32:
push offset DivrSetResult ;Returns to DivrSetResult
jmp Load32Real
EM_ENTRY eFDIV64
eFDIV64:
push offset DivSetResult
jmp Load64Real ;Returns to DivSetResult
EM_ENTRY eFDIVR64
eFDIVR64:
push offset DivrSetResult
jmp Load64Real ;Returns to DivrSetResult
EM_ENTRY eFDIVRPreg
eFDIVRPreg:
push offset PopWhenDone
EM_ENTRY eFDIVRreg
eFDIVRreg:
xchg esi,edi
EM_ENTRY eFDIVRtop
eFDIVRtop:
mov ecx,EMSEG:[esi].ExpSgn
mov ebx,EMSEG:[esi].lManHi
mov esi,EMSEG:[esi].lManLo
DivrSetResult:
;cl has tag of dividend
mov ebp,offset tFdivrDisp
mov EMSEG:[Result],edi ;Save result pointer
mov ah,cl
mov al,EMSEG:[edi].bTag
and ah,not 1 ;Ignore single vs. double on dividend
cmp ax,1
.erre bTAG_VALID eq 1
.erre bTAG_SNGL eq 0
jz DivrDouble ;Divisor was double
ja TwoOpResultSet
;.erre DivrSingle eq $ ;Fall into DivrSingle
;*********
DivrSingle:
;*********
;Computes op1/op2
;Op1 is double, op2 is single (low 32 bits are zero)
mov edx,ebx
mov eax,esi ;Mantissa in edx:eax
mov ebx,EMSEG:[edi].ExpSgn
mov edi,EMSEG:[edi].lManHi
jmp DivSingleReg
SDivBigUnderflow:
;Overflow flag set could only occur with denormals (true exp < -32768)
or EMSEG:[CURerr],Underflow
test EMSEG:[CWmask],Underflow ;Is exception masked?
jnz UnderflowZero ;Yes, return zero (in emfmul.asm)
add ecx,Underbias shl 16 ;Fix up exponent
jmp ContSdiv ;Continue with multiply
EM_ENTRY eFDIVPreg
eFDIVPreg:
push offset PopWhenDone
EM_ENTRY eFDIVreg
eFDIVreg:
xchg esi,edi
EM_ENTRY eFDIVtop
eFDIVtop:
mov ecx,EMSEG:[esi].ExpSgn
mov ebx,EMSEG:[esi].lManHi
mov esi,EMSEG:[esi].lManLo
DivSetResult:
;cl has tag of divisor
mov ebp,offset tFdivDisp
mov EMSEG:[Result],edi ;Save result pointer
mov al,cl
mov ah,EMSEG:[edi].bTag
and ah,not 1 ;Ignore single vs. double on dividend
cmp ax,1
.erre bTAG_VALID eq 1
.erre bTAG_SNGL eq 0
jz DivDouble ;Divisor was double
ja TwoOpResultSet
;.erre DivSingle eq $ ;Fall into DivSingle
;*********
DivSingle:
;*********
;Computes op2/op1
;Op2 is double, op1 is single (low 32 bits are zero)
xchg edi,ebx ;Mantissa in edi, op2 ptr to ebx
xchg ebx,ecx ;ExpSgn to ebx, op2 ptr to ecx
mov edx,EMSEG:[ecx].lManHi
mov eax,EMSEG:[ecx].lManLo
mov ecx,EMSEG:[ecx].ExpSgn ;Op2 loaded
DivSingleReg:
;dividend mantissa in edx:eax, exponent in high ecx, sign in ch bit 7
;divisor mantissa in edi, exponent in high ebx, sign in bh bit 7
xor ch,bh ;Compute result sign
xor bx,bx ;Clear out sign and tag
sub ecx,1 shl 16 ;Exponent adjustment needed
sub ecx,ebx ;Compute result exponent
.erre TexpBias eq 0 ;Exponents not biased
jo SDivBigUnderflow ;Dividing denormal by large number
ContSdiv:
;If dividend >= divisor, the DIV instruction will overflow. Check for
;this condition and shift the dividend right one bit if necessary.
;
;In previous versions of this algorithm for 24-bit and 53-bit mantissas,
;this shift was always performed without a test. This meant that a 1-bit
;normalization might be required at the end. This worked fine because
;32 or 64 bits were calculated, so extra precision was available for
;normalization. However, this version needs all 64 bits that are calculated,
;so we can't afford a normalization shift at the end. This test tells us
;up front how to align so we'll be normalized.
xor ebx,ebx ;Extend dividend
cmp edi,edx ;Will DIV overflow?
ja DoSdiv ;No, we're safe
shrd ebx,eax,1
shrd eax,edx,1
shr edx,1
add ecx,1 shl 16 ;Bump exponent to account for shift
DoSdiv:
div edi
xchg ebx,eax ;Save quotient in ebx, extend remainder
div edi
mov esi,eax
;We have a 64-bit quotient in ebx:esi. Now compare remainder*2 with divisor
;to compute round and sticky bits.
mov eax,-1 ;Set round and sticky bits
shl edx,1 ;Double remainder
jc RoundJmp ;If too big, round & sticky set
cmp edx,edi ;Is remainder*2 > divisor?
ja RoundJmp
;Observe, oh wondering one, how you can assume the result of this last
;compare is not equality. Use the following notation: n=numerator,
;d=denominator,q=quotient,r=remainder,b=base(2^64 here). If
;initially we had n < d then there was no shift and we will find q and r
;so that q*d+r=n*b, if initially we had n >= d then there was a shift and
;we will find q and r so that q*d+r=n*b/2. If we have equality here
;then r=d/2 ==> n={possibly 2*}(2*q+1)*d/(2*b), since this can only
;be integral if d is a multiple of b, but by definition b/2 <= d < b, we
;have a contradiction. Equality is thus impossible at this point.
cmp edx,1 ;Check for zero remainder
sbb eax,-2 ;eax==0 if CY, ==1 if NC (was -1)
RoundJmp:
jmp EMSEG:[RoundMode]
;*******************************************************************************
DDivBigUnderflow:
;Overflow flag set could only occur with denormals (true exp < -32768)
or EMSEG:[CURerr],Underflow
test EMSEG:[CWmask],Underflow ;Is exception masked?
jnz UnderflowZero ;Yes, return zero (in emfmul.asm)
add ecx,Underbias shl 16 ;Fix up exponent
jmp ContDdiv ;Continue with multiply
DivrDoubleSetFlag:
;Special entry point used by FPATAN to set bit 6 of flag dword pushed
;on stack before call.
or byte ptr [esp+4],40H
;*********
DivrDouble:
;*********
;Computes op1/op2
mov edx,ebx
mov eax,esi ;Mantissa in edx:eax
mov ebx,EMSEG:[edi].ExpSgn
mov esi,EMSEG:[edi].lManHi
mov edi,EMSEG:[edi].lManLo
jmp short DivDoubleReg
HighHalfEqual:
;edx:eax:ebp = dividend
;esi:edi = divisor
;ecx = exponent and sign of result
;
;High half of dividend is equal to high half of divisor. This will cause
;the DIV instruction to overflow. If whole dividend >= whole divisor, then
;we just shift the dividend right 1 bit.
cmp eax,edi ;Is dividend >= divisor?
jae ShiftDividend ;Yes, divide it by two
;DIV instruction would overflow, so skip it and calculate the effective
;result. Assume a quotient of 2^32-1 and calculate the remainder. See
;detailed comments under MaxQuo below--this is a copy of that code.
push ecx ;Save exp. and sign
mov ebx,-1 ;Max quotient digit
sub eax,edi ;Calculate correct remainder
;Currently edx == esi, but the next instruction ensures that is no longer
;true, since eax != 0. This will allow us to skip the MaxQuo check at
;DivFirstDigit.
add edx,eax ;Should set CY if quotient fit
mov eax,edi ;ecx:eax has new remainder
jc ComputeSecond ;Remainder was positive
;Quotient doesn't fit. Note that we can no longer ensure that edx != esi
;after making a correction.
mov ecx,edx ;Need remainder in ecx:eax
jmp DivCorrect1
;*********
DivDouble:
;*********
;Computes op2/op1
mov eax,edi ;Move op2 pointer
mov edi,esi
mov esi,ebx ;Mantissa in esi:edi
mov ebx,ecx ;ExpSgn to ebx
mov ecx,EMSEG:[eax].ExpSgn ;Op2 loaded
mov edx,EMSEG:[eax].lManHi
mov eax,EMSEG:[eax].lManLo
DivDoubleReg:
;dividend mantissa in edx:eax, exponent in high ecx, sign in ch bit 7
;divisor mantissa in esi:edi, exponent in high ebx, sign in bh bit 7
xor ch,bh ;Compute result sign
xor bx,bx ;Clear out sign and tag
sub ecx,1 shl 16 ;Exponent adjustment needed
sub ecx,ebx ;Compute result exponent
.erre TexpBias eq 0 ;Exponents not biased
jo DDivBigUnderflow ;Dividing denormal by large number
ContDdiv:
;If dividend >= divisor, we must shift the dividend right one bit.
;This will ensure the result is normalized.
;
;In previous versions of this algorithm for 24-bit and 53-bit mantissas,
;this shift was always performed without a test. This meant that a 1-bit
;normalization might be required at the end. This worked fine because
;32 or 64 bits were calculated, so extra precision was available for
;normalization. However, this version needs all 64 bits that are calculated,
;so we can't afford a normalization shift at the end. This test tells us
;up front how to align so we'll be normalized.
xor ebp,ebp ;Extend dividend
cmp esi,edx ;Dividend > divisor
ja DoDdiv
jz HighHalfEqual ;Go compare low halves
ShiftDividend:
shrd ebp,eax,1
shrd eax,edx,1
shr edx,1
add ecx,1 shl 16 ;Bump exponent to account for shift
DoDdiv:
push ecx ;Save exp. and sign
;edx:eax:ebp = dividend
;esi:edi = divisor
;
;Division algorithm from Knuth vol. 2, p. 237, using 32-bit "digits":
;Guess a quotient digit by dividing two MSDs of dividend by the MSD of
;divisor. If divisor is >= 1/2 the radix (radix = 2^32 in this case), then
;this guess will be no more than 2 larger than the correct value of that
;quotient digit (and never smaller). Divisor meets magnitude condition
;because it's normalized.
div esi ;Guess first quotient "digit"
;Check out our guess.
;Currently, remainder in edx = dividend - (quotient * high half divisor).
;The definition of remainder is dividend - (quotient * all divisor). So
;if we subtract (quotient * low half divisor) from edx, we'll get
;the true remainder. If it's negative, our guess was too big.
mov ebx,eax ;Save quotient
mov ecx,edx ;Save remainder
mul edi ;Quotient * low half divisor
sub ebp,eax ;Subtract from dividend extension
sbb ecx,edx ;Subtract from remainder
mov eax,ebp ;Low remainder to eax
jnc DivFirstDigit ;Was quotient OK?
DivCorrect1:
dec ebx ;Quotient was too big
add eax,edi ;Add divisor back into remainder
adc ecx,esi
jnc DivCorrect1 ;Repeat if quotient is still too big
DivFirstDigit:
cmp ecx,esi ;Would DIV instruction overflow?
jae short MaxQuo ;Yes, figure alternate quotient
mov edx,ecx ;Remainder back to edx:eax
;Compute 2nd quotient "digit"
ComputeSecond:
div esi ;Guess 2nd quotient "digit"
mov ebp,eax ;Save quotient
mov ecx,edx ;Save remainder
mul edi ;Quotient * low half divisor
neg eax ;Subtract from dividend extended with 0
sbb ecx,edx ;Subtract from remainder
jnc DivSecondDigit ;Was quotient OK?
DivCorrect2:
dec ebp ;Quotient was too big
add eax,edi ;Add divisor back into remainder
adc ecx,esi
jnc DivCorrect2 ;Repeat if quotient is still too big
DivSecondDigit:
;ebx:ebp = quotient
;ecx:eax = remainder
;esi:edi = divisor
;Now compare remainder*2 with divisor to compute round and sticky bits.
mov edx,-1 ;Set round and sticky bits
shld ecx,eax,1 ;Double remainder
jc DDivEnd ;If too big, round & sticky set
shl eax,1
sub edi,eax
sbb esi,ecx ;Subtract remainder*2 from divisor
jb DDivEnd ;If <0, use round & sticky bits set
;Observe, oh wondering one, how you can assume the result of this last
;compare is not equality. Use the following notation: n=numerator,
;d=denominator,q=quotient,r=remainder,b=base(2^64 here). If
;initially we had n < d then there was no shift and we will find q and r
;so that q*d+r=n*b, if initially we had n >= d then there was a shift and
;we will find q and r so that q*d+r=n*b/2. If we have equality here
;then r=d/2 ==> n={possibly 2*}(2*q+1)*d/(2*b), since this can only
;be integral if d is a multiple of b, but by definition b/2 <= d < b, we
;have a contradiction. Equality is thus impossible at this point.
;No round bit, but set sticky bit if remainder != 0.
or eax,ecx ;Is remainder zero?
add eax,-1 ;Set CY if non-zero
adc edx,1 ;edx==0 if NC, ==1 if CY (was -1)
DDivEnd:
mov esi,ebp ;Result in ebx:esi
mov eax,edx ;Round/sticky bits to eax
pop ecx ;Recover sign/exponent
jmp EMSEG:[RoundMode]
MaxQuo:
;ebx = first quotient "digit"
;ecx:eax = remainder
;esi:edi = divisor
;On exit, ebp = second quotient "digit"
;
;Come here if divide instruction would overflow. This must mean that ecx == esi,
;i.e., the high halves of the dividend and divisor are equal. Assume a result
;of 2^32-1, thus remainder = dividend - ( divisor * (2^32-1) )
; = dividend - divisor * 2^32 + divisor. Since the high halves of the dividend
;and divisor are equal, dividend - divisor * 2^32 can be computed by
;subtracting only the low halves. When adding divisor (in esi) to this, note
;that ecx == esi, and we want the result in ecx anyway.
;
;Note also that since the dividend is a previous remainder, the
;dividend - divisor * 2^32 calculation must always be negative. Thus the
;addition of divisor back to it should generate a carry if it goes positive.
mov ebp,-1 ;Max quotient digit
sub eax,edi ;Calculate correct remainder
add ecx,eax ;Should set CY if quotient fit
mov eax,edi ;ecx:eax has new remainder
jc DivSecondDigit ;Remainder was positive
jmp DivCorrect2
