浮点舍入误差会将报告结果移动到范围内

Posted

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了浮点舍入误差会将报告结果移动到范围内相关的知识，希望对你有一定的参考价值。

我正在开发一个函数来报告测试结果以及该特定测试结果的下限和上限。这三个值将使用指定的公式（aX + b）/ c进行转换，其中X是testResult / lowerLimit / upperLimit，a，b和c是浮点数。

如果报告的测试结果在转换之前在指定限制之内/之外，则转换后也应在限制之内/之外，以确保报告结果的有效性。

我已经确定了两种情况，其中无效的测试结果将在转换后在范围内移动但我还没有找到测试结果在转换之前的范围内并且在转换之后将超出指定限制的情况。这种情况甚至会发生吗？我不相信吗？它可以？

下面是一些代码，它产生了我提到的两个案例以及更正，以确保报告的测试结果的有效性。

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:
	Lower limit:	%d	
	Test result:	%d	
	Upper limit:	%d	
", lowerLimit, testResult, upperLimit);
    Convert(testResult, lowerLimit, upperLimit);

    lowerLimit = 500;
    testResult = 504;
    upperLimit = 503;

    printf("INPUT:
	Lower limit:	%d	
	Test result:	%d	
	Upper limit:	%d	
", 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:
	Lower limit:	%d	
	Test result:	%d	
	Upper limit:	%d	
", lowerLimit+1, testResult, upperLimit);
        }
        else /* belowLowerLimit == FALSE => testResult > upperLimit */
        {
            printf("OUTPUT:
	Lower limit:	%d	
	Test result:	%d	
	Upper limit:	%d	
", lowerLimit, testResult, upperLimit-1);
        }
    }
    else if ((TRUE == insideLimitBefore) && (FALSE == insideLimitAfter))
    {
       /* Is this case even possible? */
    }
    else
    {
        /* Do nothing */
    }
}
答案

查找测试结果在转换前的范围内并且在转换后超出指定限制的情况。这种情况甚至会发生吗？

不，给予理智的a,b,c, lowerLimit, testResult, upperLimit。

在lo,x,hi线性转换之前给定3个lo <= x <= hi值LinearMapping()，lo_new <= x_new <= hi_new将保持相同的关系，只要转换是（正）线性（没有除以0，a，b，c不是非A数）。没有转换float超出范围的int32_t。

主要原因是x的边缘情况在内部或极限，[lo...hi]，LinearMapping()可能会降低所有3的有效精度。新的x现在可能等于lo或hi和==偏好“范围内”。所以没有改变lo <= x <= hi。

OP最初发现“无效测试结果将在转换后的范围内移动”的例子，因为x正好在[lo...hi]之外，现在有效的精确度降低使得x等于lo或hi。由于==倾向于“在范围内”，所以可以看到从外到内的移动。

注意：如果LinearMapping()具有-1的负斜率，那么lo <= x <= hi很容易被破坏。作为1 <= 2 <= 3 - > -1 > -2 > -3。这使lowerLimit > upperLimit“在范围内”不能满足任何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
", lowerLimit,
      lowerLimit <= testResult ? "<=" : "> ", testResult,
      testResult <= upperLimit ? "<=" : "> ", upperLimit);
  testResult = LinearMapping(testResult);
  lowerLimit = LinearMapping(lowerLimit);
  upperLimit = LinearMapping(upperLimit);
  printf("After  %d %s %d %s %d

", 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
另一答案

...我还没有找到测试结果在转换前的范围内并且在转换后超出指定限制的情况。这种情况甚至会发生吗？我不相信吗？它可以？

是的，理论上可能会发生这种情况，尽管由于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标准说实施应该在FLT_EVAL_METHOD中定义<float.h>。如果FLT_EVAL_METHOD为0，则表示所有浮点运算和常量都使用其标称类型进行计算。在这种情况下，如果在源代码中使用单个浮点类型，则不会发生双舍入。如果FLT_EVAL_METHOD为1，则使用float评估double操作。在这种情况下，您可以通过使用double而不是float来避免双舍入。如果是2，则使用long double评估操作，并且可以通过使用long double避免双舍入。如果FLT_EVAL_METHOD为-1，则用于评估的浮点格式是不确定的，因此需要考虑双舍入。
  • 对于值的特定值或间隔和/或已知的浮点格式，可能证明不会发生双舍入。例如，假设您的输入都是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.

如果我们将它舍入为double，我们可以看到要舍入的位，1000000000001011 ...，在舍入位置大于½，所以我们向上舍入。如果我们将它舍入为long double，则要舍入的位数为01011 ....这些向下舍入，留下：

1.011100000010000000000000000000000000000000000000000010000000000
                                                      ^ Bits rounded away in double.

现在，当我们舍入到两倍时，要舍入的位是10000000000，这是中点。规则说要舍入以使低位均匀，因此结果是：

1.011100000010000000000000000000000000000000000000000010000000000

以上是关于浮点舍入误差会将报告结果移动到范围内的主要内容，如果未能解决你的问题，请参考以下文章

准确预测任意浮点格式之间转换的舍入误差

float精度问题

C语言程序设计，谭浩强老师第三版里面的一个关于浮点型数据的舍入误差问题

关于浮点型误差的解决方法

Enumerable#sum 如何避免浮点舍入错误？

NOIP2017普及组解题报告