比较C或C++中浮点值的两个和

comparing two sums of floating point values in C or C++

本文关键字：两个 C++ 比较更新时间：2023-10-16

假设您得到了两组根据IEEE754实现的浮点变量，这些变量将被视为根据标准中的公式计算的精确值。所有法律价值都是可能的。集合中变量的数量可以是任何自然数。

从数学意义上讲，什么是比较上述变量所代表的值的精确和的好方法。由于域的性质，该问题可以很容易地表示为将单个和与零进行比较。你们可以忽略NaN或无穷大存在的可能性，因为它和核心问题无关。(这些值可以简单而独立地进行检查，并以适合该问题特定应用的方式进行操作。)

一种天真的方法是简单地求和和和比较，或者求一组的值和减去另一组的数值。

bool compare(const std::vector<float>& lhs, const std::vector<float>& rhs)
{
float lSum = 0.0f;
for (auto value : lhs)
{
lSum += value;
}
float rSum = 0.0f;
for (auto value : rhs)
{
rSum += value;
}
return lSum < rSum;
}

很明显，天真的方法存在问题，正如在其他关于浮点运算的问题中提到的那样。大多数问题都与两个困难有关：

浮点值相加的结果因顺序而异

某些值集的某些加法顺序可能导致中间溢出(计算的中间结果超出了可用数据类型支持的范围)

float small = strtof("0x1.0p-126", NULL);
float big = strtof("0x1.8p126", NULL);
std::cout << std::hexfloat << small + big - big << std::endl;
std::cout << std::hexfloat << (big-2*small) + (big-small) + big - (big+small) - (big+2*small) << std::endl;

此代码将产生0和inf；这说明了排序如何影响结果。希望订购的问题也不是微不足道的。

float prev;
float curr = 0.0f;
do
{
prev = curr;
curr += strtof("0x1.0p-126", NULL);
} while (prev != curr);
std::cout << std::hexfloat << curr << std::endl;

如果有足够的时间实际完成计算，该代码将导致0x1.000000p-102，而不是像天真地预期的那样，导致0x1.fffffep127(建议实际观察将当前初始化更改为`strtof("0x1.ffff000p-103")。)；这说明了加法的中间结果和特定加数之间的比例如何影响结果。

关于获得最佳精度，已经说了很多，例如这个问题。

手头的问题不同之处在于，我们不想最大限度地提高精度，但我们有一个定义明确的函数，需要准确地实现。

虽然对一些人来说，这可能是一种有用的练习，但考虑以下场景：这些值集之间的比较可能是在各种环境中独立对整个数据集执行其他操作的基石。一些系统的同步、完美操作可能取决于这种比较是否得到了良好的定义和决定性的实现，而不考虑加数顺序和是否实现IEEE754的特定架构。

这，或者只是好奇。

在讨论中，Kahan求和算法被认为是相关的。然而，该算法是最小化误差的合理尝试。它既不保证结果的正确符号，也不独立于运算的顺序(至少保证集合排列的结果是一致的，如果错误的话)。

最明显的解决方案之一是使用/实现定点算术，使用足够数量的比特来精确地表示每个可能的操作数值，并保持精确的中间结果。

然而，也许这可以只使用浮点运算，以保证结果的正确符号。如果是这样，则需要在解决方案中解决溢出问题(如上面的一个示例所示)，因为这个问题具有特定的技术方面。

(以下是原始问题。)

我有两组多个浮点值(浮点或双精度)。我想为这个问题提供一个完美的答案，哪个集合的总和更大。由于浮点运算中存在伪影，在某些情况下，根据运算顺序的不同，朴素方法的结果可能是错误的。更不用说简单的求和会导致溢出。我无法为自己提供任何努力，因为我只有模糊的想法，所有这些想法都很复杂，没有说服力

一种可能的方法是使用超累积器计算和：这是一种计算浮点数精确和的算法。尽管这些想法已经存在了一段时间，但这个词相对较新。

在某种意义上，你可以把它看作是Kahan求和的扩展，其中顺序求和存储为一个值数组，而不仅仅是一对。然后，主要的挑战变成了弄清楚如何在各种值之间分配精度。

一些相关论文和代码：

Y。朱和海耶斯。"算法908：浮点流的在线精确求和">ACM数学软件汇刊(ACM TOMS)，37(3)：37:1-37:132010年9月。doi:10.145/1824801.1824815
- 不幸的是，论文和代码都在付费墙后面，但这似乎是C++代码
R。M.Neal，"使用小型和大型超累加器的快速精确求和"。2015.arXiv:11505.05571
- 可用的C代码
M。T.Goodrich，A.Eldawy，"浮点数求和的并行算法"。2016年。arXiv:1605.05436
- 这个和上面的Java代码

Post最初也是一个C语言，因此我的代码适用于此
我现在看到的帖子只是C++，但我在下面的文章中很少看到不适用于C++的内容。

简化为查找FP数列表之和的符号

比较两组数字就像把第二组的否定加在第一组上，然后找到联合列表的和的符号。该符号映射到2个原始集合的>、==或<。

仅执行精确的FP数学

假设：FP使用类似IEEE的数字，包括次法线，基数为2，并且对于某些操作是精确的：

具有相同二进制指数和不同符号的a +b的相加。
从0.5 <= |x| < 1.0范围内的数字减去相同符号0.5。
ldexp*()(将数字分解为有效部分和指数部分)函数返回一个精确的值。

按指数形成数组

形成一个和sums[]的数组，其值将仅为(0 or 0.5 <= |sums[i]| < 1.0)，每个可能的指数和一些大于最大值的指数一个。需要较大的和来累积超过FP_MAX的|total_sum|。这需要最多log2(SIZE_MAX)个元素。

将数字集添加到sums[]

对于数字集的每个元素，按照其二进制指数将其添加到相应的sums[]。这是关键，因为可以将相同符号和不同符号的FP数与公共FP二进制指数相加。加法可能导致具有相同符号值的进位和具有不同符号值的取消-这是可以处理的。传入的一组数字不需要排序。

归一化sum[]

对于ones[]上的每个元素，确保减少0.5、0.0或-0.5以外的任何值，将剩余部分添加到较小的ones[]。

检查sum[]的最高有效数字

最高有效(非零)one[s]是结果的符号。

以下代码使用float作为集合的FP类型执行任务。一些并行计算是使用double来检查健全性的，但对float的计算没有贡献。

最后的规范化步骤通常迭代两次。即使是最坏的情况集，我怀疑也会迭代float符号的二进制宽度，大约23次。

解决方案看起来大约是O(n)，但确实使用了一个与FP指数范围大小差不多的数组。

#include <assert.h>
#include <stdbool.h>
#include <float.h>
#include <stdio.h>
#include <time.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#if RAND_MAX/2 >= 0x7FFFFFFFFFFFFFFF
#define LOOP_COUNT 1
#elif RAND_MAX/2 >= 0x7FFFFFFF
#define LOOP_COUNT 2
#elif RAND_MAX/2 >= 0x1FFFFFF
#define LOOP_COUNT 3
#elif RAND_MAX/2 >= 0xFFFF
#define LOOP_COUNT 4
#else
#define LOOP_COUNT 5
#endif
uint64_t rand_uint64(void) {
uint64_t r = 0;
for (int i = LOOP_COUNT; i > 0; i--) {
r = r * (RAND_MAX + (uint64_t) 1u) + ((unsigned) rand());
}
return r;
}
typedef float fp1;
typedef double fp2;
fp1 rand_fp1(void) {
union {
fp1 f;
uint64_t u64;
} u;
do {
u.u64 = rand_uint64();
} while (!isfinite(u.f));
return u.f;
}
int pre = DBL_DECIMAL_DIG - 1;

void exact_add(fp1 *sums, fp1 x, int expo);
// Add x to sums[expo]
// 0.5 <= |x| < 1
// both same sign.
void exact_fract_add(fp1 *sums, fp1 x, int expo) {
assert(fabsf(x) >= 0.5 && fabsf(x) < 1.0);
assert(fabsf(sums[expo]) >= 0.5 && fabsf(sums[expo]) < 1.0);
assert((sums[expo] > 0.0) == ( x > 0.0));
fp1 half = x > 0.0 ? 0.5 : -0.5;
fp1 sum = (sums[expo] - half) + (x - half);
if (fabsf(sum) >= 0.5) {
assert(fabsf(sums[expo]) < 1.0);
sums[expo] = sum;
} else  {
sums[expo] = 0.0;
if (sum) exact_add(sums, sum, expo);
}
exact_add(sums, half, expo+1);  // carry
}
// Add  x to sums[expo]
// 0.5 <= |x| < 1
// differing sign
void exact_fract_sub(fp1 *sums, fp1 x, int expo) {
if(!(fabsf(x) >= 0.5 && fabsf(x) < 1.0)) {
printf("%d %en", __LINE__, x);
exit(-1);
}
assert(fabsf(x) >= 0.5 && fabsf(x) < 1.0);
assert((sums[expo] > 0.0) != ( x > 0.0));
fp1 dif = sums[expo] + x;
sums[expo] = 0.0;
exact_add(sums, dif, expo);
}
// Add x to sums[]
void exact_add(fp1 *sums, fp1 x, int expo) {
if (x == 0) return;
assert (x >= -FLT_MAX && x <= FLT_MAX);
//while (fabsf(x) >= 1.0) { x /= 2.0; expo++; }
while (fabsf(x) < 0.5) { x *= (fp1)2.0; expo--; }
assert(fabsf(x) >= 0.5 && fabsf(x) < 1.0);
if (sums[expo] == 0.0) {
sums[expo] = x;
return;
}
if(!(fabsf(sums[expo]) >= 0.5 && fabsf(sums[expo]) < 1.0)) {
printf("%en", sums[expo]);
printf("%d %en", expo, x);
exit(-1);
}
assert(fabsf(sums[expo]) >= 0.5 && fabsf(sums[expo]) < 1.0);
if ((sums[expo] > 0.0) == (x > 0.0)) {
exact_fract_add(sums, x, expo);
} else {
exact_fract_sub(sums, x, expo);
}
}
void exact_add_general(fp1 *sums, fp1 x) {
if (x == 0) return;
assert (x >= -FLT_MAX && x <= FLT_MAX);
int expo;
x = frexpf(x, &expo);
exact_add(sums, x, expo);
}
void sum_of_sums(const char *s, const fp1 *sums, int expo_min, int expo_max) {
fp1 sum1 = 0.0;
fp2 sum2 = 0.0;
int step = expo_max >= expo_min ? 1 : -1;
for (int expo = expo_min; expo/step <= expo_max/step; expo += step) {
sum1 += ldexpf(sums[expo], expo);
sum2 += ldexp(sums[expo], expo);
}
printf("%-20s = %+.*e %+.*en", s, pre, sum2, pre, sum1);
}

int test_sum(size_t N) {
fp1 a[N];
fp1 sum1 = 0.0;
fp2 sum2 = 0.0;
for (size_t i = 0; i < N; i++) {
a[i] = (fp1) rand_fp1();
sum1 += a[i];
sum2 += a[i];
}
printf("%-20s = %+.*e %+.*en", "initial  sums", pre, sum2, pre, sum1);
int expo_min;
int expo_max;
frexpf(FLT_TRUE_MIN, &expo_min);
frexpf(FLT_MAX, &expo_max);
size_t ln2_size = SIZE_MAX;
while (ln2_size > 0) {
ln2_size >>= 1;
expo_max++;
};
fp1 sum_memory[expo_max - expo_min + 1];
memset(sum_memory, 0, sizeof sum_memory);  // set to 0.0 cheat
fp1 *sums = &sum_memory[-expo_min];
for (size_t i = 0; i<N; i++)  {
exact_add_general(sums, a[i]);
}
sum_of_sums("post add  sums", sums, expo_min,  expo_max);
// normalize
int done;
do {
done = 1;
for (int expo = expo_max; expo >= expo_min; expo--) {
fp1 x = sums[expo];
if ((x < -0.5) || (x > 0.5)) {
//printf("xxx %4d %+.*e ", expo, 2, x);
done = 0;
if (x > 0.0) {
sums[expo] = 0.5;
exact_add(sums, x - (fp1)0.5, expo);
} else {
sums[expo] = -0.5;
exact_add(sums, x - -(fp1)0.5, expo);
}
}
}
sum_of_sums("end  sums", sums, expo_min,  expo_max);
} while (!done);
for (int expo = expo_max; expo >= expo_min; expo--) {
if (sums[expo]) {
return (sums[expo] > 0.5) ? 1 : -1;
}
}
return 0;
}
#define ITERATIONS 10000
#define MAX_NUMBERS_PER_SET 10000
int main() {
unsigned seed = (unsigned) time(NULL);
seed = 0;
printf("seed = %un", seed);
srand(seed);
for (unsigned i = 0; i < ITERATIONS; i++) {
int cmp = test_sum((size_t)rand() % MAX_NUMBERS_PER_SET + 1);
printf("Compare %dnn", cmp);
if (cmp == 0) break;
}
printf("Success");
return EXIT_SUCCESS;
}

无穷大和NaN也可以处理，在一定程度上，留待以后处理。

由2个浮点数求和得到的浮点数只是近似值。给定i₁和i₂求和，我们可以通过以下操作找到浮点求和中误差的近似值：

i₁+i₂=i₁₂
i₁₂-i-₂=i_~1//em>
1i-_~1Δ

对于n数的求和，我们能得出的最接近的近似是计算n-1加法运算的误差，然后再次取n-2的误差求和。您将重复此过程n-2次，或者直到所有错误都变为0.0

可以做几件事来将误差计算速度提高到0.0：

使用较大的浮点类型，例如long double
在求和之前对列表进行排序，以便将小数字添加到小数字，将大数字添加到大数字

现在您可以评估准确性对您的重要性。我要告诉你，在一般情况下，考虑到你得到的结果仍然是一个近似值，上述运算的计算费用是惊人的。

普遍接受的解决方案是卡汉的总结——这是速度和精度之间的幸福结合。Kahan不会将误差保持到求和的末尾，而是将其滚动到每次加法中，防止其值在最高精度浮点范围之外升级。假设我们得到vector<long double> i1，我们可以对其进行Kahan的求和，如下所示：

auto c = 0.0L;
const auto sum = accumulate(next(cbegin(i1)), cend(i1), i1.front(), [&](const auto& sum, const auto& input) {
const auto y = input - c;
const auto t = sum + y;
c = t - sum - y;
return t;
} ) - c;

确定地执行此比较的可能性之一是创建一个不动点算术类，该类的精度等于所使用的类型，并且不限制绝对值。

它可以是实现以下公共方法的类：

FixedPoint(double d);
~FixedPoint();
FixedPoint operator+(const FixedPoint& rhs);
FixedPoint operator-(const FixedPoint& rhs);
bool isPositive();

(每个支持的浮点类型都需要单独的构造函数。)

根据具体情况，实施将需要一个固定的、根据结构或动态大小决定的bool集合；可能是std::bitset、vector<bool>或静态或动态bool阵列。

为了便于实现，我建议实现2的补码编码。

这是一个明显且成本高昂的解决方案，如果这种比较是某些系统的核心，则会损害性能。希望有更好的解决方案