Python和Numpy中自定义NaN浮点数的行为

Question

我想在浮点数的NaN值中加入一些额外的信息。我在Python中使用的是单精度的IEEE 754浮点数（32位浮点数）。那么Python和NumPy是如何处理这些值的呢？

理论

根据IEEE 754-2008标准，如果指数位（23到30位）被设置，并且有效数字位中至少有一位被设置，那么这个数字就被认为不是一个数字（NaN）。所以，如果我们把浮点数转换成32位整数表示，满足以下条件的值都是可以的：

• i & 0x7f800000 == 0x7f800000
• i & 0x007fffff != 0

这样我就有很多选择了。不过，标准似乎还说有效数字的最高位是is_quiet，应该被设置，以避免在计算中出现异常。

实际测试

Python 2.7

为了确认，我做了一些测试，结果很有趣：

import math
import struct

std_nan = struct.unpack("f4", struct.pack("I", 0x7fc00000))[0]
spec_nan = struct.unpack("f4", struct.pack("I", 0x7f800001))[0]
spec2_nan = struct.unpack("f4", struct.pack("I", 0x7fc00001))[0]

print "{:08x}".format(struct.unpack("I", struct.pack("f4", std_nan))[0])
print "{:08x}".format(struct.unpack("I", struct.pack("f4", spec_nan))[0])
print "{:08x}".format(struct.unpack("I", struct.pack("f4", spec2_nan))[0])

这给出的结果是：

7fc00000
7fc00001 <<< should be 7f800001
7fc00001

这些测试似乎表明，有些东西（struct.unpack？）总是会设置is_quiet位。

NumPy

我在NumPy中也做了同样的测试，因为在这里我可以确保转换不会改变任何位：

import numpy as np

intarr = np.array([0x7f800001], dtype='uint32')
f = np.fromstring(intarr.tostring(), dtype='f4')
print np.isnan(f)

这给出的结果是：

RuntimeWarning: invalid value encountered in isnan
[True]

但是如果这个值被替换为0x7fc00001，就不会出现错误。

假设

如果我设置is_quiet并使用剩下的位来存储我自己的信息，Python和NumPy都会很高兴。Python自己处理这个位，而NumPy则依赖于底层语言的实现和/或硬件的浮点实现。

问题

我的假设是否正确？有没有什么官方文档可以证明或反驳？还是说这是一种依赖于平台的情况？

我在这里找到了一些相关内容：如何区分Python中的不同类型的NaN浮点数，但我找不到任何官方说明，关于在Python或NumPy中如何处理携带额外信息的NaN。

Answer 1

经过一段时间的思考和查看源代码后，我觉得我可以回答我自己的问题。我的假设几乎是正确的，但并不是全部。

由于NumPy和Python处理数字的方式差别很大，这个回答分为两个部分。

Python和NumPy中NaN的真实情况

NumPy

这可能与平台有关，但在大多数平台上，NumPy使用的是gcc内置的isnan，这个方法运行得很快。运行时的警告通常来自更底层的硬件。在大多数情况下，NumPy可能会用几种方法来判断一个值是否是NaN，比如用x != x，这在AMD 64平台上有效，但在gcc中，实际上是依赖gcc，它可能使用了一些非常简短的代码来实现这个功能。

所以，理论上我们无法保证NumPy如何处理NaN，但在实际中，在更常见的平台上，它会按照标准来做，因为硬件就是这样工作的。NumPy本身对NaN的类型并不在意。（除了某些NumPy特定的、不被硬件支持的数据类型和平台。）

Python

这里的情况就有趣了。如果平台支持IEEE浮点数（大多数都支持），Python会使用C库来进行浮点运算，因此在大多数情况下几乎直接调用硬件指令。所以在这方面应该和NumPy没有区别。

但是... 在Python中通常没有32位浮点数。Python的浮点对象使用的是C语言的double，这是64位格式。那在这两种格式之间如何转换特殊的NaN呢？为了看看实际情况，下面这段小C代码可以帮助我们：

/* nantest.c - Test floating point nan behaviour with type casts */

#include <stdio.h>
#include <stdint.h>

static uint32_t u1 = 0x7fc00000;
static uint32_t u2 = 0x7f800001;
static uint32_t u3 = 0x7fc00001;

int main(void)
    {
    float f1, f2, f3;
    float f1p, f2p, f3p;
    double d1, d2, d3;
    uint32_t u1p, u2p, u3p;
    uint64_t l1, l2, l3;

    // Convert uint32 -> float
    f1 = *(float *)&u1; f2 = *(float *)&u2; f3 = *(float *)&u3;

    // Convert float -> double (type cast, real conversion)
    d1 = (double)f1; d2 = (double)f2; d3 = (double)f3;

    // Convert the doubles into long ints
    l1 = *(uint64_t *)&d1; l2 = *(uint64_t *)&d2; l3 = *(uint64_t *)&d3;

    // Convert the doubles back to floats
    f1p = (float)d1; f2p = (float)d2; f3p = (float)d3;

    // Convert the floats back to uints
    u1p = *(uint32_t *)&f1p; u2p = *(uint32_t *)&f2p; u3p = *(uint32_t *)&f3p;

    printf("%f (%08x) -> %lf (%016llx) -> %f (%08x)\n", f1, u1, d1, l1, f1p, u1p);
    printf("%f (%08x) -> %lf (%016llx) -> %f (%08x)\n", f2, u2, d2, l2, f2p, u2p);
    printf("%f (%08x) -> %lf (%016llx) -> %f (%08x)\n", f3, u3, d3, l3, f3p, u3p);

    return 0;
    }

这段代码的输出是：

nan (7fc00000) -> nan (7ff8000000000000) -> nan (7fc00000)
nan (7f800001) -> nan (7ff8000020000000) -> nan (7fc00001)
nan (7fc00001) -> nan (7ff8000020000000) -> nan (7fc00001)

从第二行可以明显看出，我们遇到了和Python一样的现象。所以，是转换为double时引入了额外的is_quiet位，这个位在64位版本中紧跟在指数后面。

这听起来有点奇怪，但实际上标准是这样说的（IEEE 754-2008，第6.2.3节）：

将一个安静的NaN从较窄的格式转换为较宽的格式，然后再转换回相同的较窄格式，除了使其规范化外，不应以任何方式改变安静NaN的有效载荷。

这并没有提到信号NaN的传播。不过，这在第6.2.1节中有解释：

对于二进制格式，有效载荷编码在尾随有效数字字段的p - 2个最低有效位中。

这里的p是精度，对于32位浮点数来说是24位。所以，我的错误在于使用了信号NaN作为有效载荷。

总结

我得到了以下几点收获：

• IEEE 754-2008支持并鼓励使用qNaNs（安静的NaN）
• 奇怪的结果是因为我尝试使用sNaNs，类型转换导致is_quiet位被设置
• 在最常见的平台上，NumPy和Python都遵循IEEE 754标准
• 实现很大程度上依赖于底层的C实现，因此保证的很少（甚至有一些Python代码承认在某些平台上NaN的处理并不如预期）
• 处理这个问题的唯一安全方法是自己动手处理有效载荷

然而，有一件事情在Python、NumPy或我遇到的其他任何语言中都没有实现。第5.12.1节：

语言标准应提供将支持格式中的NaN转换为外部字符序列的可选功能，这些序列在基本NaN字符序列后附加一个后缀，以表示NaN的有效载荷（见6.2）。有效载荷后缀的形式和解释由语言定义。语言标准应要求任何此类可选输出序列在将外部字符序列转换为支持格式时被接受为输入。