Floating point rounding error moves reported result inside range

Question

I am working on a function to report test results together with lower and upper limit for that specific test result. These three values will be converted with a specified formula (aX + b)/c where X is testResult/lowerLimit/upperLimit and a,b and c are floating point numbers.

If the reported test result is inside/outside the specified limits before the conversion it shall also be inside/outside the limit after the conversion in order to ensure the validity of the reported results.

I have identified two cases where an invalid test result will move inside the range after the conversion but 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?

Below is some code that produces the two cases that I mentioned together with corrections to ensure the validity of the reported test result.

TLDR: Can the ((TRUE == insideLimitBefore) && (FALSE == insideLimitAfter)) case in the code below happen?

#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 */
    }
}

Answer 1

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?

No, given sane a,b,c, lowerLimit, testResult, upperLimit.

Given 3 values lo,x,hi that are lo <= x <= hi before the linear conversion in LinearMapping() the lo_new <= x_new <= hi_new will maintain the same relationship as long as the conversion is (positively) linear (no division by 0, a,b, c are not Not-A-Numbers). No conversion of a float that is out of range of int32_t.

The primary reason is on edge cases of x just inside or at the limit, [lo...hi] , the LinearMapping() may reduce the effective precision of all 3. The new x may now equal lo or hi and == favors "in range". So no change in lo <= x <= hi.

OP originally found examples of "invalid test result will move inside the range after the conversion" because the x was just outside [lo...hi] and the effective precision reduction now made x equal to either lo or hi. Since == favors "in range", the move from outside to inside is seen.

Note: if the LinearMapping() has a negative slope like -1, then lo <= x <= hi can easily be broken. as 1 <= 2 <= 3 --> -1 > -2 > -3. This makes lowerLimit > upperLimit which "in range" cannot be satisfied for any x.

---

For reference, OP's code simplified:

#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;
}

Output

Before 504 >  503 <= 1000
After  219 <= 219 <= 435

Before 500 <= 504 >  503
After  217 <= 219 <= 219

Answer 2

… 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?

Yes, this can occur, in theory, although due to C behavior rather than due to the underlying floating-point operations. The C standard does not guarantee that IEEE-754 floating-point arithmetic is used or even that the same precision through evaluation of an expression, and this can cause identical inputs to an expression to have different results.

Although LinearMapping is shown as a single routine, the compiler may inline it. That is, where the routine is called, the compiler may replace the call with the body of the routine. Furthermore, when it does this in different places, it may evaluate the expressions using different methods. So, in this code, LinearMapping may be evaluated using different floating-point operations in each call:

testResult = LinearMapping(testResult);
lowerLimit = LinearMapping(lowerLimit);
upperLimit = LinearMapping(upperLimit);

This means that (a * testResult + b) / c might be evaluated using only 32-bit floating-point arithmetic, while (a * upperLimit + b) / c might be evaluated using 64-bit floating-point arithmetic, converted to 32-bit after the division. (I have merged your three assignment statements into one expression here for brevity. The issue applies either way.)

One consequence of this is double rounding. When a result is calculated with one precision and then converted to another, two roundings occur, one in the initial calculation and a second in the conversion. Consider an exact mathematical result that looks like this:

    1.xxxxx01011111111xxxx1xx
             ^       ^ Start of bits to be rounded in wider format.
             | Start of bits to be rounded in narrower format.

If this were the result of a calculation in the narrower format, we would examine the bits 011111111xxxx1xx and round them down (they are less than ½ at the position where we are rounding), so the final result would be 1.xxxxx01. However, if we do the calculation in the wider format first, the bits to be rounded away are 1xxxx1xx (more than ½), and these are rounded up, making the intermediate result 1.xxxxx0110000000. When we convert to the narrower format, the bits to be rounded away are 10000000, which is exactly the midpoint (½), so the round-to-nearest-ties-to-even rule tells us to round up, which makes the final result 1.xxxxx10.

Thus, even if testResult and upperLimit are equal, the results of applying LinearMapping to them may be unequal, and testResult may appear to be outside the interval.

There may be some ways to avoid this problem:

• If your C implementation conforms to Annex F of the C standard (which essentially says that it uses IEEE-754 operations and binds them to C operators in expected ways) or at least conforms to certain parts of it, then double rounding should not occur with well written source code.

• The C standard says an implementation should define FLT_EVAL_METHOD in <float.h>. If FLT_EVAL_METHOD is 0, it means all floating-point operations and constants are evaluated with their nominal type. In this case, double rounding will not occur if you use a single floating-point type in your source code. If FLT_EVAL_METHOD is 1, float operations are evaluated with double. In this case, you could avoid double rounding by using double instead of float. If it is 2, operations are evaluated with long double, and you could avoid double rounding by using long double. If FLT_EVAL_METHOD is -1, the floating-point format used for evaluation is indeterminable, so double rounding is a concern.

• For specific values or intervals of values and/or known floating-point formats, it might be possible to prove that double rounding does not occur. E.g., given that your inputs are all int32_t, that the linear mapping parameters are certain values, and that only 32-bit or 64-bit IEEE-754 binary floating-point is in use, it might be possible to prove double rounding does not occur.

Even if your implementation conforms to Annex F or defines FLT_EVAL_METHOD to be non-negative, you must still take care in your code not to use expressions that have type double which you then assign to objects of type float. This will cause double rounding because your source code explicitly requires it, rather than because C is lax about floating-point.

As a specific example, consider (1.0 * 13546 + 1.0) / 2.3. If the floating-point constants are represented in 64-bit binary floating-point (53-bit significands) and the expression is evaluated in 64-bit binary, the result is 5890.0000000000009094947017729282379150390625. However, if the same constants are used (in 64-bit binary) but the expression is evaluated with Intel’s 80-bit binary (64-bit significands) and then converted to 64-bit, the result is 5890.

In this case, the exact mathematical quotient is:

1.01110000001000000000000000000000000000000000000000001000000000001011…
                                                          ^ Bits rounded away in double.

If we round this in double, we can see the bits to be rounded away, 1000000000001011…, are greater than ½ at the rounding position, so we round up. If we round it in long double, the bits to be rounded away are 01011…. These round down, leaving:

1.011100000010000000000000000000000000000000000000000010000000000
                                                      ^ Bits rounded away in double.

Now, when we round to double, the bits to be rounded away are 10000000000, which is the midpoint. The rule says to round to make the low bit even, so the result is:

1.011100000010000000000000000000000000000000000000000010000000000