Bug fix for cmpReals - mixed compare on negative values

[ridgeworks]

Note: unable to generate a PR on recent clone due to outstanding PR #1144.

File: pl-gmp.c, Function: cmp_i_f(int64_t i1, double d2)

currently has a bug which compare two negative values as equal when they are close in value, e.g.,

?- C is cmpr(-1,-0.999999999).
C = 0.

Buggy code:

/* See https://stackoverflow.com/questions/58734034/ */
static int cmp_i_f(int64_t i1, double d2)
{ int64_t i1_lo = i1 & 0x00000000FFFFFFFF;
  int64_t i1_hi = i1 & 0xFFFFFFFF00000000;

  return cmp_f_f((double)i1_lo, d2-(double)i1_hi);
}

should be replaced by:

/* See https://stackoverflow.com/questions/58734034/ */
static int cmp_i_f(int64_t i1, double d2)
{ int invert;                  // invert result when both args negative
  // i1 and d2 have different signs? 
  switch (signbit(d2))         // 1 if negative else 0
  { case 0: if (i1 <  0) return CMP_LESS;    else invert =  1; break;  // if i1 < 0
    case 1: if (i1 >= 0) return CMP_GREATER; else invert = -1;         // if i1 >= 0 
  }
  // i1 and d2 have same sign, compare magnitudes and invert result when both negative
  int64_t abs_i1 = abs(i1);
  return cmp_f_f((double)(abs_i1 & 0x00000000FFFFFFFF), fabs(d2)-(double)(abs_i1 & 0xFFFFFFFF00000000)) * invert;
}

or equivalent.

[JanWielemaker]

I think there is still a problem if i1 is the largest negative integer as in that case abs(i1) is wrong, no?

[ridgeworks]

I think there is still a problem if i1 is the largest negative integer as in that case abs(i1) is wrong, no?

Good catch.

Since at that magnitude of integers, floating point values are very sparse, comparing with LONG_MIN+1 is equivalent:

?- A is float(-9223372036854775808).  % LONG_MIN
A = -9.223372036854776e+18.

?- A is float(-9223372036854775807).  % LONG_MIN+1
A = -9.223372036854776e+18.

?- A is float(-9223372036854775800).  % LONG_MIN+8
A = -9.223372036854776e+18.

so how about:

/* See https://stackoverflow.com/questions/58734034/ */
static int cmp_i_f(int64_t i1, double d2)
{ int invert;                  // invert result when both args negative
  // i1 and d2 have different signs? 
  switch (signbit(d2))         // 1 if negative else 0
  { case 0: if (i1 <  0) return CMP_LESS;    else invert =  1; break;  // if i1 < 0
    case 1: if (i1 >= 0) return CMP_GREATER; else invert = -1;         // if i1 >= 0 
  }
  // i1 and d2 have same sign, compare magnitudes and invert result when both negative
  int64_t abs_i1 = llabs((i1 == LONG_MIN) ? i1+1 : i1);  // corner case when i1 == LONG MIN, cmp same as i1+1
  return cmp_f_f((double)(abs_i1 & 0x00000000FFFFFFFF), fabs(d2)-(double)(abs_i1 & 0xFFFFFFFF00000000)) * invert;
}

[JanWielemaker]

Trying to simplify things, I came to

static int
cmp_i_f(int64_t i1, double d2)
{ if ( isnan(d2) )
  { return CMP_NOTEQ;
  } else
  { double d1_lo, d1_hi;
#define TOD(i) ((double)((int64_t)(i)))
    if ( i1 >= 0 )
    { d1_lo = TOD(i1 & 0x00000000FFFFFFFF);
      d1_hi = TOD(i1 & 0xFFFFFFFF00000000);
    } else
    { d1_lo = TOD(i1 | 0xFFFFFFFF00000000);
      d1_hi = TOD(i1 | 0x00000000FFFFFFFF)+1.0;
    }
#undef TOD
    return SCALAR_TO_CMP(d1_lo, d2-d1_hi);
  }
}

Testing using this works fine until M = 53. Above that we start to see errors. That looks wrong to me. Am I right?

t(N, M) :-
    Min is -(1<<M),
    Max is 1<<M,
    forall(between(1,N,_),
           ( random_between(Min, Max, I),
             tcmpr(I)
           )).

tcmpr(Expr) :-
    I is Expr,
    PInf is float(1<<1000),
    NInf is -PInf,
    F is float(I),
    N is nexttoward(F, NInf),
    P is nexttoward(F, PInf),
    assertion(-1 =:= cmpr(I, P)),
    assertion(0 =:= cmpr(I, F)),
    assertion(1 =:= cmpr(I, N)).

[JanWielemaker]

Ok. I think this is to be expected. Floats have 54 bits precision. Above that, float(I) and I may no longer compare equal. That means we should be safe.

[ridgeworks]

A couple of issues with the test:

1. float(1<<1000) is less than floating point infinity. Why not just use N is nexttoward(D,-inf) and P is nexttoward(D,inf) respectively` (or shift by a value greater than 1023).
2. Once M is greater than a certain magnitude, assertion(0 =:= cmpr(I, F)) will start failing due to roundoff errors in float(I). (largish ints cannot be precisely represented by a floating point value). That's one issue (among others) that cmpr is trying to address.

(Oops, too late).

[JanWielemaker]

Why not just use N is nexttoward(D,-inf)

I didn't want to change the Prolog flags and I only wanted to test 64 bit ints, so these are big enough.

[ridgeworks]

I tried replacing my version with the one you committed and got the following build error:

Undefined symbols for architecture x86_64:
  "_SCALAR_TO_CMP", referenced from:
      _cmpReals in pl-gmp.c.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)
ninja: build stopped: subcommand failed.

Maybe a platform issue? Or an old clone?

[JanWielemaker]

I tried replacing my version with the one you committed and got the following build error:

Undefined symbols for architecture x86_64:
  "_SCALAR_TO_CMP", referenced from:
      _cmpReals in pl-gmp.c.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)
ninja: build stopped: subcommand failed.

Maybe a platform issue? Or an old clone?

No. It is defined in the current HEAD. The macro provides an efficient way to map the comparison of two scalars to -1, 0, 1.

[ridgeworks]

So is that something that the two functions immediately above cmp_i_f (cmp_i_i and cmp_f_f) should use as well?

[JanWielemaker]

Ideally, yes. There is not a big hurry. The final code is a little shorter and faster, but barely noticeable.

[ridgeworks]

Agreed - more of a code clean up (for ease of future maintenance).