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

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

在数学意义上，比较精确的由所述变量表示的值的总和的好方法。 由于域的性质，问题可以很容易地表示为将单个总和与零进行比较。 您可以忽略存在NaN或无限的可能性，因为它与核心问题无关。 （可以轻松独立地检查这些值，并以适合此问题的特定应用的方式采取行动。）

一种天真的方法是简单地求和和比较，或者将一组的值和另一组的值相加。

bool compare(const std::vector& lhs, const std::vector& 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 （将curr初始化更改为`strtof（“0x1.fff000p-103”）将是建议实际观察这一点。）; 这说明了添加的中间结果与特定加数之间的比例如何影响结果。

关于获得最佳精度的说法很多，例如。 这个问题 。

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

虽然对某些人来说这可能是有用的想法运动似乎充其量存在争议，但请考虑以下情形：这些值集之间的比较可能是在各种环境中独立对整个数据集执行的其他操作的基石。 一些系统的同步，完美操作可能依赖于这种比较被很好地定义和确定性地实现，而不管加数顺序和实现IEEE754的特定体系结构。

这个，或者只是好奇心。

在讨论中， Kahan求和算法已被提及为相关的。 然而，该算法是最小化误差的合理尝试。 它既不保证结果的正确符号，也不依赖于操作的顺序（至少保证一致的，如果错误的话，对于集合的排列）。

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

但是，这可以通过仅使用浮点运算来保证正确的结果符号。 如果是这样，溢出问题（如上面的一个例子中所示）需要在解决方案中解决，因为这个问题具有特定的技术方面。

（以下是原始问题。）

我有两组多个浮点（浮点或双精度）值。 我想为这个问题提供一个完美的答案。 由于浮点运算中的伪像，在某些极端情况下，天真方法的结果可能是错误的，这取决于操作的顺序。 更不用说简单的总和会导致溢出。 我不能为我提供任何努力，因为我所拥有的只是模糊的想法，所有这些都很复杂而且没有说服力。

Answer 1

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

在某种意义上，您可以将其视为Kahan求和的扩展，其中序列和存储为值数组，而不仅仅是一对。 然后，主要的挑战就是弄清楚如何在各种值之间分配精度。

一些相关的论文和代码：

• YK Zhu和WB Hayes。 “算法908：浮点流的在线精确求和”。 ACM数学软件交易 （ACM TOMS），37（3）：37：1-37：13，2010年9月.doi： 10.1145 / 1824801.1824815

  • 不幸的是，论文和代码背后是付费专区，但这似乎是C ++代码 。

• RM Neal，“使用小型和大型超级累积器的快速精确求和”。 2015. arXiv： 1505.05571

  • C代码可用

• MT Goodrich，A。Eldawy“用于求和浮点数的并行算法”。 2016. arXiv： 1605.05436

  • 这个和上面的Java代码

Answer 2

Post最初也是一个C，因此我的代码适用于此。
我现在看到post只是C ++，但我在下面看不到很容易适用于C ++。

简化为找到FP编号列表总和的符号

比较2组数字就像是将第二组的否定附加到第一组，然后找到联合列表总和的符号。 此符号映射到2个原始集的> ， ==或< 。

仅执行精确的FP数学运算

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

1. 添加具有相同二进制指数和不同符号的a +b 。

2. 从0.5 <= |x| < 1.0的数字减去相同的符号0.5 0.5 <= |x| < 1.0范围。

3. ldexp*() （将数字转换为有效和指数部分）函数返回一个确切的值。

每个指数的表单数组

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

将这组数字添加到sums[]

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

归一化sum[]

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

检查sum[]是否为最重要的数字

最重要的（非零） one[s]是结果的符号。

---

下面的代码使用float作为set的FP类型来执行任务。 一些并行计算使用double来检查是否合理，但不会对float计算做出贡献。

最后的标准化步骤通常重复两次。 即使是最坏的情况集，我怀疑会迭代float的二进制宽度，大约23次。

解决方案似乎是关于O(n) ，但确实使用了一个大约FP指数范围大小的数组。

 #include  #include  #include  #include  #include  #include  #include  #include  #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(uf)); return uf; } 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 %e\n", __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("%e\n", sums[expo]); printf("%d %e\n", 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 %+.*e\n", 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 %+.*e\n", "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= 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 = %u\n", seed); srand(seed); for (unsigned i = 0; i < ITERATIONS; i++) { int cmp = test_sum((size_t)rand() % MAX_NUMBERS_PER_SET + 1); printf("Compare %d\n\n", cmp); if (cmp == 0) break; } printf("Success"); return EXIT_SUCCESS; } 

Infinities和NaN也可以在一定程度上处理，以备日后使用。

Answer 3

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

我₁ + i ₂ = 我₁₂
我₁₂ – i ₂ = i _~1
i ₁ – i _~1 = _iΔ

对于n 个加法运算的总和，我们可以得出的最接近的近似值是计算n – 1个加法运算的误差，然后再将n – 1个 n – 1误差相加。 你将重复这个过程n – 2次或直到所有错误都变为0.0

可以采取一些措施将错误计算驱动到0.0：

1. 使用较大的浮点类型，例如long double
2. 在求和之前对列表进行排序，这样您就可以将小数字添加到小数字，将大数字添加到大数字

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

普遍接受的解决方案是Kahan的Summation，它是速度和精度之间的完美结合。 而不是将误差保持到求和的结尾，Kahan将把它滚动到每个加法中，防止它的值在最高精度浮点范围之外升级。 假设我们给了vector i1我们可以按如下方式运行Kahan的Summation：

 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; 

Answer 4

确定地进行这种比较的可能性之一是创建一个精确定点算术的类，其等级与使用的类型相同，并且不限制绝对值。

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

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

（每个受支持的浮点类型都需要单独的构造函数。）

根据具体情况，实施需要一系列固定的，决定建筑或动态尺寸的集合; 可能是std::bitset ， vector或静态或动态bool数组。

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

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