Rounding error in (**)/2

[j4n_bur53]

@abramobagnara

(see C11 7.1.3p1 and J.2 item 106).

The define is here, what does C11 7.1.3p1 mean, some link?

[abramobagnara]

@j4n_bur53

The define is here, what does C11 7.1.3p1 mean, some link?

n1570.pdf

See chapter 7.1.3 paragraph 1

[jan]

Surely parts of the arithmetic implementation are a mess. That has been discussed between various people. Part of this is hard to avoid because the dynamic typing causes quite a few code paths and the fact that certain features may be there or may be lacking (GMP) causes yet more complicates and finally there are conditions dealing with old broken C libraries as well as current bugs Still, a few helper functions could improve this

The cr_() stuff got there because @ridgeworks at some point wanted to include a “correct rounding C library” and I wanted that to touch the existing code as little as possible. It is probably better to use cr_() in the body of the arithmetic and bind these symbols depending if and how we wish to do correct rounding.

Anyway, I have more pressing things to do and anyone is invited to improve the code and send a PR

[ridgeworks]

@jan

Sounds good. Please note that I did merge the bigint-to-float branch. I think that is still ok as that includes the mpq_to_double() patches we want anyway. Your code should add a fast-path implementation for mpz_to_double() that now uses mpq_to_double().

Almost. There are a couple of additional patches to mpq_to_double; one to avoid over-rounding in the FE_UPWARD and FE_DOWNWARD rounding modes, and one to ensure that infinity is returned the result will overflow a 64 bit float.

The new function is named mpz_float which is called from mpz_to_double rather than mpz_fdiv with a divisor argument of MPZ_ONE.

Changes are confined to gmp.c as described and Tests/rational/test_ieee_754.pl which will have some additional rounding tests cases for floating various “types” of integers.

[ridgeworks]

@abramobagnara

This kind of #define should always be avoided as it is undefined behavior according to C standard (see C11 7.1.3p1 and J.2 item 106).

IMO the C standard is deficient because it doesn’t guarantee correct rounding as per IEEE754 of results from pow (and lots of other library functions). This “local” re-defintion of pow only ensures that the results are rounded (sometimes over-rounded) in the defined direction, e.g., a “round to positive” result will always be equal to or greater than the precise value.

But, as I keep claiming, I’m not proficient in C programming and don’t particularly aspire to be, so if there’s a better way of doing it, have a go.

Aside: There is at least one partial correctly rounded math library, but there are issues incorporating it into the SWIP build process IIRC.

[abramobagnara]

It is just enough to avoid to use pow name for your macro and use e.g. rounding_pow for both macro definitions and their uses.

Of course if in some configuration you want that rounding_pow uses standard library pow you can write in the related #if branch #define rounding_pow pow

[ridgeworks]

@j4n_bur53

it seems to me, there would be at least 5 test cases,
to arrive at 100% code coverage.

More, I think. Additional cases to cover overflow (to infinity) , negative integers, and the other three rounding modes. (You get “round to zero” for free because that’s what mpz_get_d already does.)

But it does depend on what your requirements are.

[ridgeworks]

@abramobagnara

It is just enough to avoid to use pow name for your macro and use e.g. rounding_pow for both macro definitions and their uses.

But that involves changing all the uses, and for all the relevant math library functions, that’s not insignificant. As structured, it’s dead simple to revert back to the original definitions just by commenting out all the macro definitions, which I think was Jan’s intent when this work was being done.

While this may be considered a bit unsanitary in general, given that the macros (and uses) are confined to a single file and it’s an all or nothing proposition, I think it’s a rather neat pragmatic solution. But I’m not a C programmer; I just have to program in C on occasion.

[j4n_bur53]

Since the new release 8.5.18 of SWI-Prolog arrived, and
to complete my creation of fuzzer test cases, I also did some
new test cases for (**)/2. I use relatively few bits in the

arguments (12 bits base and 6 bits exponent), but the result
has more bits than a usual float (12*2^6 = 768 bits), which
are then rounded. Unlike for div, the ISO core standard is less

specific what to do for pow. So its just some discrepancy
comparison, this time even more. Now that SWI-Prolog claims
doing things exact, there might be nevertheless a small bug,

when the exponent is negative, i.e. some glitch in the
reciprocal. For example I find the following failure:

/* SWI-Prolog 8.5.18 Windows */
?- case6(N, P, Q, Y), Y =\= P**Q.
N = 6,
P = 563,
Q = -56,
Y = 9.365476712718912e-155 .

SWI-Prolog doesn’t want 9.365476712718912e-155 . It prefers
9.365476712718935e-155 . But does this have some inner consistency?
Rather not I guess. The more consistent value is the one from the

test case if we do the below computation. So this could be some bug:

/* SWI-Prolog 8.5.18 Windows */
?- current_prolog_flag(iso, X).
X = true.

?- X is 563** -56.
X = 9.365476712718935e-155.

?- X is 1/(563**56).
X = 9.365476712718912e-155.

The computations seem not to depend on the iso flag value:

?- set_prolog_flag(iso, false).
true.

?- X is 563** -56.
X = 9.365476712718935e-155.

?- X is 1/(563**56).
X = 9.365476712718912e-155.

[jan]

@j4n_bur53

So this could be some bug:

The precision of the MinGW runtime isn’t very good as we have seen various times. On Linux the result is consistent (X = 9.365476712718912e-155). Maybe, after all, we should go for the alternative math library @ridgeworks once looked at. That at least would make all results consistent and platform independent.

I’m seriously wondering how many people care about the precision of floating point operations (assuming the error is not more than a couple of ulps). Is that worth the additional size of the sources and binaries and the additional complexity of building the system?

[ridgeworks]

@jan

I’m seriously wondering how many people care about the precision of floating point operations (assuming the error is not more than a couple of ulps).

I strongly suspect very few people care.Otherwise, there would be correct, precise C math libraries in common usage.

What I care about is that the float rounding modes correctly bound any precise real values. The underlying assumption is that compiled arithmetic is IEEE compliant for all rounding modes and the functions in math libraries are correct within 1 ulp in default “to_nearest” rounding mode. Largely that seems to be true, but rigorous testing is hard and there are no guarantees.

@j4n_bur53

But I am not publicitly testing intervals, I rather leave that to the interval owner.

Most likely buggy intervals will affect clpBNR, but again, also not my roadworks.

Hence my concern. But interval arithmetic is not built into SWIP, so it’s not Jan’s “roadworks” either. It does require IEEE compliance as noted above, unless you’re building a mixed Prolog and C implementation.

[ridgeworks]

@j4n_bur53

So this form of intervals reuses some of the existing Prolog arithmetic.
But if you or Jan W. don’t find floats interesting, no need to respond to this
thread. I wouldn’t care less, really. Same for other “people”,

Please don’t over-react. Interval arithmetic is entirely based on the underlying arithmetic support. Floats, and the flags that control them, matter (as do rationals). If they didn’t, I wouldn’t bother participating in this thread (as you suggest).

My only point was that intervals are constructed on top of the underlying Prolog arithmetic machinery so (perhaps) they are needn’t be part of this discussion other than making sure other rounding modes are supported.

[abramobagnara]

@jan

I’m seriously wondering how many people care about the precision of floating point operations (assuming the error is not more than a couple of ulps). Is that worth the additional size of the sources and binaries and the additional complexity of building the system?

Precision and monotonicity is fundamental for interval arithmetic, otherwise we risk that an operation on two non empty intervals gives as result an empty intervals, that is clearly unacceptable.

[j4n_bur53]

Relax, don’t worry. Somebody wrote:

@ridgeworks

Reimplemented “to_nearest” rounding and it now passes all 100
test cases in cases.p plus the failed case reported. Any additional
test data is welcome.

But maybe this somebody didn’t expect that after testing float/1, I
will also do some (/)/2 and (**)/2 testing. But this is only natural,
after unit testing float/1, to also make some integration testing

of (/)/2 and (**)/2 that depend on float/1. But thats mostlikely
alread it. I already wrote where I intend to stop:

@j4n_bur53

But I am not publicitly testing intervals, I rather leave that to the
interval owner. Most likely buggy intervals will affect clpBNR, but
again, also not my roadworks.

Hope this Helps!

Edit 15.10.2022
I made a small attempt to also test (+)/2, via this test case here:

@j4n_bur53

/* SWI-Prolog 8.5.18 */
?- X is 665723957508124892951320985600.0+26388279066624.0.
X = 6.65723957508125e+29.

It was even a test case totally inside floats. But then the discrepancy
among systems was caused by a parsing problem. If I find a way to
do something nasty, which is not parsing, I might also

do some new (+)/2 test cases and some new (*)/2 test cases. Or
somebody else does this. Or there are already some test cases.
Who knows what the future will bring?

[j4n_bur53]

Now I found a test case, where the Unix platform also fails.
Was increasing the bits in the case6/3 fuzzer:

pow2(B, E, X) :- E < 0, integer(E), !, X is 1/(B^(-E)).
pow2(B, E, X) :- X is B^E.

?- current_prolog_flag(iso, X).
X = false.

/* SWI-Prolog 8.5.20 WLS2 Nok */
?- repeat, L is -(1<<15), H is (1<<15)+1,
M is -(1<<3), J is (1<<3)+1, random(L, H, P),
 random(M, J, Q), pow2(P, Q, A), B is P**Q, A=\=B.
P = -18573,
Q = -7,
A = -1.3116730299820406e-30,
B = -1.3116730299820408e-30 .

I have SWIPL on Unix platform only since weekend, could not do this
testing beforehand, so I guess I need to update case6/3 test cases,
because they showed Unix platform completely passing, which

seems to be wrong. I used in the old fuzzer 12 bits + 6 bits, the
above uses 15 bits + 3 bits, and new nasty test cases pop up.

Edit 25.10.2022
Its a test case where JDK 19 non-strict also fails. Cool!
I couldn’t believe that it is that good. It was too good to
be true. So its also not that good!

/* Jekejeke Prolog 1.5.5, JDK 19 non-strict Nok */
?- X is -18573** -7.
X = -1.3116730299820408E-30.

[ridgeworks]

@j4n_bur53

Now I found a test case, where the Unix platform also fails.

This seems to be a test case dependant on when the conversion to float occurs. If the prefer_rationals flag is false:

?- set_prolog_flag(prefer_rationals,false).
true.

?- X is 2** -2.
X = 0.25.

?- P = -18573, Q = -7, R is P**Q.
P = -18573,
Q = -7,
R = -1.3116730299820408e-30.

?- P = -18573, Q = -7, R is 1/(P**(-Q)).
P = -18573,
Q = -7,
R = -1.3116730299820406e-30.

but if true:

?- set_prolog_flag(prefer_rationals,true).
true.

?- X is 2** -2.
X = 1r4.

?- P = -18573, Q = -7, R is P**Q.
P = -18573,
Q = -7,
R = -1.3116730299820408e-30.

?- P = -18573, Q = -7, R is 1/(P**(-Q)).
P = -18573,
Q = -7,
R = -1.3116730299820408e-30.

My interpretation: The early float conversion introduces a rounding error which is propagated during division resulting in a 1 ulp difference.

[j4n_bur53]

@j4n_bur53

maybe if explicit float/1 conversion is used? If (^)/2 is
pure, doesn’t touch the float x float signature.

I don’t know whether SWI-Prolog 8.5.20 is pure? I don’t
know what float(X)**Y does. It could do at least two things:

• Variant 1: If X happens to be a integral float
  float(X)**Y = float(integer(X)^Y)

• Variant 2: No need to check whether X is integral float
  float(X)**Y = Math.pow(float(X), float(Y))

I don’t know what it does. My test cases for (**)/2 do
all have integer arguments. I didn’t test float(X)**Y so far.
I am only at the beginning of (**)/2 testing. I didn’t test

much of (**)/2. And most likely I will also stop at integer
arguments, because with non-integer arguments I would
just open another never ending death march of test cases?

The test case e**0.8618974796837966 was out of band,
because some floating point value constants dscussion.

[jan]

@j4n_bur53

I don’t know whether SWI-Prolog 8.5.20 is pure? I don’t
know what float(X)**Y does.

The float function is defined in ISO. It converts a non-float (read rational) to a float according to the active rounding mode. **/2 does what +/2, */2, etc. do: if one of the arguments is float all the others are converted to float and the operation is executed using libm. If both are rationals (that includes integers), and the function is well defined on rationals, perform the operation using rational numbers.

At least, that is how it supposed to work. For compatibility reasons prefer_rationals is false and integer division that does not divide exactly as well as exponentiation with negative exponents and integer arguments perform float arithmetic.

I’m still wondering whether prefer_rationals should be the default. As Richard O’Keefe argued long ago, this is what ISO should have decided from the beginning. ISO arithmetic is based on a rather outdated and simplified model of C+libm practices with as only extension that functions should return an overflow exception rather than an incorrect result (which is not implemented by several Prolog’s: GNU-Prolog 1.4.5: A is 10^300 → A = 0, XSB 3.8.0 → A = -9223372036854775808).

Very little breaks by preferring rationals, but some computations may result in huge numbers and get really slow. SWI-Prolog allows for setting a limit on the (byte) size of rational numbers and define what to do when this is exceeded (error or convert to float). That solves this problem, but makes the result a little hard to predict.

[j4n_bur53]

@ridgeworks

?- P = -18573, Q = -7, R is 1/(P**(-Q)).
P = -18573,
Q = -7,
R = -1.3116730299820406e-30.

?- set_prolog_flag(prefer_rationals,true).
true.

?- P = -18573, Q = -7, R is 1/(P**(-Q)).
P = -18573,
Q = -7,
R = -1.3116730299820408e-30.

So there is a bug in long division? Did you use iso=false?
What platform did you use?

Edit 25.10.2022
I just figured out that there is indeed this bug, that (/)/2 is
buggy, and that I might need more (/)/2 test cases as well,
since they passed all, but now we have one more failing.

With my own exactpow/3, which uses an emulated exact
long division I also get a different result from SWI-Prolog (/)/2,
namely I get this result and JDK 19 non-strict is correct.

?- exactpow(-18573, -7, X).
X = -1.3116730299820408E-30.

Problem is I only tested (/)/2 with iso=true, and not with
iso=false, thats why of course I could not detect some
problems with (/)/2 that fall into iso=false.

Maybe Unix has nevertheless a quite good (**)/2?