Proper rational numbers (prototype for testing) (Discussion)

Great! Then open/4 will do automagically resolve $HOME ?
In path expressions, also the dot should regain the binary operator status, for instance ?- edit(file(my_file.pl)) would be nice. Really, just aesthetic preferences … worthless, so sorry for the noise :slight_smile: .

Wow. A lot to catch up on.

I have to admit I’m waffling. I fear that the Logtalk usage of / as a non-arithmetic separator is the tip of the iceberg. There aren’t many characters suitable for generic separators. Others include - : . | \ and most of those have limitations in SWI Prolog syntax (e.g., : would be nice for type systems but is used for module qualification, . for field expressions but used for dictionaries, | in list syntax, etc.) So I have some sympathy for not further restricting the few nice separator characters that are left.

I also tend to think that a single syntax for a rational number, used for both input and output, would be preferable, rather than a variety under configuration control (Eclipse compatibility being an exception). And I’d also like to see something minimal, like the single infix character under discussion.

With that in mind, I’ve started to appreciate the PrQ syntax used in J. J was co-designed by Kenneth Iverson as a successor to APL (without the extended font requirements), so I think it has some credibility even if it isn’t widely used.

PrQ avoids all the issues with pre-existing operators and makes the parsing of things like 1r2.0 a bit more obvious (possibly .(1r2,0) depending on the parsing of . expressions). The use of lowercase characters in numbers has been established (non-decimal radix). Like the infix separator being discussed, it avoids the multi-token lookahead that the Ruby syntax entails.

Given a clean slate, i would still prefer 1/2 but I don’t think we have that luxury.

1 Like

For sorting, you could simply use predsort/3:

predsort (+Pred, +List, -Sorted)
Sorts similar to sort/2, but determines the order of two terms by calling Pred(-Delta, +E1, +E2) . This call must unify Delta with one of < , > or = .
https://www.swi-prolog.org/pldoc/man?predicate=predsort/3

Just define this comparator first:

arithmetic_compare(<,X,Y) :- X < Y, !.
arithmetic_compare(=,X,Y) :- X =:= Y, !.
arithmetic_compare(>,_,_).

It wouldn’t need rational number literats. Works already with the status quo. Here is an example run:

Welcome to SWI-Prolog (threaded, 64 bits, version 8.1.21)

?- predsort(arithmetic_compare,[1 rdiv 3,2 rdiv 7],X).
X = [2 rdiv 7, 1 rdiv 3].

In my system I don’t have predsort, but I have recently added sort/3 with a sort option list. This sort option list is also used in bags, aggregates and tabling. This would allow adding rational numbers already to these predicates without the introduction of a special literal.

What if rdiv/2 can recognize or internally store specially already normalized rational numbers, so that we gain speed. The Prolog interpreter parser would do something towards this end. The literal would then be denoted as rdiv/2 as before.

Agree. In fact, in last hour, I completed changing all dates in the Logtalk distribution source files to ISO 8601 format, update the compiler to support this format as an alternative when type-checking dates in info/1 directives, and run all the tests (close to 6000). Around 20 tests failed due to code using A/B terms as compound terms. For example, I have an eight puzzle formulation where the code includes:

	initial_state(four_steps, [2/2,1/3,3/2,2/3,3/3,3/1,2/1,1/1,1/2]).
	initial_state(five_steps, [2/1,1/2,1/3,3/3,3/2,3/1,2/2,1/1,2/3]).
	initial_state(eighteen_steps, [2/2,2/3,1/3,3/1,1/2,2/1,3/3,1/1,3/2]).

	goal_state(goal, [2/2,1/3,2/3,3/3,3/2,3/1,2/1,1/1,1/2]).

	print_state([S0,S1,S2,S3,S4,S5,S6,S7,S8]) :-
		member(Y, [3, 2, 1]),
		nl,
		member(X, [1, 2, 3]),
		member(Tile-X/Y, [' '-S0,1-S1,2-S2,3-S3,4-S4,5-S5,6-S6,7-S7,8-S8]),
		write(Tile),
		fail.
	print_state(_) :-
		nl.

I have a port of ROK code for random number generator that includes the code:

	random(A0, A1, A2, B0, B1, B2, Random) :-
		B0 is (A0*171) mod 30269,
		B1 is (A1*172) mod 30307,
		B2 is (A2*170) mod 30323,
		Float is A0/30269 + A1/30307 + A2/30323,
		Random is Float - truncate(Float).

This particular one is easy to fix. Just change to Float is A0/30269.0 + A1/30307.0 + A2/30323.0. But illustrates what will be a common source of nuisances: (/)/2 no longer returning a float. Note that in general both arguments may only be known at runtime…

Ironically, I have an example of type parameterized sort with one of the example types being rational number represented as A/B. the tests fail as the code contains definitions as:

	N1/D1 < N2/D2 :-
		{N1*D2 < N2*D1}.

but now:

?- 2/7 = N/D.
false.

You can argue that native support for rational numbers makes the example moot but that’s besides the point: code that is perfectly valid Prolog code no longer works. Btw, ECLiPSe have no issues with the code; all tests pass. There are more failed tests due to A/B no longer resulting in a float. For example:

!     population_variance_2_03: failure 
!       test assertion failed: =~=(1156/5,231.2)
!       in file /Users/pmoura/logtalk/library/statistics/tests.lgt between lines 197-198

I also have a port of Ivan Bratko state-space search algorithms that include clauses such as:

	search(Space, State, Threshold, Solution, Cost) :-
		expand([], l(State, 0/0), Threshold, _, yes, Path, Space, Cost),
		list::reverse(Path, Solution).

	expand(Path, l(State,Cost/_), _, _, yes, [State|Path], Space, Cost) :-
		Space::goal_state(State).
	expand(Path, l(State,F/G), Threshold, Tree, Solved, Solution, Space, Cost) :-
		F =< Threshold,
		(	bagof(Next/Cost2, (Space::next_state(State, Next, Cost2), \+ Space::member_path(Next, Path)), Successors) ->
			succlist(G, Successors, Trees, Threshold, Space),
			bestf(Trees, F2, Threshold),
			expand(Path, t(State, F2/G, Trees), Threshold, Tree, Solved, Solution, Space, Cost)
		;	Solved = never
		).

	succlist(_, [], [], _, _).
	succlist(G0, [State/Cost| Rest], Trees, Threshold, Space) :-
		G is G0 + Cost,
		Space::heuristic(State, H),
		F is G + H,
		succlist(G0, Rest, Trees2, Threshold, Space),
		insert(l(State, F/G), Trees2, Trees, Threshold).

Note the State/Cost term in the last predicate.

So, in just the above examples, we haver code from three different authors (myself, ROK, and Ivan Bratko; hey! I’m in good company :smile:) using A/B as a compound term or expecting A/B to evaluate always to a float.

I don’t mind (in fact, I already did the work) of switching to the ISO 8601 format for dates. But that is not what worries me in the proposed “natural” representation for rational numbers.

I have used that but predsort is significantly slower than sort which is implemented in C. On long lists it can make a huge difference.

And I really (really) dislike rdiv as a notation for a native rational literal. Maybe it’s acceptable when rationals are actually stored as compound terms, but the whole point of the exercise is to make them first class citizens. rdiv will likely persist for backwards compatibility reasons, but as the reference manual warns:

In the long term, it is likely that rational numbers will become atomic as well as a subtype of number . User code that creates or inspects the rdiv(M,N) terms will not be portable to future versions.

2 Likes

Thanks for sharing. We get concrete input!

This is a clear case. ?- list_rationals. will find it and gives the choice between changing the operator or changing the rational_syntax flag.

This results from prefer_rationals and ECLiPSe should fall into the same trap if this flag is set. The fix is indeed trivial (besides that anyone on random number generation would consider this code completely outdated).

At least we agree this is the most natural representation :slight_smile:

Luckily 1156/5 =:= 231.2 is true. Not really sure what is going on here, but if I extract the core comparison from Logtalk as below, ?- =~=(1156/5,231.2) succeeds.
`

'=~='(Float1, Float2) :-
    (   epsilon(Epsilon), abs(Float1 - Float2) < 100*Epsilon
    ->  true
    ;   abs(Float1 - Float2) < 0.00001 * max(abs(Float1), abs(Float2))
    ).

epsilon(Epsilon) :-
    Epsilon is epsilon.

Another case of a problem detected by list_rationals/0 and really easy to fix as the literal form only appears ones and this (0)/(0) is good enough (and looks really pretty :slight_smile: )

What does all this tell us? I guess that both syntax and semantics to use rationals all along does make some code break. ECLiPSe apparently observed the same as decided to make the default for prefer_rationals false. As the 1_3 syntax is (compared to ISO) a strict extension, nothing goes wrong. This leaves two options

  • Turn this stuff off by default. That means nobody gets hurt, but we keep using the old bad habits except for a few people preferring elegance.
  • Turn this stuff on by default. This will have no implications for most of us, cause some to change the flags and some to make their code robust.

My current assessment is to turn the flags off by default for the next release. With enough people turning this on we can get this all straight and in due time we can reconsider the default.

1 Like

Do you have some benchmark results how slow predsort/2 is? What about sort/4, would this also be “slow” as predsort/3? Is it possible to plug in arithmetic comparison into sort/4? Maybe with the help of simulating some literals, by some preliminary operator in the Prolog system.

Question is whether a native sort can really speed up sorting rational numbers. Comparing two rational numbers involves two multiplications and a comparison, since. Disclaimer: I dont know whether GMP does something smarter, or whether Henrici had a trick for rational comparison:

A/B < C/D  :<=> A*D < C*B

So its a little bit more expensive to compare rational numbers than integer numbers. I guess its also more expensive than comparing float numbers. In float numbers you can use a simple shift operation, so float comparison is fast.

Also even if you restrict rational numbers to len(numerator)+len(denomintor)<max, duing comparison you still might need larger products. If you have chosen max so that the stuff fits in a register, comparison is still more expensive and needs more registers access.

I would be really currious how the optimal sort for rational numbers turns out. Maybe you need a literal task force and a sort task force.

Edit 01.02.2020:
The =:=/2 on the other hand is not a problem. If A/B and C/D are already normalized, you simply have:

A/B =:= C/D  :<=> A =:= C /\ B =:= D

On another note, maybe you get something faster for comparison if you use Henrici subtraction, and then check the resulting sign.

This is a much bigger issue than syntax IMO. If rational numbers are really supported as numbers, then I believe the precision of an arithmetic operation should be preserved. Once a value is converted to a float, you’ve lost it, so it should never be the default. In this example, Float will be a rational number unless explicit steps are taken to convert it (as you’ve done). And the underlying question is when does it matter if the value is a float or a rational - they’re both just numbers. If you want to “demote” a value to float it’s easy enough to do.

SWI-Prolog already does maintain precision for integers, e.g.,

?- X is 4/2.
X = 2.

which I think is the right thing to do.

I don’t think anyone is arguing that there will not be existing broken Prolog code in using P/Q syntax for rational literals. That has been recognized from the beginning of this whole discussion. But if that principle was sacrosanct, we wouldn’t even have gotten to this point. The open question is how much pain to balance a perceived gain. And I have some sympathy with your argument that it’s not worth it.

BTW, I suspect Eclipse would have issues if you used any integers with embedded underscores as SWIP permits. But this just proves that ‘/’ is used more than embedded underscores, which is unsurprising.

Other than the 0/0 in search/5 (oops, another use of ‘/’), I’m not sure what the problems is here. There are no other literal rational numbers, and not much arithmetic evaluation to speak of. Why wouldn’t this continue to work as before (other than replacing 0/0 with 0/(0))?

Well, just test:

trat(Len) :-
    length(List, Len),
    maplist(random_rational, List),
    time(sort(List, _)).

random_rational(R) :-
    R is random(1 000 000 000) / (random(1 000 000 000)+1).

tfloat(Len) :-
    length(List, Len),
    maplist(random_float, List),
    time(sort(List, _)).

random_float(R) :-
    R is float(random(1 000 000 000) / (random(1 000 000 000)+1)).
103 ?- trat(1 000 000).
% 1 inferences, 2.741 CPU in 2.741 seconds (100% CPU, 0 Lips)
true.

104 ?- tfloat(1 000 000).
% 0 inferences, 1.366 CPU in 1.366 seconds (100% CPU, 0 Lips)
true.

This seems no show-stopper to me.

The relative difference from 1.366 seconds to 2.741 seconds
is more than 100%. The premium fee for exactness.

Finally, cross multiply and compare.
https://gmplib.org/repo/gmp/file/tip/mpq/cmp.c#l116

It seems that GMP uses the method I wrote down already
before. At the moment I don’t know whether 100% is good

or bad, and whether Henrici is applicable or whether there
are other methods to make comparison faster.

Edit 02.02.2020 I:
Why Henrici? For usual arithmetic it gives some speed-up.
Take this here:

horner(N, R) :-
   horner(N, 0, R).

horner(1, S, R) :- !,
   R is 1+(-1 rdiv 2)*(-1 rdiv 2)*S.
horner(N, S, R) :-
   M is N-1,
   H is 1+(-1 rdiv 2)*(((-1 rdiv 2)-M) rdiv N)*S,
   horner(M, H, R).

test :-
   between(1, 500, _),
   horner(100, _),
   fail.
test.

Its a strange way to compute sqrt(2):

?- horner(100,X), Y is float(X).
X = 22506280506048041472675379598885135569
33772436904539887974113033493227786869015
3259936903 rdiv 159143435651131725489722319
40698266883214596825515126958094847260581
103904401068017057792,
Y = 1.414213562373095.

Without Henrici in my system:

?- time(test).
% Up 933 ms, GC 8 ms, Thread Cpu 922 ms (Current 08/21/18 18:52:11)
Yes

With Henrici in my system:

?- time(test).
% Up 461 ms, GC 8 ms, Thread Cpu 438 ms (Current 08/21/18 18:51:32)
Yes

Edit 02.02.2020 II:
Here its also considered an unfortunate situation:
“multiplication will take more than linear time”
https://cs.stackexchange.com/questions/6266/comparing-rational-numbers

And this guy is more optimistic:
Bill Dubuque, continued fractions
https://math.stackexchange.com/a/14314/4414

Thanks for sharing. It revealed one bug in handling rationals in optimized mode. Using rationals we get the little program below. I guess that for a computer scientist the rdiv is ok, but I bet mathematicians are more happy with the code below :slight_smile:

horner(N, R) :-
   horner(N, 0, R).

horner(1, S, R) :- !,
   R is 1+(-1/2)*(-1/2)*S.
horner(N, S, R) :-
   M is N-1,
   H is 1+(-1/2)*(((-1/2)-M)/N)*S,
   horner(M, H, R).

test(N) :-
   forall(between(1, N, _), horner(100, _)).

Timing (Intel® Core™ i7-5557U), GMP 6.

101 ?- time(test(500)).
% 51,003 inferences, 0.080 CPU in 0.080 seconds (100% CPU, 640767 Lips)

Timing using valgrind produces a picture that illustrates quite nicely that Prolog is pretty useful for arithmetic as 40% of the time is spent in GMP (the Prolog part includes program startup, so the actual percentage is even higher). I think this picture is worth sharing:

1 Like

Cool! One question, will the new rationals allow packing and unpacking of numerator and denominator? Before could use unification for that purpose.

A test case could be, write a Prolog program to greedily find by Fibonacci algorithm an Egyptian fraction expansion of a given rational number:

Example:

 ?- egyptian(4 rdiv 13,X).
 X = 1 rdiv 4+(1 rdiv 18+1 rdiv 468).

How will such Prolog code look like in the future? I used this old code:

egyptian(1 rdiv X,1 rdiv X) :- !.
egyptian(X rdiv Y,1 rdiv Z+E) :-
    Z is (Y+X-1)//X,
    T is (X rdiv Y)-(1 rdiv Z),
    egyptian(T, E).

I have merged the rational branch to master. Summary of status and open issues:

  • Rationals is an atomic type. All integers are rationals.
  • Rationals are always in their canonical form. This means their identity can be tested using Prolog term identity (=, ==). Canonical means:
    • They are integer if possible
    • Their denominator (Numerator/Deniminator) is always a positive integer
    • The numerator can be any integer (except 0 as it is represented as an integer in that case)
    • The numerator and denominator have no common divisors.
  • Rationals replace integers in the standard order of terms. Note that in SWI-Prolog numbers are ordered by value where floats are smaller than rationals if they have the same value.

Then we have two flags:

  • prefer_rationals
    Causes all arithmetic that has a precise rational value to produce a rational result. At least, that is the theory. If this is not the case it is a bug, except for the pure float functions that happen to produce a rational value for a specific rational input (e.g., cos(0) is 1.0 rather than 1).
  • rational_syntax
    This flag chooses between natural (1/3), compatibility (1R3) or none,

After Paulo’s tests I’m convinced that for now the defaults are conservative, i.e., false and compatibility. My hope would be that in due time we can properly support rationals by default and use true and natural. So, please evaluate true and natural.

At the moment, following the other Prolog flags that control syntax, each module is created with the default setting. It is probably wise to add a :- set_prolog_flag(rational_syntax, <something>). in every module that uses literal rationals in the source code.

Remaining issues and notes

  • Need to sync the 1R3 notation with ECLiPSe if they are willing to give up (or provided an alternative to) 1_3. So, the compatibility syntax may change in the coming weeks.
  • If you want to help, please add a lot of tests to ../src/Tests/core/test_rational.pl. This is now practically empty (although there are tests about rationals in several other unit tests).
  • Code using the foreign language interface to walk over arbitrary terms must be aware it can get across a rational number. PL_is_rational() tells a number is rational (or integer). Extraction and unification can only be done through the GMP exchange functions PL_get_mpq() and PL_unify_mpq().
  • QLF files and external records (PL_record_external(), fast_read/2 and friends are not compatible. For QLF that is not unusual and if the system detects an incompatible version it automatically recompiles the source. Dropped compatibility for binary terms is a pity, but was hard to avoid.
  • Several constants in SWI-Prolog.h have changed. As a result, foreign plugins (shared objects, dlls) need to be recompiled. The source code is compatible.

The next release will probably follow shortly :slight_smile:

1 Like

Like ECLiPSe, there are the functions numerator/1 and demoninator/1. In addition there is rational(@Term, -Numerator, -Demoninator) which performs both a type test and extracts the Numerator and Denominator. So, you get (a bit less pretty)

Below is the minimal change. You can change the rdiv to /, but as you want the explicit term representation here

egyptian(1 rdiv X,1 rdiv X) :- !.
egyptian(X rdiv Y,1 rdiv Z+E) :-
    Z is (Y+X-1)//X,
    T is (X rdiv Y)-(1 rdiv Z),
    rational(T, N, D),
    egyptian(N rdiv D, E).

You can also go for real rationals:

egyptian(X, X) :-
    1 =:= numerator(X),
    !.
egyptian(In,N+E) :-
    rational(In, X, Y),
    Z is (Y+X-1)//X,
    N is 1/Z,
    T is (X rdiv Y)-(1 rdiv Z),
    egyptian(T, E).
?- egyptian(4/13,X).
X = 1/4+(1/18+1/468).

So my comment was independent of the items in the list being sorted, i.e., nothing to do with rational numbers and how they are implemented. predsort is implemented in Prolog, sort is implemented as a C primitive. So for a list of 1000 variables (already sorted):

?- length(L,1000),time((between(1,1000,_),predsort(compare,L,L1),fail)).
16,868,000 inferences, 1.810 CPU in 1.837 seconds (99 CPU, 9319780 Lips)
false.

?- length(L,1000),time((between(1,1000,_),sort(L,L1),fail)).
2,000 inferences, 0.033 CPU in 0.036 seconds (92 CPU, 60037 Lips)
false.

So if I’m forced to use predsort to numerically sort a list of mixed rational and integer (or float) numbers, then I pay a factor of about 50. If rationals are true numbers such that I can use the standard order of terms to sort, then I can use the native sort predicate which is an instant win and swamps any costs due to the actual comparison.

Personally I’d be happy if it was just a factor of two given two floats can be compared in a single machine instruction.

Amazing progress given a week ago I just wanted to capture this as an issue, not really expecting much to happen.

I’m pretty happy with all of this. The default settings are fine and will make all this largely invisible to existing code. I have a slight preference for 1r3 notation (as in J) over 1R3 but that’s a minor quibble. (The “camel case” effect also makes it a bit more visually distinctive.)

There was a background discussion on tripwires and such in case the rational components got too big. Should that be added to the remaining issues?

As for the longer term, I still have some doubts about supporting multiple syntaxes for rational literals. I’d probably go for getting rid of the rational_syntax flag and just using the compatibility syntax, whatever it turns out to be. After all this discussion, having the meaning of 1/3 depend on a flag setting doesn’t seem right. It complicates reading and porting code (despite the per module setting), as well as entering queries in the “PDE”. Are there any other similar cases of altering semantics based on a flag?

2 Likes

So you loose deep indexing. When you go from:

egyptian(1 rdiv X,1 rdiv X) :- !.

To this here:

egyptian(X, X) :- 1 =:= numerator(X), !.

Well sort or predsort doesn’t change anything that rational number comparison is non-linear in size. You need to test between float and rational, inside sort and predsort. These are 4 tests. At the moment we have two tests for sort (tfloat and trat by Jan W. made before palindrom day).

I guess if number literals are present you don’t need arithmetic_compare/2 anymore, and the SWI-Prolog specific compare/3 has a similar result. Not sure. For example usually compare/3 cannot compare across different number types:

Currently:

?- compare(X,2,9 rdiv 4).
X =  (<).

?- compare(X,2,7 rdiv 4).
X =  (<).

Will this change somehow? At the moment there is the following SWI-Prolog specific incompatibility with the ISO core standard. I only state the incompatibility, I don’t make a moral judgement. But you see the implementation difference here in the case of float:

GNU Prolog 1.4.5 (64 bits):

?- compare(X,2,2.25).
X = (>)
?- compare(X,2,1.75).
X = (>)

Welcome to SWI-Prolog (threaded, 64 bits, version 8.1.15):

?- compare(X,2,2.25).
X =  (<).
?- compare(X,2,1.75).
X =  (>).

I guess analogously SWI-Prolog will now do something for rational numbers and change what happens currently with rdiv/2 compounds. Prolog systems that don’t follow the ideas of SWI-Prolog have still the option to implement arithmetic_compare/3 (nothing from ISO, my idea).

ISO core standard would require a list of such incompatibilities for each Prolog system to be published. I am not expert on these lists. But ECLiPSe Prolog is one of the few Prolog systems that a year ago or so published such lists for his system:

Specification of implementation defined features
http://eclipseclp.org/doc/bips/lib/iso_strict/

Specification of implementation defined features
http://eclipseclp.org/doc/bips/lib/iso/

Because 9 rdiv 4 is a term. If we read rationals as such we get

101 ?- compare(X,2,9/4).
X =  (<).

102 ?- compare(X,2,7/4).
X =  (>).
1 Like

True. This is still a bit open to debate, but knowing Joachim a bit I think he goes for R and I would like to sync this with ECLiPSe.

Probably. A tripwire system can help preparing code to deal with rationals as default.

As usual, one clear standard without flags is best. I think that should be simply 1/3 and always prefer rationals, but that requires people to update their software to make it compatible with this syntax (not hard) and actively decide to convert to floats (often by using a float constant in the expression) there were rationals are not suitable (read: too expensive and we are happy with the float approximation). I think this is a fair compromise. Note that different modules can do as they please and as the thing exchanged is the same rational it will all happily work together, unlike similar existing practice with chars vs. codes.