?- set_prolog_flag(prefer_rationals, true),
X is 10r1^(-5).
X = 1r100000.
Please close the same bug report on Github.
?- set_prolog_flag(prefer_rationals, true),
X is 10r1^(-5).
X = 1r100000.
Please close the same bug report on Github.
I’m not quite sure what you would expect. It seems to me that :
[eclipse 2]: X is 2^(-5).
X = 0.03125
is incorrect. The answer should be a rational (1r32 like SWIP). If you want to make it a float, use the float function.
The float rounding flag only affects floating point operations. Rationals are precise so rounding is never required. Interval arithmetic is the same; if the bounds are precise, no rounding needed.
From rationalize/1 doc:
…the result is only accurate within the rounding error of floating point numbers …
But:
?- X is rational(pi-3).
X = 39854788871587r281474976710656.
?- X is (pi-3) - 39854788871587r281474976710656.
X = 0.0.
Remember that π is an irrational number, so can’t be precisely represented by either a rational or a float. But pi, the arithmetic function returns the floating point value closest to π, and that’s what is used in arithmetic evaluation. Neither SWIP or Eclipse, or pretty much any other (non-interval) arithmetic system, is correct in the mathematical sense.
Yeah. This should imply that for every (normal) float F the following should hold (I think):
F =:= rationalize(F).
Which for pi-3 is not the case. So, there is indeed something fishy. As @anon95304481 correctly states, the implementation is an iteration that stops if abs(P/Q-F) < episilon. That is highly dubious.
I tried P/Q == F using C float comparison as termination condition. For pi-3 that results in the same result as ECLiPSe, which seems promising. Unfortunately, some quick tests prove this doesn’t always terminate. I tried to remedy this by making the loop terminate on fixed point, but it seems possible that it never reaches fixed point ![]()
Now luckily ECLiPSe is also open source (using MPL) and although there are some indirections it wasn’t too hard to find the function that does the actual work. I didn’t have the time to figure out what really happens. It seems more complicated and doing a lot of work using GMP functions. Most likely it is a lot slower
Correctness seems important here though.
It seems to me that that difference is due to rational vs. rationalize:
?- X is pi-3, X =:= rationalize(X).
false.
?- X is pi-3, X =:= rational(X).
X = 0.14159265358979312.
Again, rationalize does not guarantee to produce an precisely equivalent rational number from a float. What am I missing?
Now that that’s clear (which it wasn’t to me before), I respectfully suggest you open another discussion rather than hijacking this one (see topic).
I think rationalize should produce the smallest N/D for which float(N/D) compares equal to the original float. Our implementation does not pass that test. The ECLiPSe one seems to.
I’m not a float expert, but it seems there are two issues. One of them is generating this rational number, where the usage of DBL_EPSILON as stop condition is (I think) incorrect, but then still,
?- 0.1 =:= 1r10.
Is true in ECLiPSe and false in SWI-Prolog. I guess this could have to do with the rounding mode?
Then doesn’t ECLiPSe have it wrong? 0.1 cannot be precisely represented by an IEEE floating point number due to the decimal to binary conversion. So this query should fail. But:
?- 0.5 =:= 1r2.
true.
It only fails if the rounding done by read (0.1) is different from the rounding done by converting 1r10. I guess it is desirable that these are the same and it seems ECLiPSe got this right. I may be wrong though, as I said before, I’m not an expert in this.
True enough, I guess. If both happen to generate the same incorrect float, the compare will work. And you’re probably right that the issue is caused by GMP’s mpq_set_d() using to_zero for conversion while the C round function (used by read?) obeys the current rounding mode (usually to_nearest). If Eclipse gets it right all the time, I’m guessing it doesn’t use GMP. Or at least not for this conversion.
I was already convinced that SWI-Prolog’s conversion between float and rational could be more precise. The current master has a copy of the ECLiPSe functions that do this. That seems to work fine in the sense that the results seem to be fully equivalent to ECLiPSe and in general better behaved.
I don’t think this is the end of the whole story. Generating random floats we see violations of
X =:= rational(X) as well as X =:= rationalize(X). I played with variations of this program:
rfloat(X) :-
% random_between(-308,308,E),
random_between(-10,10,E),
X is random_float * 10**E.
t(X) :-
assertion(X =:= rational(X)),
assertion(X =:= rationalize(X)),
% assertion(rational(X) =:= rationalize(X)),
true.
r(N) :-
forall(between(1,N,_),
( rfloat(F),
t(F))).
This seems to hold most of the time, as in ?- r(10 000). often succeeds. Going to a million you still
find violations that you can copy to ECLiPSe which also fails. As you go to the more extremely small or large integers the number of violations rapidly increases.
The whole story is not really clear to me at the moment. I don’t think performance is the main issue as I don’t see applications converting a lot of floats to rationals. After all, if you are in float space things are imprecise and fast. Why move to rationals then? The result is already imprecise and that doesn’t get fixed. ECliPSe rational to double conversion should respect float rounding and, as I read the comments, should perform better than GMP for really large numerator and denominator values. Also here, it isn’t really clear why one would do this conversion as it looses the carefully maintained exact result.
I think the above calls for as-precise-as-possible as primary target, offering performance.
But still, neither X =:= rational(X), nor X =:= rationalize(X) is true for all floats. It is now true for many more floats, but mathematically that is no progress
The issue is now whether it is possible to guarantee these relations, at which price and how important it is?
See this post
Example run:
102 ?- r(1 000 000).
ERROR: Assertion failed: user:(4.1200187083136876e-13=:=rationalize(4.1200187083136876e-13))
This one works for rational/1, but notably as you increase the range and run more tests you also find failing rational/1 tests. The failure is consistent with ECLiPSe (well, the about 5 I checked) using the current git version of SWI-Prolog.
No. X =:= rational(X) where X is a float compares a float to a rational, which converts the rational into a float before doing C == on the two double values.
The picture is getting a bit blurred. Quite quickly my test finds:
5.750009747987844e-308 =\= rational(5.750009747987844e-308)
The errors only seem to occur quite close to min/max float though.
The ECLiPSe 7.0 binary succeeds here, but when I compile ECLiPSe from source it behaves as SWI-Prolog (which is to be expected as they use now the same implementation for the relevant conversions).
You test on 8.1.22, while I test on the git version which uses the ECLiPSe rational->double conversion. This implementation is claimed to have two advantages over the GMP version: it is more efficient on really huge numerator/denominator and it respects the float rounding mode.
And I believe copy/pasted text just as well as screenshots.
Yes. The problem is with float/1 (or more precisely with the GNU MPQ --> double conversion). The ECLiPSe algorithm basically takes the two integers, if both are really long if strips bits from the right that exceed the double representation anyway, converts both to a double and divides.
This number has a denominator of 322 digits, so this creates Inf and the result is 0.0, which doesn’t match the input. The odd part is the the ECLiPSe binary I downloaded works fine, but the version I compiled from their source doesn’t (yes, both 7.0_53). Both versions agree that float(denominator) is Inf.
Right now I understand the source, but not why the distributed binary gets it right. Understanding the source also hints at the solution: use mpz_get_d_2exp() and I now have consistent results 
And what would ?- rational(X,2,4). do? As is, the system only has canonical rational numbers and thus == and = work as expected on them. It seems a bad idea to break this rule and I don’t like relying on the user to provide a canonical numerator/denoninator pair.
I’m not against making rational/3 multi-moded, although it was added as a type test providing some additional information similar to is_dict/2 in the days that rational numbers were terms. I can also imagine we deprecate rational/3. With atomic rationals the value seems small.
Thanks to some creative work from Joachim Schimpf, in the current GIT version we see for every float X that is not NaN or Inf:
X =:= rational(X)
X =:= rationalize(X).
This is also true for subnormal X, but if the flag float_underflow is error subnormal floats generate an evaluation_error(float_underflow) exception. The latter is still not set in stone.
Applaus for Joachim!
We could, but I’m more inclined to deprecate rational/3. It had its use as a quick way to decompose both integers and rationals, but with atomic rationals that is probably less needed. The numerator and denominator functions are pretty efficient and rdiv or / are pretty good at creating rationals. I have little clue whether there are applications would really benefit from a rational/3. I guess one of the few arguments might be to allow for non-canonical rationals. So far I’m not inclined to support these.
This at least suggests differently. It is actually a bit odd as it does move the sign, but doesn’t divide by 2.
Above that there is more explicitly:
True. You have to use the GIT source. That also applies to ECLiPSe. Well, the binary release is surely not good enough. I don’t know how quickly ECLiPSe updates the public repository and provides releases. You know for a long time that you get the latest source for SWI-Prolog from GIT.
A simple random test generator generates zillions of counter examples on the old versions. It passed 10M random tests on the current sources quite a couple of times, so I think about 100M numbers have been tried.