安全浮点分部
我的代码中有一些地方，我想确保2个任意浮点数（32位单精度）的除法不会溢出。 目标/编译器不保证（明确地）对-INF / INF的良好处理和（不完全保证IEEE 754的exception值 – （可能未定义） – 并且目标可能会改变）。 此外，我无法对这几个特殊地点的输入进行保存，我必须使用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; } 编辑：根据Pascal Cuoq的建议，将((abs_denum * FLT_MAX) < (float32_t)fabs(num))更改为((abs_denum * FLT_MAX) <= (float32_t)fabs(num)) 。
在((abs_denum * FLT_MAX) < (float32_t)fabs(num) ，乘积abs_denum * FLT_MAX可以向下舍入并且最终等于fabs(num) 。这并不意味着num / denum在数学上不大于FLT_MAX ，并且您应该担心它可能会导致您想要避免的溢出。您最好将此替换为< by <= 。
对于替代解决方案，如果double类型可用并且比float更宽，则计算(double)num/(double)denum可能更经济。 如果float是binary32ish而double是binary64ish，那么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。
当商接近FLT_MAX时num, denom小心使用num, denom FLT_MAX 。
以下使用受OP启发的测试但远离FLT_MAX附近的FLT_MAX 。 正如@Pascal Cuoq指出的那样，舍入可能只会将结果推到边缘。 相反，它使用FLT_MAX/FLT_RADIX FLT_MAX*FLT_RADIX和FLT_MAX*FLT_RADIX FLT_MAX/FLT_RADIX阈值。
通过使用FLT_RADIX （通常为2）进行缩放，代码应始终获得精确结果。 在任何舍入模式下舍入都不会感染结果。
就速度而言，“快乐路径”，即结果肯定不会溢出应该是一个快速的计算。 仍然需要进行unit testing，但评论应该提供这种方法的要点。
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在评论中提到的各种方法：
-
copysign()或signbit() - 使用工会
另外，一些函数，如果实现良好，可以在+/- 0区分，如atan2f(zero, -1.0)和sprintf(buffer, "%+f", zero) 。
注意：为简单起见，使用float与float32_t 。
注意：也许使用fabsf()而不是fabs() 。
轻微：建议denom （分母）代替denum 。
为了避免带有舍入的拐角情况以及什么不是，你可以按下除数上的指数 – 使用frexp（）和ldexp（） – 并担心结果是否可以缩放而不会溢出。 或者frexp（）这两个参数，并且手工完成exponenent工作。