How to divide an unsigned 8-bit integer by 3 without divide or multiply instructions (or lookup tables)

Question

On CPUs without divide or multiply instructions like the Zilog Z80 and the MOS 6502, how could division by three be efficiently implemented?

I know in practice in game code you'd usually use a 256-byte lookup table, but I want to know the mathematical way because I'm curious. Also, the same technique would probably be useful in FPGA logic design, which I don't know much about, but where I think you wouldn't use a lookup table. (Well that intuition was wrong!)

When I look this up online I mostly find either

1. General division routines which iterate and allow any arbitrary denominator.
2. Summing of several divisions by powers of two, which are done using right-shift instructions: n >> 2 + n >> 4 + n >> 6

An unrolled version for a specific hard-coded denominator should presumably always be faster than a general method using a loop.

But my attempts at the "sum of shifts" method looks close for small numbers but not for larger numbers:

numerator
|           sum-of-shifts
|           |           true answer
↓           ↓           ↓

3           0           1
6           1           2
9           2           3
12          3           4
15          3           5
... snip ...
243         78          81
246         79          82
249         80          83
252         81          84
255         81          85

Have I got something wrong in my sum-of-shifts implementation? Or is the real best method something different? Or am I maybe missing some adjustment to make up for this tending to round down each step?

(If this belongs on StackOverflow (or math.SE) feel free to move it. My thinking was this used to be used in games and demos but nowadays everything has fast multiply and probably also fast divide, except maybe FPGA but FPGA is relevant for retro as well anyway.)

Answer 1

You are so close to a solution: the maximum error (remainder) of the shift-and-add method is only 9, so you could simply use a 10-element table to fix it:

int div3(int x){
    const int t[]={0,0,0,1,1,1,2,2,2,3};
    //int y = (x>>2) + (x>>4) + (x>>6); // max error is 12 if done one-line
    int y = (x>>2) + (x>>4);
    y+=y>>4;
    int b = x - (y<<1) - y;
    y+=(t[b]);
    return y;
}

ideone here

Update:

It turns out that 9 is small enough of an error that we can accurately calculate [0, 9]/3 without using LUT, and within the limit of 8-bit register. So here's a version that strictly doesn't use any LUT.

int div3(int x){
    int y = (x>>2) + (x>>4);
    y+=y>>4;        // slightly better than (x>>2) + (x>>4) + (x>>6)
    x=(x-y-(y<<1)); // x is remainder of x-3*y, maximum x is 9
    y+=((x<<3)+(x<<1)+x)>>5; // 11/32 is good enough approximation of 1/3 for [0, 9]
    return y;
}

ideone here

Update 2:

It turns out 11/32 is good up until 31, and there happens to be an easy way to narrow down the remainder to [0, 30].

By dividing x into x = 16*a + b where a is upper 4 bits and b is lower 4 bits, we have x = 3*5*a + (a+b), where a+b is between [0, 30]. So x/3 = 5*a + (a+b)*11/32.

To be honest, I am quite surprised this actually works!

int div3(int x) {
    int a = x>>4;
    int b = x&0x0f;
    int y = a + (a<<2);
    b = a+b;
    y+=((b<<2)+b+(b>>1))>>4; // use 5.5/16 to avoid overflow, result is the same
    return y;
}

ideone here

Update 3:

Please forgive me for keep adding to this answer, but I came across some interesting insights when commenting on this answer:

You can take advantage of the fact that 4=3+1, i.e. let x = a* 4+b = a*3+(a+b) then let (a+b) be the new x until x is less than 4.

This leads to a recursive solution like this (test here):

// recursion naturally ends when x<=3 but we can finish early at x<=5
// for any x between [0,5], x/3 could be computed by comparing with 2

int div3(int x){
    if(x>5){
        int a = x>>2;
        int b = x&0x03;
        return a + div3(a+b);
    } else {
        return x>2;
    }
}

The recursion itself isn't particularly interesting, except for its simplicity. But what if we fully unroll this recursion?

For an integer x we rewrite it as base-4, i.e. x = a*4^3 + b*4^2 + c*4 + d.

Then by unrolling the recursion once, we have x = (a*4^2*3 + b*4*3 + c*3) + (a*4*2 + b*4 + c) + d = (x>>2)*3 + (x>>2) + d .

Keep unrolling the recursion, eventually, we have x = (x>>2)*3 + (x>>4)*3 + (x>>6)*3 + (a+b+c+d).

Equivalently x = 3*((a>>2)+(x>>4)+(x>>6)) + (a+b+c+d).

Here we are back to the origin: this is how the shift-add-add approximation works. It also predicts the maximum error to be 12 (because the maximum of a+b+c+d is 12 when a b c d all equal to 3). Since this equation is exact, the remainder x%3 is exactly (a+b+c+d)%3 as well.

In general, we can easily expand this to all x/(2^n-1) and x%(2^n-1) problems of all x and n.

Answer 2

The question isn't particularly retrocomputing related, but can come up a bit more in that context so I'll answer it here rather than punt to another StackExchange site.

Division is equivalent to multiplication by the reciprocal. Computing the reciprocal does of course itself involve division, and chaining the operations loses precision, so this is less than useful for general-purpose division. However, it does win for division by a compile-time constant. Modern compilers do this all the time, for example LLVM produces this for a trivial function which does an 8 bit division by 3 and returns the result:

        uxtb    r0, r0
        movs    r1, #171
        muls    r0, r1, r0
        lsrs    r0, r0, #9
        bx      lr

171 is ceil(512/3). What you see in that code is zero-extension to 32 bits, a 32 bit multiplication by the precomputed 512/3, and then a division by 512 using a shift. Although it doesn't apply in this case, some divisors require a fixup step, and this is well outside retrocomputing and very much into heavy computer science territory, so I'll just point you at e.g. http://homepage.cs.uiowa.edu/~jones/bcd/divide.html for the fine detail.

On a retrocomputer, you would pick a different code sequence. If you do not require a precise answer, you can skip the fixup and also multiply by 86 (ceil(256/3)) so that it requires a division by 256 instead, which does not involve a shift at all as you'd compute a result in e.g. BC and document the function as returning its result in B.

What's left is multiplication by a constant. You can do this with quarter-squares (except that uses a lookup table, albeit one which is readily reused elsewhere in your program and is effectively free), shift-and-add in a loop as a general multiplication function—this is ultimately much the same thing as your sum-of-shifts mentioned in the question—or unrolled shift-and-add in a bespoke divide-by-three function. Unrolled would be fastest, but since 86 is 0b01010110 (which has 4 bits set) then that would require four iterations of substantially-identical code and it's debatable whether this is a good use of precious memory.

Answer 3

I am not really familiar with the Z80 or MOS6502, but generally on any 8-bit CPU you can use shifting and subtracting to resolve each result bit.

For starters, maybe best would be to avoid cases where you have to divide by 3 to begin with so you don't have to do it.

For this specific case of dividing by 3, the largest result fits into 7 bits, and the largest value of 3 shifted left so that fits in 8 bits is 192, or 0xC0.

The algorithm¹ for dividing X by 3 is:

0. If 0xC0 can be subtracted from X then do it and add 0x40 to result.
1. If 0x60 can be subtracted from X then do it and add 0x20 to result.
2. If 0x30 can be subtracted from X then do it and add 0x10 to result.
3. If 0x18 can be subtracted from X then do it and add 0x08 to result.
4. If 0x0C can be subtracted from X then do it and add 0x04 to result.
5. If 0x06 can be subtracted from X then do it and add 0x02 to result.
6. If 0x03 can be subtracted from X then do it and add 0x01 to result.

This can obviously be a loop, as each value is just shifted by one bit right per loop.

And can easily be extended into a generic function, or kept inline with fixed values for maximum performance.

There could be other ways and variations, and how to do it most efficiently depends on available CPU opcodes, but this is a general idea.

Another idea is to multiply with suitable number that represents one third and then take the upper bits. This may work as Z80 has 16-bit register pairs. I don't have exact values at hand, but if you multiply using bitshifts by 85 and then take upper 8 bit register, you have effectively multiplied with 85/256 which is 0.332. As a correction, the value needs to be corrected by adding a suitable number (85) before taking the high byte. This works well with 8-bit MCUs that have multiply instruction.

¹ In Python (with godbolt):

def divide_by_3(x):
    def apply(n, x, r):
        if x >= 3 * n:
            return x - 3 * n, r + n
        else:
            return x, r

    x, r = apply(0x40, x, 0)
    x, r = apply(0x20, x, r)
    x, r = apply(0x10, x, r)
    x, r = apply(0x08, x, r)
    x, r = apply(0x04, x, r)
    x, r = apply(0x02, x, r)
    _, r = apply(0x01, x, r)

    return r

for i in range(256):
    n = divide_by_3(i)

    if n != i // 3:
        print(f"Failed to divide {i} by 3, got {n} expected {i//3}")

Answer 4

In comments, WimC pointed to working 6502 code that divides an 8-bit unsigned integer by 3, with quotient rounded towards zero as required by C and C++. This code computes x / 3 as (4/3 * x) / 4, where the factor 4/3 is approximated as 1 + 1/4 + 1/16 + 1/64, which is an underestimate:

res = (arg + (arg / 4) + (arg / 16) + (arg / 64)) / 4

Clearly this computation requires 9 bits throughout, which is awkward with 8-bit processors. But we can re-arrange the computation as

res = (((((arg) / 4) + arg) / 4 + arg) / 4 + arg) / 4

and further expand this into

res = ((((((((arg) / 2) / 2 + arg) / 2) / 2 + arg) / 2) / 2 + arg) / 2) / 2

The advantage of this is that many 8-bit processors allowed for the efficient computation of (a + b) / 2, in that the addition would result in a 9-bit result, the lower 8 bit of which were stored in the accumulator while the most significant bit was stored in the carry flag. A subsequent rotate-right-through-carry instruction could then perform a 9-bit right rotation, effecting the division of the 9-bit sum by two. In the case of the 6502, the instruction for this is ROR.

So far we have glossed over the fact that our chosen approximation to 4/3 is an underestimate and that integer division truncates. This suggests that our computed quotient will frequently be too low, and a quick test verifies this. We can try to counteract this by adding a correction constant. I do not know of any particular algorithm to derive such corrections, I usually search for these by trial and error. In this case, we get the correct quotient for all 256 possible dividends by computing:

res = ((((((((arg + 85) / 2) / 2 + arg) / 2) / 2 + arg) / 2) / 2 + arg) / 2) / 2

One immediately notices that the linked 6502 code uses a different correction constant. This is due to two factors: (1) The addition instruction used is not a plain addition but an ADC (add with carry) which incorporates the carry set by a preceding LSR (logical shift right). This tends to make the result of the addition larger than in the pseudocode above. (2) To minimize cycle count, the code adds the correction after the first division by two. Instead of starting off with CLC / ADC #43 / ROR / LSR, it only needs three instructions: LSR / ADC #21 / LSR, with the value of the carry flag well defined after the initial LSR.

Best I recall the Z80 instruction set, the algorithm presented for the 6502 could be translated one to one for the Z80: LSR becomes SRL, ROR turns into RRA, and ADC remains ADC.

Below are tested 8086 implementations of the variants mentioned above. On an 8086 one would want to use AAM 3 to divide a byte-size operand in AL by 3, so this serves merely as a reference that may be more familiar to a modern audience.

__asm mov al, byte ptr [arg];
__asm mov cl, al;
__asm shr al, 1;
__asm adc al, 21;
__asm shr al, 1;
__asm adc al, cl;
__asm rcr al, 1;
__asm shr al, 1;
__asm adc al, cl;
__asm rcr al, 1;
__asm shr al, 1;
__asm adc al, cl;
__asm rcr al, 1;
__asm shr al, 1;
__asm mov byte ptr[res], al;

__asm mov al, byte ptr [arg];
__asm mov cl, al;
__asm clc;
__asm adc al, 43;
__asm rcr al, 1;
__asm shr al, 1;
__asm adc al, cl;
__asm rcr al, 1;
__asm shr al, 1;
__asm adc al, cl;
__asm rcr al, 1;
__asm shr al, 1;
__asm adc al, cl;
__asm rcr al, 1;
__asm shr al, 1;
__asm mov byte ptr[res], al;

__asm mov al, byte ptr [arg];
__asm mov cl, al;
__asm add al, 85;
__asm rcr al, 1;
__asm shr al, 1;
__asm add al, cl;
__asm rcr al, 1;
__asm shr al, 1;
__asm add al, cl;
__asm rcr al, 1;
__asm shr al, 1;
__asm add al, cl;
__asm rcr al, 1;
__asm shr al, 1;
__asm mov byte ptr[res], al;

Answer 5

Well one of the mysteries of life is that often the answer pops into your head minutes or seconds after posting the question online.

(I puzzled over this for four or fives days, along with another old Z80 coder I met online.)

As I alluded to, it was the loss of precision and there were two ways to fix it.

1. Do the adding before part of the shifting.
2. Add 1. I remember having to do this a lot in integer stuff to do with rounding. I have a bit of a feeling of why but can't put it into words.

Here is JavaScript code:

let foo = n => 1 + ((n + (n>>2) + (n>>4) + (n>>6))>>2)

And here is some JS to print out a table with just the 85 8-bit values that are exact multiples of 3:

[...Array(85).keys()].map(n => (n+1)*3).map(n => n + " " + foo(n) + " " + n/3)

One thing I'm not sure about though. JavaScript integers are not 8-bit so there could be an overflow doing the adds before the >> 2. I'll have a think about that.

(Actually I just noticed that my code is wrong when the numerator is 0. It returns 1 instead of 0. I tried adding ½ instead of 1 by doing >> 1 then + 1 then >> 1 but then the answers were off for n ≥ 231)

Answer 6

If one feeds the bits of an arbitrary-length number into a finite state transducer in left-to-right order, it's possible to output that number divided by three. In each state indicated by an uppercase letter, the machine may receive a 0 or 1, output the indicated bits.

    0    1
A - 0A   0B
B - 0C   1A
C - 1B   1C

I'm not sure how to best realize such a state machine in code, but perhaps having an unrolled loop for state A, an unrolled loop which uses state B when an even number of bits remain and C when an odd number of bits remain, and an unrolled loop which uses state C when an even number remain and B when an odd number remain, might be a good approach depending upon what shift and rotate operations are available on the target platform.

Answer 7

Here's my two penneth...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

int main(void)
{
    unsigned char x;
    unsigned char xover3;
    unsigned char xover3approx;
    int i, errorcount = 0;
    int diff, maxdiff = 0, maxdiffcount = 0;
    for (i=0; i<256; i++) {
      x = i;
      xover3 = x/3;
      xover3approx = 1 + (x >> 2) + (x >> 4) + (x >> 6) + (x >> 7);
      diff = abs(xover3 - xover3approx);
      if (diff != 0) {
        if (diff > maxdiff) {
          maxdiff = diff;
          maxdiffcount = 0;
        }
        if (diff == maxdiff) maxdiffcount++;
        if (diff > 1) printf("%i %i %i \n",i,  (int)xover3, (int)xover3approx);
        errorcount++;
      }
    }
    printf("No. of errors = %i\n", errorcount);
    printf("Max. error = %i (for %i values)\n", maxdiff, maxdiffcount);
}

produces...

63 21 19 
111 37 35 
123 41 39 
126 42 40 
127 42 40 
128 42 44 
255 85 83 
No. of errors = 126
Max. error = 2 (for 7 values)

So 1 + (x >> 2) + (x >> 4) + (x >> 6) + (x >> 7) is correct more than half the time, and in error by more than 1 for only 7 values.

Answer 8

You can rephrase it by saying how to multiply by the reciprocal of 3 which is binary 0.01010101(..)

You can calculate a*01010110 (rounding up to handle exact multiples) and then return the high byte.

Z80 code to divide a by 3:

      ld b, 0     ; 
      ld c, $56   ; initialize with reciprocal = ceil((1/3)*256)
      ld hl, 0    ; accumulate in BC
loop: rra         ; divide A by 2, putting low bit in carry
      jr nc, skip ; if low bit was 0, don't accumulate
      add hl, bc  ; low bit was 1, accumulate the reciprocal
skip: sla c       ; multiply reciprocal by 2
      rl b        ; "
      or a        ; any bits left?
      jr nz, loop ; yes? - then do next bit
      ld a, h     ; copy high byte of result to A     

Answer 9

Consider the following x86 code (used as inline-assembly in Delphi code):

function DIVByteBy3Fast(B:Byte):WORD;
asm
                    xor     AH, AH          // Convert BYTE in AL -> WORD in AX.
                    mov     CX, AX          // copy of B -> CX

                    // call    MulWordBy10  // AX := 10  * B
                    mov     DX, AX          // MulWordBy10 - inlined here.
                    add     AX, AX
                    add     AX, AX
                    add     AX, DX
                    add     AX, AX

                    mov     DX, AX          // DX := 10  * B
                    //call    MulWordBy16   // AX := 160 * B
                    add     AX, AX          // MulWordBy16 - inlined here. 
                    add     AX, AX
                    add     AX, AX
                    add     AX, AX

                    ADD     AX, DX          // AX := 170 * B
                    ADD     AX, CX          // AX := 171 * B

                    SHR     AX, 9           // AX := (171/512) * B = B DIV 3
end;

Basically, it divides an unsigned BYTE by 3 by using "multiply by reciprocal" -principle.

Especially on the Z80 CPU you could use the same principle, since although a Z80 is a 8-bit CPU, it still supports 16-bit additions:

ADD HL, DE ADD HL, BC

Answer 10

Multiplying by the reciprocal then shifting right works fine. 256*2/3 ≈ 171 = 10101011b. That takes one shift and one addition per each 1 in the multiplier.

A quick test in Octave to show that this produces correct results:

>> any(floor([0:255]*171/512) - floor([0:255]/3))
ans = 0

We can save one addition, though:

1010'1011b * x = 
              = 1010'1100b * x - x * 1b
              = 1011'0000b * x - 100b * x - 1b * x
              = 1100'0000b * x - 1'0000b * x - 100b * x - 1b * x
              = 1'0000'0000b * x - 100'0000b * x - 1'0000b - 100b * x - 1b * x

So, we've turned 5 additions into 4 subtractions. The first and last multiplications are free, since they are "shifts" by either 0 or 1 full bytes and on 8-bit systems we got access to byte registers directly.

An implementation in C looks as follows:

#include <assert.h>
#include <stdint.h>
#include <stdio.h>

uint8_t div3(uint16_t dividend) {
    uint16_t result = dividend << 8;
    result -= dividend;
    dividend <<= 2;
    result -= dividend;
    dividend <<= 2;
    result -= dividend;
    dividend <<= 2;
    result -= dividend;
    return result >> 9;
}

int main() {
    for (int i = 0; i < 256; i++) {
        assert(div3(i) == i/3);
        printf("%3d%4d\n", div3(i), i/3);
    }
}

The Z80 version could look like this:

; 8-bit integer divide by 3 function.
; Inputs:    H=dividend
; Outputs:   H=floor(quotient)
; Clobbered: B,L,D,E

INTDIV3:
     LD   L, 0         ; HL = dividend << 8
     LD   D, L
     LD   E, H         ; DE = 16-bit dividend
     AND  A            ; C=0 (CCF idiom)
     SBC  HL, DE       ; result(HL) -= dividend(DE)
     LD   B, 3         ; do below 3 times
LOOP:
     SHL  E
     RL   D
     SHL  E
     RL   D            ; dividend(DE) <<= 2, C=0
     SBC  HL, DE       ; result(HL) -= dividend(DE)
     DJNZ LOOP
;
     SRL  H            ; result(HL) >>= 9
     RET

The loop could be unrolled, saving a few cycles and preserving the B register's value.

I've also "emulated" the above in C. You can run it on godbolt.

#include <assert.h>
#include <stdint.h>

union {
  struct { uint8_t C[1], B[1], E[1], D[1], L[1], H[1], A[1]; _Bool CF; };
  struct { uint16_t BC[1], DE[1], HL[1]; };
} r;

void LD(uint8_t *d, uint8_t s) { *d = s; }
void AND(uint8_t s) { *r.A &= s, r.CF = 0; }
void SBC(uint16_t *d, uint16_t s) { *d -= s + r.CF; }
void SHL(uint8_t *d) { r.CF = *d & 0x80, *d <<= 1; }
void RL(uint8_t *d) {
    uint8_t CF = *d & 0x80;
    *d = (*d << 1) | r.CF, r.CF = CF;
}
void SRL(uint8_t *d) { r.CF = *d & 1; *d >>= 1; }
_Bool DJNZ() { return (--*r.B); }

uint8_t div3(uint8_t dividend)
{
    *r.H = dividend;
    LD(r.L, 0);             // HL = dividend << 8
    LD(r.D, *r.L);
    LD(r.E, *r.H);          // DE = 16-bit dividend
    AND(*r.A);              // C=0 (CCF idiom)
    SBC(r.HL, *r.DE);       // result(HL) -= dividend(DE)
    LD(r.B, 3);
    do {
        SHL(r.E);
        RL(r.D);
        SHL(r.E);
        RL(r.D);            // dividend(DE) <<= 2, C=0
        SBC(r.HL, *r.DE);   // result(HL) -= dividend(DE)
    } while (DJNZ());
    SRL(r.H);               // result(HL) >>= 9
    return *r.H;
}

int main() {
    for (int i = 0; i < 256; i++)
        assert(div3(i) == i/3);
}

Answer 11

A silly way to divide by 3 is to change base to 3 and then shift right. Example:

127(10) = 11201(3)
1120(3) = 0 + 2*3 + 1*9 + 1 * 27  = 6 + 9 + 27 = 15 + 27 = 42

No LUT needed. This is much like user3528438's answer, only he switches to base 16. Try it online.

// x=a*4+b=a*3+(a+b), all credit is due to user3528438
int div3(int const x)
{
  if (x > 3)
  {
    int a = x >> 2;

    return a + div3(x - (a << 1) - a);
  }
  else if (x == 3)
  {
    return 1;
  }
  else
  {
    return 0;
  }
}

a non-recursive solution:

// x=a*4+b=a*3+(a+b)
int div3(int x)
{
  int q = 0;

  while (x > 3) // x > 5
  {
    int a = x >> 2;

    q += a;
    //x -= (a << 1) + a; // reduce x
    x = a + (x & 0x3); // x = a + b
  }

  // return q + (x > 2);    
  return q + (3 == x);
}

--

// x=2a+b=3a+(b-a)
int div3(int x)
{
  int q = 0;

  while ((x < 0) || (x > 5))
  {
    int a = x >> 1;

    q += a;
    x = (x & 1) - a; // possible approximation: x = -a;
  }

  return q + (x > 2);
}

this is usually approximated as:

int div3(int x)
{
  int q = 0;

  do
  {
    q += x >>= 1;
  } while (x = -x);

  return q;
}

--

// fixed-point iteration:
// x_(n+1) = a - x_n/2
// you can also try:
// x_(n+1) = (a - x_n) / 2
// can be improved, if you're math-savvy
int div3(int x)
{
  int a = x;

  x >>= 1;

  x = a - (x>>1);
  x = a - (x>>1);
  x = a - (x>>1);
  x = a - (x>>1);
  x = a - (x>>1);
  x = a - (x>>1);
  x = a - (x>>1);
  x = a - (x>>1);

  return x >> 1;
}

--

uint8_t div3(uint8_t x)
{
  uint8_t a = x;
  x >>= 1;

  x = (a + (x >> 1)) >> 1;
  x = (a + (x >> 1)) >> 1;
  x = (a + (x >> 1)) >> 1;
  x = (a + (x >> 1) + 1) >> 1;

  return x >> 1;
}

-- this one just approximates

uint8_t div3(uint8_t x)
{
  uint8_t a = x - (x >> 2);

  x = a - (x >> 3);
  x = a - (x >> 3);
  x = a - (x >> 3);

  return x >> 1;
}