struct (un)packing of half-precision nan floats is non-invertible

[seh-dev]

Bug report

Bug description:

I noticed that chaining struct.unpack() and struct.pack() for IEEE 754 Half Precision floats (e) is non-invertible for nan. E.g.:

import struct

original_bytes = b'\xff\xff'

unpacked_float = struct.unpack('e', original_bytes)[0]  # nan
repacked_bytes = struct.pack('e', unpacked_float)  # b'\x00\xfe'  != b'\xff\xff'

IEEE nans aren't unique, so this isn't that surprising... However I found it curious that the same behavior is not exhibited for float (f) or double (d) format, where every original bit pattern I tested could be recovered from the unpacked nan object.

Is this by design?

Here's a quick pytest script that tests over a broad range of nan/inf/-inf cases for each encoding format.

# /// script
# requires-python = ">=3.11"
# dependencies = ["pytest"]
# ///
import struct
import pytest


# Floating Point Encodings Based on IEEE 754 per https://en.wikipedia.org/wiki/IEEE_754#Basic_and_interchange_formats
# binary 16 (half precision) - 1 bit sign, 5 bit exponent, 11 bit significand
# binary 32 (single precision) - 1 bit sign, 8 bit exponent, 23 bit significand
# binary 64 (double precision) - 1 bit sign, 11 bit exponent, 52 bit significand


MAX_TEST_CASES = 100000  # limit number of bit patterns being sampled so we aren't waiting too long


@pytest.mark.parametrize(["precision_format", "precision", "exponent_bits"], [("f", 32, 8), ("d", 64, 11), ("e", 16, 5)])
@pytest.mark.parametrize("sign_bit", [0, 1])
@pytest.mark.parametrize("endianness", ["little", "big"])
def test_struct_floats(precision_format: str, precision: int, exponent_bits: int, sign_bit: int, endianness: str):
    significand_bits = precision - exponent_bits - 1

    n_tests = min(MAX_TEST_CASES, 2**significand_bits)

    significand_patterns = [significand_bits * "0", significand_bits * "1"] + [
        bin(i + 1)[2:] for i in range(1, 2**significand_bits, 2**significand_bits // n_tests)
    ]

    for i in range(n_tests):
        binary = str(sign_bit) + "1" * exponent_bits + significand_patterns[i]
        if endianness == "big":
            format = ">" + precision_format
        elif endianness == "little":
            format = "<" + precision_format
        else:
            raise NotImplementedError()

        test_bytes = int(binary, base=2).to_bytes(precision // 8, endianness)

        unpacked = struct.unpack(format, test_bytes)
        assert len(unpacked) == 1

        repacked = struct.pack(format, unpacked[0])

        assert (
            repacked == test_bytes
        ), f"struct pack/unpack was not invertible for format {format} with raw value: {test_bytes} -> unpacks to {unpacked[0]}, repacks to {repacked}"

if __name__ == "__main__":
    pytest.main([__file__])

CPython versions tested on:

3.13, 3.11, 3.12

Operating systems tested on:

Linux, Windows

Linked PRs

• gh-130317: fix PyFloat_Pack/Unpack[24] for NaN's with payload #130452
• gh-130317: Fix strict aliasing in PyFloat_Pack8() #133150
• gh-130317: Skip test_pack_unpack_roundtrip_for_nans() on x86 #133155
• gh-130317: fix test_pack_unpack_roundtrip() and add docs #133204

[skirpichev]

It seems you are on IEEE-platform, or unpacking special values will fail for float and double formats. So for those formats, pack/unpack functions work by copying bits.

But not PyFloat_Pack2() and PyFloat_Unpack2(). E.g. the later just ignores all payload in the nan value and maps data to one or another quiet nan:

cpython/Objects/floatobject.c

Lines 2402 to 2405 in 12e1d30

	else {
	/* NaN */
	return sign ? -fabs(Py_NAN) : fabs(Py_NAN);
	}

The PyFloat_Pack2() also ignores all payload from double nan:

cpython/Objects/floatobject.c

Lines 2050 to 2059 in 12e1d30

	else if (isnan(x)) {
	/* There are 2046 distinct half-precision NaNs (1022 signaling and
	1024 quiet), but there are only two quiet NaNs that don't arise by
	quieting a signaling NaN; we get those by setting the topmost bit
	of the fraction field and clearing all other fraction bits. We
	choose the one with the appropriate sign. */
	sign = (copysign(1.0, x) == -1.0);
	e = 0x1f;
	bits = 512;
	}

Is this by design?

Looks as a bug for me.

CC @mdickinson

Edit: assuming doubles are binary64, following patch fix your tests:

diff --git a/Objects/floatobject.c b/Objects/floatobject.c
index 3b72a1e7c3..e473fb72fe 100644
--- a/Objects/floatobject.c
+++ b/Objects/floatobject.c
@@ -2048,14 +2048,16 @@ PyFloat_Pack2(double x, char *data, int le)
         bits = 0;
     }
     else if (isnan(x)) {
-        /* There are 2046 distinct half-precision NaNs (1022 signaling and
-           1024 quiet), but there are only two quiet NaNs that don't arise by
-           quieting a signaling NaN; we get those by setting the topmost bit
-           of the fraction field and clearing all other fraction bits. We
-           choose the one with the appropriate sign. */
         sign = (copysign(1.0, x) == -1.0);
         e = 0x1f;
-        bits = 512;
+
+        uint64_t v;
+
+        memcpy(&v, &x, sizeof(v));
+        bits = v & 0x1ff;
+        if (v & 0x800000000000) {
+            bits += 0x200;
+        }
     }
     else {
         sign = (x < 0.0);
@@ -2401,7 +2403,16 @@ PyFloat_Unpack2(const char *data, int le)
         }
         else {
             /* NaN */
-            return sign ? -fabs(Py_NAN) : fabs(Py_NAN);
+            uint64_t v = ((sign? 0xff00000000000000 : 0x7f00000000000000)
+                          + 0xf0000000000000);
+
+            if (f & 0x200) {
+                v += 0x800000000000;
+                f -= 0x200;
+            }
+            v += f;
+            memcpy(&x, &v, sizeof(v));
+            return x;
         }
     }
 

FYI: #55943. Probably the reason why payload was ignored is that the patch was adapted from numpy sources.

[skirpichev]

@tim-one, does it looks as an issue for you?

[skirpichev]

PR is ready for review: #130452

[vstinner]

Fixed in the main branch (future Python 3.14) by change 6157135. The change is not backported to 3.13 branch since it's a minor issue.

[vstinner]

I reopen the issue, there are failures on x86 (32-bit): #130452 (comment)

[skirpichev]

Hmm, I think that underlying reason is same as for 32-bit Windows.

Here reproducer. I create qNaN (1), add payload (2), change it to sNaN (3), assign to a temporary variable (4) and, finally - trying to return a sNaN from a function (5).

In 64-bit:

$ gcc a.c -Wall -Wpedantic -Werror -lm && ./a.out
+nan = 0 11111111111 1000000000000000000000000000000000000000000000000000
positive quiet nan with payload: 0
+nan = 0 11111111111 1000000000000000001011011101110011001010010111000111
positive quiet nan with payload: 12311111111
+nan = 0 11111111111 0000000000000000001011011101110011001010010111000111
positive signaling nan with payload: 12311111111
+nan = 0 11111111111 0000000000000000001011011101110011001010010111000111
positive signaling nan with payload: 12311111111
+nan = 0 11111111111 0000000000000000001011011101110011001010010111000111
positive signaling nan with payload: 12311111111

In 32-bit:

$ gcc -m32 a.c -Wall -Wpedantic -Werror -lm && ./a.out
+nan = 0 11111111111 1000000000000000000000000000000000000000000000000000
positive quiet nan with payload: 0
+nan = 0 11111111111 1000000000000000001011011101110011001010010111000111
positive quiet nan with payload: 12311111111
+nan = 0 11111111111 0000000000000000001011011101110011001010010111000111
positive signaling nan with payload: 12311111111
+nan = 0 11111111111 1000000000000000001011011101110011001010010111000111
positive quiet nan with payload: 12311111111
+nan = 0 11111111111 1000000000000000001011011101110011001010010111000111
positive quiet nan with payload: 12311111111

It seems, in 32-bit - double loaded to FPU:

$ objdump -d a.out | grep fl
    131b:       dd 00                   fldl   (%eax)
    1469:       dd 80 8c e0 ff ff       fldl   -0x1f74(%eax)
    14d3:       dd 45 e8                fldl   -0x18(%ebp)
    1524:       dd 45 e0                fldl   -0x20(%ebp)

#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <inttypes.h>
const unsigned long long float_layout = (1 << 30) + (1 << 22);
const unsigned long long double_layout = (1ULL << 62) + (1ULL << 51);

void print_bits(const size_t size, void const * const ptr,
                const unsigned long long layout)
{
    unsigned char *b = (unsigned char*) ptr;
    unsigned char byte;
    long long i, j;
    for (i = size - 1; i >= 0; i--) {
        for (j = 7; j >= 0; j--) {
            if (layout & 1ULL<<(8*i+j)) {
                printf(" ");
            }
            byte = (b[i] >> j) & 1;
            printf("%u", byte);
        }
    }
    printf("\n");
}
void print_double(const double *x)
{
    uint64_t v = *((uint64_t *)x);
    uint64_t man = (v - (v & (0xffffffffULL << 52)));
    int sign = (v >> 63) & 1, bexp = ((v & (0x7ffULL << 52)) >> 52);
    int exp = bexp;
    if (exp && exp != 2047) {
        exp -= 1023;
    }
    printf("%+a = ", *x);
    print_bits(8, &v, double_layout);
    if (bexp != 2047) {
        printf(("%s0x%d.%" PRIx64 " p%+d\n"), sign ? "-" : "+",
               (bexp != 0), man, bexp ? exp : -1023);
    }
    else {
        if (man) {
            printf("%s %s nan with payload: %llu\n",
                   sign ? "negative" : "positive",
                   man & (1ULL << 51) ? "quiet" : "signaling",
                   man & ~(1ULL << 51));
        }
        else {
            printf("%sinf\n", sign ? "-" : "+");
        }
    }
}
double test(){
    double snan = NAN;
    print_double(&snan);
    uint64_t *v = ((uint64_t *)&snan);
    *v |= 12311111111;
    print_double(&snan);
    *v -= (1ULL << 51);
    print_double(&snan);
    double v2 = snan;
    print_double(&v2);
    v = ((uint64_t *)&v2);
    if (*v & (1ULL << 51)) {
        *v -= (1ULL << 51);
    }
    return v2;
}
int main(int argc, char* argv[])
{
    double v3 = test();
    print_double(&v3);
    return 0;
}

[skirpichev]

No, please don't close this yet.

I think the best option is just skip testing of just snan's round-trip in 32-bit mode. Unfortunately we can't do here anything else with the current PyFloat_Pack/Unpack API. (Passing by reference - works.)

Maybe it worth documenting?

[vstinner]

Maybe it worth documenting?

Documentating the issue sounds like a good option.

(Passing by reference - works.)

I don't think that it's worth it it to add a new API just for sNaN.

[skirpichev]

Documenting the issue sounds like a good option.

PR is ready: #133204

I don't think that it's worth it it to add a new API just for sNaN.

Sure. Though, maybe it's something we could keep in mind. Using PyObject * for argument and as return value will solve this issue:

unsigned char* PyFloat_Pack(PyObject *x, size_t size, int le);
PyObject* PyFloat_Unpack(const unsigned char *p, size_t size, int le);

Two functions vs 6, we also can support someday IEEE sizes>=128.

Currently, we can workaround this problem in struct.pack/unpack(), at cost of code complexity. But I doubt it worth.

[vstinner]

Currently, we can workaround this problem in struct.pack/unpack(), at cost of code complexity. But I doubt it worth.

Most platforms are now 64-bit and don't seem to be affected by the issue. I don't think that it's worth it to invest time on fixing struct.pack/unpack(). And I would prefer that PyFloat_Pack/Unpack*() would remain consistent with struct.pack/unpack().