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.

2

Happens not always:

/* SWI-Prolog 8.5.8 */
?- X is nexttoward(0.5,0), Y is X+0.5.
X = 0.49999999999999994,
Y = 1.0.

?- X is nexttoward(1.5,0), Y is X+0.5.
X = 1.4999999999999998,
Y = 1.9999999999999998.

It has to do that when you add 0.5 to X in the first query, it
needs shift alignment, the floating point number cannot keep
the same exponent to perform the addition, and then it has to

round looking at the last place.

Edit 23.03.2022:
But its quite amazing that Math.round() does indeed do it more
precisely, thats the result from Java, JavaScript and Python:

/* Other Systems */
?- X is 0.49999999999999994, Y is X+0.5.
X = 0.49999999999999994, Y = 1.0.

?- X is 0.49999999999999994, Y is round(X).
X = 0.49999999999999994, Y = 0.0.

ECLiPSe Prolog and TauProlog also agree. But TauProlog
gives the result as integer instead as float.

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.

4

Thats a again a can of worms. Just out of curiousity I was
checking out the ISO core standard, and it tells me:

/* 9.1.6.1 Floating point to integer rounding functions */
round(x) := floor(x + 1/2)

So I made some additional tests, first Scryer Prolog v0.8.123:

?- X is round(0.5).
   X = 1.
?- X is round(-0.5).
   X = 0.

And then SWI-Prolog 8.5.8:

?- X is round(0.5).
X = 1.
?- X is round(-0.5).
X = -1.

Do not misunderstand me. I am not preaching some orthodoxy.
I just wonder again whether there is a flag or something so that
this can be queried or controlled.

For example my definition of round function usually says,
and what SWI-Prolog possibly implemented so far:

/* Alternative to 9.1.6.1 ISO core standard */
round(x) = truncate(x + sign(x)*1/2)

And not what the ISO core standard gives.

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: