Do Developers Understand IEEE Floating Point?

According to the authors of the eponymous paper (http://pdinda.org/Papers/ipdps18.pdf), apparently not. They presented a quiz that was supposed to test basic understanding of FP to a bunch of software engineers and academics. The results were, to put it mildly, not great (skip to section II.B to see the types of questions they used and to IV.A for the punchline). Curiously, "Formal Training" in FP didn't help much (see Fig. 19).

The study might be flawed or biased for all sorts of reasons, but their findings definitely match my personal experience. To many otherwise very qualified software engineers, FP code is dark magic that you should avoid if you can and liberally splinkle with random fabs(x - answer) < EPS everywhere if you can't. Which is pretty sad considering that IEEE 754 actually gives you very strong guarantees on your computations (but you have to understand what they are).

Now, FP (and numerical analysis in general) is a very complex and nuanced topic, and in no way I'm an expert in this area (refer to https://www.springer.com/gp/book/9783319765259 if you want qualified advice). But at least I think I know a nice way to get a basic grasp on IEEE FP internals and hopefully demystify it: try and implement a FP32 (https://en.wikipedia.org/wiki/Single-precision_floating-point_format) to FP16 (https://en.wikipedia.org/wiki/Half-precision_floating-point_format) conversion routine.

Why?

* It doesn't require a lot of code. An unoptimized implementation easily fits in 25...30 lines.
* It's trivial to verify your routine is correct: there are only so much 32-bit floats, so you can just test them all (https://randomascii.wordpress.com/2014/01/27/theres-only-four-billion-floatsso-test-them-all/) and compare against a reference implementation (e.g., https://en.wikipedia.org/wiki/F16C).
* You will have to learn how exactly the bits of a FP number are arranged in memory (keywords: biased exponent, normalized significand, signed zeroes, positive/negative infinities).
* More importantly, because the conversion is lossy, you will have to learn how the slightly more obscure parts of the IEEE standard work (keywords: denormal numbers, round-to-nearest-even rounding mode, signaling/quiet NaN's).

By the way, if you want to know what an optimized version looks like, here's the code I eventually shipped in our FP16 library: https://github.com/catboost/catboost/blob/master/library/float16/float16.cpp#L14 All credit goes to Fabien Gessen and his ingenious public domain implementation; all I did was modify it to properly handle NaN payloads and add a copious number of comments.

When memory-safe code isn't

This short ~30-line Rust program [1] does not contain any unsafe keywords. It slices an array of length 3 and tries to access the 16777216'th element of the slice. Obviously, the resulting slice has length 2**64 - 1, so the program tries to write to an invalid memory location, segfaults and dies.

Oh wait.

This is not supposed to happen. In fact, the whole selling point of Rust is that it can't happen, provided that your code doesn't use unsafe, the unsafe code from the standard library is free from UB (see [2] for an attempt to formalize this), and there are no bugs in the Rust compiler. Any illegal memory access you make should terminate your program quickly and ruthlessly before you actually read or write data, unless you want your system to be completely compromised. If this sounds overly dramatic to you, read this OpenBSD [3] advisory. The OpenBSD developers take security extremely seriously (I mean this in a literal way, not in the "oops, your data has just been leaked" way), and yet they just found a remote code execution hole in their software, all due to a single out-of-bounds read.

It turns out that in this case both the Rust compiler and the Rust standard library have no bugs. The culprit is LLVM -- the code generation framework Rust uses internally. Fortunately, this particular issue is incredibly hard to trigger accidentally (see [4] for the whole investigation with all the gory details), but still, I had a couple of teachable moments when I was trying to debug this.

First, contrary to what I believed, Rust-related LLVM bugs happen all the time. In retrospect, this should be pretty obvious, since LLVM is at least twice as large as the Rust compiler and the standard library combined. An ongoing effort to upgrade Rust to LLVM 10 [5] lists lots of examples; as far as I can see, none of those are memory safety issues, but there are a lot of scary stuff there -- infinite loops in LLVM passes, compile-time regressions, performance regressions, etc.

And second, the complexity is there for a reason. The stuff LLVM does to your code is mind-boggling; already after a dozen of passes I can barely recognize the initial program I wrote, because all abstractions are stripped away and the control flow is aggressively restructured. I can't say that I understand even 10% of what's going on in LLVM passes, but the part I really liked (incidentally, the one that caused the LLVM bug in question) is the ScalarEvolution framework [6], [7], which allows you to infer all sorts of info about the loops in your program.

[1]: https://play.rust-lang.org/?version=stable&mode=release&edition=2018&gist=0e6866dbb74fd88796c73b218bffcbf6
[2]: https://plv.mpi-sws.org/rustbelt/stacked-borrows/paper.pdf
[3]: https://www.openwall.com/lists/oss-security/2020/02/24/5
[4]: https://github.com/rust-lang/rust/issues/69225
[5]: https://github.com/rust-lang/rust/pull/67759
[6]: https://www.youtube.com/watch?v=AmjliNp0_00
[7]: https://github.com/llvm/llvm-project/blob/master/llvm/lib/Analysis/ScalarEvolution.cpp

PSA: I'm too lazy to figure out how to enable comments in Telegram channels, but if you have any sort of feedback on anything I wrote, please send me a Telegram DM (@ivan_komarov). Or drop me a line at Ivan.Komarov@dfyz.info. Either way works.

Adventures in quantum factoring

It all started this summer when I finally got round to participating in this year's edition of Q# Coding Contest (Q# is a pretty slick F#-inspired programming language for developing quantum algorithms by Microsoft). Most of the challenges they wanted us to solve were entertaining, but kind of anticlimactic -- I mean, checking if a number is divisible by 3 is fun, but that's not exactly the problem we really need quantum computers for. Sometime after the contest, my curiosity got the best of me, and I started to think if I could find a small but useful algorithm in the Q# documentation that I could implement on my own to get a sense of what real quantum programming looks like.

This is where it hit me -- I may not know or care much about quantum chemistry and Hartree–Fock theory, but I'm fond of cryptography, and post-quantum crypto is all the rage these days. To understand why exactly we suddenly need to deprecate basically all currently deployed public-key crypto, I set out to implement Shor's algorithm for the quantum simulator built into Q# to try and break RSA factor small integers.

When I was doing preliminary research, I stumbled upon a paper with a poignant title that showed how easy it is to (intentionally or accidentally) cheat by creating a specialized quantum circuit for factoring a particular integer if you already know its factors. Naturally, I wanted to avoid that and implement a fully general algorithm. I was also determined to use as few qubits as possible because I hoped to eventually factor 15 (or 21 if I get fortunate) on real quantum hardware. With that in mind, I found a beginner-friendly paper that looked sufficiently close to SOTA in terms of the number of qubits used and got to work.

In contrast to the Q# Coding Contest challenges, Shor's algorithm has a non-trivial classical part that has to run on a regular computer. In fact, the algorithm is not really about factoring; it's about finding a period of f(x) = a^x mod n. So first you have to do a bunch of work to recast factoring in terms of period finding, and then use the period p you found to identify the factors you are after. It gets a little involved because what you actually get from the quantum part of the algorithm is not p proper but an approximation to x/p (where x is a number less than p), but all in all, the classical part is not terribly exciting.

Not surprisingly, the exciting part is the quantum circuit for period finding. There's a really beautiful algorithm called Quantum Phase Estimation (see this excellent book for more details) that says you can accurately estimate the eigenvalue of a quantum operator by preparing an eigenstate, applying the operator a few times controlled on a qubit that is in an equal superposition of basis states, applying an inverse Fourier transform, and then measuring. For period finding, the quantum operator is simply the modular multiplication modulo n. However, it turns out that, while quantum computers excel at computing Fourier transforms, they are astonishingly bad at such complex operations as checks notes adding two numbers together modulo another number, since all quantum operations have to be reversible, which means you have resort to tricks such as dummy quantum registers that are always set to zero (and sure enough, more Fourier transforms).

Debugging this is even harder than debugging CUDA code, but it feels like magic when it finally works. My initial implementation would give me periods that kinda sorta looked right (hard to tell, it's an approximation after all), but then I looked through my code and found that I messed up the order of bits when rotating qubits. When I fixed that, the periods suddenly started to look exactly right; at that moment, I gasped quite literally and audibly.

This is the code I ended up with. Bugs and performance problems of course are possible (I would appreciate a code review from my quantum-savvy readers!), but at least this program can factor all integers not greater than 100 in a few minutes on my laptop. The largest integer I tried was 209, which took 141 seconds to factor.

It would be only logical to conclude with a report on how I ran my implementation on a real quantum computer, but sadly this is the part where I need to rant a little. The amount of hype around quantum computing is incredible; it's worse than Deep Learning and Bitcoin combined. Everyone (Microsoft included) keeps blabbing endlessly about Harnessing the Quantum Power to Transform Your Digital Business, Getting Ready To Tap Into Quantum Potential, and Quantum-Inspired Optimization (whatever that means). It's next to impossible even to know my program will run at all. Will there be enough qubits available to factor at least 15? (probably yes) Will I exceed the circuit size limit? (probably not) The algorithm I use implies it is possible to measure a qubit twice, can I do that? (who the hell knows) At the time of writing, as far as I can tell, the only company with something vaguely resembling a specification of their hardware is IBM, and, as expected, they don't support Q#.

I'm still moderately optimistic because Azure Quantum "will move to public preview later this year", so worst-case scenario, I'll just have to debug my code by trial and error later this year. Stay tuned for updates!

Feedback needed

I'm trying to decide on what the next post in this channel should be about. I basically have two options on my mind at this moment.

First, I just completed the Flare-On 2020 challenge. This was a brutal CTF-like reverse engineering marathon that lasted one month for me, and I'm thinking of creating a write-up like this once it's officially over. I'm not remotely as good at write-ups as MalwareUnicorn is, but I hope mine can be at least moderately entertaining (I definitely learned a lot this year).

Second, I have a DEC Alpha multiprocessor system (I bought one when I was playing around with CPU memory models), and a partially completed project to restore the LLVM backend for Alpha. They used to have a backend in LLVM quite a few years ago, but it was eventually ditched as no one in their right mind uses Alpha anymore. Except me, I guess. I'm not sure I have what it takes to finish this project any time soon, but, again, I learned a lot while working on it, so why not write a few paragraphs about weird CPU instructions, selection DAGs, and LLVM backends?

Please help me choose!

Understanding Software Dynamics, part 1

I'm sorry for the, uh, slight delay with that DEC Alpha post I promised. I still hope I'm going to finish it someday.

Instead, I want to break radio silence with a series of posts about the recent Dick Sites' book called Understanding Software Dynamics, which is one of the best books about performance engineering I've ever seen. This certainly is the first technical book I read from cover to cover in quite a while.

In the first post of the series, I want to explore one of the key takeaways from the book: "when estimates and reality differ, there is always good learning". This statement sounds extremely trite, but it perfectly describes what happened to me when I was doing one of the earlier exercises (yep, it comes with the exercises!) from the book. You are supposed to mark the incr variable as volatile in the loop that looks like this:

static const int kIterations = 1000 * 1000000;
uint64_t sum = 0;
/*volatile added here*/ volatile int incr = ...; // the value of incr doesn't matter

for (int i = 0; i < kIterations; ++i) {
 sum += incr;
}

Then Dick wants you to re-compile the program with both -O0 and -O2, and explain any performance changes (or lack of changes) compared to the volatile-less version.

The running time of the -O2 version after adding volatile matched my expectations pretty closely (the compiler stopped optimizing the loop away).

I didn't expect any major performance changes in the -O0 version, though. Surprisingly, compiling this with -O0 resulted in a ≈2.5x increase in the number of CPU cycles consumed by the loop!

There could be dozens of plausible reasons behind this, so initially I didn't think much of it. However, inspecting the disassembly revealed that both versions compiled to almost identical code. The only change was the order of i and incr on the stack.

More precisely, when volatile is absent, i is stored at [rbp-0x30], and incr is at [rbp-0x2c]. After adding volatile, the locations of i and incr are swapped.

(if you wonder why the variables are spilled to the stack in the first place, remember it's an -O0 build)

The slowdown was easily reproducible, so I couldn't shrug this off as some fluke. I started ruling out the possible suspects:
* both programs computed the same result
* they both ran at the same CPU frequency, and the wall times were consistent with the number of cycles
* no other process was using the CPU core (or its logical sibling) the programs were running on
* both incr and i comfortably fit in the L1 data cache in both programs, so the number of L1d accesses was essentially equal to the number of memory accesses, and the number of L1d misses was essentially zero
* the programs were of exactly the same size and the layout of code/data was exactly the same, which excluded weird loop alignment issues

None of this made any sense to me. Both programs executed the same straightforward single-threaded loop with one addition, two loads, and one store, but one of them was running circles around the other.

Since my CPU was an AMD Ryzen 3970X, I eventually resorted to running AMD μProf and looking at every possible performance counter I could find. I failed to find any discernible difference. This is a small subset of the counters I tried that, as far as I can see, says "these programs behave identically, except that the one at the bottom wastes significantly more CPU cycles".

In desperation, I distilled the full-blown C code into a small, self-contained assembly program and started tinkering randomly, trying to figure out what makes slowdown go away. I observed the following facts, which are not that useful but moderately interesting:
* The C code stored i and incr as 32-bit variables, extending them with cdqe when appropriate. For simplicity, I switched to 64-bit variables. The problem persisted.
* The sum variable was also spilled to the stack. I moved it to a register, the problem persisted.
* Converting add [...], 1 into inc [...] also

didn't change anything.
* I started accessing the variables directly with [rsp + X] instead of propertly tracking rbp and using [rbp - X]. The problem persisted.

At this point, I began playing with stack offsets, and this is where some curious patterns finally started to emerge:
* The slowdown happens only when then the stack offset (i.e., 0x30 in [rbp-0x30]) of i is not divisible by 8. The absolute value of the offset and the actual address of i can be arbitrary.
* The offset and address of incr are irrelevant. That is, if the stack offset of i is divisible by 8, there is no slowdown.
* On every loop iteration, i is first modified in memory with ++i (inc [...]), and then loaded from memory again for i < kIterations (cmp [...], 999999999).

Finally, with a little help from my friends and their hardware, I discovered something amazing: the slowdown only happens on AMD Zen 2 processors! Both programs have the same speed when run on AMD Zen 3 or Intel CPUs (unfortunately, I couldn't get my hands on a Zen 1).

Example output on AMD Zen 2:

$ ./test_amd_puzzler.sh
Mode
l name:                      AMD Ryzen Threadripper 3970X 32-Core Processor
FAST VERSION (bin size = 9200)
cycles = 2392617802
result = 42000000000
SLOW VERSION (bin size = 9200)
cycles = 6092499439
result = 42000000000

Example output on AMD Zen 3:

$ ./test_amd_puzzler.sh
Model name:                      AMD Ryzen 7 5800HS with Radeon Graphics
FAST VERSION (bin size = 9248)
cycles = 1695570624
result = 42000000000
SLOW VERSION (bin size = 9248)
cycles = 1692961866
result = 42000000000

Example Intel output:

$ ./test_amd_puzzler.sh
Model name:          Westmere E56xx/L56xx/X56xx (Nehalem-C)
FAST VERSION (bin size = 9352)
cycles = 6752197154
result = 42000000000
SLOW VERSION (bin size = 9352)
cycles = 6722025261
result = 42000000000

Based on this, I strongly suspect that the fast version of the program hits the obscure undocumented Zen-2-only feature Anger Fog found a couple of years ago. He calls it "Mirroring memory operands" (see section 21.18). Essentially, under some conditions, the CPU allows you to perform successive store/load/RMW instructions operating on the same memory address with zero latency. One of these conditions goes like this:

> The memory address may have an index register, a scale factor, and an offset. The offset must be divisible by 4 in 32 bit mode, and by 8 in 64 bit mode.

Sounds really familiar! In other words, it's not that the slow version is unreasonably slow. It's more like the fast version is unreasonably fast.

I'm still not 100% sure my explanation is correct, for two reasons:
* Agner Fog says that another condition is that "the offset must be in the range -128 to +127 to fit into an 8-bit signed integer", but in my experiments, offsets as high as 1024 worked just fine.
* He also says that the memory operand mirroring machinery is completely absent from Zen 3. However, the number of cycles in the output above for Zen 3 is so low there is a possibility that whatever makes the fast version go fast on Zen 2 might also apply to Zen 3 (although it doesn't depend on offsets being divisible by 8 anymore).

If anyone knows how to definitively find out what's going on, please do let me know! You can run your own experiments by modifying the assembly and then running the same shell script I used to get the example outputs above.