奇怪uint32_t到float数组的转换
Strange uint32_t to float array conversion
我有以下代码片段：
#include <cstdio>
#include <cstdint>
static const size_t ARR_SIZE = 129;
int main()
{
uint32_t value = 2570980487;
uint32_t arr[ARR_SIZE];
for (int x = 0; x < ARR_SIZE; ++x)
arr[x] = value;
float arr_dst[ARR_SIZE];
for (int x = 0; x < ARR_SIZE; ++x)
{
arr_dst[x] = static_cast<float>(arr[x]);
}
printf("%s\n", arr_dst[ARR_SIZE - 1] == arr_dst[ARR_SIZE - 2] ? "OK" : "WTF??!!");
printf("magic = %0.10f\n", arr_dst[ARR_SIZE - 2]);
printf("magic = %0.10f\n", arr_dst[ARR_SIZE - 1]);
return 0;
}
如果我在 MS Visual Studio 2015 下编译它，我可以看到输出是：
WTF??!!
magic = 2570980352.0000000000
magic = 2570980608.0000000000
所以最后一个arr_dst元素和前面的不一样，但是这两个值是通过相同的值转换得到的，填充了arr数组！
这是一个错误吗？
我注意到如果我按以下方式修改转换循环，我会得到 "OK" 结果：
for (int x = 0; x < ARR_SIZE; ++x)
{
if (x == 0)
x = 0;
arr_dst[x] = static_cast<float>(arr[x]);
}
所以这可能是矢量化优化的一些问题。
此行为不会在 gcc 4.8 上重现。有什么想法吗？
MSVC++ 使用的 32 位 IEEE-754 二进制浮点数仅提供 6-7 位十进制数字的精度。您的起始值完全在该类型的 范围 之内，但它似乎不能完全由该类型表示，大多数 uint32_t 类型的值确实如此。
同时，x86 或 x86_64 处理器的浮点单元使用比 MSVC++ 的 64 位 double 更广泛的表示。似乎在循环退出后，最后计算的数组元素以其扩展精度形式保留在 FPU 寄存器中。然后程序可以直接从寄存器中使用该值，而不是从内存中读回它，它必须对先前的元素执行此操作。
如果程序通过将较窄的表示提升为较宽的表示而不是相反的方式来执行 == 比较，那么这两个值可能确实比较不相等，因为从扩展精度到 float 和返回失去精度。在任何情况下，当传递给 printf() 时，这两个值都将转换为类型 double；如果它们确实比较不相等，那么这些转换的结果很可能也不同。
我没有使用 MSVC++ 编译选项，但很可能有一个可以消除这种行为。此类选项有时会以 "strict math" 或 "strict fp" 等名称命名。但是请注意，在 FP-heavy 程序中打开这样的选项（或关闭它的相反选项）可能会非常昂贵。
我对 PowerPC 实现 (Freescale MCP7450) 进行了调查，因为恕我直言，它们的记录比英特尔提出的任何 voodoo 都要好得多。
事实证明，浮点单元、FPU 和向量单元对于浮点运算可能有不同的舍入。 FPU 可以配置为使用四种舍入模式之一；四舍五入到最接近的（默认），截断，朝向正无穷大和朝向负无穷大。然而，矢量单元只能四舍五入到最接近的值，少数 select 指令具有特定的四舍五入规则。 FPU 的内部精度为 106 位。矢量单元满足 IEEE-754 但文档没有说明更多。
查看您的结果，转换 2570980608 更接近原始整数，表明 FPU 具有比向量单元或不同舍入模式更好的内部精度。
unsigned 和 float 之间的转换在 x86 上并不简单；没有针对它的单一指令（直到 AVX512）。一种常见的技术是转换为带符号的，然后修复结果。有多种方法可以做到这一点。 （参见 ，并非所有结果都完美无缺。）
MSVC 使用一种策略对前 128 个元素进行矢量化，然后对最后一个标量元素使用不同的策略（不会进行矢量化），这涉及转换为 double，然后从 double至 float.
gcc 和 clang 从它们的向量化和标量方法中产生 2570980608.0 结果。 2570980608 - 2570980487 = 121 和 2570980487 - 2570980352 = 135（inputs/outputs 没有四舍五入），因此 gcc 和 clang 在这种情况下会产生正确的四舍五入结果（误差小于 0.5ulp）。 IDK 如果这对每个可能的 uint32_t 都是正确的（但只有 2^32 个，we could exhaustively check）。 MSVC 的矢量化循环的最终结果有略多于 0.5ulp 的误差，但标量方法为此输入正确舍入。
IEEE 数学要求 + - * / 和 sqrt 产生正确的舍入结果（误差小于 0.5ulp），但其他函数（比如 log）没有这么严格的要求。 IDK 对 int->float 转换的舍入有什么要求，所以 IDK 如果 MSVC 所做的是严格合法的（如果你没有使用 /fp:fast 或任何东西）。
另见 Bruce Dawson 的 Floating-Point Determinism blog post（他关于 FP 数学的优秀系列的一部分），尽管他没有提到整数<->FP 转换。
我们可以在 OP 链接的 asm 中看到 MSVC 做了什么（简化为只有有趣的指令并手动注释）:
; Function compile flags: /Ogtp
# assembler macro constants
_arr_dst$ = -1040 ; size = 516
_arr$ = -520 ; size = 516
_main PROC ; COMDAT
00013 mov edx, 129
00018 mov eax, -1723986809 ; this is your unsigned 2570980487
0001d mov ecx, edx
00023 lea edi, DWORD PTR _arr$[esp+1088] ; edi=arr
0002a rep stosd ; memset in chunks of 4B
# arr[0..128] = 2570980487 at this point
0002c xor ecx, ecx ; i = 0
# xmm2 = 0.0 in each element (i.e. all-zero)
# xmm3 = __xmm@4f8000004f8000004f8000004f800000 (a constant repeated in each of 4 float elements)
####### The vectorized unsigned->float conversion strategy:
$LL7@main: ; do{
00030 movups xmm0, XMMWORD PTR _arr$[esp+ecx*4+1088] ; load 4 uint32_t
00038 cvtdq2ps xmm1, xmm0 ; SIGNED int to Single-precision float
0003b movaps xmm0, xmm1
0003e cmpltps xmm0, xmm2 ; xmm0 = (xmm0 < 0.0)
00042 andps xmm0, xmm3 ; mask the magic constant
00045 addps xmm0, xmm1 ; x += (x<0.0) ? magic_constant : 0.0f;
# There's no instruction for converting from unsigned to float, so compilers use inconvenient techniques like this to correct the result of converting as signed.
00048 movups XMMWORD PTR _arr_dst$[esp+ecx*4+1088], xmm0 ; store 4 floats to arr_dst
; and repeat the same thing again, with addresses that are 16B higher (+1104)
; i.e. this loop is unrolled by two
0006a add ecx, 8 ; i+=8 (two vectors of 4 elements)
0006d cmp ecx, 128
00073 jb SHORT $LL7@main ; }while(i<128)
#### End of vectorized loop
# and then IDK what MSVC smoking; both these values are known at compile time. Is /Ogtp not full optimization?
# I don't see a branch target that would let execution reach this code
# other than by falling out of the loop that ends with ecx=128
00075 cmp ecx, edx
00077 jae $LN21@main ; if(i>=129): always false
0007d sub edx, ecx ; edx = 129-128 = 1
...一些更荒谬的编译时已知跳跃稍后...
######## The scalar unsigned->float conversion strategy for the last element
$LC15@main:
00140 mov eax, DWORD PTR _arr$[esp+ecx*4+1088]
00147 movd xmm0, eax
# eax = xmm0[0] = arr[128]
0014b cvtdq2pd xmm0, xmm0 ; convert the last element TO DOUBLE
0014f shr eax, 31 ; shift the sign bit to bit 1, so eax = 0 or 1
; then eax indexes a 16B constant, selecting either 0 or 0x41f0... (as whatever double that represents)
00152 addsd xmm0, QWORD PTR __xmm@41f00000000000000000000000000000[eax*8]
0015b cvtpd2ps xmm0, xmm0 ; double -> float
0015f movss DWORD PTR _arr_dst$[esp+ecx*4+1088], xmm0 ; and store it
00165 inc ecx ; ++i;
00166 cmp ecx, 129 ; } while(i<129)
0016c jb SHORT $LC15@main
# Yes, this is a loop, which always runs exactly once for the last element
相比之下，clang 和 gcc 也不会在编译时优化整个事情，但它们确实意识到它们不需要清理 loop，并在各自的循环之后只做一个标量存储或转换。 （clang 实际上会完全展开所有内容，除非你告诉它不要这样做。）
参见Godbolt compiler explorer上的代码。
gcc 只是将 16b 的上半部分和下半部分分别转换为浮点数，然后将它们与乘以 65536 相加。
Clang 的 unsigned -> float 转换策略很有趣：它根本不使用 cvt 指令。我认为它将无符号整数的两个 16 位半部分直接填充到两个浮点数的尾数中（使用一些技巧来设置指数（按位布尔值和 ADDPS），然后像 gcc 一样将低半部分和高半部分加在一起。
当然，如果您编译为 64 位代码，标量转换只需将 uint32_t 零扩展为 64 位，并将其作为带符号的 int64_t 转换为浮点数。 Signed int64_t 可以表示 uint32_t 的每一个值，x86 可以高效地将 64 位 signed int 转换为 float。但这并没有矢量化。
我有以下代码片段：
#include <cstdio>
#include <cstdint>
static const size_t ARR_SIZE = 129;
int main()
{
uint32_t value = 2570980487;
uint32_t arr[ARR_SIZE];
for (int x = 0; x < ARR_SIZE; ++x)
arr[x] = value;
float arr_dst[ARR_SIZE];
for (int x = 0; x < ARR_SIZE; ++x)
{
arr_dst[x] = static_cast<float>(arr[x]);
}
printf("%s\n", arr_dst[ARR_SIZE - 1] == arr_dst[ARR_SIZE - 2] ? "OK" : "WTF??!!");
printf("magic = %0.10f\n", arr_dst[ARR_SIZE - 2]);
printf("magic = %0.10f\n", arr_dst[ARR_SIZE - 1]);
return 0;
}
如果我在 MS Visual Studio 2015 下编译它，我可以看到输出是：
WTF??!!
magic = 2570980352.0000000000
magic = 2570980608.0000000000
所以最后一个arr_dst元素和前面的不一样，但是这两个值是通过相同的值转换得到的，填充了arr数组！
这是一个错误吗？
我注意到如果我按以下方式修改转换循环，我会得到 "OK" 结果：
for (int x = 0; x < ARR_SIZE; ++x)
{
if (x == 0)
x = 0;
arr_dst[x] = static_cast<float>(arr[x]);
}
所以这可能是矢量化优化的一些问题。
此行为不会在 gcc 4.8 上重现。有什么想法吗？
MSVC++ 使用的 32 位 IEEE-754 二进制浮点数仅提供 6-7 位十进制数字的精度。您的起始值完全在该类型的 范围 之内，但它似乎不能完全由该类型表示，大多数 uint32_t 类型的值确实如此。
同时，x86 或 x86_64 处理器的浮点单元使用比 MSVC++ 的 64 位 double 更广泛的表示。似乎在循环退出后，最后计算的数组元素以其扩展精度形式保留在 FPU 寄存器中。然后程序可以直接从寄存器中使用该值，而不是从内存中读回它，它必须对先前的元素执行此操作。
如果程序通过将较窄的表示提升为较宽的表示而不是相反的方式来执行 == 比较，那么这两个值可能确实比较不相等，因为从扩展精度到 float 和返回失去精度。在任何情况下，当传递给 printf() 时，这两个值都将转换为类型 double；如果它们确实比较不相等，那么这些转换的结果很可能也不同。
我没有使用 MSVC++ 编译选项，但很可能有一个可以消除这种行为。此类选项有时会以 "strict math" 或 "strict fp" 等名称命名。但是请注意，在 FP-heavy 程序中打开这样的选项（或关闭它的相反选项）可能会非常昂贵。
我对 PowerPC 实现 (Freescale MCP7450) 进行了调查，因为恕我直言，它们的记录比英特尔提出的任何 voodoo 都要好得多。
事实证明，浮点单元、FPU 和向量单元对于浮点运算可能有不同的舍入。 FPU 可以配置为使用四种舍入模式之一；四舍五入到最接近的（默认），截断，朝向正无穷大和朝向负无穷大。然而，矢量单元只能四舍五入到最接近的值，少数 select 指令具有特定的四舍五入规则。 FPU 的内部精度为 106 位。矢量单元满足 IEEE-754 但文档没有说明更多。
查看您的结果，转换 2570980608 更接近原始整数，表明 FPU 具有比向量单元或不同舍入模式更好的内部精度。
unsigned 和 float 之间的转换在 x86 上并不简单；没有针对它的单一指令（直到 AVX512）。一种常见的技术是转换为带符号的，然后修复结果。有多种方法可以做到这一点。 （参见
MSVC 使用一种策略对前 128 个元素进行矢量化，然后对最后一个标量元素使用不同的策略（不会进行矢量化），这涉及转换为 double，然后从 double至 float.
gcc 和 clang 从它们的向量化和标量方法中产生 2570980608.0 结果。 2570980608 - 2570980487 = 121 和 2570980487 - 2570980352 = 135（inputs/outputs 没有四舍五入），因此 gcc 和 clang 在这种情况下会产生正确的四舍五入结果（误差小于 0.5ulp）。 IDK 如果这对每个可能的 uint32_t 都是正确的（但只有 2^32 个，we could exhaustively check）。 MSVC 的矢量化循环的最终结果有略多于 0.5ulp 的误差，但标量方法为此输入正确舍入。
IEEE 数学要求 + - * / 和 sqrt 产生正确的舍入结果（误差小于 0.5ulp），但其他函数（比如 log）没有这么严格的要求。 IDK 对 int->float 转换的舍入有什么要求，所以 IDK 如果 MSVC 所做的是严格合法的（如果你没有使用 /fp:fast 或任何东西）。
另见 Bruce Dawson 的 Floating-Point Determinism blog post（他关于 FP 数学的优秀系列的一部分），尽管他没有提到整数<->FP 转换。
我们可以在 OP 链接的 asm 中看到 MSVC 做了什么（简化为只有有趣的指令并手动注释）:
; Function compile flags: /Ogtp
# assembler macro constants
_arr_dst$ = -1040 ; size = 516
_arr$ = -520 ; size = 516
_main PROC ; COMDAT
00013 mov edx, 129
00018 mov eax, -1723986809 ; this is your unsigned 2570980487
0001d mov ecx, edx
00023 lea edi, DWORD PTR _arr$[esp+1088] ; edi=arr
0002a rep stosd ; memset in chunks of 4B
# arr[0..128] = 2570980487 at this point
0002c xor ecx, ecx ; i = 0
# xmm2 = 0.0 in each element (i.e. all-zero)
# xmm3 = __xmm@4f8000004f8000004f8000004f800000 (a constant repeated in each of 4 float elements)
####### The vectorized unsigned->float conversion strategy:
$LL7@main: ; do{
00030 movups xmm0, XMMWORD PTR _arr$[esp+ecx*4+1088] ; load 4 uint32_t
00038 cvtdq2ps xmm1, xmm0 ; SIGNED int to Single-precision float
0003b movaps xmm0, xmm1
0003e cmpltps xmm0, xmm2 ; xmm0 = (xmm0 < 0.0)
00042 andps xmm0, xmm3 ; mask the magic constant
00045 addps xmm0, xmm1 ; x += (x<0.0) ? magic_constant : 0.0f;
# There's no instruction for converting from unsigned to float, so compilers use inconvenient techniques like this to correct the result of converting as signed.
00048 movups XMMWORD PTR _arr_dst$[esp+ecx*4+1088], xmm0 ; store 4 floats to arr_dst
; and repeat the same thing again, with addresses that are 16B higher (+1104)
; i.e. this loop is unrolled by two
0006a add ecx, 8 ; i+=8 (two vectors of 4 elements)
0006d cmp ecx, 128
00073 jb SHORT $LL7@main ; }while(i<128)
#### End of vectorized loop
# and then IDK what MSVC smoking; both these values are known at compile time. Is /Ogtp not full optimization?
# I don't see a branch target that would let execution reach this code
# other than by falling out of the loop that ends with ecx=128
00075 cmp ecx, edx
00077 jae $LN21@main ; if(i>=129): always false
0007d sub edx, ecx ; edx = 129-128 = 1
...一些更荒谬的编译时已知跳跃稍后...
######## The scalar unsigned->float conversion strategy for the last element
$LC15@main:
00140 mov eax, DWORD PTR _arr$[esp+ecx*4+1088]
00147 movd xmm0, eax
# eax = xmm0[0] = arr[128]
0014b cvtdq2pd xmm0, xmm0 ; convert the last element TO DOUBLE
0014f shr eax, 31 ; shift the sign bit to bit 1, so eax = 0 or 1
; then eax indexes a 16B constant, selecting either 0 or 0x41f0... (as whatever double that represents)
00152 addsd xmm0, QWORD PTR __xmm@41f00000000000000000000000000000[eax*8]
0015b cvtpd2ps xmm0, xmm0 ; double -> float
0015f movss DWORD PTR _arr_dst$[esp+ecx*4+1088], xmm0 ; and store it
00165 inc ecx ; ++i;
00166 cmp ecx, 129 ; } while(i<129)
0016c jb SHORT $LC15@main
# Yes, this is a loop, which always runs exactly once for the last element
相比之下，clang 和 gcc 也不会在编译时优化整个事情，但它们确实意识到它们不需要清理 loop，并在各自的循环之后只做一个标量存储或转换。 （clang 实际上会完全展开所有内容，除非你告诉它不要这样做。）
参见Godbolt compiler explorer上的代码。
gcc 只是将 16b 的上半部分和下半部分分别转换为浮点数，然后将它们与乘以 65536 相加。
Clang 的 unsigned -> float 转换策略很有趣：它根本不使用 cvt 指令。我认为它将无符号整数的两个 16 位半部分直接填充到两个浮点数的尾数中（使用一些技巧来设置指数（按位布尔值和 ADDPS），然后像 gcc 一样将低半部分和高半部分加在一起。
当然，如果您编译为 64 位代码，标量转换只需将 uint32_t 零扩展为 64 位，并将其作为带符号的 int64_t 转换为浮点数。 Signed int64_t 可以表示 uint32_t 的每一个值，x86 可以高效地将 64 位 signed int 转换为 float。但这并没有矢量化。