Безопасное подразделение с плавающей точкой

В моем коде есть места, где я хочу убедиться, что деление двух произвольных чисел с плавающей запятой (32-битная одинарная точность) не будет переполнено. Цель / компилятор не гарантирует (достаточно явно) удобную обработку -INF/INF и (не полностью гарантирует IEEE 754 для исключительных значений - (возможно, неопределенных) - и цель может измениться). Также я не могу сделать предположения о сохранении входных данных для этих нескольких особых мест, и я привязан к стандартным библиотекам C90.

Я прочитал " Что должен знать каждый компьютерщик" об арифметике с плавающей точкой, но, если честно, я немного растерялся.

Итак... Я хочу спросить сообщество, поможет ли следующий фрагмент кода, и есть ли лучшие / более быстрые / точные / корректирующие способы сделать это:

#define SIGN_F(val) ((val >= 0.0f)? 1.0f : -1.0f)

float32_t safedivf(float32_t num, float32_t denum)
{
   const float32_t abs_denum = fabs(denum);
   if((abs_denum < 1.0f) && ((abs_denum * FLT_MAX) <= (float32_t)fabs(num))
       return SIGN_F(denum) * SIGN_F(num) * FLT_MAX;
   else
       return num / denum;
}

Изменить: Изменено ((abs_denum * FLT_MAX) < (float32_t)fabs(num)) в ((abs_denum * FLT_MAX) <= (float32_t)fabs(num)) как рекомендует Паскаль Куок.

Источник

4 ответа

Решение

В ((abs_denum * FLT_MAX) < (float32_t)fabs(num), продукт abs_denum * FLT_MAX может округляться и заканчиваться равным fabs(num), Это не значит что num / denum не математически больше, чем FLT_MAXи вы должны быть обеспокоены тем, что это может вызвать переполнение, которого вы хотите избежать. Вам лучше заменить это < от <=,

Для альтернативного решения, если double тип доступен и шире чем float, это может быть более экономичным, чтобы вычислить (double)num/(double)denum, Если float является бинарным и double двоичный 64, единственный способ double деление может переполниться, если denum равно (а) нулю (и если denum это ноль, твой код тоже проблематичен).

double dbl_res = (double)num/(double)denum;
float res = dbl_res < -FLT_MAX ? -FLT_MAX : dbl_res > FLT_MAX ? FLT_MAX : (float)dbl_res;
Источник

Извините, я не могу комментировать из-за моей низкой репутации, но вы можете попытаться извлечь показатели и мантиссы чисел и чисел и убедиться, что это условие:

((exp(num) - exp (denum)) > max_exp) &&  (mantissa(num) >= mantissa(denum))

И в соответствии со знаком входных данных генерировать соответствующий INF.

Источник

Тщательно работать с num, denom когда частное близко FLT_MAX,

Следующее использует тесты, вдохновленные OP, но держится подальше от результатов рядом FLT_MAX, Как отмечает @Pascal Cuoq, округление может просто сдвинуть результат с края. Вместо этого он использует пороги FLT_MAX/FLT_RADIX а также FLT_MAX*FLT_RADIX,

Масштабируя с FLT_RADIXобычно 2, код всегда должен получать точные результаты. Ожидается, что округление в любом режиме округления не повлияет на результат.

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

static int SD_Sign(float x) {
  if (x > 0.0f)
    return 1;
  if (x < 0.0f)
    return -1;
  if (atan2f(x, -1.0f) > 0.0f)
    return 1;
  return -1;
}

static float SD_Overflow(float num, float denom) {
  return SD_Sign(num) * SD_Sign(denom) * FLT_MAX;
}

float safedivf(float num, float denom) {
  float abs_denom = fabsf(denom);
  // If |quotient| > |num|
  if (abs_denom < 1.0f) {
    float abs_num = fabsf(num);
    // If |num/denom| > FLT_MAX/2 --> quotient is very large or overflows
    // This computation is safe from rounding and overflow.
    if (abs_num > FLT_MAX / FLT_RADIX * abs_denom) {
      // If |num/denom| >= FLT_MAX*2 --> overflow
      // This also catches denom == 0.0
      if (abs_num / FLT_RADIX >= FLT_MAX * abs_denom) {
        return SD_Overflow(num, denom);
      }
      // At this point, quotient must be in or near range FLT_MAX/2 to FLT_MAX*2
      // Scale parameters so quotient is a FLT_RADIX * FLT_RADIX factor smaller.
      if (abs_num > 1.0) {
        abs_num /= FLT_RADIX * FLT_RADIX;
      } else {
        abs_denom *= FLT_RADIX * FLT_RADIX;
      }
      float quotient = abs_num / abs_denom;
      if (quotient > FLT_MAX / (FLT_RADIX * FLT_RADIX)) {
        return SD_Overflow(num, denom);
      }
    }
  }
  return num / denom;
}

SIGN_F() необходимо учитывать в denum является +0.0 или же -0.0, Различные методы, упомянутые @Pascal Cuoq в комментарии:

  1. copysign() или же signbit()
  2. Используйте союз

Дополнительно, некоторые функции, когда они хорошо реализованы, дифференцируются на +/- ноль, как atan2f(zero, -1.0) а также sprintf(buffer, "%+f", zero),

Примечание: используется float против float32_t для простоты.
Примечание: возможно использовать fabsf() скорее, чем fabs(),
Незначительный: Предложить denom (знаменатель) вместо denum,

Источник

Чтобы избежать угловых случаев с округлением, а что нет, вы можете помассировать показатель степени делителя - с помощью frexp() и ldexp() - и беспокоиться о том, можно ли уменьшить результат без переполнения. Или frexp() оба аргумента, и выполняйте экспоненциальную работу вручную.

Источник
Другие вопросы по тегам