сравнение двух сумм значений с плавающей запятой в C или переполнении стека

Сообщение изначально было также написано на C, и мой код относится к этому.
Теперь я вижу, что пост предназначен только для C ++, но в следующем я не вижу ничего такого, что было бы неприменимо к C ++.

Упростить поиск знака суммы списка чисел ФП

Сравнение двух наборов чисел аналогично добавлению отрицания второго набора к первому и затем нахождению знака суммы объединенного списка. Этот знак отображается на >, == или же < из 2 оригинальных наборов.

Выполнять только точную математику FP

предположение: FP использует IEEE-подобные числа, включая субнормалы, основание 2, и является точным для определенных операций:

1. Добавление a +b с одинаковым двоичным показателем и отличающимся знаком.

2. Вычитание того же знака 0.5 из числа в 0.5 <= |x| < 1.0 спектр.

3. ldexp*() (разбить число на значимые и экспоненциальные части) функция возвращает точное значение.

Массив формы на экспоненту

Формируем массив сумм sums[] чьи ценности будут только когда-либо (0 or 0.5 <= |sums[i]| < 1.0)по одному для каждого возможного показателя и для некоторых показателей, превышающих макс. Более крупные необходимы для накопления |total_sum| это превышает FP_MAX, Это нужно до log2(SIZE_MAX) больше элементов.

Добавить набор чисел в sums[]

Для каждого элемента набора номера, добавьте его в соответствующий sums[] согласно его двоичному показателю. Это является ключевым моментом, так как может быть добавлено одинаковое число знаков и числа FP с разными знаками с общим двоичным показателем FP именно так. Добавление может привести к переносу с одинаковыми значениями знака и отмене с различными значениями знака — это обрабатывается. Входящий набор номеров не нужно сортировать.

нормировать sum[]

Для каждого элемента на ones[], убедитесь, что любые значения не 0,5, 0,0 или -0,5 уменьшены, оставшаяся часть добавлена ​​к меньшему ones[],

Осмотреть sum[] для самой значимой цифры

Наиболее значимый (ненулевой) one[s] является признаком результата.

Приведенный ниже код выполняет задачу, используя float как тип FP набора. Некоторые параллельные вычисления выполняются с использованием double проверить на здравомыслие, но не способствует float расчет.

Шаг нормализации в конце повторяется, как правило, дважды. Я подозреваю, что даже в худшем случае будет повторяться двоичная ширина float signicand, около 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 также могут быть обработаны, до некоторой степени, оставить это на потом.

Число с плавающей запятой, полученное в результате суммирования двух чисел с плавающей запятой, является только приближение. Дано я₁ а также я₂ Подводя итог, мы можем найти приближение ошибки в суммировании с плавающей точкой, выполнив это:

я₁ + я₂ знак равно я₁₂
я₁₂ — я₂ знак равно я_{~ 1}
я₁ — я_{~ 1} знак равно я_Δ

Ближайший приближение мы могли бы придумать для суммирования N числа будут вычислять ошибку для N — 1 операции сложения, а затем суммировать те N — 1 ошибки снова принимают N — 2. И вы повторите этот процесс N — 2 раз или пока все ошибки не ушли в 0.0

Есть пара вещей, которые можно сделать, чтобы ускорить вычисления ошибок до 0,0:

1. Используйте больший тип с плавающей точкой, например long double
2. Сортируйте список перед суммированием, чтобы вы добавляли маленькие числа к маленьким и большие числа к большим

Теперь вы можете оценить, насколько важна для вас точность. Я скажу вам, что в общем случае вычислительные затраты на вышеуказанную операцию возмутительны, учитывая результат, который вы получите все еще быть приближением.

Общепринятое решение Суммирование Кахана это счастливый брак между скоростью и точностью. Вместо того, чтобы удерживать ошибку до конца суммирования, Кахан будет прокручивать ее в каждом добавлении, предотвращая эскалацию ее значения за пределы диапазона с плавающей запятой с наивысшей точностью. Скажи, что нам дали vector<long double> i1 мы могли бы запустить Суммирование Кахана на нем следующим образом:

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.

Это очевидное и очень дорогое решение, которое ухудшило бы производительность, если бы это сравнение было ядром какой-либо системы. Надеюсь, есть лучшее решение.