Python 和 C++ long double 给出不同的结果

Question

我正在尝试调试一段 C++ 代码和另一段返回不同结果的 Python 代码之间的不一致。

一般来说，假设初始浮点值相同，我一直期望两个环境之间的浮点运算能够产生相同的结果。例如，如果我在双精度值之间执行简单的除法，例如：

d1 = np.double(4e2)
d2 = np.double(4e4)
d3 = d1 / d2
------------
double x1 = 4e2;
double x2 = 4e4;
double x3 = x1 / x2;

结果在各个位上都是相同的：

Python:
4.0000000000000000000000000e+02
 00000000 00000000 00000000 00000000 00000000 00000000 01111001 01000000
4.0000000000000000000000000e+04
 00000000 00000000 00000000 00000000 00000000 10001000 11100011 01000000
1.0000000000000000208166817e-02
 01111011 00010100 10101110 01000111 11100001 01111010 10000100 00111111
-----------
Cpp:
4.00000000000000000000000e+02
00000000 00000000 00000000 00000000 00000000 00000000 01111001 01000000 
4.00000000000000000000000e+04
00000000 00000000 00000000 00000000 00000000 10001000 11100011 01000000 
1.00000000000000002081668e-02
01111011 00010100 10101110 01000111 11100001 01111010 10000100 00111111 

但是，分别使用

np.longdouble

和

long double

时似乎不会发生这种情况。初始值不仅包含一些随机初始化的位（对于这两种实现，重新运行会导致输出位的后面部分有点随机设置......为什么？），而且所得除法始终会产生不同的结果：

d1 = np.longdouble(4e2)
d2 = np.longdouble(4e4)
d3 = d1 / d2
-----
4.0000000000000000000000000e+02
 00000000 00000000 00000000 00000000 00000000 00000000 00000000 11001000
 00000111 01000000 00110110 00100011 11111101 01111111 00000000 00000000
4.0000000000000000000000000e+04
 00000000 00000000 00000000 00000000 00000000 00000000 01000000 10011100
 00001110 01000000 00110110 00100011 11111101 01111111 00000000 00000000
1.0000000000000000208166817e-02
 00001010 11010111 10100011 01110000 00111101 00001010 11010111 10100011
 11111000 00111111 00110110 00100011 11111101 01111111 00000000 00000000

而对于C++

long double x1 = 4e2l;
long double x2 = 4e4l;
long double x3 = x1 / x2;
------
4.00000000000000000000000e+02
00000000 00000000 00000000 00000000 00000000 00000000 00000000 11001000 
00000111 01000000 01101010 10111110 11100000 01111111 00000000 00000000 
4.00000000000000000000000e+04
00000000 00000000 00000000 00000000 00000000 00000000 01000000 10011100 
00001110 01000000 01100111 10111110 11100000 01111111 00000000 00000000 
9.99999999999999999979671e-03
00001010 11010111 10100011 01110000 00111101 00001010 11010111 10100011 
11111000 00111111 01111000 10111110 11100000 01111111 00000000 00000000 

就我而言，这会产生不幸的影响，因为对

round

的调用最终导致两个实现的方向不同。所以我的问题是，到底发生了什么？长双特殊在什么情况下会导致两个结果不匹配？

如果相关，以下是我输出相应二进制表示的方法：

def binar(num):
    num = np.array([num])
    s = ""
    ss = []
    count = 0
    for v in num.view(np.int8):
        s += " " + str(np.binary_repr(v, width=8))
        count += 1
        if count == 8:
            ss.append(s)
            s = ""
    if s != "":
        ss.append(s)
    return ss

和

void printBinary(const char* prefix, auto num)
{
    unsigned char* it = reinterpret_cast<unsigned char*>(&num);
    for (std::size_t i = 0; i < sizeof(num); i+=8) {
        std::cout << prefix;
        for (std::size_t j = 0; j < 8; j++)
            std::cout << std::bitset<8>(it[i+j]) << ' ';
        std::cout << '\n';
    }
}

编辑：根据要求，以下是使用

printf("%La")

打印的 Python 和 C++ 的值（对于 Python，我使用 Pybind11 来传递长双精度值）：

Python:
d1 = 0xc.8p+5
d2 = 0x9.c4p+12
d3 = 0xa.3d70a3d70a3d8p-10
----------
C++
x1 = 0xc.8p+5
x2 = 0x9.c4p+12
x3 = 0xa.3d70a3d70a3d70ap-10

Answer 1

检查 numpy 针对您的特定安装使用的指令（见下文）。

在我的特定机器上，我无法重现您的差异。我总是使用 C++（使用 godbolt）、pip 安装的 Numpy 副本（在我的机器上）以及从头编译的 Numpy 版本（在我的机器上）获得“C++”结果。

为了测试 C++，我使用了 godbolt。我在支持最新版本的 clang、gcc 和 icc 时得到了相同的结果。

• 网址：https://godbolt.org/z/bhPo9M43c
• 使用的关键指令：FDIVP（x87指令）

对于Python，我尝试了pip安装版本和自己从头编译。我使用 Cython 在 Jupyter 中运行了以下单元格（因为这对我来说比使用 pybind11 更容易）。

# ---- CELL 1 ----------------
%load_ext Cython

# ---- CELL 2 ----------------
%%cython
# distutils: define_macros=NPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION

import numpy as np
cimport numpy as cnp
import cython
from libc.stdio cimport snprintf

def to_str(cnp.ndarray x):
    cdef char buf[256]
    cdef long double v = (<long double *>x.data)[0]
    snprintf(buf, 256, "%0.20Le    %0.16La\n", v, v)
    return bytes(buf).decode('ascii')
    
# ---- CELL 3 ----------------
d1 = np.longdouble(4e2)
d2 = np.longdouble(4e4)
d3 = d1 / d2
print(to_str(np.asarray(d3)))

输出相当于你的“C++”结果。

9.99999999999999999980e-03    0xa.3d70a3d70a3d70a0p-10

numpy 使用的关键指令是 FDIVRP（另一条 x87 指令，只是处理堆栈与 FDIVP 略有不同）。

objdump -d ~/.local/lib/python3.10/site-packages/numpy/core/_multiarray_umath.cpython-310-x86_64-linux-gnu.so | grep -A50 '<LONGDOUBLE_divide>:'
00000000000e55b0 <LONGDOUBLE_divide>:
...snip...
   e55ee:       de f9                   fdivrp %st,%st(1)
...snip...

尝试在您的特定 numpy 库上运行

objdump

，看看它的编译方式是否不同。