Collective groan from the 6502.org community.

But wait! I'm just wondering if I can improve on what I've got, and there's no better collection of minds than here.

So, I have two unsigned 16-bit values A and B, where A<B, and I want to compute A/B to 8-bit fractional precision.

My starting point was to look at a simpler case: dividing two 8 bit values, or, in this case, 16-bit divided by 8-bit, yielding an 8-bit value. To my moderate surprise, such a routine is already presented in the 6502 wiki:

Given that the LSB of the numerator is always zero, this can be subtly optimized to the following:

That actually looks pretty efficient, and I can use it if I know that both my dividends are <256. But I needed a 16-bit case, so I came up with the following:

This normally completes (depending on the inputs) in between 250 and 300 cycles, which is OK, but I'm wondering if I can do better.

I wondered if there was something in doing the subtract, caching the result, and only committing it if it didn't overflow, but I feel like the routine above can short-circuit that work when not required, and probably gives a small performance gain (I guess I would need to measure it).

But the main reason I'm posting this is to ask if anyone knows a better way, or a totally different way! I'm interested in performance here, so if there's a table-based approach which uses a modest amount of RAM for tables, I'm open to ideas.

Note that I don't need an enormous amount of precision: 1/256 is fine. The use case here is as an input to an arctan table: an index from 0-255 (representing 0...1) yielding an angle between 0 and 45 degrees, where 45 degrees is represented by 128. With this kind of resolution, I don't need any more precision in the division results, as every angle is covered in the table. What I do need is to calculate hundreds of these things in real time, so, the quicker the better!

Over to you Hope everyone's doing OK during these strange times.

Have you looked into non-restoring division? That's often faster. Instead of comparing and then subtracting, you add or subtract depending on the previous result. (I'm not sure if it will be faster in this case, as with a bytewise comparison you can get away with doing half the work most of the time.)

Not sure really what you mean there Ed - do you have an example?

The next thing I tried was to try forcing the dividends to be 8-bit in length by shifting them both until the highest set bit of the denominator is at bit 7 of the 8-bit input.

This finishes, more typically, in around 170 cycles, which is a distinct improvement, but the results are less accurate, occasionally coming out 1 less than the correct result. I think this is within an acceptable error tolerance, and it could possibly be corrected by rounding the dividends up according to the underflow bit. Any other ideas?

_________________
http://WilsonMinesCo.com/ lots of 6502 resources
The "second front page" is http://wilsonminesco.com/links.html .
What's an additional VIA among friends, anyhow?

Thanks Garth!

Yes, of course I went to the usual places to see what was already on the 6502.org site, and I'm painfully aware that 6502 division is a FAQ, but I figured my needs were sufficiently distinct from the "generic" 32-bit/16-bit division code that it was worth starting a topic on it.

I've always found the 6502 division algorithm surprisingly arcane, and I always need to re-derive it each time I need a variant on it, and I guess it's no wonder that it's a topic which comes up time and time again.

In my case, I'm happy with an approximation (already the required precision isn't really that great!) so I guess that's more the direction I'd like discussion to go. Just as there's the famous "difference of two squares" multiplication hack, I was wondering if a table-based approach to division was viable. I guess my next test can be to try a reciprocal multiply with the f(x)=x*x/4 tables, after having coerced the dividends into 8-bit approximate quantities. Anticipating better performance and worse accuracy (as is the way).

_________________
http://WilsonMinesCo.com/ lots of 6502 resources
The "second front page" is http://wilsonminesco.com/links.html .
What's an additional VIA among friends, anyhow?

I haven't been here for a while! Just amusing myself with a little 6502 prototype during these strange lockdown days.

I hadn't come across your large tables page before, and, while I don't have much more than about 4k to spend on tables, I've convinced myself that "difference of squares" multiplication tables, combined with a LSB/MSB reciprocal multiplication table for 256 values is probably enough for what I need. I'm going to keep accumulating truncation errors though, so I'm hoping it can all be ironed out by rounding up here and there. If errors are consistent, it won't be as bad as if they vary with inputs and lead to 'wobbling' outputs!

I might have a go at this, as I have an idea for how to go about it, but it'll have to wait until some urgent day-job shenanigans gets out of the way. It's something that I unfortunately have to spend even my Sunday preparing for, and it doesn't leave me with the mental bandwidth necessary to really get to grips with the 6502.

Tomorrow I'm going to try the following:

A reciprocal table of 16-bit precision: f(x)=65535/x for x=0..255.
The usual difference of squares multiplication table: f(x)=x*x/4 for x=0..511

To divide unsigned A/B, first shift down both until the highest bit in B is at bit 7. Since typically the values are <4000, it's more efficient to shift down. Loop is unrolled, of course No obvious accuracy win in rounding up (ADC #0 using the final carry shifted out).

Then, using the reciprocal table and the difference of squares multiplication, calculate:
lsb(A * recip_msb(B)) + msb(A * recip_lsb(B))

Quick tests on various pairs of values suggest that it won't give results any worse than those I was getting before, but I need to implement it and test some corner cases.

Final step is to put that result into an arctan table: f(x)=arctan(x/256)*512/pi for x=0..255 to get an approximate bearing in an octant of a circle. Phew!

Okay, so here's a quick and dirty version of compare-and-conditional-subtract, in which the subtraction itself is also the comparison: This is actually pretty straightforward and compact, but I calculate it as taking about 304 cycles on average (ie. when half the result bits are set). Minimum runtime is 284, excluding JSR/RTS, when the result is zero.

OTOH, I don't see how your original 16-bit routine gets to under 300 cycles on average, as the "typical" code path I can see is 45 cycles times 8 iterations, which is 360. Is there something non-obvious about the input data which makes the fast path happen more often than the slow and medium ones?

Not exactly but I found a previous thread where Bruce chimes in:

integer division and multiplication

Edit: if you have fast multiplication sorted out, maybe try Goldschmidt division.

That's a nice way to do it - thanks! I've used a similar trick when consuming data bit-by-bit (to determine when you need to buffer your next byte of data without having to count bits). I unrolled the 8-bit version, though it probably wouldn't be particularly worth doing it here. It might be worth trying to keep some value in the accumulator to save a memory shift each time round.


To be honest, I didn't take a particularly rigorous approach to timing it. I just used BBC BASIC to set up random inputs, and used a VIA timer to time the routine (minus the JSR/RTS) with interrupts off. Here's an example of the results (input, correct result, calculated result, cycle count):

This is with B from 0-65535, and A<B.

Interestingly, if I reduce the range of B to 0-4000, the average timing is generally higher:

Some more detailed analysis could undoubtedly reveal why, but I don't have the time or energy to really look at it at the moment!

I wanted to try your code and drop some timing results in, but I was finding I was getting the wrong results about 50% of the time - going to take a look to see if I've made a mistake somewhere, and will report back! Timing-wise, it was looking to be far more consistent, taking around 300 cycles each time.

Hmm, I don't think your code is equivalent - specifically, the algorithm needs to always subtract if left-shifting the numerator at the top of the loop shifted a set bit out. Unfortunately I can't think of a neat way to express that in your code without ruining its sleekness!

Example output:

Thanks for links Ed! Will try digesting that with a coffee! A skim read suggests that Goldschmidt division works with reals rather than integers, but I haven't really taken it in properly.

I think if I can get the required precision with a table-based approach, I'll go with that. The problem is that tables generally require reducing the domain of the inputs, but I've found so far that reducing the dividend and divisor to 8-bit values reduces the accuracy of the results (up to +/-2 of the correct result). I will probably try various different implementations (reciprocal multiply vs 8-bit division vs 16-bit division) and see how bad I can bear the inaccuracies to be, compared to the performance advantage.

Here are the results of the table-based approach.

First of all, the tables. They are page-aligned and need 2.5k:

Next some one-time initialization to set up some zero page locations:

And now the division code:

This trick for multiplication is well-known, so I won't go over it again unless anyone would like an explanation. The nice thing is that it's extremely quick to get a result for A*B (A, B both unsigned 8-bit). What's clear is that the most work is now done normalizing the inputs so they fit in 8 bits. To optimize this better, it would check the top 4 bits of the denominator, and, if non-zero, shift in the other direction into the high byte, so that no more than 4 shifts ever need to be performed. For the purposes of this test, I'll keep the inputs to 12-bit values.

The final part of the routine is multiplying the numerator, A, by the reciprocal of the denominator, B. Since the reciprocal is a 16-bit quantity, we actually have to do:

A * recip_msb * 256 + A * recip_lsb

We are only concerned with bits 8-15 of the result, so we add the LSB of (A * recip_msb) to the MSB of (A * recip_lsb).

Results time (inputs, actual value, computed value, cycles taken):

It's so much faster than any other approach, but there are plenty of off-by-one errors (which is more due to the truncation of the inputs to 8-bit quantities than any errors from the reciprocal table, although the latter doesn't help). Whether this is sufficiently accurate remains to be seen, but it's nearly three times quicker than the full division approach, so it's definitely a contender.