How ClickHouse Rebuilds a Vector From Its Bit Planes

TL;DR: A QBit column stores a vector as one bit plane per bit position. Read the top p planes and you do not have a vector — you have p bits per dimension, and you must decide what number those bits stand for. ClickHouse’s answer is to round to the centre of the coarse cell (set the most significant dropped bit), to carve the IEEE zero and infinity cells out of that rule, and — at p = 1 — to throw the arithmetic away entirely and count the differing sign bits with XOR and popcount.

The storage side of a bit-transposed vector column is the easy half. The hard half happens after the read: you asked for eight bit planes, you got eight bits per dimension, and now you have to turn that into a number you can subtract from another number.

The obvious reconstruction is to zero-fill the bits you did not read. That is what a truncating cast does, and it is wrong in a way that quietly destroys the search. This article is about the reconstruction — the rounding rule, the IEEE corner cases it has to dodge, and the special case where the whole arithmetic collapses into a popcount.

The bits come back as a bit pattern

The stored form of a QBit(Float32, N) is 32 FixedString streams, one per bit position of the float. Reading p of them and scattering them back into words gives you, for each dimension, a word with p significant bits and 32 − p zeros at the bottom.

The scatter is the inverse of the transposition, and it is written to be branch-free:

cpp
// src/DataTypes/Serializations/SerializationQBit.h:152
/// Generic kernel that untransposes one bit plane. Defined in the header s.t. it inlines into hot loops.
DECLARE_DEFAULT_CODE(
    template <typename T>
    ALWAYS_INLINE void untransposeBitPlaneImpl(const UInt8 * __restrict src, T * __restrict dst, size_t stride_len, T bit_mask)
    {
        const size_t bytes_per_fs = stride_len / 8;
        ssize_t row_base = stride_len - 1;

        for (size_t b = 0; b < bytes_per_fs; ++b, row_base -= 8)
        {
            const uint8_t v = src[b];

            /// Fast out on common all-zeros case
            if (!v)
                continue;

            for (int i = 0; i < 8; ++i)
            {
                /// Mask is 0...0 if current bit is 0, 1...1 if it is 1. Use it to avoid a branch
                T mask = static_cast<T>(-T((v >> i) & 1));
                dst[row_base - 7 + i] |= (mask & bit_mask);
            }
        }
    })

Two details: the |= lets several planes accumulate into the same destination words (which is exactly what happens when you read eight planes and want an eight-bit value), and -T(bit) turns a bit into an all-ones or all-zeros word without a branch.

The distance function processes 32 rows at a time, untransposing the planes into a contiguous block:

cpp
// src/Functions/array/distanceTransposed.cpp:1146
        for (size_t base_row = 0; base_row < input_rows_count; base_row += block_size)
        {
            const size_t rows_in_block = std::min(block_size, input_rows_count - base_row);
            const size_t words_in_block = rows_in_block * padded_array_size;

            memset(block.data(), 0, words_in_block * sizeof(CalcT));

            /// Untranspose, for each stride group, its `precision` bit planes into that group's slice of every row of the block
            for (size_t group = 0; group < num_groups; ++group)
            {
                for (size_t bit = 0; bit < precision; ++bit)
                {
                    const auto & col = assert_cast<const ColumnFixedString &>(*planes[group * precision + bit].column);
                    Word bit_mask = static_cast<Word>(Word(1) << (sizeof(Word) * 8 - 1 - bit));

                    for (size_t r = 0; r < rows_in_block; ++r)
                    {
                        const UInt8 * src = reinterpret_cast<const UInt8 *>(col.getChars().data()) + (base_row + r) * bytes_per_group;
                        untranspose_kernel(
                            src, reinterpret_cast<Word *>(block_row(r) + group * padded_group), padded_group, bit_mask);
                    }
                }
            }

At this point the block holds precision significant bits per element. It is not yet a vector of floats. It is a vector of truncated floats.

The centre-fill

Zero-filling the dropped bits reconstructs the lower edge of the coarse cell that the truncated value belongs to. That is a systematically downward biased estimate, and at low precision it degenerates completely:

cpp
// src/Functions/array/distanceTransposed.cpp:1121
        /// A value truncated to `precision` bit planes is rounded to the centre of its coarse cell by setting the most significant
        /// dropped bit. Zero-filling the dropped bits would instead reconstruct the cell's lower edge, which biases every value
        /// towards zero and degenerates at low precision: e.g. for BFloat16 at precision 1 only the sign bit survives, so every
        /// value would be reconstructed as +-0.0 and the distance would be the same for every row. When `precision` covers the
        /// whole word, the fill is zero; the dispatch in
        /// `executeImpl` downcasts to a narrower CalcT only while `precision` is strictly below the narrow width, so this happens
        /// exactly when `precision` covers the whole *element* and no bits are dropped at all. `centre_fill` is a single bit strictly
        /// below the `precision` kept planes.
        const Word centre_fill = precision < sizeof(Word) * 8 ? static_cast<Word>(Word(1) << (sizeof(Word) * 8 - 1 - precision)) : Word(0);

One bit. centre_fill is the most significant dropped bit — the bit immediately below the kept planes — and OR-ing it into every reconstructed word moves the value from the bottom of its cell to the middle.

The degenerate case is spelled out: for BFloat16 at precision 1, only the sign bit survives, and zero-filling would reconstruct every element as ±0.0. Every distance would be zero. With the centre-fill, the sign bit survives and the next bit is set, so the element becomes ±2.0 — the top exponent bit — and the sign is all that is left, which is exactly what one bit of precision is supposed to give you.

The centre-fill is applied in one of three shapes, depending on the element type and the precision:

cpp
// src/Functions/array/distanceTransposed.cpp:1206
            constexpr bool is_int8 = std::is_same_v<CalcT, Int8>;
            /// BFloat16 and Float32 both carry an 8-bit exponent; Float64 carries 11. Unused for Int8.
            constexpr size_t exponent_bits = std::is_same_v<CalcT, Float64> ? 11 : 8;

            const bool apply_centre = is_int8 || precision == 1 || precision > exponent_bits;
            const bool collapse_zero_cell = !is_int8 && precision > exponent_bits;

Read it as three cases:

Int8 — the code is a linear integer, so its centre is the arithmetic midpoint of the code range. Applied unconditionally.
precision == 1 — sign-only truncation. The single dropped bit that gets set is the top exponent bit, which reconstructs a bounded ±2.0. Every value keeps just its sign.
precision > exponent_bits — the whole exponent survives and only mantissa bits are dropped, so the centre is a meaningful midpoint within the value’s own binade. This is the only case where the “zero cell” and “infinity cell” carve-outs are needed.

And there is a fourth case that is deliberately not centred: when precision truncates exponent bits, the most significant dropped bit is an exponent bit, and setting it is a multiplicative jump across many binades — not an approximation in value space at all:

cpp
// src/Functions/array/distanceTransposed.cpp:1201
            /// For a float at `2 <= precision <= exponent_bits` the most significant dropped bit is an *exponent* bit, so setting
            /// it is a multiplicative jump across many binades - not a usable approximation in value space, and (once squared in
            /// the kernel) architecture-sensitive. There the bounded lower edge of the coarse exponent cell is kept instead (no
            /// centre), as before this reconstruction change: a smaller `precision` then trades accuracy for speed without blowing
            /// magnitudes up by orders of magnitude, as the function contract requires.

So the rule is not “always centre” — it is “centre when the centre is a bounded, value-space-meaningful approximation, and leave the lower edge otherwise”. That distinction is the difference between a quantisation scheme that is approximately the original vector and one that multiplies it by 2^k.

The IEEE carve-outs

OR-ing a bit into a float word is not a rounding operation — it is a bit operation, and bit operations do not know about IEEE. Two categories of float must be left alone.

The zero cell. +0.0 is an all-zero word and -0.0 is a word with only the sign bit set. If you OR the centre-fill into them, +0.0 becomes a tiny positive subnormal and -0.0 becomes a tiny negative subnormal. The function must detect them and leave them alone:

cpp
// src/Functions/array/distanceTransposed.cpp:1232
                else if (collapse_zero_cell)
                {
                    /// Two IEEE categories are carved out of the midpoint rule (the quantized LUT likewise only reconstructs
                    /// finite cells), both detected after masking off the sign bit:
                    ///  - the zero cell: `+0.0` (an all-zero word) and `-0.0` (only the sign bit kept) must stay zero rather than
                    ///    gain a spurious tiny magnitude from `centre_fill`, which would otherwise turn a stored `-0.0` into a
                    ///    non-zero negative subnormal;
                    ///  - the non-finite cell (all exponent bits set): `+-inf` has a zero kept mantissa, so OR-ing `centre_fill`
                    ///    would flip it to a `NaN` and change the IEEE category of a legitimate input. Leaving it untouched keeps
                    ///    `+-inf` exactly infinite, and keeps a `NaN` a `NaN` whenever a set mantissa bit survives the truncation;
                    ///    a `NaN` whose payload lies entirely in the dropped bits truncates to the `+-inf` encoding and therefore
                    ///    reconstructs to `+-inf` (see the reconstruction-policy comment above - this cell is inherently ambiguous
                    ///    at reduced precision, and preserving the canonical infinity avoids corrupting a genuine `+-inf`).
                    /// Only the finite, non-zero words in between are centred. `exponent_bits` is only meaningful for a float, so
                    /// this branch (and its `exponent_mask`) is guarded from the `Int8` instantiation by `is_int8` above.
                    constexpr Word non_sign_mask = static_cast<Word>(~(Word(1) << (sizeof(Word) * 8 - 1)));
                    constexpr Word exponent_mask
                        = static_cast<Word>(((Word(1) << exponent_bits) - 1) << (sizeof(Word) * 8 - 1 - exponent_bits));
                    for (size_t r = 0; r < rows_in_block; ++r)
                    {
                        Word * row = words + r * padded_array_size;
                        for (size_t d = 0; d < used_dims; ++d)
                        {
                            const Word magnitude = row[d] & non_sign_mask;
                            if (magnitude != 0 && (magnitude & exponent_mask) != exponent_mask)
                                row[d] |= centre_fill;
                        }
                    }
                }

The detection is two bit tests: mask off the sign bit, and if the rest is zero it is the zero cell; if the exponent bits are all ones it is the non-finite cell. Everything in between gets the centre-fill.

Why it matters. Consider the case the comment calls out: a stored -0.0 that gains a tiny negative magnitude. In a cosine distance, a vector of stored zeros would no longer have a zero dot product with anything, and the distance would stop being 1. Worse, the comment notes that “a reduced-precision cosine distance report identical zero vectors as maximally dissimilar” — a single stored zero can change the ranking of the whole result set.

The infinity cell. ±inf has a zero mantissa, so OR-ing the centre-fill into it sets a mantissa bit and turns it into a NaN. That changes the IEEE category of a legitimate input, which is not a rounding question at all. The rule here is explicitly lossy and says so: a NaN whose payload lives entirely in the dropped mantissa bits truncates to the ±inf encoding and reconstructs as ±inf. Preserving the canonical infinity, the comment says, is preferred over fabricating a NaN, “because doing so would corrupt a genuinely stored ±inf”.

Narrowing the arithmetic to match

Reconstruction and arithmetic width are the same decision. If only p bits survive, the reference vector can be narrowed to a type that holds p bits, and the SIMD kernel can be one that processes twice as many lanes:

cpp
// src/Functions/array/distanceTransposed.cpp:646
        /// We need to find two types: the type of the reference vector and the type of the calculation.
        /// The type of calculation is determined by the value of `precision. For example, if col_x is Float32 and p = 16, we will only have
        /// 16 meaningful bits to calculate the distance. So we can downcast the reference vector to BFloat16 and do calculations faster.
        auto dispatch_by_accum_type = [&]<typename RefT>(auto func)
        {
            /// Int8 has only 8 bits and cannot be downcasted to a narrower type, so calculations are always done in Int8.
            if constexpr (std::is_same_v<RefT, Int8>)
            {
                return func.template operator()<Int8, Int8>();
            }
            else
            {
                /// Downcast only while the centre-fill midpoint (the most significant dropped bit) stays representable in the
                /// narrower calculation word, i.e. `precision` is strictly below the narrow width. At `precision == 16` a Float32
                /// element still drops 16 mantissa bits and its midpoint is the bit just below a BFloat16 word, so that boundary
                /// case must keep the full Float32 width; a BFloat16 element at `precision == 16` drops nothing and stays BFloat16.
                auto calc_type
                    = (precision < 16 || std::is_same_v<RefT, BFloat16>
                           ? TypeToTypeIndex<BFloat16>
                           : (precision <= 32 ? TypeToTypeIndex<Float32> : TypeToTypeIndex<Float64>));

                /// Float64 cannot be downcasted to Float32 or BFloat16 in an easy way by reordering bits. That is why with it we always do
                /// calculations in full width. Alternatively, we could static_cast each element when calculating, but it is slower.
                if (std::is_same_v<RefT, Float64>)
                    return func.template operator()<RefT, Float64>();

The condition is not “precision ≤ 16” but “precision < 16”, and the reason is the centre-fill again: at exactly 16 bits, a Float32 element still drops 16 mantissa bits, and its midpoint is the bit just below a BFloat16 word — which a BFloat16 cannot represent. So the boundary case keeps the full Float32 width. The reconstruction rule and the arithmetic width are coupled, and the code says so.

Precision 1: the sign vector

At one bit of precision, the reconstruction is degenerate in a useful way: only the sign of each element survives. The general path would reconstruct every element as ±2.0 (for a float) or ±64 (for a raw Int8):

cpp
// src/Functions/array/distanceTransposed.cpp:684
            /// A single bit plane keeps only the sign, and the general path would reconstruct every element to the `centre_fill`
            /// of its sign cell: +-64 (0x40) for a raw Int8 and +-2.0 (the top exponent bit) for a float. The one-bit paths
            /// reproduce these magnitudes without untransposing the plane.
            if (precision == 1)
            {
                constexpr Float64 scale = std::is_same_v<T, Int8> ? 64.0 : 2.0;

Once both sides are sign vectors of the same magnitude, the distance depends only on how many dimensions disagree in sign — the Hamming distance. And for the three metrics, the answer is a closed form in the Hamming distance:

cpp
// src/Functions/array/distanceTransposed.cpp:121
    /// Every dimension whose signs differ contributes (2 * scale)^2 to the squared distance, the rest contribute nothing.
    static Float64 fromHammingDistance(size_t hamming, size_t /*dims*/, Float64 scale)
    {
        return 2.0 * scale * std::sqrt(static_cast<Float64>(hamming));
    }

cpp
// src/Functions/array/distanceTransposed.cpp:187
    /// Both sign vectors have the norm scale * sqrt(dims) and their dot product is scale^2 * (dims - 2 * hamming), so the
    /// cosine distance is 2 * hamming / dims whatever the scale.
    static Float64 fromHammingDistance(size_t hamming, size_t dims, Float64 /*scale*/)
    {
        return 2.0 * static_cast<Float64>(hamming) / static_cast<Float64>(dims);
    }

cpp
// src/Functions/array/distanceTransposed.cpp:251
    /// Each dimension with agreeing signs contributes scale^2, each dimension with differing signs contributes -scale^2.
    static Float64 fromHammingDistance(size_t hamming, size_t dims, Float64 scale)
    {
        return scale * scale * (static_cast<Float64>(dims) - 2.0 * static_cast<Float64>(hamming));
    }

The symmetric version: XOR and popcount

The symmetric path (enabled by qbit_one_bit_symmetric_distance) reduces the reference vector to its signs too, and then the whole distance is a popcount over the packed plane bytes:

cpp
// src/Functions/array/distanceTransposed.cpp:964
    /// Symmetric one-bit (`precision == 1`) path, enabled by the setting `qbit_one_bit_symmetric_distance`.
    ///
    /// With a single bit plane only the sign of every stored element is known, and the general path would reconstruct each
    /// element as +-scale (`centre_fill` gives +-2.0 for a float and +-64 for a raw Int8; the Lloyd-Max path gives the +-m prefix
    /// centroid of row 1 of its LUT). The reference vector is reduced to its signs the same way, so both sides are sign vectors
    /// of magnitude `scale` and the distance depends only on the number of dimensions whose signs differ: the Hamming distance
    /// between the stored sign bit plane and the sign bit plane of the reference. It is computed with XOR and popcount directly
    /// on the packed plane bytes, without untransposing anything, and the kernel's distance for every possible Hamming
    /// distance in `[0, used_dims]` is precomputed once into a lookup table (`Kernel::fromHammingDistance`).

The popcount itself is a hand-written loop over 64-bit words:

cpp
// src/Functions/array/distanceTransposed.cpp:57
/// Number of differing bits between two byte strings of `bytes` bytes.
ALWAYS_INLINE size_t hammingDistance(const UInt8 * __restrict a, const UInt8 * __restrict b, size_t bytes)
{
    size_t result = 0;
    size_t i = 0;
    for (; i + sizeof(UInt64) <= bytes; i += sizeof(UInt64))
        result += static_cast<size_t>(std::popcount(unalignedLoad<UInt64>(a + i) ^ unalignedLoad<UInt64>(b + i)));
    for (; i < bytes; ++i)
        result += static_cast<size_t>(std::popcount(static_cast<uint8_t>(a[i] ^ b[i])));
    return result;
}

64 dimensions per popcount, and no untransposing at all — the stored plane and the reference’s sign plane are both already in the packed layout.

The asymmetric version: a 256-entry table

The default one-bit path keeps the reference at full precision, so the distance is not a pure Hamming count. But it is still a function of a single scalar: the signed sum S = Σ sᵢ · yᵢ, where sᵢ is the stored sign and yᵢ the reference element.

cpp
// src/Functions/array/distanceTransposed.cpp:845
    /// Default one-bit (`precision == 1`) path: the reference vector keeps its full precision (asymmetric distance), so the
    /// results are those of the general path for the sign-only reconstruction, computed without untransposing anything.
    ///
    /// Every stored element is known only as a sign s_i and reconstructed as scale * s_i, so all three kernels are functions of
    /// the signed sum S = sum s_i * y_i, of the squared norm of the reference and of the dimension (`Kernel::fromSignedSum`).
    /// S is computed straight from the packed sign plane: for a constant reference, a 256-entry table per plane byte holds the
    /// partial sum of +-y over the 8 dimensions of that byte for every bit pattern, so a row costs one table lookup per byte of
    /// its sign plane instead of an untranspose of the plane into `used_dims` words followed by a distance kernel.

So: for each byte of the sign plane (8 dimensions), and for each of the 256 possible bit patterns, precompute the signed contribution of that byte to S. Then a row costs one table lookup per byte of the plane:

cpp
// src/Functions/array/distanceTransposed.cpp:900
            /// table[(group * bytes_per_group + byte) * 256 + pattern] = sum over the 8 dimensions of that byte of +-y_d, where the
            /// sign is negative for the dimensions whose bit is set in `pattern`.
            PODArray<TableT> table(plane_bytes * 256);
            for (size_t group = 0; group < num_groups; ++group)
            {
                for (size_t byte = 0; byte < bytes_per_group; ++byte)
                {
                    TableT weights[8];
                    for (size_t bit = 0; bit < 8; ++bit)
                    {
                        const auto d = dimension_of(byte, bit);
                        weights[bit] = d ? static_cast<TableT>(y[group * stride + *d]) : TableT(0);
                    }
                    TableT * entries = table.data() + (group * bytes_per_group + byte) * 256;
                    for (size_t pattern = 0; pattern < 256; ++pattern)
                    {
                        TableT sum = 0;
                        for (size_t bit = 0; bit < 8; ++bit)
                            sum += ((pattern >> bit) & 1) ? -weights[bit] : weights[bit];
                        entries[pattern] = sum;
                    }
                }
            }

256 entries × 8 bits of work each = 2048 operations to build the table, once. After that, each row is a sequence of lookups:

cpp
// src/Functions/array/distanceTransposed.cpp:924
            for (size_t row = 0; row < input_rows_count; ++row)
            {
                /// Independent partial sums, so the lookups are not one dependent add chain.
                AccT partial[4]{};
                for (size_t group = 0; group < num_groups; ++group)
                {
                    const UInt8 * stored = stored_plane(group) + row * bytes_per_group;
                    const TableT * entries = table.data() + group * bytes_per_group * 256;
                    for (size_t byte = 0; byte < bytes_per_group; ++byte)
                        partial[byte % 4] += entries[byte * 256 + stored[byte]];
                }
                const Float64 signed_sum = static_cast<Float64>(partial[0] + partial[1] + partial[2] + partial[3]);
                result_data[row] = Kernel::fromSignedSum(signed_sum, y2, used_dims, scale);
            }

Four independent accumulators so the lookups are not one dependent add chain — the same trick used in the dot-product kernel.

The lookup only works because the reference vector is constant, which is the normal case (a query vector). With a per-row reference, the table would have to be rebuilt per row, and the code falls back to adding the ±y terms directly:

cpp
// src/Functions/array/distanceTransposed.cpp:939
        else
        {
            /// One reference vector per row: no table pays off, so add the +-y_d up directly from the bits.
            for (size_t row = 0; row < input_rows_count; ++row)
            {
                const YT * y = ref_values.data() + (row == 0 ? 0 : ref_offsets[row - 1]);
                AccT sum = 0;
                for (size_t group = 0; group < num_groups; ++group)
                {
                    const UInt8 * stored = stored_plane(group) + row * bytes_per_group;
                    for (size_t byte = 0; byte < bytes_per_group; ++byte)
                        for (size_t bit = 0; bit < 8; ++bit)
                            if (const auto d = dimension_of(byte, bit))
                            {
                                const AccT value = static_cast<AccT>(y[group * stride + *d]);
                                sum += ((stored[byte] >> bit) & 1) ? -value : value;
                            }
                }
                result_data[row] = Kernel::fromSignedSum(static_cast<Float64>(sum), squared_norm(y), used_dims, scale);
            }
        }

The quantised sibling: a different reconstruction entirely

There is a second family of QBit columns, QBit(Int8), where the stored codes are not raw floats but quantiser indices — 256 levels optimised for a standard-normal source, which is what a unit-normalised embedding coordinate approximately is:

cpp
// src/Common/LloydMaxQuantizer.h:11
/// 256-level Gaussian Lloyd-Max scalar quantizer (the MSE-optimal scalar quantizer for a standard-normal source).
/// It is intended for values that are approximately N(0, 1) - e.g. a unit-norm embedding coordinate after a random
/// orthogonal (randomized Hadamard) rotation, scaled by sqrt(dimensions) to reach unit variance.
///
/// The code is the index 0..255 of the Lloyd-Max cell, stored as Int8 (index - 128), so the sign bit equals the sign of
/// the value (+0/-0 map to the positive/negative central cells) and the top `b` bits form a valid embedded 2^b-level
/// quantizer: bit-truncation of the code yields coarser Int4/Int2/binary codes.

The last sentence is the key property: because the code is a table index and the table is monotone, truncating the code’s bits yields a valid coarser quantiser. So a QBit(Int8) searched at 4-bit precision is a genuine 16-level quantiser, not a mangled 8-bit one.

But the reconstruction is not a centre-fill. For a quantiser the right estimate of a truncated code is the conditional mean of the cells that the retained prefix covers — the minimum-error choice, and a completely different rule:

cpp
// src/Functions/array/distanceTransposed.cpp:1300
    /// Quantized variant: the QBit holds Int8 Lloyd-Max codes and the reference vector is either the full-precision Float32 query or
    /// a quantized Array(Int8) query.
    /// `planes` holds `num_groups * precision` FixedString bit-plane columns in group-major order (plane[g * precision + b]).
    /// Each stride group is untransposed into a contiguous slice of a `used_dims`-element buffer of raw code bytes, which is then
    /// dequantized to Float32 prefix centroids before the Float32 distance kernel runs. For `p < 8`, each centroid is the
    /// Gaussian conditional mean over the union of fine Lloyd-Max cells represented by the retained prefix; at `p = 8` the
    /// existing `dequantizeInt8ToBFloat16` reconstruction is preserved exactly.
    /// An Array(Int8) reference is a complete (not partially read) quantized query, so it is dequantized at full precision (row 8 of
    /// the LUT, i.e. `dequantizeInt8ToBFloat16`) before the kernel runs; a Float32 reference is used verbatim.

and the LUT that holds those centroids is indexed by precision:

cpp
// src/Functions/array/distanceTransposed.cpp:1350
        /// Reconstruction table for this precision: maps a raw code byte to its Float32 prefix centroid (or the exact
        /// BFloat16 Lloyd-Max reconstruction at precision 8).
        const std::array<Float32, 256> & dequant = LloydMax::transposedDequantLUT()[precision];

One row of a 256 × 9 table per precision level: 256 code values, each mapped to the centroid of the cell the retained prefix covers.

Trade-offs

What the design buys:

The reconstruction is a rounding rule, not a truncation. Centre-filling makes the reduced-precision vector an unbiased estimate of the original instead of a systematically smaller one, which is what makes the reduced precision usable rather than merely fast.
The IEEE categories survive. Zero stays zero, infinity stays infinity, and NaN stays NaN whenever any of its payload survives the truncation. Without that, a stored zero vector could come out as a tiny subnormal and change the ranking of a cosine search.
The arithmetic width follows the precision. Eight bits of precision means BFloat16 arithmetic and a kernel that processes twice as many lanes, with the boundary case (precision == 16) handled explicitly because the centre-fill no longer fits in the narrower word.
Precision 1 is not a slow path. It is a popcount over the packed planes (symmetric) or a 256-entry table lookup per plane byte (asymmetric) — both of which skip the untranspose entirely.
The quantised path uses the statistically right estimator. Conditional means over the retained cells, not a bit trick.

What it costs:

The rule has cases. apply_centre has three disjuncts and collapse_zero_cell exists only to undo part of the rule. Every new element type or precision regime means another case, and the comment block explaining them is longer than the code.
Some information is irretrievably lost. A NaN whose payload lives in the dropped bits comes back as ±inf; the code says so and calls the cell “inherently ambiguous at reduced precision”.
The centre-fill is only valid in value space for some precisions. When the dropped bit is an exponent bit, the centre is a multiplicative jump, so the rule is disabled and the biased lower edge is used instead. So the reconstruction quality is not monotone in precision.
The fast paths need a constant reference. The 256-entry table is built once per query; with a per-row reference vector the code falls back to bit-by-bit accumulation and loses most of the advantage.
The symmetric one-bit distance is a different function. The setting’s own documentation says it “ignores the magnitudes of the reference elements: a dimension where the reference is large counts as much as one where it is close to zero”, which is why it is off by default.

The bargain, in one sentence: a bit-transposed column is only useful if the reader has a principled way to turn “the top eight bits of every dimension” back into numbers, and ClickHouse’s principled way is a rounding rule that knows about IEEE — with a completely different, statistically motivated rule for the quantised variant, and a shortcut that bypasses reconstruction entirely when only one bit is needed.