31

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.)

14
  • 7
    On an 8-bit cpu you do divisions either efficiently or without lookup tables.
    – lvd
    Commented Nov 9, 2022 at 11:07
  • 4
    Also sum of shifts must be infinite: x/3 = x/(4*3/4) = (x/4)(1/(1-1/4))=(x/4)(1+1/4+1/16+1/64+...)=x(1/4+1/16+1/64+1/256+...). Which amount of shifts is enough and whether the rounding is needed depends on the range and bitness of your divident.
    – lvd
    Commented Nov 9, 2022 at 11:11
  • 10
    Sidenote: "FPGA programming, [...] where I think you wouldn't use a lookup table." FPGA is all about lookup - all logic is made from lookup tables. SCNR. Then again and more to the topic, using a 8 (9) bit lookup table would as well be a common solution as they contain lots of RAM. Preloading one such block as ROM is pretty standard.
    – Raffzahn
    Commented Nov 9, 2022 at 13:31
  • 2
    On the division part, implementation is always a trade off in resources (registers, memory, size and speed). A lookup will need 256 bytes plus maybe 3-10 bytes code, zero memory, 1..2 registers and a few cycles execution. In contrast a loop, like @justme's example will go with zero table but 50-100 bytes of code and 2..4 registers, up to 5 bytes of memory and several hundred cycles. So the table is more often than not the right solution as well - unless one is extreme memory strapped, or willing to accept the slight distortion of the shift method shown :)) (same for what lvd suggests)
    – Raffzahn
    Commented Nov 9, 2022 at 13:40
  • 8
    There’s a 6502 “divide by three” here. Did not try it, but it seems not obvious on first glance.
    – WimC
    Commented Nov 9, 2022 at 21:30

11 Answers 11

28

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.

2
  • I really like your update #2 version! And it still works if you change int to unsigned char. Commented Nov 10, 2022 at 16:37
  • 1
    what you actually did is: you changed base to 16 and then did your magic. Commented Dec 11, 2022 at 10:41
15

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.

4
  • 1
    You really couldn't assume the CPU here to have a multiplier or a 32 bit register. Register size is the bigger difficulty here: if the maximum input is n bits, then to get a bit-exact solution, the reciprocal must be more precise than 1/(2^n), hence the reciprocal must be n+1 bits(as is the case in your LLVM example). This method isn't even feasible on a 16-bit machine. Commented Nov 9, 2022 at 12:41
  • 4
    @user3528438 The question was very specifically to divide an unsigned 8-bit number by 3; so this method — multiply by 171, shift right by 9 — works with 16-bit arithmetic. It's just ((x << 7) + (x << 5) + (x << 3) + (x << 1) + x) >> 9.
    – Tommy
    Commented Nov 9, 2022 at 14:37
  • 3
    @Tommy or the equivalent ((x << 6) + (x << 4) + (x << 2) + x + (x >> 1)) >> 8 which brings back this answer’s idea of just using the upper byte after the 16 bit operation (but without the errors you get when just multiplying with 86).
    – Holger
    Commented Nov 9, 2022 at 17:03
  • 3
    If the result is only required to be correct for multiples of 3 (as in the test cases in the question), you can just multiply by 171 in an 8-bit register. If you don't have a multiplication instruction, this can be done in 5 instructions by t = x + (x << 4); return -((t << 2) + t). Commented Nov 10, 2022 at 16:00
13

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 algorithm1 for dividing X by 3 is:

  1. If 0xC0 can be subtracted from X then do it and add 0x40 to result.
  2. If 0x60 can be subtracted from X then do it and add 0x20 to result.
  3. If 0x30 can be subtracted from X then do it and add 0x10 to result.
  4. If 0x18 can be subtracted from X then do it and add 0x08 to result.
  5. If 0x0C can be subtracted from X then do it and add 0x04 to result.
  6. If 0x06 can be subtracted from X then do it and add 0x02 to result.
  7. 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.

1 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}")
11
  • 14
    @OmarL It's the same algorithm you do division of any two numbers with pen and paper. Try it.
    – Justme
    Commented Nov 9, 2022 at 8:58
  • 6
    @OmarL, when we did it in base 10 in elementary school, we called it long division. It's simpler in base 2 because there's only two possible choices for each next digit of the quotient. Commented Nov 9, 2022 at 15:47
  • 4
    This is basically binary search. (All digit-by-digit algorithm degrades to binary search if the digits are binary). A while ago I wrote an algorithm to calculate integer square root bit-by-bit, and only to find out I have just implemented binary search. Commented Nov 9, 2022 at 16:48
  • 4
    I was initially confused by the algorithm you presented, because absent indication to the contrary I assumed that X was immutable and all steps could be executed in parallel, which clearly is not the case. Maybe "then do it" should be reworded as "then assign X - 0XC0 to X" to make it clear X is mutated at each step. Commented Nov 10, 2022 at 15:41
  • 3
    @stanri Of course it can be done like that, but question asked for efficient ways. You can use a simple algorithm which takes long to run for large numbers, so loop will run up to 85 times for input of 255. The long division can be done in a loop too, and it only needs to run the loop 7 times.
    – Justme
    Commented Nov 11, 2022 at 17:01
8

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;
2
  • 85 is very interesting! I noticed mid-way through my experiments before posting here that for some numbers doing a shift right one bit on each pair of bits worked. That was equivalent to shifting the whole byte one bit to the right followed by an AND with 0x55 which is 01010101. But it didn't work for numbers with bits carried between the pairs so I gave up in, not being clever enough (-: Commented Nov 10, 2022 at 14:46
  • I think it was well after midnight when I tried to emulate the 6502 code in JS. As an old Z80 guy I was used to the mnemonics indicating which opcodes made use of the carry. But even when I realized some of the 6502 ones involved carry without a mnemonic hint and altered my JS it still didn't work. So I knew it was time to let my brain sleep (-: Commented Nov 10, 2022 at 14:49
6

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)

6
  • So e.g. 1/3 = 1?
    – Tommy
    Commented Nov 9, 2022 at 11:56
  • 2
    @Tommy no, 0/3=1 too. Commented Nov 9, 2022 at 12:44
  • 3
    "Well one of the mysteries of life is that often the answer pops into your head minutes or seconds after posting the question online." - My old coworkers and I used to call this the 'sounding board' principle. Frequently, the simple act of explaining a problem to someone else will help you to grok the answer your looking for - but the person can be replaced by a 'sounding-board' just as effectively. It works. Just explain your dilemma to an empty room and amaze yourself with the answer appearing as if by magic!
    – Geo...
    Commented Nov 9, 2022 at 20:45
  • 4
    it's also called rubberducking Commented Nov 9, 2022 at 21:19
  • @Geo... The funny thing this time was I had already explained the problem to an old Z80 programmer I met in a TRS-80 Discord server so I was surprised this only happened on explaining it the second time, a first for me I think. Commented Nov 10, 2022 at 0:26
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.

2
  • 2
    That state machine is a neat representation of the algorithm how to do long division on binary numbers. My answer has a specific pseudocode implementation of long division for dividing a byte and it keeps the state differently, by looking at two bits at a time, modifying the state by subtraction, and moving on to look at the next input bit.
    – Justme
    Commented Nov 10, 2022 at 11:34
  • These are (binary) Genaille rods for division by three. Took me a while to make the connection.
    – WimC
    Commented Nov 11, 2022 at 20:46
3

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.

3

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     
1
  • 1
    I just checked that algorithm (not the assembler code) in C: For numbers x>=128 this will round up the result of numbers with (x%3)=2. However, multiplying with 0.010101011 seems to work. This operation could be done by multiplying with 0xAB instead of 0x56 and shifting right by 9 bits. Or by initializing L with a/2 instead of 0 and multiplying with 0x55. Commented Dec 11, 2022 at 14:55
2

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

1
  • Instead of shr ax,9 you would use shr ah,1 on 8086 (or an adjusted version that only needs a right-shift by 8 so wouldn't need a shift at all.) Actual 8086 doesn't have immediate shifts, only by 1 or by a count in cl, and cycle cost is linear with the shift count, no barrel shifter until 386. But if shr reg, imm8 is available (186 and later), you'd use shl ax, 4, not a chain of 4x 2-byte add instructions. Commented Nov 11, 2022 at 7:41
2

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);
}
-2

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;
}
12
  • 3
    How do you convert something to base 3 without lookup tables or heavy math operations, before you can simply shift right a base-3 number? And maybe convert it back to something usable?
    – Justme
    Commented Dec 11, 2022 at 10:59
  • 1
    This is not an answer the question, which is "On CPUs without divide or multiply instructions like the Zilog Z80 and the MOS 6502, how could division by three be efficiently implemented?".
    – TonyM
    Commented Dec 11, 2022 at 11:20
  • 1
    Shift to right is not inefficient itself. But converting 8-bit integer from base 2 to base 3 and back to base 2 is, unless you can prove otherwise. That is the missing key part in your answer. And the CPU does not work in base-3 numbers, it works in base-2 numbers, so you can only shift base-2 numbers efficiently.
    – Justme
    Commented Dec 11, 2022 at 12:15
  • 1
    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. E.g. 255 = 63*4+3 = 63*3 + 66 = 63*3 + 16*4 + 2 = 63*3 + 16*3 + 18 = 63*3 + 16*3 + 4*4 +2 = 63*3 + 16*3 + 4*3 + 6 = 63*3 + 16*3 + 4*3 + 1*4 + 2 = 63*3 + 16*3 + 4*3 + 1*3 + 3 = 85*3. Commented Dec 11, 2022 at 13:03
  • 1
    You can easily form a recursive solution that way. I don't think I like it because it's no faster and only marginally smaller than existing answers. However it's also the only answer that can handle arbitrarily large inputs. (try adapting any of the other solutions to uint64). Commented Dec 11, 2022 at 13:28

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .