Dividing unsigned 8-bit numbers

They are actually useful for optimising compilers too! Mul or mul+shifts is often faster than div

Yes, if you know something is an exact multiple of n = r*2^k where r is odd, you can divide out the multiple by right-shifting k followed by modular multiplication by the modular multiplicative inverse of r.

       00011011 (27)
     x 01010101 (-0.333...)
    ------------------------
       11110111 (-9)
    
    invert and add one:
    
       00001001 (9, like we wanted)

Can you provide an example with details? Thanks!

Sort of. (After all, a "reciprocal" for our purposes is just a multiplicative inverse in the reals, so it makes sense that it would be related to the multiplicative inverse in other domains.)

3 times 171 is 513. So to divide by 3, we could multiply by 171 and then divide by 513. Dividing by 513... isn't any easier, but 513 is close to 512, so we hope that dividing by 512 (which is trivial - just a right-shift) gets us close enough.

Suppose for a moment we try dividing 3 by 3 using this trick. First we'll multiply by 171 to get 513. Consider that value in binary:

    1000000001
     ^^^^^^^^^
      ~~~~~~~~

Dividing by 512 will shift away the ^ bits. For floor division, we therefore want the ^ underlined value to be close to zero. (That way, when we divide 255 (say) by 3, the error won't be big enough to overflow into the result bits.)

The multiplicative inverse fact is equivalent to telling us that the ~ underlined bits are exactly 1. Conveniently, that's close to 0 - but we didn't account for all the ^ underlined bits. For example, the multiplicative inverse of 7 (mod 256) is 183, but 7 times 183 is 1281. That's close to 1280, but that doesn't really help us - we could right-shift by 8 but then we still have to divide by 5. If we just ignore the problem and divide by 1024 (right-shift by 10), of course we get a lot of wrong results. (Even 6 / 7 would give us 1 instead of 0.)

It also turns out that we'll need more bits for accurate results in the general case. I think it's possible without overflowing 16-bit numbers, but it definitely requires a bit more trickery for problematic divisors like (from my testing) 195. I thought I remembered the details here but proved myself wrong at the Python REPL :(

    uint8_t div195(uint8_t x) { return x >= 195; }

It's not entirely a coincidence but also not a general result that one should use the modular inverse as multiplier.

171 * 3 = 2^9 + 1, which is not surprising as we know that 171 * 3 = 1 (mod 2^8). So rearranged we have 171 / 2^9 = 1/3 + 1/(2^9*3) which shows it's a close approximation of 1/3.

Worst-case, you just multiply the result by the divisor and subtract that from the dividend.

Put in a couple hours to see how well it worked. The approximately correct form is slightly faster than just doing an FP division without proper rounding.

    uint8_t magic_divide(uint8_t a, uint8_t b) {
        float x = static_cast<float>(b);
        int i = 0x7eb504f3 - std::bit_cast<int>(x);
        float reciprocal = std::bit_cast<float>(i);
        return a*1.94285123f*reciprocal*(-x*reciprocal+1.43566f);
    }

No errors, ~20 cycles latency on skylake by uiCA. You can optimize it a bit further by approximating the newton raphson step though:

    uint8_t approx_magic_divide(uint8_t a, uint8_t b) {
        float x = static_cast<float>(b);
        int i = 0x7eb504f3 - std::bit_cast<int>(x);
        float reciprocal = std::bit_cast<float>(i);
        reciprocal = 1.385356f*reciprocal;
        return a*reciprocal;
    }

Slightly off when a is a multiple of b, but approximately correct. 15.5 cycles on skylake because of the dependency chain.

Don't have time to figure out the actual numbers myself, but here's an example doing the same for various transcendental functions: https://github.com/J-Montgomery/hacky_float_math/blob/623ee9...

a*(1.xxx/b) = a*(-1*1.xxx*log2(b)), which means you should be able to do a*(fma(b, magic, constant)) with appropriate conversions on either side. Should work in 32 bits for u8s.

[deleted]

Alderlake supports AVX512-FP16. Still only 9.6x faster than div. Most likely reciprocal is just too slow.

I've never thought of a micro-lookup-table!

That's cool. What have you used those for?

pretty sure a memory access is faster than the methods presented in the article.

If you scale up the multipliers, you should be able to eliminate the variable shift, which would reduce the lookup table to 512 bytes.

This seems like something that could be in hardware to allow native execution of the instruction. Is that not the case anywhere?

If a SIMD ISA exists, someone must know that it exists, because definitionally we only apply the term "SIMD ISA" to things that were consciously created to be such. So we could simply check every such example. Saying "only known example" is indeed silly.

But e.g. in mathematics, if we say that "almost every member of set X has property Y; the only known exception is Z" then there absolutely could be more exceptions, even if we pool the knowledge of every mathematician. It isn't necessary that X is finite, or even enumerable. It could be possible for exceptions other than Z to exist even though every other member of the set that we know about has the property. It could be possible to prove that there are at most finitely many exceptions in an infinite set, and only know of Z but not be able to rule out the possibility of more exceptions than that.

We don't even need to appeal to infinities. For example, there are problems in discrete math where nobody has found the exact answer (which necessarily is integer, by the construction of the problem) but we can prove upper and lower bounds. Suppose we find a problem where the known bounds are very tight (but not exact) and the bounded value is positive. Now, construct a set of integers ranging from 0 up to (proven upper bound + 1) inclusive... you can probably see where this is going.

The latter doesn't apply to SIMD ISAs, because we know all the interesting (big hand-wave!) properties of all of them rather precisely - since they're designed to have those properties.

Yeah but my point is that as a reader I'm trying to figure out which ISAs actually don't provide this (vs. which ISAs the author lacks knowledge of), and I still don't know what those are. The sentence looks like it's supposed to tell me, but it doesn't.