emu8051/tests/sqroot.asm
2014-01-26 00:25:48 -05:00

375 lines
8.1 KiB
NASM

; Test program to verify correct emu8051 operation
; Taken from http://www.vzsite.us/8051/
;
; Test desc: 32-bit division
; Test output: PC = $FFF0
; Test output: SP = $30
; Test output: PSW = $04
; Test output: A = $44
; Test output: B = $44
;17 Jan 00 re-written for consistency with assembler/compiler byte order
;21 May 99 8051 source code
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;SQROOT1 is called to calculate the square root of a 32-bit number.
;
;Call:
;
; r0 => MSB of 32-bit input
; lcall SQROOT1
;
;Return:
;
; computed square root is in acc(LSB's) and b(MSB's)
; (root is also in est(MSB's) and est+1(LSB's))
;
;SQROOT1 uses the formula
;
; estimate = (last_estimate + input/last_estimate)/2
;
;Method is described as Newton's, Newton-Raphson, and Babylonian method.
;DIV32U is called to do the 32/16-bit division.
;SQROOT1 perfoms a fixed number of iterations to converge to the root.
;
;SQROOT1 will destroy all registers except r0, dptr.
;
;Original Author: John Veazey, Ridgecrest, CA, 18 Apr 99
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
q0 EQU 70h
q1 EQU 71h
est EQU 72h ; 2 bytes: The estimated value being developed
sqrcnt EQU 74h ; Counts iterations
square EQU 75h
TOS EQU 30h ; Adresse du dessus de la pile.
ORG 0000h ; Reset vector
MOV SP,#TOS ; Init stack pointer
;; Set 32-bit input value
MOV square+0,#012h
MOV square+1,#034h
MOV square+2,#032h
MOV square+3,#010h
;; Set pointer to square
mov r0,#square
LCALL SQROOT
LJMP 0FFF0h
SQROOT:
;
;Save r0 on the stack
;
mov a,r0
push acc
;
;Initialize the iteration counter
;
mov sqrcnt,#7 ;There will be seven calculations
;
;Find the MSB set in input, call it n as in 2**n. Calculate
; r = ceiling(n+1)/2, then set up the first estimate to
; be (2**r)-1.
;
mov a,@r0 ;Is MS byte of input not-zero?
jz sqr112
mov r1,#33 ; Yes, set r1 to the bit number+2
sjmp sqr140
sqr112:
inc r0 ; No,
mov a,@r0 ;Is 2nd MS byte of input not-zero?
jz sqr114
mov r1,#25 ; Yes
sjmp sqr140
sqr114:
inc r0
mov a,@r0 ;Is 3rd MS byte of input not-zero?
jz sqr116
mov r1,#17 ; Yes
sjmp sqr140
sqr116:
inc r0
mov est,#0
mov a,@r0 ;Is LS byte of input zero?
jnz sqr122
mov est+1,#0 ; Yes, return with a zero
ljmp sqr900 ; because input is zero
sqr122:
mov est+1,#1
dec a ;Is LS byte of input = 1?
jz sqr900 ; Yes, return with a one
inc a ; No
mov r1,#9
sqr140:
dec r1
rlc a
jnc sqr140
mov a,r1 ;Form ceiling[(n+1)/2]
clr c
rrc a ; (divide by 2, add remainder)
jnc sqr142
inc a
sqr142:
mov r1,a
clr a
mov est,a
sqr144:
setb c ;Get the 1 to shift in
rlc a
mov b,a
mov a,est
rlc a
mov est,a
mov a,b
dec r1
cjne r1,#0,sqr144
mov est+1,a
;
;Load input into DIV32U dividend register
;
sqr200:
pop acc
push acc
mov r0,a
mov a,@r0 ;Set MSB's
mov r7,a
inc r0
mov a,@r0
mov r6,a
inc r0
mov a,@r0
mov r5,a
inc r0
mov a,@r0
mov r4,a
;
;Load last estimate into DIV32U divisor registers
;
mov a,est+1 ;Set LSB's
mov r0,a
mov a,est ;Set MSB's
mov r1,a
;
;Call DIV32U to do the 32/16-bit division (input/est)
;
lcall DIV32U ;(r7,r6,r5,r4)/(r1,r0) = (r5,r4)
jnc sqr310 ;If overflow, set to maximum number
mov r4,#-1
mov r5,#-1
;
;Add the last estimate to the quotient
;
sqr310:
mov a,r4
add a,est+1 ;sets c
mov est+1,a
mov a,est
addc a,r5 ;sets c
;
;Divide sum by 2 and save as new estimate
;
rrc a ;sets c
mov est,a
mov a,est+1
rrc a ;c discarded
mov est+1,a
;
;Decrement the iteration counter and repeat if not zero
;
djnz sqrcnt,sqr200
;
;Make the normal return
;
sqr900:
pop acc ;Restore caller's pointer
mov r0,a
mov b,est ;Get answer in (b,a)
mov a,est+1
ret
DIV32U:
;
;Clear the working quotient
;
clr a
mov q1,a
mov q0,a
;
;Clear the msb's of a 32-bit working divisor (r3,r2,r1,r0)
;
mov r3,a
mov r2,a
;
;b counts the number of places+1 the divisor was initially
; shifted left to align its ms bit set with the ms bit set
; in the dividend
;
mov b,#1
;
;Make an error return if trying to divide by zero
;
mov a,r1
orl a,r0
jnz du100
ljmp du920 ;Make the error return
;
;Just return with quotient and remainder zero if dividend is zero
;
du100:
mov a,r7
orl a,r6
orl a,r5
orl a,r4
jnz du200
mov r7,a
mov r6,a
ljmp du910 ;Make a normal return
;
;Align the msb set in the demoninator with the msb set in the
; numerator. Increment the shift count in b each time a shift left
; is performed.
;
du200:
mov a,r3 ;Stop if msb set
clr c
rlc a
jc du600
subb a,r7 ;Compare r3 & r7, (c clear)
jz du210 ; jump if r3=r7
jnc du600 ; jump if r3>r7
sjmp du240 ; r3<r7
du210:
mov a,r6 ;r3=r7, so compare r2 & r6
subb a,r2
jc du600 ; jump if r2>r6
jnz du240 ; jump if r2<r6
mov a,r5 ;r2=r6, so compare r1 & r5
subb a,r1
jc du600 ; jump if r1>r5
jnz du240 ; jump if r1<r5
mov a,r4 ;r1=r5, so compare r0 & r4
subb a,r0
jc du600 ; jump if r0>r4
du240:
clr c ;Now shift the denominator
mov a,r0 ; left 1 bit position
rlc a
mov r0,a
mov a,r1
rlc a
mov r1,a
mov a,r2
rlc a
mov r2,a
mov a,r3
rlc a
mov r3,a
inc b ;Increment b counter and
sjmp du200 ; continue
;
;Compare the shifted divisor with the remainder (what's
; left of the dividend)
;
du600:
mov a,r7
clr c
subb a,r3
jc du720 ;jump if r3>r7
jnz du700 ;jump if r3<r7
mov a,r6
subb a,r2
jc du720 ;jump if r2>r6
jnz du700 ;jump if r2<r6
mov a,r5
subb a,r1
jc du720 ;jump if r1>r5
jnz du700 ;jump if r1<r5
mov a,r4
subb a,r0
jc du720 ;jump if r0>r4
;
;Divisor is equal or smaller, so subtract it off and
; get a 1 for the quotient
;
du700:
mov a,r4
clr c
subb a,r0
mov r4,a
mov a,r5
subb a,r1
mov r5,a
mov a,r6
subb a,r2
mov r6,a
mov a,r7
subb a,r3
mov r7,a
clr c
cpl c ;Get a 1 for the quotient
sjmp du730
;
;Divisor is greater, get a 0 for the quotient
;
du720:
clr c
;
;Shift 0 or 1 into quotient
;
du730:
mov a,q0
rlc a
mov q0,a
mov a,q1
rlc a
mov q1,a
jc du920 ;overflow - make the error return
;
;Now shift the denominator right 1, decrement the counter
; in b until b = 0
;
du740:
clr c
mov a,r3
rrc a
mov r3,a
mov a,r2
rrc a
mov r2,a
mov a,r1
rrc a
mov r1,a
mov a,r0
rrc a
mov r0,a
djnz b,du600
;
;Move quotient and remainder
;
mov a,r5
mov r7,a
mov a,r4
mov r6,a
mov a,q1
mov r5,a
mov a,q0
mov r4,a
;
;Make the normal return
;
du910:
clr c
ret
;
;Make the error return
;
du920:
clr c
cpl c
ret