Reciprocal Approximation with 1 Subtraction

Yep. Along the same lines. This one is even simpler, though, as it requires only a single integer CPU instruction (and the simplest of all instructions too).

If you want full precision, you need to do three Newton-Raphson iterations after the initial approximation. One iteration is:

    y = y * (2.0F - x * y);

For some more explanation: the main idea behind both tricks is that the IEEE floating point formats are designed so that the most significant bits represent its exponent, that is, floor(log2(x)). Hence reinterpret-casting a float x to an unsigned integer uint(x) approximates a multiple of log2(x). So these kinds of approximations work like logarithm arithmetic: float(C + a*uint(x)) approximates x^a, for a suitable constant C. Quake's invsqrt is a=-1/2, this post is a=-1.

More detail on https://en.wikipedia.org/wiki/Fast_inverse_square_root#Alias... .

IEEE754 floats are a very well-designed binary format, and the fact that these approximations are possible is part of this design; indeed, the first known instance of this trick is for a=1/2 by W. Kahan, the main designer of IEE754.

There are basic integer operations in the FP/SIMD units on most CPUs, so there’s no generally need to “move back and forth” unless you need to branch on the result of a comparison, use a value as an address, or do some more specialized arithmetic.

These kinds of tricks are still used today. They're not so useful if you need a reciprocal or square root, since CPUs now have dedicated hardware for that, but it's different if you need a _cube_ root or x^(1/2.4).

There is also the case with machines that lack FP support, like some ARM Cortex M variants.

EDIT: after double-checking my work, I realized I have a better bound on maximum error, but not a better average error. So, the magic number depends on the goal or metric, but mean relative error seems reasonable. Leaving my original comment here, but note the big caveat that I’m half wrong.

One can do better still - 0x7EF311BC is a near optimal value at least for inputs in the range of [0.001 .. 1000].

The simple explanation here is:

The post’s number 0x7EFFFFFF results in an approximation that is always equal to or greater than 1/x. The value 0x7EEEEBB3 is better, but it’s less than 1/x around 2/3rds of the time. My number 0x7EF311BC appears to be as well balanced as you can get, half the time greater and half the time less than 1/x.

To find this number, I have a Jupyter notebook that plots the maximum absolute value of relative error over a range of inputs, for a range of magic constants. Once it’s setup, it’s pretty easy to manually binary search and find the minimum. The plot of max error looks like a big “V”. (Edit while the plot of mean error looks like a big “U” near the minimum.

The optimal number does depend on the input range, and using a different range or allowing all finite floats will change where the optimal magic value is. The optimal magic number will also change if you add one or more Newton iterations, like in that github snippet (and also seen in the ‘quake trick’ code).

PPS maybe 0x7EF0F7D0 is a pretty good candidate for minimizing the average relative error…?

  Initial approximation:
    Good bits min: 4
    Good bits avg: 5.242649912834
    Error max: 0.0505102872849 (4.30728 bits)
    Error avg: 0.0327344845327 (4.93304 bits)

  1 NR step:
    Good bits min: 8
    Good bits avg: 10.642581939697
    Error max: 0.00255139507338 (8.61450 bits)
    Error avg: 0.00132373889641 (9.56117 bits)

  2 NR steps:
    Good bits min: 17
    Good bits avg: 19.922843217850
    Error max: 6.62494557693e-06 (17.20366 bits)
    Error avg: 2.62858584054e-06 (18.53728 bits)

  3 NR steps:
    Good bits min: 23
    Good bits avg: 23.674004554749
    Error max: 1.19249960972e-07 (22.99951 bits)
    Error avg: 3.44158509521e-08 (24.79235 bits)

We used memcpy everywhere in our runtime and after the 10th or so time doing it, it becomes less awkward.

Yes, most casts via pointers are formally UB in both C and C++, and you can induce weird behavior by compilers if you transmit casted pointers through a function boundary (see [0] for a standard example: notice how the write to *pf vanishes into thin air). Since people like doing it in practice, GCC and Clang have an -fno-strict-aliasing flag to disable this optimization, and the MSVC compiler doesn't use it in the first place (except for restrict pointers in C). They don't go too far with it regardless, since lots of code using the POSIX socket API casts around struct pointers as a matter of course.

Apart from memcpy(), the 'allowed' methods include unions in C (writing to one member and reading from another), and bit_cast<T>() and std::start_lifetime_as<T>() in C++.

[0] https://godbolt.org/z/dxMMfazoq

It would be better if compilers were fixed. Current treatment of UB is insane. The default should be no UB and a pragma that lets you turn it on for specific parts. That used to be the case and some of the world's most performant software was written for compilers that didn't assume no pointer aliasing nor no integer overflow (case in point: Quake).

The big problem, apart from "memcpy is a bit tedious to write", is that it is impossible to guarantee absence of UB in a large program. You will get inlined functions from random headers and suddenly someone somewhere is "dereferencing" a pointer by storing it as a reference, and then your null pointer checks are deleted. Or an integer overflows and your "i>=0" array bounds assert is deleted. I have seen that happen several times and each time the various bits of code all look innocent enough until you sit down and thoroughly read each function and the functions it calls and the functions they call. So we compile with the worst UB turned off (null is a valid address, integers can overflow, pointers can alias, enums can store any integer) and honestly, we cannot see a measurable difference in performance.

UB as it is currently implemented is simply an enormous footgun. It would be a lot more useful if there were a profiler that could find parts of the code, which would benefit from UB optimisations. Then we could focus on making those parts UB-safe, add asserts for the debug/testing builds, and turn on more optimisations for them. I am quite certain nobody can write a fully UB-free program. For example, did you know that multiplying two "unsigned short" variables is UB? Can you guarantee that across all the template instantiations for custom vector types you've written you never multiply unsigned shorts? I'll leave it as an exercise to the reader as to why that is.

A few years ago with GCC 10, I had some strange UD2 instructions (I think) inserted around some inline assembly that was fixed by replaced my `reinterpret_cast`s with `bit_cast`.

Yes, if implemented as shown.

You could use intrinsics for the bit casting & so on and it would be well-defined to the extent that those intrinsics are.

(I understand some people consider SIMD intrinsics in general to be UB at language level, in part because they let you switch between floats & ints in a hardware-specific way.)

The proper solution is to use std::bit_cast in modern C++ or otherwise use memcpy, and of course know what you're doing.

Some things that could mess with you:

* Floating-point endianity is not the same as integer endianity.

* Floating-point alignment requirements differ from integer alignment requirements.

* The compuler is configured to use something else than 32-bit binary32 IEEE 754 for the type "float".

* The computer does not use two's complement arithmetic for integers.

In practice, these are not real problems.

Not sure I understand your number 0x7eb504f3 - does it require using the 1.385 factor? Is it possible there is a typo in the number? That value doesn’t make sense to me. I’m measuring error here as the absolute value of relative error, e.g., | (v - v_approx) / v |. With that constant alone plugged into the poster’s code, I get a much less accurate approximation, with a minimum error of at least ~0.27 (meaning it’s always far away from the target) and worst case error of ~0.293 over the input range [1/1000…1000]. For comparison, the original 0x7effffff error is min:0 max:0.125. With 0x7eeeebb3, the error I get is min:0 max:0.0667. With 0x7ef311bc, error is min:0 max:0.0505. It’s important that the min error over a large interval is zero, because it means the approximation actually touches the target at least once.

Engineering mathematical calculations always looks like "jank." Take a look at Plauger's The C Standard Library. It is full of "magic" numbers.

But then again, representing irrational numbers in binary is inherently janky.

Yes of course there’s a more graceful approach. You can use an actual divide instruction or routine instead of the bit-casting subtraction trick.

It’s not about interfaces. The point is that a single subtract is simpler math, and (probably) faster than a divide, and it gets you surprisingly close to the right answer. The reason it works is subtracting two exponents is a subtraction in log space which is equivalent to a divide in linear space.

People don’t use this trick in practice that much if at all, especially these days with GPUs that do reciprocals and square roots in hardware, it’s just an interesting thing to know about that helps solidify our understanding of floating point numbers and logarithms.

If you're not constrained to software solutions you have a whole world of opportunities. E.g. if it's a graphics or neural net pipeline you can pour tricks like this (or better) onto it. If it's a CPU then you can add special instructions that do exponent manipulation and the likes.

For more general algorithms along these lines, see exercises 25–28 in section 1.2.2 of Knuth[1] (and note that the printed solution to exercise 28 in early printings has an error[2]).

[1] D.E. Knuth, The art of computer programming. Vol. 1, third edition, Addison-Wesley, Reading, MA, 1997.

[2] https://www.werkema.com/2021/10/14/my-knuth-check/

Excellent! Will have a look.