浮点舍入错误将报告的结果移动到范围内
Floating point rounding error moves reported result inside range
我正在开发一个功能来报告测试结果以及该特定测试结果的下限和上限。这三个值将使用指定的公式 (aX + b)/c 进行转换,其中 X 为 testResult/lowerLimit/upperLimit,a、b 和 c 为浮点数。
如果报告的测试结果是inside/outside转换前的指定限值,它也应是inside/outside转换后的限值,以确保报告结果的有效性。
我已经确定了两种情况,其中无效测试结果在转换后会移动到范围内,但我还没有找到测试结果在转换前在范围内而在转换后超出指定限制的情况转换。这种情况还能发生吗?我不相信?可以吗?
下面是生成我提到的两个案例的一些代码以及更正以确保报告的测试结果的有效性。
TLDR:下面代码中的 ((TRUE == insideLimitBefore) && (FALSE == insideLimitAfter)) 情况会发生吗?
#include <stdio.h>
#include <stdint.h>
#define TRUE (uint8_t)0x01
#define FALSE (uint8_t)0x00
int32_t LinearMapping(const int32_t input);
void Convert(int32_t testResult, int32_t lowerLimit, int32_t upperLimit);
int main(void)
{
int32_t lowerLimit = 504;
int32_t testResult = 503;
int32_t upperLimit = 1000;
printf("INPUT:\n\tLower limit:\t%d\t\n\tTest result:\t%d\t\n\tUpper limit:\t%d\t\n", lowerLimit, testResult, upperLimit);
Convert(testResult, lowerLimit, upperLimit);
lowerLimit = 500;
testResult = 504;
upperLimit = 503;
printf("INPUT:\n\tLower limit:\t%d\t\n\tTest result:\t%d\t\n\tUpper limit:\t%d\t\n", lowerLimit, testResult, upperLimit);
Convert(testResult, lowerLimit, upperLimit);
return 0;
}
int32_t LinearMapping(const int32_t input)
{
float retVal;
const float a = 1.0;
const float b = 1.0;
const float c = 2.3;
retVal = a * input;
retVal += b;
retVal /= c;
return (int32_t)retVal;
}
void Convert(int32_t testResult, int32_t lowerLimit, int32_t upperLimit)
{
uint8_t insideLimitAfter;
uint8_t belowLowerLimit;
uint8_t insideLimitBefore = ((lowerLimit <= testResult) && (testResult <= upperLimit)) ? TRUE : FALSE;
if (FALSE == insideLimitBefore)
{
/* testResult is either below or above lowerLimit/upperLimit respectively */
if (testResult < lowerLimit)
{
belowLowerLimit = TRUE;
}
else /* testResult > upperLimit */
{
belowLowerLimit = FALSE;
}
}
testResult = LinearMapping(testResult);
lowerLimit = LinearMapping(lowerLimit);
upperLimit = LinearMapping(upperLimit);
insideLimitAfter = ((lowerLimit <= testResult) && (testResult <= upperLimit)) ? TRUE : FALSE;
if ((FALSE == insideLimitBefore) && (TRUE == insideLimitAfter))
{
if (TRUE == belowLowerLimit)
{
printf("OUTPUT:\n\tLower limit:\t%d\t\n\tTest result:\t%d\t\n\tUpper limit:\t%d\t\n", lowerLimit+1, testResult, upperLimit);
}
else /* belowLowerLimit == FALSE => testResult > upperLimit */
{
printf("OUTPUT:\n\tLower limit:\t%d\t\n\tTest result:\t%d\t\n\tUpper limit:\t%d\t\n", lowerLimit, testResult, upperLimit-1);
}
}
else if ((TRUE == insideLimitBefore) && (FALSE == insideLimitAfter))
{
/* Is this case even possible? */
}
else
{
/* Do nothing */
}
}
to find a case where the test result is inside the range before the conversion and will be outside the specified limits after conversion. Can this case even occur?
不,考虑到理智 a,b,c, lowerLimit, testResult, upperLimit
。
给定 3 个值 lo,x,hi
lo <= x <= hi
在 LinearMapping()
中的线性转换之前 lo_new <= x_new <= hi_new
将保持相同的关系,只要转换是(正)线性(不除以 0,a
、b
、c
不是非数字)。不转换超出 int32_t
.
范围的 float
主要原因是在 x
的边缘情况下 [lo...hi]
, x
可能会降低所有 3 个的有效精度。新的 x
现在可能等于 lo
或 hi
并且 ==
有利于 "in range"。所以 lo <= x <= hi
.
没有变化
OP 最初找到了 "invalid test result will move inside the range after the conversion" 的示例,因为 x
就在 [lo...hi]
之外,现在有效的精度降低使得 x
等于 lo
或 hi
。由于==
偏爱"in range",可见由外向内的走势。
注意:如果LinearMapping()
有负斜率,比如-1,那么lo <= x <= hi
很容易被打破。作为 1 <= 2 <= 3
--> -1 > -2 > -3
。这使得 lowerLimit > upperLimit
其中 "in range" 不能满足任何 x
.
作为参考,OP 的代码已简化:
#include <stdio.h>
#include <stdint.h>
int LinearMapping(const int input) {
const float a = 1.0;
const float b = 1.0;
const float c = 2.3;
float retVal = a * input;
retVal += b;
retVal /= c;
return (int) retVal;
}
void Convert(int testResult, int lowerLimit, int upperLimit) {
printf("Before %d %s %d %s %d\n", lowerLimit,
lowerLimit <= testResult ? "<=" : "> ", testResult,
testResult <= upperLimit ? "<=" : "> ", upperLimit);
testResult = LinearMapping(testResult);
lowerLimit = LinearMapping(lowerLimit);
upperLimit = LinearMapping(upperLimit);
printf("After %d %s %d %s %d\n\n", lowerLimit,
lowerLimit <= testResult ? "<=" : "> ", testResult,
testResult <= upperLimit ? "<=" : "> ", upperLimit);
}
int main(void) {
Convert(503, 504, 1000);
Convert(504, 500, 503);
return 0;
}
输出
Before 504 > 503 <= 1000
After 219 <= 219 <= 435
Before 500 <= 504 > 503
After 217 <= 219 <= 219
… I have yet to find a case where the test result is inside the range before the conversion and will be outside the specified limits after conversion. Can this case even occur? I don't believe so? Can it?
是的,这在理论上是可能发生的,尽管是由于 C 行为而不是由于底层浮点运算。 C 标准不保证使用 IEEE-754 浮点运算,甚至不保证通过计算表达式具有相同的精度,这可能导致表达式的相同输入有不同的结果。
虽然 LinearMapping
显示为单个例程,但编译器可能会将其内联。也就是说,在调用例程的地方,编译器可以用例程的主体替换调用。此外,当它在不同的地方执行此操作时,它可能会使用不同的方法评估表达式。因此,在此代码中,可以在每次调用中使用不同的浮点运算来计算 LinearMapping
:
testResult = LinearMapping(testResult);
lowerLimit = LinearMapping(lowerLimit);
upperLimit = LinearMapping(upperLimit);
这意味着 (a * testResult + b) / c
可能仅使用 32 位浮点运算进行计算,而 (a * upperLimit + b) / c
可能使用 64 位浮点运算进行计算,之后转换为 32 位师。 (为简洁起见,我已将您的三个赋值语句合并为一个表达式。该问题适用于任何一种方式。)
这样做的一个结果是双舍入。当以一种精度计算结果然后转换为另一种精度时,会发生两次舍入,一次在初始计算中,第二次在转换中。考虑如下所示的精确数学结果:
1.xxxxx01011111111xxxx1xx
^ ^ Start of bits to be rounded in wider format.
| Start of bits to be rounded in narrower format.
如果这是更窄格式的计算结果,我们将检查位 011111111xxxx1xx 并将它们向下舍入(在我们舍入的位置它们小于 ½),因此最终结果将是 1 .xxxxx01。但是,如果我们首先以更宽的格式进行计算,则要四舍五入的位为 1xxxx1xx(大于 ½),并且将这些四舍五入,使中间结果为 1.xxxxx0110000000。当我们转换为更窄的格式时,要四舍五入的位是 10000000,正好是中点 (½),所以四舍五入到最接近偶数规则告诉我们四舍五入,这使得最终结果 1.xxxxx10.
因此,即使 testResult
和 upperLimit
相等,对它们应用 LinearMapping
的结果也可能不相等,并且 testResult
可能看起来超出了区间.
可能有一些方法可以避免这个问题:
如果您的 C 实现符合 C 标准的附件 F(本质上说它使用 IEEE-754 操作并以预期的方式将它们绑定到 C 运算符)或至少符合某些部分,那么双舍入不应该出现在编写良好的源代码中。
C 标准说实现应该在 <float.h>
中定义 FLT_EVAL_METHOD
。如果 FLT_EVAL_METHOD
为 0,则表示所有浮点运算和常量都使用其标称类型进行评估。在这种情况下,如果您在源代码中使用单个浮点类型,则不会发生双舍入。如果 FLT_EVAL_METHOD
为 1,则 float
操作使用 double
计算。在这种情况下,您可以使用 double
而不是 float
来避免双重舍入。如果它是 2,则使用 long double
计算操作,并且可以使用 long double
避免双舍入。如果 FLT_EVAL_METHOD
为 -1,用于评估的浮点格式是不确定的,因此需要考虑双舍入。
对于特定值或值区间 and/or 已知的浮点格式,可能可以证明不会发生双舍入。例如,假设您的输入都是 int32_t
,线性映射参数是特定值,并且仅使用 32 位或 64 位 IEEE-754 二进制浮点数,则可能可以证明不会发生双舍入。
即使您的实现符合附件 F 或将 FLT_EVAL_METHOD
定义为非负数,您仍然必须在代码中注意不要使用类型为 double
的表达式,然后将其赋值float
类型的对象。这将导致双舍入,因为您的源代码明确要求它,而不是因为 C 对浮点数不严格。
作为具体示例,请考虑 (1.0 * 13546 + 1.0) / 2.3
。如果浮点常量以 64 位二进制浮点数(53 位有效数)表示,并且以 64 位二进制计算表达式,则结果为 5890.0000000000009094947017729282379150390625。但是,如果使用相同的常量(64 位二进制),但使用 Intel 的 80 位二进制(64 位有效数)计算表达式,然后转换为 64 位,则结果为 5890。
在这种情况下,精确的数学商为:
1.01110000001000000000000000000000000000000000000000001000000000001011…
^ Bits rounded away in double.
如果我们将其四舍五入,我们可以看到要四舍五入的位,1000000000001011…,在四舍五入的位置大于½,所以我们四舍五入。如果我们用 long double 对其进行四舍五入,则要四舍五入的位是 01011……。这些四舍五入,剩下:
1.011100000010000000000000000000000000000000000000000010000000000
^ Bits rounded away in double.
现在,当我们四舍五入时,四舍五入的位数是 10000000000,也就是中点。规则说四舍五入使低位均匀,所以结果是:
1.011100000010000000000000000000000000000000000000000010000000000
我正在开发一个功能来报告测试结果以及该特定测试结果的下限和上限。这三个值将使用指定的公式 (aX + b)/c 进行转换,其中 X 为 testResult/lowerLimit/upperLimit,a、b 和 c 为浮点数。
如果报告的测试结果是inside/outside转换前的指定限值,它也应是inside/outside转换后的限值,以确保报告结果的有效性。
我已经确定了两种情况,其中无效测试结果在转换后会移动到范围内,但我还没有找到测试结果在转换前在范围内而在转换后超出指定限制的情况转换。这种情况还能发生吗?我不相信?可以吗?
下面是生成我提到的两个案例的一些代码以及更正以确保报告的测试结果的有效性。
TLDR:下面代码中的 ((TRUE == insideLimitBefore) && (FALSE == insideLimitAfter)) 情况会发生吗?
#include <stdio.h>
#include <stdint.h>
#define TRUE (uint8_t)0x01
#define FALSE (uint8_t)0x00
int32_t LinearMapping(const int32_t input);
void Convert(int32_t testResult, int32_t lowerLimit, int32_t upperLimit);
int main(void)
{
int32_t lowerLimit = 504;
int32_t testResult = 503;
int32_t upperLimit = 1000;
printf("INPUT:\n\tLower limit:\t%d\t\n\tTest result:\t%d\t\n\tUpper limit:\t%d\t\n", lowerLimit, testResult, upperLimit);
Convert(testResult, lowerLimit, upperLimit);
lowerLimit = 500;
testResult = 504;
upperLimit = 503;
printf("INPUT:\n\tLower limit:\t%d\t\n\tTest result:\t%d\t\n\tUpper limit:\t%d\t\n", lowerLimit, testResult, upperLimit);
Convert(testResult, lowerLimit, upperLimit);
return 0;
}
int32_t LinearMapping(const int32_t input)
{
float retVal;
const float a = 1.0;
const float b = 1.0;
const float c = 2.3;
retVal = a * input;
retVal += b;
retVal /= c;
return (int32_t)retVal;
}
void Convert(int32_t testResult, int32_t lowerLimit, int32_t upperLimit)
{
uint8_t insideLimitAfter;
uint8_t belowLowerLimit;
uint8_t insideLimitBefore = ((lowerLimit <= testResult) && (testResult <= upperLimit)) ? TRUE : FALSE;
if (FALSE == insideLimitBefore)
{
/* testResult is either below or above lowerLimit/upperLimit respectively */
if (testResult < lowerLimit)
{
belowLowerLimit = TRUE;
}
else /* testResult > upperLimit */
{
belowLowerLimit = FALSE;
}
}
testResult = LinearMapping(testResult);
lowerLimit = LinearMapping(lowerLimit);
upperLimit = LinearMapping(upperLimit);
insideLimitAfter = ((lowerLimit <= testResult) && (testResult <= upperLimit)) ? TRUE : FALSE;
if ((FALSE == insideLimitBefore) && (TRUE == insideLimitAfter))
{
if (TRUE == belowLowerLimit)
{
printf("OUTPUT:\n\tLower limit:\t%d\t\n\tTest result:\t%d\t\n\tUpper limit:\t%d\t\n", lowerLimit+1, testResult, upperLimit);
}
else /* belowLowerLimit == FALSE => testResult > upperLimit */
{
printf("OUTPUT:\n\tLower limit:\t%d\t\n\tTest result:\t%d\t\n\tUpper limit:\t%d\t\n", lowerLimit, testResult, upperLimit-1);
}
}
else if ((TRUE == insideLimitBefore) && (FALSE == insideLimitAfter))
{
/* Is this case even possible? */
}
else
{
/* Do nothing */
}
}
to find a case where the test result is inside the range before the conversion and will be outside the specified limits after conversion. Can this case even occur?
不,考虑到理智 a,b,c, lowerLimit, testResult, upperLimit
。
给定 3 个值 lo,x,hi
lo <= x <= hi
在 LinearMapping()
中的线性转换之前 lo_new <= x_new <= hi_new
将保持相同的关系,只要转换是(正)线性(不除以 0,a
、b
、c
不是非数字)。不转换超出 int32_t
.
float
主要原因是在 x
的边缘情况下 [lo...hi]
, x
可能会降低所有 3 个的有效精度。新的 x
现在可能等于 lo
或 hi
并且 ==
有利于 "in range"。所以 lo <= x <= hi
.
OP 最初找到了 "invalid test result will move inside the range after the conversion" 的示例,因为 x
就在 [lo...hi]
之外,现在有效的精度降低使得 x
等于 lo
或 hi
。由于==
偏爱"in range",可见由外向内的走势。
注意:如果LinearMapping()
有负斜率,比如-1,那么lo <= x <= hi
很容易被打破。作为 1 <= 2 <= 3
--> -1 > -2 > -3
。这使得 lowerLimit > upperLimit
其中 "in range" 不能满足任何 x
.
作为参考,OP 的代码已简化:
#include <stdio.h>
#include <stdint.h>
int LinearMapping(const int input) {
const float a = 1.0;
const float b = 1.0;
const float c = 2.3;
float retVal = a * input;
retVal += b;
retVal /= c;
return (int) retVal;
}
void Convert(int testResult, int lowerLimit, int upperLimit) {
printf("Before %d %s %d %s %d\n", lowerLimit,
lowerLimit <= testResult ? "<=" : "> ", testResult,
testResult <= upperLimit ? "<=" : "> ", upperLimit);
testResult = LinearMapping(testResult);
lowerLimit = LinearMapping(lowerLimit);
upperLimit = LinearMapping(upperLimit);
printf("After %d %s %d %s %d\n\n", lowerLimit,
lowerLimit <= testResult ? "<=" : "> ", testResult,
testResult <= upperLimit ? "<=" : "> ", upperLimit);
}
int main(void) {
Convert(503, 504, 1000);
Convert(504, 500, 503);
return 0;
}
输出
Before 504 > 503 <= 1000
After 219 <= 219 <= 435
Before 500 <= 504 > 503
After 217 <= 219 <= 219
… I have yet to find a case where the test result is inside the range before the conversion and will be outside the specified limits after conversion. Can this case even occur? I don't believe so? Can it?
是的,这在理论上是可能发生的,尽管是由于 C 行为而不是由于底层浮点运算。 C 标准不保证使用 IEEE-754 浮点运算,甚至不保证通过计算表达式具有相同的精度,这可能导致表达式的相同输入有不同的结果。
虽然 LinearMapping
显示为单个例程,但编译器可能会将其内联。也就是说,在调用例程的地方,编译器可以用例程的主体替换调用。此外,当它在不同的地方执行此操作时,它可能会使用不同的方法评估表达式。因此,在此代码中,可以在每次调用中使用不同的浮点运算来计算 LinearMapping
:
testResult = LinearMapping(testResult);
lowerLimit = LinearMapping(lowerLimit);
upperLimit = LinearMapping(upperLimit);
这意味着 (a * testResult + b) / c
可能仅使用 32 位浮点运算进行计算,而 (a * upperLimit + b) / c
可能使用 64 位浮点运算进行计算,之后转换为 32 位师。 (为简洁起见,我已将您的三个赋值语句合并为一个表达式。该问题适用于任何一种方式。)
这样做的一个结果是双舍入。当以一种精度计算结果然后转换为另一种精度时,会发生两次舍入,一次在初始计算中,第二次在转换中。考虑如下所示的精确数学结果:
1.xxxxx01011111111xxxx1xx ^ ^ Start of bits to be rounded in wider format. | Start of bits to be rounded in narrower format.
如果这是更窄格式的计算结果,我们将检查位 011111111xxxx1xx 并将它们向下舍入(在我们舍入的位置它们小于 ½),因此最终结果将是 1 .xxxxx01。但是,如果我们首先以更宽的格式进行计算,则要四舍五入的位为 1xxxx1xx(大于 ½),并且将这些四舍五入,使中间结果为 1.xxxxx0110000000。当我们转换为更窄的格式时,要四舍五入的位是 10000000,正好是中点 (½),所以四舍五入到最接近偶数规则告诉我们四舍五入,这使得最终结果 1.xxxxx10.
因此,即使 testResult
和 upperLimit
相等,对它们应用 LinearMapping
的结果也可能不相等,并且 testResult
可能看起来超出了区间.
可能有一些方法可以避免这个问题:
如果您的 C 实现符合 C 标准的附件 F(本质上说它使用 IEEE-754 操作并以预期的方式将它们绑定到 C 运算符)或至少符合某些部分,那么双舍入不应该出现在编写良好的源代码中。
C 标准说实现应该在
<float.h>
中定义FLT_EVAL_METHOD
。如果FLT_EVAL_METHOD
为 0,则表示所有浮点运算和常量都使用其标称类型进行评估。在这种情况下,如果您在源代码中使用单个浮点类型,则不会发生双舍入。如果FLT_EVAL_METHOD
为 1,则float
操作使用double
计算。在这种情况下,您可以使用double
而不是float
来避免双重舍入。如果它是 2,则使用long double
计算操作,并且可以使用long double
避免双舍入。如果FLT_EVAL_METHOD
为 -1,用于评估的浮点格式是不确定的,因此需要考虑双舍入。对于特定值或值区间 and/or 已知的浮点格式,可能可以证明不会发生双舍入。例如,假设您的输入都是
int32_t
,线性映射参数是特定值,并且仅使用 32 位或 64 位 IEEE-754 二进制浮点数,则可能可以证明不会发生双舍入。
即使您的实现符合附件 F 或将 FLT_EVAL_METHOD
定义为非负数,您仍然必须在代码中注意不要使用类型为 double
的表达式,然后将其赋值float
类型的对象。这将导致双舍入,因为您的源代码明确要求它,而不是因为 C 对浮点数不严格。
作为具体示例,请考虑 (1.0 * 13546 + 1.0) / 2.3
。如果浮点常量以 64 位二进制浮点数(53 位有效数)表示,并且以 64 位二进制计算表达式,则结果为 5890.0000000000009094947017729282379150390625。但是,如果使用相同的常量(64 位二进制),但使用 Intel 的 80 位二进制(64 位有效数)计算表达式,然后转换为 64 位,则结果为 5890。
在这种情况下,精确的数学商为:
1.01110000001000000000000000000000000000000000000000001000000000001011… ^ Bits rounded away in double.
如果我们将其四舍五入,我们可以看到要四舍五入的位,1000000000001011…,在四舍五入的位置大于½,所以我们四舍五入。如果我们用 long double 对其进行四舍五入,则要四舍五入的位是 01011……。这些四舍五入,剩下:
1.011100000010000000000000000000000000000000000000000010000000000 ^ Bits rounded away in double.
现在,当我们四舍五入时,四舍五入的位数是 10000000000,也就是中点。规则说四舍五入使低位均匀,所以结果是:
1.011100000010000000000000000000000000000000000000000010000000000