Ideal divisors: when a division compiles down to just a multiplication

A few months back I wrote a blog post very much inspired by Dan’s post, though a little more academic: https://kevinventullo.com/2020/12/21/2-adic-logarithms-and-f...

The idea is to exploit deeper structures in p-adic number theory to develop useful and novel algorithms with unsigned integer arithmetic, e.g. fast integer exponentiation. Feedback welcome!

(Warning: math heavy)

> I don't know how gcc generates its approximate 2-adic reciprocals. Possibly it uses something based on the Euclidean GCD algorithm. I wasn't able to find the precise line of source in a reasonable time.

It's 0x1FFF...F / 7 = 0x249...249249249.

Just about the oddest thing I've seen; I guess, if I were the author, it'd be my "1 in 10000" day?

There's no specific optimisation for these specific numbers. There's a more general optimisation to turn division by a fixed number into a multiplication and shift where possible (because division is a lot slower than multiplication), and there is another optimisation that eliminates useless shifts. Combined they yield this result for these specific numbers.

Division is orders of magnitude slower than multiplication.

They're often special cases of general algorithms. You can do many divisions as combinations of multiplications & shifts; in this case, the resulting code simplifies down to the smallest possible.

I'll run through a base 2 example and you should be able to extrapolate for any number

eg. find the equivalent multiplier of 1/31 in 8bits. This must exist because 31 has co-prime factors to 2^8 (the rules of a solvable diophantine require co-prime factors and we use that below).

Now the goal is to find a number that when multiplied by 31 and mod 2^8 = 1. This will give us a multiplicative equivalent to 1/31.

We can use the diophantine equation. 31x - 256y = 1. This equation asks 'what multiplied by 31 and subtracted by any number of 2^8 gives us 1?'. ie. It's asking what multiplier = 1/31 % 2^8

(I'm not going to work through the diophantine equation here but it's a log(n) equation and easy to do but I'm lazy so using wolfram alpha which has the diophantine equation built in)

Solving diophantine(31x - 256y = 1) tells us 'x must be of the form 256n + 233'. So (256n + 223)*31 = 1 under modulo 256

So now we know whenever we see a division by 31 we can multiply by 223 (or 256n + 223 if you really want a larger multiplier) and we'll get a number that's equivalent to 1/31 under modulo 2^8. You can do the above process for any number under any base as long as the number is co-prime to the base.

It's typically faster to calculate a reciprocal (especially at compile time) and then multiply it, but this is often at the expense of precision. In C++, for example, the -O3 flag on gcc will still emit a division (divss) to preserve compliant floating point precision in this example: https://godbolt.org/z/PY3Wjno4P

However, when we enable `-ffast-math`, we permit the compiler to do a bit more aggressive optimization in this domain (trading off precision), and we see multiplication (mulss) is emitted: https://godbolt.org/z/KEb1nc43z

This is computer science, not algebra.

The linked article is about integer division and has nothing to do with floating point, and in fact I can't see how this idea could be applied to floating-point arithmetic. Can it?

I have a 27" monitor also, and it looks fine to me.