6

I am porting C code to Delphi and find an issue in the way the compilers (Delphi 10.4.1 and MSVC2019, both targeting x32 platform) handle comparison of +NAN to zero. Both compilers use IEEE754 representation for double floating point values. I found the issue because the C-Code I port to Delphi is delivered with a bunch of data to validate the code correctness.

The original source code is complex but I was able to produce a minimal reproducible application in both Delphi and C.

C-Code:

#include <stdio.h>
#include <math.h>

double AngRound(double x) {
    const double z = 1 / (double)(16);
    volatile double y;
    if (x == 0) 
        return 0;
    y = fabs(x);
    /* The compiler mustn't "simplify" z - (z - y) to y */
    if (y < z)
        y = z - (z - y);      // <= This line is *NOT* executed
    if (x < 0)
        return -y;
    else
        return y;             // <= This line is executed
}

union {
    double d;
    int bits[2];
} u;


int main()
{
    double lon12;
    double ar;
    int    lonsign;

    // Create a +NAN number IEEE754
    u.bits[0] = 0;
    u.bits[1] = 0x7ff80000;

    lon12    = u.d;                // Debugger shows lon12 is +nan
    if (lon12 >= 0)
        lonsign = 1;
    else
        lonsign = -1;              // <= This line is executed
    // Now lonsign is -1

    ar = AngRound(lon12);
    // Now ar is +nan

    lon12 = lonsign * ar;
    // Now lon12 is +nan
}

Delphi code:

program NotANumberTest;

{$APPTYPE CONSOLE}

{$R *.res}

uses
  System.SysUtils;

type
    TRec = record
       case t : Boolean of
       TRUE:  (d    : Double);
       FALSE: (bits : array [0..1] of Cardinal);
    end;

function AngRound(x : Double) : Double;
const
    z : Double = 1 / Double(16);
var
    y : Double;
begin
    if x = 0 then
        Result := 0
    else begin
        y := abs(x);
        if y < z then
            // The compiler mustn't "simplify" z - (z - y) to y
            y := z - (z - y);           // <= This line is executed
        if x < 0 then
            Result := -y                // <= This line is executed
        else
            Result := y;
    end;
end;

var
    u       : TRec;
    lon12   : Double;
    lonsign : Integer;
    ar      : Double;
begin
    // Create a +NAN number IEEE754
    u.bits[0] := 0;
    u.bits[1] := $7ff80000;

    lon12 := u.d;                       // Debugger shows lon12 is +NAN
    if lon12 >= 0 then
        lonsign := 1                    // <= This line is executed
    else
        lonsign := -1;
    // Now lonsign is +1

    ar := AngRound(lon12);
    // Now ar is -NAN

    lon12 := lonsign * ar;
    // Now lon12 is -NAN
end.

I have marked the lines with are executed after a comparison. Delphi evaluate (lon12 >= 0) to TRUE when lon12 variable equal +NAN. MSVC evaluate (lon12 >= 0) to FALSE when lon12 variable equal +NAN.

lonsign has different values in C and Delphi.

AngRound receiving +NAN as argument return different values.

Final value of lon12 is (fatally) different.

Machine code generated by the compilers are different:

Delphi generated machine code:

Delphi generated machinecode

MSVC2019 generated machine code:

MSVC2019 generated machine code

The comparison result seems more logical in Delphi: (lon12 >= 0) is TRUE when lon12 is +NAN. Does this means the bug is in MSVC2019 compiler? Should I consider the test data set of the original C-Code carry in error?

fpiette
  • 11,983
  • 1
  • 24
  • 46
  • See [What would cause the C/C++ <, <=, and == operators to return true if either argument is NaN?](https://stackoverflow.com/a/23642563/576719) for why the MSVC compiler can use other operators with a compiler switch working with NaN. – LU RD Feb 01 '21 at 11:50
  • 1
    A note about the comment that the compiler must not simplify `z - (z - y)` to `y`. To do that in C, use `z - (double) (z - y)`. Per the C standard, the cast forces the compiler to “discard” excess precision, meaning it cannot treat the expression as an infinitely precise mathematical `z - (z - y)` that equals `y`; it must produce a precision-limited `double` result for `z - y` and then subtract that from `z`. (However, there could still be double-rounding effects if the C implementation has some `long double` type with more precision.) – Eric Postpischil Feb 01 '21 at 11:53
  • @EricPostpischil Thanks for the note. This C-Code is not mine. FYI, the original code is there: https://geographiclib.sourceforge.io/html/C/. The author provides test data and when his code is compiled with MSVC2019, all test are passed. Of course this doesn't mean there is no bugs. – fpiette Feb 01 '21 at 13:01

1 Answers1

11

First of all, your Delphi program does not behave as you describe, at least on the Delphi version readily available to me, XE7. When your program is run, an invalid operation floating point exception is raised. I'm going to assume that you have actually masked floating point exceptions.

Update: It turns out that at some time between XE7 and 10.3, Delphi 32 bit codegen switched from fcom to fucom which explains why XE7 sets the IA floating point exception, but 10.3 does not.

Your Delphi code is very far from minimal. Let's try to make a truly minimal example. And let's look at other comparison operators.

{$APPTYPE CONSOLE}

uses
  System.Math;

var
  d: Double;
begin
  SetFPUExceptionMask(exAllArithmeticExceptions);
  SetSSEExceptionMask(exAllArithmeticExceptions);
  d := NaN;
  Writeln(d > 0);
  Writeln(d >= 0);
  Writeln(d < 0);
  Writeln(d <= 0);
  Writeln(d = d);
  Writeln(d <> d);
end.

Under 32 bit in XE7, this outputs

TRUE
TRUE
FALSE
FALSE
TRUE
FALSE

Under 32 bit in 10.3.3 (and 10.4.1 as you report in a comment below), this outputs

TRUE
TRUE
TRUE
TRUE
FALSE
TRUE

Under 64 bit in XE7 and 10.3.3 (and 10.4.1 as your report), this outputs

FALSE
FALSE
FALSE
FALSE
FALSE
TRUE

The 64 bit output is correct. The 32 bit output for both variants are incorrect. This we we can see by referring to What is the rationale for all comparisons returning false for IEEE754 NaN values?

all comparisons with the operators ==, <=, >=, <, > where one or both values is NaN returns false, contrary to the behaviour of all other values.

For your 32 bit Delphi code, you will need to workaround this bug and include special case code whenever it needs to handle such comparisons. Unless of course, by some happy chance, that you are using 10.4 and it already fixes the issue.

David Heffernan
  • 601,492
  • 42
  • 1,072
  • 1,490
  • My Delphi code is written to be as close as possible to the C-Code and show what is going wrong. An no, I didn't turned off any floating point exception and it does NOT trigger any exception (Delphi 10.4.1 and compile for 32 bits). And your code is NOT the same case: Delphi "NaN" is NOT +NaN, it is -NaN. Delphi NaN binary representation in hex is $FFF8000000000000 and my issue is related to -NaN wich is $7FF8000000000000. – fpiette Feb 01 '21 at 12:46
  • NaNs are not signed. Since you've got 10.4, what output does the program in my answer produce? – David Heffernan Feb 01 '21 at 13:00
  • Your code compiled as 32 bits produce TRUE TRUE TRUE TRUE FALSE TRUE and do not trigger any exception (without disabling it). For x64, it produce the exception so disabling it are required. The result is FALSE FALSE FALSE FALSE, TRUE. – fpiette Feb 01 '21 at 13:05
  • So now you believe me that 32 bit Delphi is the issue here? – David Heffernan Feb 01 '21 at 13:09
  • 1
    Note that NaN can have several internal representations (as explained [here](https://en.wikipedia.org/wiki/NaN#Floating_point)). However *"applications are not required to provide distinct semantics for those distinct NaN values"*. The sign bit is *"most often ignored in applications"*. – Olivier Feb 01 '21 at 13:39
  • 1
    Re “NaNs are not signed”: Do you mean in Delphi or more generally in IEEE 754? IEEE 754 does have a sign bit for NaNs (and there is actually some discussion going on right now in the IEEE-754 standard mailing list about its handling). – Eric Postpischil Feb 01 '21 at 14:25
  • @EricPostpischil Well sure they have a sign bit. But are there any operations where the sign bit is used? It's not relevant for comparison operators for sure. – David Heffernan Feb 01 '21 at 14:35
  • 4
    @DavidHeffernan: The standard does not specify the interpretation of the sign bit of a NaN, but it does specify its handling even for NaNs. The `copy`, `negate`, and `abs` operations copy, flip, and clear the sign bit. There is even a `copySign` bit to copy it from one operand and paste it onto another and an `isSignMinus` to test it. The `totalOrder` predicate is affected by the sign bit of a NaN. – Eric Postpischil Feb 01 '21 at 14:44