Another floating point casualty: round/integer arithmetic functions

Help! Predicate
1

SWIP Version 8.4.1:

?- X is nexttoward(0.5,-inf), N is round(X).
X = 0.49999999999999994,
N = 1.

?- X is nexttoward(0.5,-inf), N is integer(X).
X = 0.49999999999999994,
N = 1.

?- X is nexttoward(-0.5,inf), N is round(X).
X = -0.49999999999999994,
N = -1.

Looking at pl-arith.c source, function ar_integer, the problem seems to lie here:

	{ r->value.i = (int64_t)(n1->value.f + 0.5);
	  if ( r->value.i < 0 )		/* Why can this happen? */
	    r->value.i = PLMAXINT;
	} else
	{ r->value.i = (int64_t)(n1->value.f - 0.5);
	  if ( r->value.i > 0 )
	    r->value.i = PLMININT;
	}

If I execute the following bit of C code in an online service (C Online Compiler & Interpreter - Replit):

    double xpos = nexttoward( 0.5,-1);
    double xneg = nexttoward(-0.5, 1);
    printf("xpos = %24.20e   xpos+0.5 = %24.20e\n",xpos, xpos+0.5);
    printf("xneg = %24.20e   xneg-0.5 = %24.20e\n",xneg, xneg-0.5);
    printf("xpos = %24.20e   lround(xpos) = %li\n",xpos,lround(xpos));
    printf("xneg = %24.20e   lround(xneg) = %li\n",xneg,lround(xneg));

it prints

xpos = 4.99999999999999944489e-01   xpos+0.5 = 1.00000000000000000000e+00
xneg = -4.99999999999999944489e-01   xneg-0.5 = -1.00000000000000000000e+00
xpos = 4.99999999999999944489e-01   lround(xpos) = 0
xneg = -4.99999999999999944489e-01   lround(xneg) = 0

so it looks like an issue with the C level arithmetic. Function lround seems to get it right so a quick fix might be to use that instead adding/subtracting 0.5.

3

Updated that. Just, lround(9223372036854775807.0) produces a large negative integer :frowning: So, I had to use the same kludge as before, if rounding a positive float returns a negative integer, use PLMAXINT and similar for the negative side. Why is this the case? lround() is supposed to raise a float exception if it cannot do its work. Is this a bug?

Added some more tests.

5

As the docs say, SWI-Prolog rounds outwards, as most languages seem to do. Possibly the ISO flag should change this. GNU-Prolog (1.4.5) is rather strange:

?- A is round(0.5).
A = 0
?- A is round(-0.5).
A = 0
?- A is round(1.5).
A is 2.

The first seems a bit strange to me as binary floats can represent 0.5 precisely, no? ECLiPSe returns round(0.5) as 0.0, which may give a clue why this happens?

7

Looks like GNU-Prolog is rounding to even, which I thought is the most common default rounding mode.

8

The code says it does the equivalent of (long)rint(val). Linux docs for rint() say:

… using the current rounding direction (see fesetround(3)) and without raising the inexact exception. When the current rounding direction is to nearest, these functions round halfway cases to the even integer in accordance with IEEE-754.

It seems there are many opinions :slight_smile:

9

It appears MinGW (9.2.1) has a broken llround(), so using llround() is still not perfect :frowning: See MinGW - Minimalist GNU for Windows / Patches / #370 lround and llround functions are incorrect This was reported in 2008 … Anyway, a number of the IEEE-754 tests fail as well :frowning: