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

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.

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
;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
;N1 = (D * Rx) >> 8
mul Rx , NH              ;2
mov N1L, PRODL           ;1
mov N1H, PRODH           ;1
mul Rx , NL              ;2
;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
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
;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
rjmp idivu_16x16_end   ;2
;Adjust result up or down by 1 if needed.
;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
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

• What is the specific assembler you're using? Commented Dec 30, 2020 at 15:51
• Also, why is the division being performed on this device? Commented Dec 30, 2020 at 16:05
• Finally, can you give some domain information on the variables being divided? Often, this information can help optimize algebraic operations or even avoid them. Commented Dec 30, 2020 at 16:09
• @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. Commented Jan 1, 2021 at 3:21
• @Reinderien I am using AVR Assembler in Atmel Studio 7. Commented Jan 1, 2021 at 3:23

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

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

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

;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
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
brcs idivu_16x16_mofl   ;2   6  43

;Adjust result up or down by 1 if needed.
;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
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
ret                     ;1   8  53

• Putting off a stab at decent naming in hopes of fully understanding the algorithm - don't hold your breath. Commented Jan 17, 2021 at 21:29
• 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. Commented Jan 18, 2021 at 7:25
• 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. Commented Jan 18, 2021 at 16:38
• 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). Commented Jan 18, 2021 at 16:47
• 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). Commented Jan 18, 2021 at 16:48

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
;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
;N1 = (D * Rx) >> 8
mul Rx , NH              ;1
movw N1L, PRODL          ;1
mul Rx , NL              ;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
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

;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
rjmp idivu_16x16_end   ;2

;Adjust result up or down by 1 if needed.
;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
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.

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.

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
mul POW2NORM, NUM1
movw NUM1NORM, r0
//mov NUM1NORM, r0
//mov NUM2NORM, r1
mul POW2NORM, NUM0
#undef POW2NORM
//mov      NUM0NORM, r0

// 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
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
// 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
nofixup:
movw r22, Q0
ret

//overflow:
//sub  NUM0, r26
//sbc  NUM1, r27
decquot:
dec Q0
increm :
//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
mul DINV1, NUM0

mul Q0, DEN0
movw r26, r0
mul Q1, DEN0

// 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
subi Q0, 1
sbci Q1, 0
nofix:
movw r22, Q0
ret

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

• Assuming it works, I like this solution. Could you please include the lookup table values, and anything else you referenced but didn't define. Commented Dec 30, 2023 at 15:16
• @user4574 I added the tables and more ways to probably save memory. I should have thought of the second one earlier (I considered shifting, but not shifting large dens very small lookup for high den < 16). Note that lo8 and hi8 are from me using gnu as. Apart from that, I don't think I have other undefined symbols, except highbyte0 which is just the code for < 256, tweaked for the new reciprocal. I'll hopefully write that later. Commented Dec 30, 2023 at 23:50
• After a bit of testing, it seems like unconditionnally adding 1 to the low byte of the reciprocal used here makes the approximation of the quotient be correct up to +/- 1, which shouldn't cause problem with the current way it's done. The high byte can be modified in the table. Commented Dec 31, 2023 at 2:45
• @user4574 After trying out the divnby1, I managed 60 cycles worst case for the entire algorithm. This removes the high rec. table and half the low one, but as I said, worsens 8 bit dens by a lot. Commented Jan 1 at 6:17
• Actually what I said about unconditionnally adding 1 was wrong. You have to add 1 when the denominator is bigger than 1. That makes the result be as expected. Note that you can also remove the need to check for the approximation being less than the quotient (special case 1 at the beginning for 2 cycles, add 1 to the low reciprocal, then add 1 if the denom == 0xAC, for 3 cycles (or 2, if you want to load the reciprocal of 0xAC manually)), but I don't know if that's interesting for you. Commented Jan 2 at 3:43

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)

• Is this really an answer though?
– Mast
Commented Jan 14, 2021 at 16:00
• @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. Commented Jan 14, 2021 at 17:02