Безопасный отдел операций с плавающей запятой

В ((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 является двоичным 32, а double двоичным 64, единственный способ переполнения деления double - это если denum равно (a) нулю (и если 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;

Вы можете попытаться извлечь показатели и мантиссы num и denum и убедиться, что условие:

((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() or signbit()
2. Используйте союз

Кроме того, некоторые функции при правильной реализации различаются по +/- нулю, например atan2f(zero, -1.0) и sprintf(buffer, "%+f", zero).

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

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