Reading a file format in 300 lines

Ferro tests itself against reference data generated by NumPy, which means something has to read NumPy's file format. That format is called .npy, and writing a reader for it is the first real code in the project.

It is also a good tour of the small, sharp details that make numerical code different from ordinary code. Almost every one of the complications below looks like a corner case and is actually a thing that will happen.

#What is in the file

The .npy format is refreshingly simple. A file has three parts, in order:

1. A six-byte magic string, \x93NUMPY, followed by two version bytes.
2. A header length, then the header itself — which is a Python dictionary written out as plain text.
3. The raw bytes of the array, with nothing in between them.

The header of a small array looks like this, padded with spaces:

{'descr': '<f4', 'fortran_order': False, 'shape': (2, 3, 4), }

That is genuinely a Python dict literal sitting in the middle of a binary file. It tells you the element type (<f4), whether the data is stored in a particular transposed order, and the dimensions. After the header come 24 elements of 4 bytes each — 96 bytes of nothing but numbers.

Parsing this needs no libraries. Find the magic string, read the length, pull out three values from the text, then interpret the remaining bytes. About three hundred lines with tests. Here is what made up most of them.

#Endianness, or which end of a number comes first

'<f4' means "little-endian 4-byte float". The f4 part is the type. The < is the interesting character.

A 32-bit number occupies four bytes, and there are two conventions for their order. Little-endian puts the least significant byte first; big-endian puts the most significant byte first. The number 1 is 01 00 00 00 in one convention and 00 00 00 01 in the other. Get it backwards and you do not get an error — you get a wildly different number. Reading little-endian 1.0 as big-endian gives you roughly $4.6 \times 10^{-41}$.

Essentially every machine you will encounter is little-endian, so this rarely comes up. The reader still handles it, because a fixture read with the wrong byte order produces confidently wrong numbers rather than a failure, and that is precisely the category of bug worth spending twenty lines to make impossible.

The prefix can be one of four characters, and Ferro's reader treats each one deliberately:

Single-byte types get | because a single byte has no internal order to argue about. Big-endian multi-byte data is refused outright rather than silently misread.

#A short story about the linter being wrong

Rust has an excellent linter called Clippy, which suggests improvements. It looked at the endianness check, which was originally written as:

let order_ok = match order {
    "<" | "|" => true,
    "=" => cfg!(target_endian = "little"),
    _ => false,
};

and suggested collapsing it to a one-liner:

let order_ok = matches!(order, "<" | "|" | "=");

Clippy's reasoning is sound as far as it goes. cfg!(target_endian = "little") is resolved at compile time, so on a little-endian machine that middle arm becomes the literal true, and a match where every arm returns a literal boolean really is better written as matches!.

But the rewrite is wrong. It hardcodes the answer that happens to be correct on the machine doing the compiling. Compile the same code for a big-endian target and the original returns false for = while the suggestion returns true — so the "improved" version would accept native-order data on a machine where native order is exactly the thing being guarded against.

The fix was to restructure so that both the linter and the logic are satisfied, with a comment explaining why it is not written the obvious way:

// Written as a boolean expression rather than a `match` because on a
// little-endian host every arm folds to a literal and clippy suggests
// collapsing it to `matches!(order, "<" | "|" | "=")` — which would
// silently accept native-order data on a big-endian host.
const NATIVE_IS_LITTLE: bool = cfg!(target_endian = "little");
let order_ok = matches!(order, "<" | "|") || (order == "=" && NATIVE_IS_LITTLE);

#Row-major and column-major

The 'fortran_order': False field is about a genuine ambiguity in how a rectangular grid becomes a flat sequence of bytes.

Take this matrix:

Memory is one-dimensional, so those six numbers have to be laid out in a line. There are two sensible ways. Row-major (also called C order, after the C language) walks across each row: 1 2 3 4 5 6. Column-major (Fortran order, after Fortran) walks down each column: 1 4 2 5 3 6.

Both are correct. Both are in wide use — NumPy defaults to row-major, MATLAB and Fortran and most linear algebra libraries default to column-major. And if you read column-major data as though it were row-major, you get a matrix that is transposed and reshaped into nonsense.

NumPy records which one it used, so the reader can handle it. Ferro's takes the simple path: it always hands back row-major order, reindexing the data if the file was column-major. That way no test ever has to know or care which layout the generator happened to produce. One fixture, smoke/fortran_f64.npy, is deliberately written column-major to prove the reindexing works.

#The shapes that are not shapes

Two edge cases in the shape field consistently break array code, so both got fixtures.

A zero-dimensional array has the shape () — an empty tuple. It is not empty; it holds exactly one element. It is a scalar that has been wrapped in array clothing, and it comes up constantly in real use, because summing an entire array gives you back a zero-dimensional array rather than a bare number. Code that assumes shapes have at least one dimension breaks here.

An empty array can have the shape (0, 3) — zero elements, but rank 2, and a meaningful second dimension. This shows up whenever you filter an array down to nothing. Code that treats "no elements" as "no shape" breaks here, and so does code that divides by the number of elements to compute an average.

Both are one-line additions to the fixture generator. Without them, both would be discovered months later, in the middle of an unrelated debugging session.

#NaN, infinity, and negative zero

Floating-point numbers include several values that are not really numbers, and all of them are legitimate results of legitimate operations.

Infinity is what you get from overflow, or from dividing a positive number by zero. It behaves sensibly under arithmetic: infinity plus one is infinity.

NaN stands for "not a number" and comes from genuinely undefined operations: zero divided by zero, the logarithm of a negative number, infinity minus infinity. NaN has a property that breaks naive code: it is not equal to itself. NaN == NaN is false. That is deliberate and standardized, and it means that comparing two arrays element by element with == will report a difference between a NaN and itself.

For a test harness this matters. If a fixture deliberately covers $\log(-1)$, then the expected value is NaN, the computed value is NaN, and the test must pass. So Ferro's comparison treats two NaNs as equal and requires infinities to match in sign, while still rejecting a NaN compared against a real number:

if actual.is_nan() || expected.is_nan() {
    return actual.is_nan() && expected.is_nan();
}
if actual.is_infinite() || expected.is_infinite() {
    return actual == expected;
}

Negative zero is the last oddity. Floating point has two zeros, 0.0 and -0.0, and they compare as equal while being distinguishable in other ways. The fixture includes both, and the test checks the sign bit specifically, because a round-trip through a file format is exactly the kind of thing that can quietly lose it.

#Where 2⁵³ comes in

The reader converts everything to a 64-bit float for comparison. One comparison path for every type is a real simplification — otherwise every test needs to know the storage type of its fixture.

It is also lossy in one place, and pretending otherwise would be a trap.

A 64-bit float has 53 bits of precision available for the integer part of its value. Every whole number up to $2^{53} = 9{,}007{,}199{,}254{,}740{,}992$ can be represented exactly. Past that, they cannot. The integer 9,007,199,254,740,993 is not representable, and converting it to a float silently gives you 9,007,199,254,740,992 instead.

Ferro will absolutely handle 64-bit integers, and quietly rounding them during a test comparison would mean a test that passes when it should fail. So there is a separate integer accessor, and a fixture holding values on both sides of the boundary, and a test that asserts the lossy path is lossy:

let widened = array.to_f64();
assert_ne!(
    widened[4] as i64, exact[4],
    "if this ever holds, f64 gained precision and this test can go"
);

An assertion that a limitation still exists looks strange. Its job is to make sure nobody later "simplifies" the integer accessor away on the reasonable-sounding grounds that the float path gives the same answers. It does not, and this test says so.

#No unsafe code, on purpose

Rust normally guarantees memory safety, and offers an unsafe escape hatch for the cases where you need to do something the compiler cannot verify. Inside an unsafe block you get the sharp tools, along with the responsibility for not cutting yourself.

There is an obvious temptation here. The file contains bytes; you want floats; those bytes already are floats. Rather than decoding each one, you could reinterpret the whole buffer as an array of floats at zero cost. Real numerical libraries do this.

Ferro's reader does not. It decodes each value individually:

NpyDType::F32 => self.decode(|c| f64::from(f32::from_le_bytes(chunk4(c)))),

This copies, and is measurably slower, and is the right call anyway. The reinterpretation trick requires the buffer to be correctly aligned in memory — a four-byte float generally must start at an address divisible by four — and a buffer read from a file has no such guarantee. Getting it wrong is undefined behavior: it might work on your machine, work in tests, and fail somewhere else.

And the performance cost is irrelevant. This code reads a handful of small files during a test run. Nothing in Ferro's runtime touches it. Trading a speedup no one can perceive for a category of bug that is genuinely difficult to diagnose would be a bad deal even if the speedup were large.

There is a second benefit. Rust has a tool called Miri that detects undefined behavior in unsafe code, and Ferro runs it in CI from the very first commit. A reader with no unsafe in it passes trivially — which means that when hand-written SIMD kernels arrive later and Miri goes red, it is reporting a real new problem rather than something that was always there.

#The result

Three hundred lines, twenty-one tests, and no dependencies beyond the error handling library. The tests build .npy files byte by byte in memory rather than reading from disk, so the parser is verified without any fixture files existing — which matters, because it is the thing every other test in the project depends on.

Next: what happens when something goes wrong, and why assertion failed is not an acceptable answer.