Fastest (in clock cycles) 16-bit x 16-bit unsigned integer division algorithm for ATMEGA1284?

Question

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

Answer 1

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

Answer 2

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.

Answer 3

I don't know if this is going to be useful to OP, but it may be useful to others.

At first, I thought actually normalising the numerator/denominator would be a terrible idea. But it turned out to be not so bad, quite good actually.

This algorithm is based on what libgcc does for double word by double word division, read the comments for slightly more info.

I count 61 cycles if using a reduced power of 2 table, 50 if aligning and using a full power table. Also a free space reduction is to use 128 bytes of power table, and jumping if the denominator is already normalised (costs 2 cycles for a skip). This is an improvement even when applying a straightening of the slowest path in the algorithm + removing unnecessary quotient operations (I count 54).

Something that could save more cycles is to see if the correction of the quotient can be improved. I also didn't bother re-doing the high byte 0 case, but there may be a way to make that be not too bad using the same algorithm, but that may make that the worse case.

Tables content: POW2 contains 1 << clz(x) (which can be obviously shortened), R1L_TBL now contains (2**16 - 1) / x - (x == 0xAC) for x >= 128. Should contain that for all x's. Probably needs adjustments in the high byte 0 case.

Edit: I missed one obvious thing: for big denominators (>= 2**14) you could simply roll with a small approximation and/or branches. The test for it would cost 2 cycles. Thus the power of 2 table can be shorted to 64, with the same cycle penalty. Other possibilities are the usage of actual shifting ( with all constant shift tricks) for only the denominator (and probs only for large ones) (probably not the numerator) <- this may save more memory (if den1 < 16 -> do lookup, otherwise shift/do other stuff described before).

Also saved a cycle from straightening the last check.

Tables:

.align 256
R1L_TBL:
.byte 255, 255, 255, 85, 255, 51, 170, 146, 255, 113, 153, 69, 85, 177, 73, 17, 255, 15, 56, 121, 204, 48, 162, 33, 170, 61, 216, 123, 36, 211, 136, 66, 255, 193, 135, 80, 28, 235, 188, 144, 102, 62, 24, 244, 209, 176, 144, 114, 85, 57, 30, 5, 236, 212, 189, 167, 146, 125, 105, 86, 68, 50, 33, 16, 255, 240, 224, 210, 195, 181, 168, 155, 142, 129, 117, 105, 94, 83, 72, 61, 51, 41, 31, 21, 12, 3, 250, 241, 232, 224, 216, 208, 200, 192, 185, 177, 170, 163, 156, 149, 143, 136, 130, 124, 118, 112, 106, 100, 94, 89, 83, 78, 73, 67, 62, 57, 52, 48, 43, 38, 34, 29, 25, 20, 16, 12, 8, 4, 255, 252, 248, 244, 240, 236, 233, 229, 225, 222, 218, 215, 212, 208, 205, 202, 199, 195, 192, 189, 186, 183, 180, 178, 175, 172, 169, 166, 164, 161, 158, 156, 153, 151, 148, 146, 143, 141, 138, 136, 134, 131, 129, 127, 124, 122, 120, 118, 116, 114, 112, 110, 108, 106, 104, 102, 100, 98, 96, 94, 92, 90, 88, 87, 85, 83, 81, 80, 78, 76, 74, 73, 71, 70, 68, 66, 65, 63, 62, 60, 59, 57, 56, 54, 53, 51, 50, 48, 47, 46, 44, 43, 41, 40, 39, 37, 36, 35, 33, 32, 31, 30, 28, 27, 26, 25, 24, 22, 21, 20, 19, 18, 17, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1
R1H_TBL:
.byte 255, 255, 128, 85, 64, 51, 42, 36, 32, 28, 25, 23, 21, 19, 18, 17, 16, 15, 14, 13, 12, 12, 11, 11, 10, 10, 9, 9, 9, 8, 8, 8, 8, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
POW2:
.byte 0, 128, 64, 64, 32, 32, 32, 32, 16, 16, 16, 16, 16, 16, 16, 16, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1

Code:

.macro SIZE name
.size \name, .-\name
.endm

.macro GFUNC name
.global \name
.func \name         
.type \name , %function
\name:
.endm

.macro ENDF name
SIZE \name
.endfunc
.endm

#define DEN0NORM r18
#define DEN1NORM r19
#define DINV  r20
#define POW2NORM r20
#define NUM0NORM  r21
#define DEN0  r22
#define DEN1  r23
#define NUM0  r24
#define NUM1  r25
#define REM0  r24
#define REM1  r25
#define NUM1NORM  r26
#define NUM2NORM  r27
#define Q0 r30
#define Q1 r31
#define DINV0 r26
#define DINV1 r27
#define TEMP r20
#define ZERO DEN1

// This function follows the double word by double word
// division algorithm present at: 
// https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libgcc/libgcc2.c
// in the function udivmoddi4
// it ignores some special cases
GFUNC div16by16
cpi DEN1, 0
breq highbyte0
noshift8:
//#define SPACE_CONTIOUS_NORM
#ifdef SPACE_CONTIOUS_NORM
mov r30, DEN0;
cpse DEN1, r1;
mov    r30, DEN1
noshift8:
clr    r31
cpi    r30, 16
brsh   hibits
subi   r30, lo8(-(POW2))
sbci   r31, hi8(-(POW2))
lpm    POW2NORM, Z
swap   POW2NORM
rjmp   normalize
hibits :
swap   r30
andi   r30, 15
subi   r30, lo8(-(POW2))
sbci   r31, hi8(-(POW2))
lpm    POW2NORM, Z
#else
mov    r30, DEN1
noshift8:
ldi    r31, hi8(POW2)
lpm    POW2NORM, Z
#endif
normalize:
mul POW2NORM, DEN1
mov DEN1NORM, r0
mul POW2NORM, DEN0
//mov DEN0NORM, r0
add DEN1NORM, r1
mul POW2NORM, NUM1
movw NUM1NORM, r0
//mov NUM1NORM, r0
//mov NUM2NORM, r1
mul POW2NORM, NUM0
#undef POW2NORM
//mov      NUM0NORM, r0
add NUM1NORM, r1

// this reciprocal is the low byte of ((1 << 16) - 1) / d with d's highest bit set
// so it can (must, if not wanting to clr the upper bit) be reused if the denominator has a zero high byte
// this is not the same reciprocal as OP's for powers of 2
mov    r30, DEN1NORM
ldi    r31, hi8(R1L_TBL)
lpm    DINV, Z

// in libgcc, this is the beginning of the line udiv_qrnnd (q0, n1, n2, n1, d1)
// implemented with the beginning of DIV2BY1 in the gmp paper:
// https://gmplib.org/~tege/division-paper.pdf
mul DINV, NUM2NORM
mov Q0, r1
add r0, NUM1NORM
adc Q0, NUM2NORM
clr Q1
// at this point, the q0 (as q1) result can be used for 8 bit dens

// find a faster way to correct the quotient?
// this seems to save approx 5 cycles over continuing the truncating div
mul Q0, DEN0
movw r26, r0
mul Q0, DEN1
add r27, r0
// here r1 is 0, no need to clear it
// for no overflow: sub 1 at 0xAC of the table
// brcs overflow
// here we simply substract from the remainder repeatedly
// if we ever go negative, we substracted one too much
fixquotglobal:
sub NUM0, r26
sbc NUM1, r27
brcs decquot
sub REM0, DEN0
sbc REM1, DEN1
brcs increm
inc Q0
sub REM0, DEN0
sbc REM1, DEN1
brcc incquot
add REM0, DEN0
adc REM1, DEN1
nofixup:
movw r22, Q0
ret 

//overflow:
//sub  NUM0, r26
//sbc  NUM1, r27
decquot:
dec Q0
increm :
add REM0, DEN0
adc REM1, DEN1
//rjmp final
movw r22, Q0
ret 

incquot:
inc Q0
movw r22, Q0
ret 

highbyte0:
mov r30, DEN0
ldi r31, hi8(R1L_TBL)
lpm DINV0, Z
ldi r31, hi8(R1H_TBL)
lpm DINV1, Z

// this is to make sure the approximation is >= quotient
cpi  DEN0, 2
brlo handle1den
inc  DINV0
not1:
cpi  DEN0, 0xAC
brne mulrec
inc  DINV0
mulrec:
mul DINV1, NUM1
movw Q0, r0
mul DINV0, NUM0
mov TEMP, r1
mul DINV0, NUM1
add TEMP, r0
adc Q0, r1
adc Q1, ZERO
mul DINV1, NUM0
add TEMP, r0
adc Q0, r1
adc Q1, ZERO

mul Q0, DEN0
movw r26, r0
mul Q1, DEN0
add r27, r0
            
// note: here this can be handled by
// the same code as 16 byte the same as original
// that costs no cycles for the worst case on the 16 / 16 path
// but it costs a jump here, which makes 8 bit the worst case
sub NUM0, r26
sbc NUM1, r27
brcc nofix
add NUM0, DEN0
adc NUM1, DEN1
subi Q0, 1
sbci Q1, 0
nofix:
movw r22, Q0
ret

handle1den:
movw r22, NUM0
clr NUM0
clr NUM1
ret 
ENDF div16by16

Answer 4

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)