'{$STAMP BS2pe}
'{$PBASIC 2.5}
' divide40x24.bse  (c) 2003 Tracy Allen http://www.owlogic.com
' computes long division, 40 bit numerator by 24 bit denominator
' N = D * Q + R, find Q and R given N and D
' This BASIC Stamp algorithm uses byte by byte method, instead of word by word
' (compare with divide32x32.bse)    Watch out for carrys!
' The byte by byte algorithm is easily ported to pic or other byte oriented micro.
' debug code commented back in will show intermediate results of calculations

n0       VAR  Byte		' numerator entered by user
n1       VAR  Byte
n2       VAR  Byte
n3       VAR  Byte
n4       VAR  Byte

d0       VAR  Byte		' denominator entered by user
d1       VAR  Byte
d2       VAR  Byte

q0       VAR  Byte		' quotient to be computed
q1       VAR  Byte
q2       VAR  Byte
q3       VAR  Byte
q4       VAR  Byte

r0       VAR  Byte		' remainder to be computed
r1       VAR  Byte
r2       VAR  Byte
r3       VAR  Byte		' this could be a bit

x0       VAR  Byte		' working variable
x1       VAR  Byte
x2       VAR  Byte
x3       VAR  Byte		' this could be a bit

cb       VAR  Bit		' carry/borrow

kn       VAR  Byte		' number of significant bits
kw       CON	40		' word size of numerator
idx      VAR  Byte		' index
ix       VAR  Byte		' for debug display

top:
' demo asks user to input two numbers to test
SEROUT 16,$54,[CR,"enter 10 hex digits for numerator, 6 for denominator (0..9 a..f)"]

SEROUT 16,$54,[CR,"numerator 0000000000: "]
SERIN 16,$54,[HEX2 n4,HEX2 n3,HEX2 n2,HEX2 n1,HEX2 n0]

SEROUT 16,$54,[CR,"denominator 000000: "]
SERIN 16,$54,[HEX2 d2,HEX2 d1,HEX2 d0]

SEROUT 16,$54,[CR,"thank you..."]

test_divide_by_zero:
IF d0|d1|d2=0 THEN
  SEROUT 16,$54,[CR,"divide by zero" ,CR]
  GOTO top
  ENDIF

find_N_bits:
GOSUB find_nbits	' find number of significant bits in N
' DEBUG 32,DEC kn  ,32			' show number

GOSUB zero_accumulators

show_N_zero:
IF kn=0 THEN
  SEROUT 16,$54,[CR,"numerator is zero"]
  GOTO display
  ENDIF

test_N_less_than_D:
' pretest n<d
' not necessary to do this
' because the main algorithm will resolve the quotient and remainder correctly
' but this can save time
FOR idx=0 TO 3
  r0(idx)=n0(idx)
  NEXT
  GOSUB R_minus_D
  IF cb=1 THEN
    SEROUT 16,$54,[CR,"numerator<denominator"]
    FOR idx=0 TO 2
      r0(idx)=n0(idx)
    NEXT
    GOTO display
  ENDIF

main_initialize:
GOSUB zero_accumulators
IF kn<kw THEN
  FOR idx=1 TO kw-kn	' left justify N
    GOSUB Nshift
  NEXT
ENDIF
' FOR idx=4 TO 0		' show left justify
'   DEBUG HEX2 n0(idx)
' NEXT

main_calculation:
  FOR idx=kn TO 1
    GOSUB Nshift
    GOSUB Rshift
    GOSUB R_minus_D
    cb=~cb
    GOSUB Qshift
    IF cb=1 THEN	' borrow did not occur; R>=D; R:=R-D=X
      r0=x0
      r1=x1
      r2=x2
    ENDIF
'      DEBUG CR,TAB
'      FOR ix=3 TO 0	' show remainder
'        DEBUG HEX2 r0(ix)
'      NEXT
'       DEBUG TAB
'      FOR ix=3 TO 0	' show workspace
'        DEBUG HEX2 x0(ix)
'      NEXT
'      DEBUG TAB,BIN1 cb,TAB
'      FOR ix=4 TO 0	' show quotient
'        DEBUG HEX2 q0(ix)
'      NEXT
  NEXT

display:
  SEROUT 16,$54,[CR,"  quotient: ",HEX2 q4,HEX2 q3,HEX2 q2,HEX2 q1,HEX2 q0]
  SEROUT 16,$54,[CR,"  remainder: ",HEX2 r2,HEX2 r1, HEX2 r0,CR,CR]

GOTO top


' -subroutines-------------------------------------------

zero_accumulators:
FOR idx=0 TO 8	' zero accumulators, quotient & remainder
  q0(idx)=0
NEXT
RETURN

' find number of significant bits in numerator n4:n3:n2:n1:n0 up to 40 bits
' a result kn=0 means numerator is zero
find_nbits:
  kn=0
  idx=4
find_nbits1:
  kn=NCD n0(idx)
  IF kn=0 AND idx>0 THEN
    idx=idx-1
    GOTO find_nbits1
  ENDIF
    kn=idx*8+kn
RETURN

' shift numerator n4:n3:n2:n1:n0 one bit left into cb, zero into lsbit
Nshift:
  cb=n4.BIT7
  n4=n4<<1|n3.BIT7
  n3=n3<<1|n2.BIT7
  n2=n2<<1|n1.BIT7
  n1=n1<<1|n0.BIT7
  n0=n0<<1
RETURN

' shift quotient q4:q3:q2:q1:q0 one bit left, and cb into lsbit
Qshift:
  q4=q4<<1|q3.BIT7
  q3=q3<<1|q2.BIT7
  q2=q2<<1|q1.BIT7
  q1=q1<<1|q0.BIT7
  q0=q0<<1|cb
RETURN

' shift remainder r2:r1:r0 left, and cb into lsbit
' shift can move one more bit into r3
' but subtraction of d2:d1:d0 will always keep result to two bytes
' workspace is at least one bit larger than D
Rshift:
  r3=r2.BIT7
  r2=r2<<1|r1.BIT7
  r1=r1<<1|r0.BIT7
  r0=r0<<1|cb
RETURN

' compute R-D and set cb=1 if a borrow occurs (D>R)
' x3:x2:x1:x0 is a workspace variable, as with R, at least one bit larger than D.
R_minus_D:
  x0=r0-d0
  IF d0>r0 THEN cb=1 ELSE cb=0   ' borrow- (x0 min r0 - r0 max 1)

  x1=d1+cb
  IF x1<d1 THEN cb=1 ELSE cb=0
  x1=r1-x1
  IF x1>r1 THEN cb=1

  x2=d2+cb
  IF x2<d2 THEN cb=1 ELSE cb=0
  x2=r2-x2
  IF x2>r2 THEN cb=1

  x3=r3-cb
  IF x3>r3 THEN cb=1 ELSE cb=0 	' could use cb=x3.bit7 so long as x3 is a byte
RETURN

' alternative form R-D using computed carry instead of IF-THEN-ELSE
' note Stamp does not have a carry/borrow bit
' remember when three terms are in a sum, either operation can generate a carry
' if one of the terms is 1, then there can only be one carry
R_minus_D_alternate:
  x0=r0-d0
  cb=d0 MIN r0 - r0 MAX 1

  x1=d1+cb
  cb=x1 MAX d1 - d1 MAX 1
  x1=r1-x1
  cb=x1 MIN r1 - r1 MAX 1 | cb

  x2=d2+cb
  cb=x2 MAX d2 - d2 MAX 1
  x2=r2-x2
  cb=x2 MIN r2 - r2 MAX 1 | cb

  x3=r3-cb
  cb=cb MIN r3 -r3 MAX 1
RETURN

DEBUG HOME