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 %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<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 = %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;
}
在一定程度上也可以处理无穷大和 NaN。