6
\$\begingroup\$

I am trying to create an optimized 16-bit division algorithm for the AVR ATMEGA1284. The goal is to reduce the number of clock cycles as much as possible.

AVR INSTRUCTION SET MANUAL: https://ww1.microchip.com/downloads/en/devicedoc/atmel-0856-avr-instruction-set-manual.pdf

AVR200: Multiply and Divide Routines: https://ww1.microchip.com/downloads/en/AppNotes/doc0936.pdf

A standard shift and subtract type division algorithm suggested by Atmel/Microchip takes between 173 clock cycles and 243 clock cycles, depending on if you unroll the loop or not.

What I have so far takes a maximum of 68 clock cycles. I ran an exhaustive test 16-hour test proving that the algorithm returns the correct result for all 2^32 combinations of inputs and outputs. So I am not looking for any validation that the algorithm returns the correct results. I am looking for ways to reduce either the code size, lookup table size, or number of clock cycles.

The constraints are as follows.

  • The dividend is an unsigned 16-bit number passed into the algorithm in a pair of 8-bit registers.
  • The divisor is an unsigned 16-bit number passed into the algorithm in a pair of 8-bit registers.
  • The algorithm returns the quotient in a pair of 8-bit registers.
  • The algorithm also returns the remainder in a pair of 8-bit registers.
  • The algorithm code (plus any lookup tables) must consume less than 5K bytes of flash memory.
  • The algorithm may return any values for division by 0.

Here is what I have so far.

.align 256
;Recipricol table #1, high byte.
;R1H_TBL[x] = min( high(2^16/x) / 256 , 255)
R1H_TBL:
.db 0xFF, 0xFF, 0x80, 0x55, 0x40, 0x33, 0x2A, 0x24, 0x20, 0x1C, 0x19, 0x17, 0x15, 0x13, 0x12, 0x11
.db 0x10, 0x0F, 0x0E, 0x0D, 0x0C, 0x0C, 0x0B, 0x0B, 0x0A, 0x0A, 0x09, 0x09, 0x09, 0x08, 0x08, 0x08
.db 0x08, 0x07, 0x07, 0x07, 0x07, 0x06, 0x06, 0x06, 0x06, 0x06, 0x06, 0x05, 0x05, 0x05, 0x05, 0x05
.db 0x05, 0x05, 0x05, 0x05, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04
.db 0x04, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03
.db 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02
.db 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02
.db 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02
.db 0x02, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
;Recipricol table #1, low byte.
;R1L_TBL[x] = min( low(2^16/x) mod 256 , 255)
R1L_TBL:
.db 0xFF, 0xFF, 0x00, 0x55, 0x00, 0x33, 0xAA, 0x92, 0x00, 0x71, 0x99, 0x45, 0x55, 0xB1, 0x49, 0x11
.db 0x00, 0x0F, 0x38, 0x79, 0xCC, 0x30, 0xA2, 0x21, 0xAA, 0x3D, 0xD8, 0x7B, 0x24, 0xD3, 0x88, 0x42
.db 0x00, 0xC1, 0x87, 0x50, 0x1C, 0xEB, 0xBC, 0x90, 0x66, 0x3E, 0x18, 0xF4, 0xD1, 0xB0, 0x90, 0x72
.db 0x55, 0x39, 0x1E, 0x05, 0xEC, 0xD4, 0xBD, 0xA7, 0x92, 0x7D, 0x69, 0x56, 0x44, 0x32, 0x21, 0x10
.db 0x00, 0xF0, 0xE0, 0xD2, 0xC3, 0xB5, 0xA8, 0x9B, 0x8E, 0x81, 0x75, 0x69, 0x5E, 0x53, 0x48, 0x3D
.db 0x33, 0x29, 0x1F, 0x15, 0x0C, 0x03, 0xFA, 0xF1, 0xE8, 0xE0, 0xD8, 0xD0, 0xC8, 0xC0, 0xB9, 0xB1
.db 0xAA, 0xA3, 0x9C, 0x95, 0x8F, 0x88, 0x82, 0x7C, 0x76, 0x70, 0x6A, 0x64, 0x5E, 0x59, 0x53, 0x4E
.db 0x49, 0x43, 0x3E, 0x39, 0x34, 0x30, 0x2B, 0x26, 0x22, 0x1D, 0x19, 0x14, 0x10, 0x0C, 0x08, 0x04
.db 0x00, 0xFC, 0xF8, 0xF4, 0xF0, 0xEC, 0xE9, 0xE5, 0xE1, 0xDE, 0xDA, 0xD7, 0xD4, 0xD0, 0xCD, 0xCA
.db 0xC7, 0xC3, 0xC0, 0xBD, 0xBA, 0xB7, 0xB4, 0xB2, 0xAF, 0xAC, 0xA9, 0xA6, 0xA4, 0xA1, 0x9E, 0x9C
.db 0x99, 0x97, 0x94, 0x92, 0x8F, 0x8D, 0x8A, 0x88, 0x86, 0x83, 0x81, 0x7F, 0x7D, 0x7A, 0x78, 0x76
.db 0x74, 0x72, 0x70, 0x6E, 0x6C, 0x6A, 0x68, 0x66, 0x64, 0x62, 0x60, 0x5E, 0x5C, 0x5A, 0x58, 0x57
.db 0x55, 0x53, 0x51, 0x50, 0x4E, 0x4C, 0x4A, 0x49, 0x47, 0x46, 0x44, 0x42, 0x41, 0x3F, 0x3E, 0x3C
.db 0x3B, 0x39, 0x38, 0x36, 0x35, 0x33, 0x32, 0x30, 0x2F, 0x2E, 0x2C, 0x2B, 0x29, 0x28, 0x27, 0x25
.db 0x24, 0x23, 0x21, 0x20, 0x1F, 0x1E, 0x1C, 0x1B, 0x1A, 0x19, 0x18, 0x16, 0x15, 0x14, 0x13, 0x12
.db 0x11, 0x0F, 0x0E, 0x0D, 0x0C, 0x0B, 0x0A, 0x09, 0x08, 0x07, 0x06, 0x05, 0x04, 0x03, 0x02, 0x01
;Recipricol table #2
;R2_TBL[x] = min( 2^16/(x+256), 255)
R2_TBL:
.db 0xFF, 0xFF, 0xFE, 0xFD, 0xFC, 0xFB, 0xFA, 0xF9, 0xF8, 0xF7, 0xF6, 0xF5, 0xF4, 0xF3, 0xF2, 0xF1
.db 0xF0, 0xF0, 0xEF, 0xEE, 0xED, 0xEC, 0xEB, 0xEA, 0xEA, 0xE9, 0xE8, 0xE7, 0xE6, 0xE5, 0xE5, 0xE4
.db 0xE3, 0xE2, 0xE1, 0xE1, 0xE0, 0xDF, 0xDE, 0xDE, 0xDD, 0xDC, 0xDB, 0xDB, 0xDA, 0xD9, 0xD9, 0xD8
.db 0xD7, 0xD6, 0xD6, 0xD5, 0xD4, 0xD4, 0xD3, 0xD2, 0xD2, 0xD1, 0xD0, 0xD0, 0xCF, 0xCE, 0xCE, 0xCD
.db 0xCC, 0xCC, 0xCB, 0xCA, 0xCA, 0xC9, 0xC9, 0xC8, 0xC7, 0xC7, 0xC6, 0xC5, 0xC5, 0xC4, 0xC4, 0xC3
.db 0xC3, 0xC2, 0xC1, 0xC1, 0xC0, 0xC0, 0xBF, 0xBF, 0xBE, 0xBD, 0xBD, 0xBC, 0xBC, 0xBB, 0xBB, 0xBA
.db 0xBA, 0xB9, 0xB9, 0xB8, 0xB8, 0xB7, 0xB7, 0xB6, 0xB6, 0xB5, 0xB5, 0xB4, 0xB4, 0xB3, 0xB3, 0xB2
.db 0xB2, 0xB1, 0xB1, 0xB0, 0xB0, 0xAF, 0xAF, 0xAE, 0xAE, 0xAD, 0xAD, 0xAC, 0xAC, 0xAC, 0xAB, 0xAB
.db 0xAA, 0xAA, 0xA9, 0xA9, 0xA8, 0xA8, 0xA8, 0xA7, 0xA7, 0xA6, 0xA6, 0xA5, 0xA5, 0xA5, 0xA4, 0xA4
.db 0xA3, 0xA3, 0xA3, 0xA2, 0xA2, 0xA1, 0xA1, 0xA1, 0xA0, 0xA0, 0x9F, 0x9F, 0x9F, 0x9E, 0x9E, 0x9D
.db 0x9D, 0x9D, 0x9C, 0x9C, 0x9C, 0x9B, 0x9B, 0x9A, 0x9A, 0x9A, 0x99, 0x99, 0x99, 0x98, 0x98, 0x98
.db 0x97, 0x97, 0x97, 0x96, 0x96, 0x95, 0x95, 0x95, 0x94, 0x94, 0x94, 0x93, 0x93, 0x93, 0x92, 0x92
.db 0x92, 0x91, 0x91, 0x91, 0x90, 0x90, 0x90, 0x90, 0x8F, 0x8F, 0x8F, 0x8E, 0x8E, 0x8E, 0x8D, 0x8D
.db 0x8D, 0x8C, 0x8C, 0x8C, 0x8C, 0x8B, 0x8B, 0x8B, 0x8A, 0x8A, 0x8A, 0x89, 0x89, 0x89, 0x89, 0x88
.db 0x88, 0x88, 0x87, 0x87, 0x87, 0x87, 0x86, 0x86, 0x86, 0x86, 0x85, 0x85, 0x85, 0x84, 0x84, 0x84
.db 0x84, 0x83, 0x83, 0x83, 0x83, 0x82, 0x82, 0x82, 0x82, 0x81, 0x81, 0x81, 0x81, 0x80, 0x80, 0x80
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;ARGUMENTS:  r16, r17, r18, r19
;  r16:r17 = N (numerator)
;  r18:r19 = D (divisor)
;RETURNS:    r20, r21
;  r20:r21 (quotient)
;  r22:r23 (remainder)
;
;DESCRIPTION:  divides an unsigned 16 bit number N by unsigned 16 bit divisor D
;  Typical run time is 57 to 68 clock cycles.
;
;ALGORITHM OVERVIEW
;
;RZERO = 0;
;if(D < 256){
;  N2 = (N * ((R1H_TBL[D] << 8) + R1L_TBL[D])) >> 16;
;  P  = N2 * D
;}else{
;  D1 = (R1H_TBL[D] * D) >> 8
;  N1 = (R1H_TBL[D] * N) >> 8
;  if(D1 < 256){
;    N2 = N1 >> 8;
;  }else{
;    N2 = N1 * R2_TBL[D1 & 0xFF];
;  }
;  P = N2 * D;
;  if(P > 65535){
;    N2 = N2 - 1    ;//Decrement quotient
;    N1 = N2 - P + D;//Calculate remainder
;    return;//return quotient in N2, remainder in N1
;  }
;}
;N1 = N - P;
;if(P > N){
;  N2 = N2 - 1;//decrease quotient
;  N1 = N1 + D;//increase reamainder
;  return;//return quotient in N2, remainder in N1
;}
;if(N1 > D){
;  N2 = N2 + 1;
;  N1 = N1 - D;
;  return;//return quotient in N2, remainder in N1
;}
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
.def NH    = r16 .def NL    = r17
.def DH    = r18 .def DL    = r19
.def N2H   = r20 .def N2L   = r21
.def N1H   = r22 .def N1L   = r23
.def PRODL = r0  .def PRODH = r1
.def PH    = r2  .def PL    = r3
.def D1H   = r4  .def D1L   = r5
.def RZERO = r6
.def Rx    = r7 

idivu_16x16:  
    clr RZERO                 ;1
    ;Check that DH is not zero
    tst DH                    ;1
    brne idivu_16x16_dhne   ;2
    ;code for D < 256   
idivu_16x8:
    ;lookup low byte of recipricol into P.
    ;where P = min(2^16 / D,2^16-1)
    mov zl, DL               ;1
    ldi zh, high(R1L_TBL*2)  ;1 
    lpm PL, Z                ;3
    ldi zh, high(R1H_TBL*2)  ;1 
    lpm PH, Z                ;3
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;calculate N2 = (P * N) >> 16
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;     NH:NL
    ;  X  RH:RL
    ;------------------------------------------
    ;   N2H    |   N2L    |  N1H     | dropped
    ;----------+----------+----------+---------
    ;          |          | H(PL*NL) | L(PL*NL)
    ;          | H(PL*NH) | L(PL*NH) |
    ;          | H(PH*NL) | L(PH*NL) |
    ; H(PH*NH) | L(PH*NH) |          |
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
    mul NL , PL     ;2  PL*NL
    mov N1H, PRODH  ;1  N1H <= H(PL*NL)
    mul NH , PH     ;2  PH*NH
    mov N2H, PRODH  ;1  N2H <= H(PH*NH)
    mov N2L, PRODL  ;1  N2L <= L(PH*NH)
    mul NH , PL     ;2  PL*NH
    add N1H, PRODL  ;1  N1H <= H(PL*NL) + L(PL*NH) 
    adc N2L, PRODH  ;1  N2L <= L(PH*NH) + H(PL*NH)
    adc N2H, RZERO  ;1  propagate carry to N2H      
    mul NL , PH     ;2  PH*NL
    add N1H, PRODL  ;1  N1H <= H(PL*NL) + L(PL*NH) + L(PH*NL)
    adc N2L, PRODH  ;1  N2L <= H(PH*NL) + L(PH*NH) + H(PL*NH)
    adc N2H, RZERO  ;1  propagate carry to N2H  
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;calculate P = N2 * DL ,note DH=0
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;    
    mul N2L, DL              ;2
    mov PL, PRODL            ;1
    mov PH, PRODH            ;1
    mul N2H, DL              ;2
    add PH, PRODL            ;1
    rjmp idivu_16x16_adj_n ;2
    ;code for D >= 256
idivu_16x16_dhne:          
    ;Lookup Rx = min(256 / DH, 255)     
    mov zl, DH               ;1
    ldi zh, high(R1H_TBL*2)  ;1 
    lpm Rx, Z                ;3
    ;D1 = (D * Rx) >> 8          
    mul Rx , DH              ;2
    mov D1L, PRODL           ;1
    mov D1H, PRODH           ;1
    mul Rx , DL              ;2
    add D1L, PRODH           ;1
    adc D1H, RZERO           ;1
    ;N1 = (D * Rx) >> 8          
    mul Rx , NH              ;2
    mov N1L, PRODL           ;1
    mov N1H, PRODH           ;1
    mul Rx , NL              ;2
    add N1L, PRODH           ;1
    adc N1H, RZERO           ;1
    ;if D1H = 0 then use Rx = 256, otherwise use table   
    tst D1H                  ;1
    brne idivu_16x16_dxhne ;2
    mov N2L, N1H             ;1
    clr N2H                  ;1
    rjmp idivu_16x16_checkn;2
    idivu_16x16_dxhne:
    ;Lookup Rx = (2 ^ 16) \ (256 + D1H)
    mov zl, D1L              ;1
    ldi zh, high(R2_TBL*2)   ;1
    lpm Rx, Z                ;3
    ;N2 = (N1 * R2) >> 16
    mul Rx  , N1H            ;2
    mov PL  , PRODL          ;1
    mov N2L , PRODH          ;1
    mul Rx , N1L             ;2
    add PL , PRODH           ;1
    adc N2L, RZERO           ;1
    clr N2H                  ;1
    idivu_16x16_checkn:
    ;Check result (it may be off by +/- 1)
    ;P = N2 * D
    ;NOTE:  N2 <=255 so NxH = 0, also P < 2^16 so we can discard upper byte of DH * NxL
    mul DL , N2L             ;2
    mov PL, PRODL            ;1
    mov PH, PRODH            ;1
    mul DH , N2L             ;2
    add PH , PRODL           ;1 
    brcc idivu_16x16_adj_n ;2
    ;if multiply overflowed then...
    ;decrement quotient
    ;calculate remainder as N - P + D
    subi N2L, 0x01           ;1
    sbci N2H, 0x00           ;1
    mov N1L, NL              ;1
    mov N1H, NH              ;1
    sub N1L, PL              ;1
    sbc N1H, PH              ;1
    add  N1L, DL             ;1
    adc  N1H, DH             ;1
    rjmp idivu_16x16_end   ;2
    ;Adjust result up or down by 1 if needed.
    idivu_16x16_adj_n:
    ;Add -P to N, with result in P
    mov N1L, NL              ;1
    mov N1H, NH              ;1
    sub N1L, PL              ;1
    sbc N1H, PH              ;1
    brsh idivu_16x16_pltn  ;2
    idivu_16x16_decn2:
    ;if P > N then decrement quotient, add to remainder
    subi N2L, 0x01           ;1
    sbci N2H, 0x00           ;1
    add  N1L, DL             ;1
    adc  N1H, DH             ;1
    rjmp idivu_16x16_end   ;2
    idivu_16x16_pltn:
    ;test remainder to D 
    cp  N1L, DL              ;1
    cpc N1H, DH              ;1
    ;if remainder < D then goto end
    brlo idivu_16x16_end   ;2
    ;if remainder >= D then increment quotient, reduce remainder
    subi N2L, 0xFF           ;1
    sbci N2H, 0xFF           ;1
    sub N1L, DL              ;1
    sbc N1H, DH              ;1
    idivu_16x16_end:
    ret
    .undef NH    .undef NL   
    .undef DH    .undef DL   
    .undef N2H   .undef N2L  
    .undef N1H   .undef N1L  
    .undef PRODL .undef PRODH
    .undef PH    .undef PL   
    .undef D1H   .undef D1L  
    .undef RZERO 
    .undef Rx
\$\endgroup\$
11
  • \$\begingroup\$ What is the specific assembler you're using? \$\endgroup\$ – Reinderien Dec 30 '20 at 15:51
  • \$\begingroup\$ Also, why is the division being performed on this device? \$\endgroup\$ – Reinderien Dec 30 '20 at 16:05
  • \$\begingroup\$ Finally, can you give some domain information on the variables being divided? Often, this information can help optimize algebraic operations or even avoid them. \$\endgroup\$ – Reinderien Dec 30 '20 at 16:09
  • \$\begingroup\$ @Reinderien The division is being used to determine the intersection points of lines in a 2D space. This is used as part of some 2D graphics rendering routines. Specifically I have some data points from an external source and am scaling the data points in both X & Y for visual display. A 2D line is drawn between adjacent data-points to create a line graph. It may happen that one or both points are totally outside the bounds of the screen once scaled. In that case I need to determine the intersections of the lines with the display area before drawing them. \$\endgroup\$ – user4574 Jan 1 at 3:21
  • \$\begingroup\$ @Reinderien I am using AVR Assembler in Atmel Studio 7. \$\endgroup\$ – user4574 Jan 1 at 3:23
3
\$\begingroup\$

Have mercy on the maintenance programmer - It may be your older self.

Separate documentation tends to become

  • separated, as in not always readily accessible
      (Murphy: when direly needed)
  • out of sync - well, that's a problem even with in-line comments
      → document, in the code (for every part created for a separate reason)
    • what it is good for
    • where it got inspired (there may be easily accessible explanations enlightening to someone unfamiliar with the problem at hand or the approach used: even name dropping may help find such)
    • what has been the incentive to write it

Adopting good practices running counter to adverse customs isn't easy and fast.
Much material about assembly programming is pre-1980s, when there was some reason to have short mnemonics for instructions and operands. (No matter pointing pen or finger at a (printed…if you were lucky) program listing: no pop-up. So better keep things all in one line…).

Please use telling names. Coding in assembly is no licence not to.
In a division implementation, I'd not imagine problems with R for remainder or Q for quotient. Resist any impulse to outsmart everyone with the likes of DVsor. N for numerator wouldn't be bad if talking about fractions, but if N2 and N1 in addition - all three in H and L flavours - weren't bad enough, along comes

;NOTE:  N2 <=255 so NxH = 0, also P < 2^16 so we can discard upper byte of DH * NxL

P is mentioned in the ALGORITHM OVERVIEW.

In one comment, you switched from

sum = sum + term2
sum = sum + term2 + term3

to

sum = sum + term2
sum = term3 + sum + term2

Even then, I'd prefer

sum += term2 
sum += term3

I am looking for ways to reduce either

  • the code size,
  • lookup table size,
  • or number of clock cycles

One source of inspiration on how to code integer arithmetic is libgcc:
A "non-performing" division would be slightly faster than a non-restoring one, but hardly faster than about 120 cycles.

Rather than trying to understand the algorithm you sketch in OVERVIEW and thinking up shortcuts myself, I scrutinised the code presented. Did you write it from scratch, or did you take some compiler output for inspiration?

Catching my eye:

  • "register order" differs from the one implied by mul or the GCC calling convention, preventing the 1 cycle&word advantage each movw offers. As this is not included in The constraints, change either one.
  • The critical ("normal"?) path turns out to be taken branches mostly. With AVR, non-taken branches are faster.
    As it turns out, table access is on the critical path. While it would seem possible to save at least 127 bytes of R1H_TBL, it would cost speed.

With an eye on recognisability rather than mnemonic naming:

; not sure about the correct "movw syntax with defs", anyway)

idivu_16x8:                 ; 8 cycles less than idivu_16x16?
; stab at catching all the micro efficiencies:
; - AVR taken branches take longer: keep the critical path straight
; - "movw" takes one cycle where two "mov"s take two
; - where possible, arrange order of computation to render "tst" redundant
; to leverage movw, low and high registers are exchanged with respect to
;  <https://codereview.stackexchange.com/revisions/254077/4>
; digits starting an aligned in-line comment are cycle counts:
;  per instruction, cumulative in basic block,
;  and worst case cumulative from idivu_16x16

;code for D < 256   
    clr RZERO               ;1       3
;lookup low and high byte of reciprocal into P.
;where P = min(2^16 / D, 2^16-1)
    mov zl, DL              ;1       4
    ldi zh, high(R1L_TBL*2) ;1   1
    lpm PL, Z               ;3   2
    ldi zh, high(R1H_TBL*2) ;1   5
    lpm PH, Z               ;3   6  10
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;calculate N2 = (P * N) >> 16
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;     NH:NL
    ;  X  RH:RL
    ;------------------------------------------
    ;   N2H    |   N2L    |  N1H     | dropped
    ;----------+----------+----------+---------
    ; H(PH*NH) | L(PH*NH) | H(PL*NL) | L(PL*NL)
    ;          | H(PL*NH) | L(PL*NH) |
    ;          | H(PH*NL) | L(PH*NL) |
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    mul NL , PL     ;2      13  PL*NL
    mov N1H, PRODH  ;1   2      N1H <= H(PL*NL)
    mul NH , PH     ;2   3      PH*NH
    movw N2H,PRODH  ;1   5      N2H <= H(PH*NH)
    ;mov N2L,PRODL  ;1   6      N2L <= L(PH*NH)
    mul NH , PL     ;2   6      PL*NH
    add N1H, PRODL  ;1   8      N1H <= H(PL*NL) + L(PL*NH) 
    adc N2L, PRODH  ;1   9      N2L <= L(PH*NH) + H(PL*NH)
    adc N2H, RZERO  ;1  10      propagate carry to N2H      
    mul NL , PH     ;2  11      PH*NL
    add N1H, PRODL  ;1  13      N1H <= H(PL*NL) + L(PL*NH) + L(PH*NL)
    adc N2L, PRODH  ;1  14      N2L <= L(PH*NH) + H(PL*NH) + H(PH*NL)
    adc N2H, RZERO  ;1  15      propagate carry to N2H  
    ;calculate P = N2 * DL, note DH=0
    mul N2L, DL     ;2  16
    movw P, PROD    ;1  18
    ;mov PH, PRODH  ;1  19
    mul N2H, DL     ;2  20
    add PH, PRODL   ;1  22
    rjmp idivu_16x16_adj_n  ;2  23  36

d1Heq:
    mov N2L, N1H            ;1      25
    clr N2H                 ;1   1
    rjmp idivu_16x16_checkn ;2   2  27

idivu_16x16:
    ;Check that DH is not zero
    tst DH                  ;1   0   0
    breq idivu_16x8         ;2   1

;code for D >= 256
;idivu_16x16_dhne:
    clr RZERO               ;1       2 *
;Lookup Rx = min(256 / DH, 255)
    mov zl, DH              ;1       3 *
    ldi zh, high(R1H_TBL*2) ;1   1
    lpm Rx, Z               ;3   2

    ;N1 = (N? * Rx) >> 8
    mul Rx , NH             ;2   5
    movw N1L,PRODL          ;1   7
    ;mov N1H,PRODH          ;1   8
    mul Rx , NL             ;2   8
    add N1L, PRODH          ;1  10
    adc N1H, RZERO          ;1  11

    ;D1 = (D * Rx) >> 8
    mul Rx , DH             ;2  12
    movw D1L,PRODL          ;1  14
    ;mov D1H,PRODH          ;1  15
    mul Rx , DL             ;2  15
    add D1L, PRODH          ;1  17
    adc D1H, RZERO          ;1  18

    ;if D1H = 0 then use Rx = 256, otherwise use table
    ;tst D1H                ;1  19
    brne d1Heq              ;2  19  22

idivu_16x16_dxhne:
    ;Lookup Rx = (2 ^ 16) \ (256 + D1H)
    mov zl, D1L             ;1      23 *
    ldi zh, high(R2_TBL*2)  ;1   1
    lpm Rx, Z               ;3   2
    ;N2 = (N1 * R2) >> 16
    mul Rx , N1H            ;2   5
    mov PL , PRODL          ;1   7
    mov N2L, PRODH          ;1   8
    mul Rx , N1L            ;2   9
    add PL , PRODH          ;1  11
    adc N2L, RZERO          ;1  12
    clr N2H                 ;1  13  36

idivu_16x16_checkn:
    ;Check result (it may be off by +/- 1)
    ;P = N2 * D
    ;NOTE:  N2 <=255 so NxH = 0,
    ;       also P < 2^16 so we can discard upper byte of DH * NxL
    mul DL, N2L             ;2      37 *
    movw PL, PRODL          ;1   2
    ;mov PH, PRODH          ;1   3
    mul DH, N2L             ;2   3
    add PH, PRODL           ;1   5
    ;brcc idivu_16x16_adj_n ;2   6  43
    brcs idivu_16x16_mofl   ;2   6  43

;Adjust result up or down by 1 if needed.
idivu_16x16_adj_n:
    ;Add -P to N, with result in P
    ;mov N1L, NL            ;1      44 *
    movw N1H, NH            ;1
    sub N1L, PL             ;1   1
    sbc N1H, PH             ;1   2
    ;brsh idivu_16x16_pltn  ;2   3  47
    brlo idivu_16x16_decn2  ;2   3  47

idivu_16x16_pltn:
    ;test remainder to D 
    cp  N1L, DL             ;1      49 *
    cpc N1H, DH             ;1   1
    ;if remainder < D then goto end
    brlo idivu_16x16_end    ;2   2  51

    ;if remainder >= D then increment quotient, reduce remainder
    subi N2L, 0xFF          ;1   3
    sbci N2H, 0xFF          ;1   4
    sub N1L, DL             ;1   5
    sbc N1H, DH             ;1   6  55
idivu_16x16_end:
    ret                     ;       56 **

idivu_16x16_decn2:
    ;if P > N then decrement quotient, add to remainder
    subi N2L, 1             ;1      49
    sbci N2H, 0             ;1   1
    add  N1L, DL            ;1   2
    adc  N1H, DH            ;1   3
    ret                     ;    4  53


idivu_16x16_mofl:
    ;if multiply overflowed then...
    ;decrement quotient
    ;calculate remainder as N - P + D
    subi N2L, 0x01          ;1      45
    sbci N2H, 0x00          ;1   1
    mov N1L, NL             ;1   2
    mov N1H, NH             ;1   3
    sub N1L, PL             ;1   4
    sbc N1H, PH             ;1   5
    add N1L, DL             ;1   6
    adc N1H, DH             ;1   7
    ret                     ;1   8  53
\$\endgroup\$
7
  • \$\begingroup\$ Putting off a stab at decent naming in hopes of fully understanding the algorithm - don't hold your breath. \$\endgroup\$ – greybeard Jan 17 at 21:29
  • \$\begingroup\$ With the byte/register order changed, all 2-instruction adjust quotient sequences could be replaced by the respective word instruction, saving one word each, but no cycle. \$\endgroup\$ – greybeard Jan 18 at 7:25
  • \$\begingroup\$ I can adjust the naming a bit, but the issue is that I don't want to use more registers than necessary, so the same registers get used for different things at different parts of the algorithm. \$\endgroup\$ – user4574 Jan 18 at 16:38
  • \$\begingroup\$ The best way to understand the algorithm is to look at the "ALGORITHM OVERVIEW" part of the comments. Basically for a divisor that is 8 bits or less, look up the16-bit reciprocal from the table and multiply by that (pretty standard). For a divisor that is more than 8 bits we look up a pair of values R1*R2 whose product makes the correct reciprocal. This allows me to use two tables of 256 entries rather than one table of 65536 entries. R1 is based only on a lookup of the upper byte. R2 is based only on a lookup of the lower byte (after multiplying by R1). \$\endgroup\$ – user4574 Jan 18 at 16:47
  • \$\begingroup\$ The first value R1 is normalizes the divisor to a form that looks like 0b1_xxxx_xxxx. The second value R2 completes the division. The other parts of the code are just to correct for errors where the result can be off by +/- 1 count. The algorithm is my own original work, but is similar in concept to Goldschmidt division. It is not at all inspired by what GCC does, because GCCs division algorithm on this part is pretty slow (over 200 cycles). \$\endgroup\$ – user4574 Jan 18 at 16:48
0
\$\begingroup\$

For assessibility and tinkering, I tried to transliterate into Python:

''' all-16-bit numerator / divisor = (quotient, remainder) ALGORITHM OVERVIEW '''
reciprocals_h = R1H_TBL = [ 255 ] + [ min(((1 << 16) // x) >> 8, 255) for x in range(1, 256) ]
reciprocals_l = R1L_TBL = [ 255 ] + [ min(((1 << 16) // x)&0xFF, 255) for x in range(1, 256) ]
R2_TBL = [min((1 << 16) // (x + 256), 255) for x in range(0, 256)]
divisor, numerator = D, N = 5, 42

if divisor < 256:
    quotient = N2 = (numerator * ((reciprocals_h[divisor] << 8)
                                 + reciprocals_l[divisor])) >> 16
    product  = P1 = quotient * divisor
else:
    D1 = (reciprocals_h[divisor] * divisor) >> 8
    remainder = N1 = (reciprocals_h[divisor] * numerator) >> 8
    quotient = (remainder >> 8 if D1 < 256
           else remainder * R2_TBL[D1 & 0xFF])
    product = quotient * divisor
    if product > 65535:
        quotient = quotient - 1                   # Decrement quotient
        remainder = quotient - product + divisor  # Calculate remainder
        return (quotient, remainder)

remainder = numerator - product
if product > numerator:
    quotient = quotient - 1                       # decrement quotient
    remainder = remainder + divisor               # increase remainder
    return (quotient, remainder)
if remainder > divisor:
    quotient = quotient + 1
    remainder = remainder - divisor
    return (quotient, remainder)
\$\endgroup\$
2
  • \$\begingroup\$ Is this really an answer though? \$\endgroup\$ – Mast Jan 14 at 16:00
  • 1
    \$\begingroup\$ @Mast definitely not an answer: one of those extended comments SE doesn't provide for. Room for user4574 to improve on the readability of the embedded ALGORITHM OVERVIEW, starting with variable naming, before this gets deleted or rated down into oblivion. A code review from me would be an instance of Don't do as I do: Do as I say. \$\endgroup\$ – greybeard Jan 14 at 17:02
0
\$\begingroup\$

Here are the updates after adding RAM based tables and using MOVW where appropriate. The RAM based tables save a few clock cycles compared to reading flash based tables.

The high/low positions of any multi-register arguments were swapped to facilitate use of MOVW, which saves a clock cycle each time its used compared to using two separate byte moves.

The max run time is now 62 clock cycles.

A preprocessor symbol (RAM_DIVIDE_TABLE) was added to select use of tables in RAM or in ROM.

There is an init_math routine that gets called once at startup to copy the tables to RAM.

I re-ran my exhaustive 16-hour test (on the actual chip) to make sure every one of the 2^32 combinations of inputs yielded the correct output when compared to a standard shift and subtract type divide routine. There were no failures.

The tables now look like this...

#define RAM_DIVIDE_TABLE
#ifdef RAM_DIVIDE_TABLE
    .dseg
    .align 256
    R1H_TBL: .byte 256
    R1L_TBL: .byte 256
    R2_TBL:  .byte 256
    .cseg
init_math:
;Copy ROM divide tables to RAM
    clr r1;counter
    ldi zl, low(R1H_TBL_ROM*2)  ;1    
    ldi zh, high(R1H_TBL_ROM*2) ;1
    ldi yh, low(R1H_TBL)
    ldi yh, high(R1H_TBL)
init_math_loop_1:
    lpm r0, Z+                  ;3
    st Y+, r0                  ;2
    inc r1                      ;1
    brne init_math_loop_1       ;2
init_math_loop_2:
    lpm r0, Z+                  ;3
    st Y+, r0                  ;2
    inc r1                      ;1
    brne init_math_loop_2       ;2
init_math_loop_3:
    lpm r0, Z+                  ;3
    st Y+, r0                  ;2
    inc r1                      ;1
    brne init_math_loop_3       ;2
    ret
#else    
    .cseg
    .align 256
#endif
;Recipricol table #1, high byte.
;R1H_TBL[x] = min( high(2^16/x) / 256 , 255)
R1H_TBL_ROM:
.db 0xFF, 0xFF, 0x80, 0x55, 0x40, 0x33, 0x2A, 0x24, 0x20, 0x1C, 0x19, 0x17, 0x15, 0x13, 0x12, 0x11
.db 0x10, 0x0F, 0x0E, 0x0D, 0x0C, 0x0C, 0x0B, 0x0B, 0x0A, 0x0A, 0x09, 0x09, 0x09, 0x08, 0x08, 0x08
.db 0x08, 0x07, 0x07, 0x07, 0x07, 0x06, 0x06, 0x06, 0x06, 0x06, 0x06, 0x05, 0x05, 0x05, 0x05, 0x05
.db 0x05, 0x05, 0x05, 0x05, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04
.db 0x04, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x03
.db 0x03, 0x03, 0x03, 0x03, 0x03, 0x03, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02
.db 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02
.db 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02, 0x02
.db 0x02, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
.db 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01
;Recipricol table #1, low byte.
;R1L_TBL[x] = min( low(2^16/x) mod 256 , 255)
R1L_TBL_ROM:
.db 0xFF, 0xFF, 0x00, 0x55, 0x00, 0x33, 0xAA, 0x92, 0x00, 0x71, 0x99, 0x45, 0x55, 0xB1, 0x49, 0x11
.db 0x00, 0x0F, 0x38, 0x79, 0xCC, 0x30, 0xA2, 0x21, 0xAA, 0x3D, 0xD8, 0x7B, 0x24, 0xD3, 0x88, 0x42
.db 0x00, 0xC1, 0x87, 0x50, 0x1C, 0xEB, 0xBC, 0x90, 0x66, 0x3E, 0x18, 0xF4, 0xD1, 0xB0, 0x90, 0x72
.db 0x55, 0x39, 0x1E, 0x05, 0xEC, 0xD4, 0xBD, 0xA7, 0x92, 0x7D, 0x69, 0x56, 0x44, 0x32, 0x21, 0x10
.db 0x00, 0xF0, 0xE0, 0xD2, 0xC3, 0xB5, 0xA8, 0x9B, 0x8E, 0x81, 0x75, 0x69, 0x5E, 0x53, 0x48, 0x3D
.db 0x33, 0x29, 0x1F, 0x15, 0x0C, 0x03, 0xFA, 0xF1, 0xE8, 0xE0, 0xD8, 0xD0, 0xC8, 0xC0, 0xB9, 0xB1
.db 0xAA, 0xA3, 0x9C, 0x95, 0x8F, 0x88, 0x82, 0x7C, 0x76, 0x70, 0x6A, 0x64, 0x5E, 0x59, 0x53, 0x4E
.db 0x49, 0x43, 0x3E, 0x39, 0x34, 0x30, 0x2B, 0x26, 0x22, 0x1D, 0x19, 0x14, 0x10, 0x0C, 0x08, 0x04
.db 0x00, 0xFC, 0xF8, 0xF4, 0xF0, 0xEC, 0xE9, 0xE5, 0xE1, 0xDE, 0xDA, 0xD7, 0xD4, 0xD0, 0xCD, 0xCA
.db 0xC7, 0xC3, 0xC0, 0xBD, 0xBA, 0xB7, 0xB4, 0xB2, 0xAF, 0xAC, 0xA9, 0xA6, 0xA4, 0xA1, 0x9E, 0x9C
.db 0x99, 0x97, 0x94, 0x92, 0x8F, 0x8D, 0x8A, 0x88, 0x86, 0x83, 0x81, 0x7F, 0x7D, 0x7A, 0x78, 0x76
.db 0x74, 0x72, 0x70, 0x6E, 0x6C, 0x6A, 0x68, 0x66, 0x64, 0x62, 0x60, 0x5E, 0x5C, 0x5A, 0x58, 0x57
.db 0x55, 0x53, 0x51, 0x50, 0x4E, 0x4C, 0x4A, 0x49, 0x47, 0x46, 0x44, 0x42, 0x41, 0x3F, 0x3E, 0x3C
.db 0x3B, 0x39, 0x38, 0x36, 0x35, 0x33, 0x32, 0x30, 0x2F, 0x2E, 0x2C, 0x2B, 0x29, 0x28, 0x27, 0x25
.db 0x24, 0x23, 0x21, 0x20, 0x1F, 0x1E, 0x1C, 0x1B, 0x1A, 0x19, 0x18, 0x16, 0x15, 0x14, 0x13, 0x12
.db 0x11, 0x0F, 0x0E, 0x0D, 0x0C, 0x0B, 0x0A, 0x09, 0x08, 0x07, 0x06, 0x05, 0x04, 0x03, 0x02, 0x01
;Recipricol table #2
;R2_TBL[x] = min( 2^16/(x+256), 255)
R2_TBL_ROM:
.db 0xFF, 0xFF, 0xFE, 0xFD, 0xFC, 0xFB, 0xFA, 0xF9, 0xF8, 0xF7, 0xF6, 0xF5, 0xF4, 0xF3, 0xF2, 0xF1
.db 0xF0, 0xF0, 0xEF, 0xEE, 0xED, 0xEC, 0xEB, 0xEA, 0xEA, 0xE9, 0xE8, 0xE7, 0xE6, 0xE5, 0xE5, 0xE4
.db 0xE3, 0xE2, 0xE1, 0xE1, 0xE0, 0xDF, 0xDE, 0xDE, 0xDD, 0xDC, 0xDB, 0xDB, 0xDA, 0xD9, 0xD9, 0xD8
.db 0xD7, 0xD6, 0xD6, 0xD5, 0xD4, 0xD4, 0xD3, 0xD2, 0xD2, 0xD1, 0xD0, 0xD0, 0xCF, 0xCE, 0xCE, 0xCD
.db 0xCC, 0xCC, 0xCB, 0xCA, 0xCA, 0xC9, 0xC9, 0xC8, 0xC7, 0xC7, 0xC6, 0xC5, 0xC5, 0xC4, 0xC4, 0xC3
.db 0xC3, 0xC2, 0xC1, 0xC1, 0xC0, 0xC0, 0xBF, 0xBF, 0xBE, 0xBD, 0xBD, 0xBC, 0xBC, 0xBB, 0xBB, 0xBA
.db 0xBA, 0xB9, 0xB9, 0xB8, 0xB8, 0xB7, 0xB7, 0xB6, 0xB6, 0xB5, 0xB5, 0xB4, 0xB4, 0xB3, 0xB3, 0xB2
.db 0xB2, 0xB1, 0xB1, 0xB0, 0xB0, 0xAF, 0xAF, 0xAE, 0xAE, 0xAD, 0xAD, 0xAC, 0xAC, 0xAC, 0xAB, 0xAB
.db 0xAA, 0xAA, 0xA9, 0xA9, 0xA8, 0xA8, 0xA8, 0xA7, 0xA7, 0xA6, 0xA6, 0xA5, 0xA5, 0xA5, 0xA4, 0xA4
.db 0xA3, 0xA3, 0xA3, 0xA2, 0xA2, 0xA1, 0xA1, 0xA1, 0xA0, 0xA0, 0x9F, 0x9F, 0x9F, 0x9E, 0x9E, 0x9D
.db 0x9D, 0x9D, 0x9C, 0x9C, 0x9C, 0x9B, 0x9B, 0x9A, 0x9A, 0x9A, 0x99, 0x99, 0x99, 0x98, 0x98, 0x98
.db 0x97, 0x97, 0x97, 0x96, 0x96, 0x95, 0x95, 0x95, 0x94, 0x94, 0x94, 0x93, 0x93, 0x93, 0x92, 0x92
.db 0x92, 0x91, 0x91, 0x91, 0x90, 0x90, 0x90, 0x90, 0x8F, 0x8F, 0x8F, 0x8E, 0x8E, 0x8E, 0x8D, 0x8D
.db 0x8D, 0x8C, 0x8C, 0x8C, 0x8C, 0x8B, 0x8B, 0x8B, 0x8A, 0x8A, 0x8A, 0x89, 0x89, 0x89, 0x89, 0x88
.db 0x88, 0x88, 0x87, 0x87, 0x87, 0x87, 0x86, 0x86, 0x86, 0x86, 0x85, 0x85, 0x85, 0x84, 0x84, 0x84
.db 0x84, 0x83, 0x83, 0x83, 0x83, 0x82, 0x82, 0x82, 0x82, 0x81, 0x81, 0x81, 0x81, 0x80, 0x80, 0x80

The divide routine looks like this...

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;ARGUMENTS:  r16, r17, r18, r19
;  r17:r16 = N (numerator)
;  r19:r18 = D (divisor)
;RETURNS:    r20, r21
;  r21:r20 (quotient)
;  r23:r22 (remainder)
;
;DESCRIPTION:  divides an unsigned 16 bit number N by unsigned 16 bit divisor D
;  Max run time is 62 clock cycles.
;
;ALGORITHM OVERVIEW
;
;RZERO = 0;
;if(D < 256){
;  N2 = (N * ((R1H_TBL[D] << 8) + R1L_TBL[D])) >> 16;
;  P  = N2 * D
;}else{
;  D1 = (R1H_TBL[D] * D) >> 8
;  N1 = (R1H_TBL[D] * N) >> 8
;  if(D1 < 256){
;    N2 = N1 >> 8;
;  }else{
;    N2 = N2 * R2_TBL[D1 & 0xFF];
;  }
;  P = N2 * D;
;  if(P > 65535){
;    N2 = N2 - 1    ;//Decrement quotient
;    N1 = N2 - P + D;//Calculate remainder
;    return;//return quotient in N2, remainder in N1
;  }
;}
;N1 = N - P;
;if(P > N){
;  N2 = N2 - 1;//decrease quotient
;  N1 = N1 + D;//increase reamainder
;  return;//return quotient in N2, remainder in N1
;}
;if(N1 > D){
;  N2 = N2 + 1;
;  N1 = N1 - D;
;  return;//return quotient in N2, remainder in N1
;}
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
.def NL    = r16 .def NH    = r17 ;numerator
.def DL    = r18 .def DH    = r19 ;divisor
.def N2L   = r20 .def N2H   = r21 ;temp variables, becomes quotient.
.def N1L   = r22 .def N1H   = r23 ;temp variables, becomes remainder.
.def PRODL = r0  .def PRODH = r1  ;hardware multiply product
.def PL    = r2  .def PH    = r3  ;product
.def D1L   = r4  .def D1H   = r5
.def RZERO = r6                   ;zero value
.def Rx    = r7 

idivu_16x16:
    clr RZERO                 ;1
    ;Check that DH is not zero
    tst DH                    ;1
    brne idivu_16x16_dhne   ;2
    ;code for D < 256   
idivu_16x8:
    ;lookup low byte of recipricol into P.
    ;where P = min(2^16 / D,2^16-1)
    mov zl, DL               ;1
    #ifdef RAM_DIVIDE_TABLE
        ldi zh, high(R1L_TBL)    ;1 
        ld PL, Z                 ;2
        ldi zh, high(R1H_TBL)    ;1
        ld PH, Z                 ;2
    #else
        ldi zh, high(R1L_TBL*2)  ;1
        lpm PL, Z                ;3
        ldi zh, high(R1H_TBL*2)  ;1
        lpm PH, Z                ;3
    #endif  
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;calculate N2 = (P * N) >> 16
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;     NH:NL
    ;  X  RH:RL
    ;------------------------------------------
    ;   N2H    |   N2L    |  N1H     | dropped
    ;----------+----------+----------+---------
    ;          |          | H(PL*NL) | L(PL*NL)
    ;          | H(PL*NH) | L(PL*NH) |
    ;          | H(PH*NL) | L(PH*NL) |
    ; H(PH*NH) | L(PH*NH) |          |
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
    
    mul NL , PL     ;1  PL*NL
    mov N1H, PRODH  ;1  N1H <= H(PL*NL)
    mul NH , PH     ;1  PH*NH
    movw N2L, PRODL
    mul NH , PL     ;1  PL*NH
    add N1H, PRODL  ;1  N1H <= H(PL*NL) + L(PL*NH) 
    adc N2L, PRODH  ;1  N2L <= L(PH*NH) + H(PL*NH)
    adc N2H, RZERO  ;1  propagate carry to N2H      
    mul NL , PH     ;1  PH*NL
    add N1H, PRODL  ;1  N1H <= H(PL*NL) + L(PL*NH) + L(PH*NL)
    adc N2L, PRODH  ;1  N2L <= H(PH*NL) + L(PH*NH) + H(PL*NH)
    adc N2H, RZERO  ;1  propagate carry to N2H  
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;calculate P = N2 * DL ,note DH=0
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;    
    mul N2L, DL              ;1
    movw PL, PRODL           ;1
    mul N2H, DL              ;1
    add PH, PRODL            ;1
    rjmp idivu_16x16_adj_n ;2
    ;code for D >= 256
idivu_16x16_dhne:          
    ;Lookup Rx = min(256 / DH, 255)     
    mov zl, DH               ;1
    #ifdef RAM_DIVIDE_TABLE
        ldi zh, high(R1H_TBL)    ;1
        ld Rx, Z                 ;2
    #else
        ldi zh, high(R1H_TBL*2)  ;1
        lpm Rx, Z                ;3
    #endif
    ;D1 = (D * Rx) >> 8          
    mul Rx , DH              ;1
    movw D1L, PRODL          ;1
    mul Rx , DL              ;1
    add D1L, PRODH           ;1
    adc D1H, RZERO           ;1
    ;N1 = (D * Rx) >> 8          
    mul Rx , NH              ;1
    movw N1L, PRODL          ;1
    mul Rx , NL              ;1
    add N1L, PRODH           ;1
    adc N1H, RZERO           ;1
    ;if D1H = 0 then use Rx = 256, otherwise use table   
    tst D1H                  ;1
    brne idivu_16x16_dxhne ;2
    
        mov N2L, N1H             ;1
    clr N2H                  ;1
    rjmp idivu_16x16_checkn;2

    idivu_16x16_dxhne:
    ;Lookup Rx = (2 ^ 16) \ (256 + D1H)
    mov zl, D1L              ;1
    #ifdef RAM_DIVIDE_TABLE
        ldi zh, high(R2_TBL)     ;1
        ld Rx, Z                 ;2
    #else
        ldi zh, high(R2_TBL*2)   ;1
        lpm Rx, Z                ;3
    #endif
    ;N2 = (N1 * R2) >> 16
    mul Rx  , N1H            ;1
    mov PL  , PRODL          ;1
    mov N2L , PRODH          ;1
    mul Rx , N1L             ;1
    add PL , PRODH           ;1
    adc N2L, RZERO           ;1
    clr N2H                  ;1

    idivu_16x16_checkn:
    ;Check result (it may be off by +/- 1)
    ;P = N2 * D
    ;NOTE:  N2 <=255 so NxH = 0, also P < 2^16 so we can discard upper byte of DH * NxL
    mul DL , N2L             ;1
    movw PL, PRODL           ;1
    mul DH , N2L             ;1
    add PH , PRODL           ;1 
    brcc idivu_16x16_adj_n ;2

    ;if multiply overflowed then...
    ;decrement quotient
    ;calculate remainder as N - P + D
    subi N2L, 0x01           ;1
    sbci N2H, 0x00           ;1
    movw N1L, NL             ;1
    sub N1L, PL              ;1
    sbc N1H, PH              ;1
    add  N1L, DL             ;1
    adc  N1H, DH             ;1
    rjmp idivu_16x16_end   ;2

    ;Adjust result up or down by 1 if needed.
    idivu_16x16_adj_n:
    ;Add -P to N, with result in P
    movw N1L, NL             ;1
    sub N1L, PL              ;1
    sbc N1H, PH              ;1
    brsh idivu_16x16_pltn  ;2

    idivu_16x16_decn2:
    ;if P > N then decrement quotient, add to remainder
    subi N2L, 0x01           ;1
    sbci N2H, 0x00           ;1
    add  N1L, DL             ;1
    adc  N1H, DH             ;1
    rjmp idivu_16x16_end   ;2

    idivu_16x16_pltn:
    ;test remainder to D 
    cp  N1L, DL              ;1
    cpc N1H, DH              ;1
    ;if remainder < D then goto end
    brlo idivu_16x16_end   ;2

    ;if remainder >= D then increment quotient, reduce remainder
    subi N2L, 0xFF           ;1
    sbci N2H, 0xFF           ;1
    sub N1L, DL              ;1
    sbc N1H, DH              ;1
    idivu_16x16_end:
    ret
    .undef NH    .undef NL   
    .undef DH    .undef DL   
    .undef N2H   .undef N2L  
    .undef N1H   .undef N1L  
    .undef PRODL .undef PRODH
    .undef PH    .undef PL   
    .undef D1H   .undef D1L  
    .undef RZERO 
    .undef Rx  

One optimization that could have a big impact on the code run time would be to see if I can eliminate overflow error checking by rounding the table values up or down by one count. Right now some combinations of table values and inputs can result in overflow, and I need to check for that. I might try running a variant of this algorithm on a PC and find those cases and see if rounding the table value up or down by one count fixes the problem without breaking anything else.

In any case I am not aware of any other 16-bit x 16-bit divide algorithms that will run on this processor in less than 62 cycles, so we did good so far.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.