TIL: CLP (Q) is a natural consequence of CLP (FD)

I was thinking about CLP (FD) and the limits of how it can be used (I.E., why can’t it be used on floats?) and I realised something I didn’t know before: obviously constraint processing over rationals and complex rationals exists, since these are just pairs of integers (or rationals) (this is how they teach it in real analysis), so you get them for free with the integers. You can’t really do the same efficiently for irrationals since whilst you can identify limits with dedekind cuts and you might try to do interval sets with these, you need to basically know them with infinite precision to prevent overlap, and there are a whole boatload of irrational numbers (more than rationals) rendering that method quite inefficient. Not sure how CLP QR attempts to get around that problem, but very interesting thought! Moral of the story: You can manage to encode most things with integers.

Between rational numbes and irrational numbers,
there is for example the class of algebraic numbers.

If a number is irrational but not algebraic, the
number is called transcendental. Examples:

  • x = sqrt(2): Is an algebraic number , since
    it is the root of the polynomial x^2 - 2 .

  • x = pi: Is a transcendental number , since it
    it is not the root of any polynomial.

Many computer algebra systems (CAS) can
compute with algebraic numbers. Cantor already
showed that there are countably infinite many

such numbers, they are not uncountable, unlike
the transcendental numbers. CAS systems can
compute with them since there exist

algorithms for the basic operations and equality,
like for example using sylvester matrixes. There
do not exist Prolog implementations for CLP(A)

where A would be the algebraic numbers, because
Prolog never went that far as CAS did. Quite unlike
LISP, which was more active in this domain:

See also:

Macsyma - 5.1 Algebraic Integers
https://people.eecs.berkeley.edu/~fateman/macsyma/docs/refman16.pdf

Algebraic numbers allow arbitrary roots, like for
example cube root 2^(1/3) , since this is again a root
of a polynomial x^3 – 2. Inside algebraic numbers

there is an interesting subclass of numbers built by
only applying square root, which occupied already
the Greek minds. These have a geometric equivalent

and are numbers constructible by ruler and compass.

2 Likes

Thanks for your insight! Wow, I didn’t know those old LISP algebra systems were still around and kicking. I wish I were less busy with other projects, or I’d be diving down that rabbit hole right now. I did wonder if maybe there was some way to represent the algebraics, being able to do that in SWI-PL would be really nice for all kinds of applications.

By the way, I was having a play around with implementing CLP(Q) myself from CLP(FD) because I think it would be neat to show how trivial it is, but since mod is seemingly not perfectly bidirectional (nor is gcd) I’m having issues with bidirectional rational multiplication/division. Any recommendations on this? Am I missing another nuance? Or simply not approaching the issue correctly?

Currently Python is keeping up with these systems:

https://12000.org/my_notes/CAS_ode_tests/index.htm

But the trend was possibly that CAS was first takled by
proof verification, many algorithms are now verified via Isabelle/HOL,
Coq, etc.. and a new trend is basically that LLMs can do all

the CAS. Also bandwagons that Prolog totally missed. The
Python thingy is interesting, since recently ChatGPT provided
me a solution in Python, run the Python code for me on their

server and then gave me back the result from Python.

1 Like

I support add-on pack(clpBNR) which is a CLP over the extended set of reals, including rationals and irrationals (which I define as including algebraics). Also now available on SWISH (see clpBNR_quickstart).

If you think it’s trivial, I don’t understand what you’re trying to do. While rationals are “just pairs of integers” (numerator and denominator), I don’t think they come “for free”. That SWIP additionally includes the GMP library to support arithmetic for rationals (and extended integers) would suggest otherwise. Similarly I don’t think a CLP over rationals comes for free because you have CLP(FD), but perhaps I’m wrong.

As an aside, floats are a subset of rationals with some interesting properties. They are compact (SWIP floats are 64 bits + cell overhead) but can represent a large range of values. And they are efficient as all commonly used CPU hardware support the IEEE standard. The downside is that arithmetic operations on floats may be imprecise, so you need rounding control to ensure overall correctness of interval enclosures (sets of reals).

Somehow I have the feeling they follow a different philosophy than
most constraint logic programming CLP(X). The later often sees an
application as doing two steps, first building a model via attribute variables

and then solving the model. This approach breaks if I want to have a
knowledge base of constraints, so that I have a kind of facts being
constraints. Sympy can do that somehow:

from sympy.physics.units import *

(5000 * gram) / (2 * gram)
# 2500

The import makes available a knowledgebase of units and their
relationship. I tried to do the same and came up with Disjoint-set data
structure
. It should work whenever you have a set of relations of the form

xj = fij(xi)

Where fij is invertible. The easiest is scaling:

?- init.

?- find(hr, sec, X).
X = 3600/1.

?- find(sec, min, X).
X = 1/60.

?- find(meters, sec, X).
false.

Source code:

find(U, V, R) :-
   deref(U, 1/1, H, R1),
   deref(V, 1/1, H, R2),
   div(R1, R2, R).

deref(U, R, V, S) :-
   link(U, R2, U2), !,
   mul(R, R2, H),
   deref(U2, H, V, S).
deref(U, R, U, R).

conversion(meters, 82/25, feet).  /* Better 30000/9144 ? */
conversion(feet,   12/1,  in  ).
conversion(hr,     60/1,  min ).
conversion(min,    60/1,  sec ).

:- dynamic link/3.

add(U, R, V) :-
   deref(U, 1/1, U2, R1),
   mul(R, R1, H),
   deref(V, 1/1, V2, R2),
   div(H, R2, J),
   assertz(link(U2, J, V2)).

init :- retractall(link(_,_,_)),
   conversion(U, R, V), add(U, R, V), fail; true.

div(A, B/C, R) :-
   mul(A, C/B, R).

mul(X/Y, Z/T, C/D) :-
   A is X*Z,
   B is Y*T,
   gcd(A, B, H),
   C is A//H,
   D is B//H.

gcd(0, A, R) :- !, R = A.
gcd(B, A, R) :-
   H is A mod B,
   gcd(H, B, R).

Sounds like you implemented some ideas from here:

Generalised Union-Find
Thom Frühwirth - 2008
https://exia.informatik.uni-ulm.de/fruehwirth//Papers/linear-lnai.pdf

It shows how to combine these functions:

Using the Linear Polynomials from the same CHR paper, one could
convert celsius and fahrenheit etc.. , which have a reference point.
But your Prolog code shows, we don’t need CHR either.

Sympy suggests something even more, to use exponent vectors of
monomials, and then convert via matrix inversion. Not sure how
this would fit into a Disjoint-set data structure ?

Sorry, it’s not really relevant to your goal, but I have recently made a pack for unit aware computation that is compatible with clpBNR, meaning you can do this:

?- qeval({A == 5000*gram / 2 * gram}). % wrong equation, see below
A = 2500.

Here is the original post: Units: a new pack for units and quantities

edit:
the correct equation was this:

?- qeval({A == 5000*gram / (2 * gram)}).
A = 2500.

By the way, I don’t really understand that sentence.
Isn’t the model building in clp* explicitly the encoding of facts and relations into constraints ?

Doesn’t work on my side. Maybe I dont have the newest releases?

?- pack_list_installed.
% Contacting server at https://www.swi-prolog.org/pack/query ... ok
i clpBNR@0.12.1                    - CLP over Reals using Interval Arithmetic - includes Rational, Integer and Boolean domains as subsets.
i units@0.3                        - unit and quantity arithmetic for swi-prolog
true.

?- use_module(library(units)).
% *** clpBNR v0.12.1 ***.
%   Arithmetic global flags will be set to prefer rationals and IEEE continuation values.
true.

?- qeval({A == 5000*gram / 2 * gram}).
ERROR: Domain error: `q{q:1,u:1,v:_15236}' expected, found `5000/2 * kind(isq:mass)^2[si:gram^2]'
Etc..

Sorry, I have recently made a lot of change to the library.
I have pushed a new version 0.4. I believe you can update and test again ?

Ah, I made a classic syntax error and I wrongly retyped the equation in my previous message.
It should have been:

?- qeval({A == 5000*gram / (2 * gram)}).

multiplication operator * has lower priority than the division / operation.
Although the error is a bit cryptic, you can see that the quantity equation is not right:
5000/2 * kind(isq:mass)^2[si:gram^2]. this means that you have the number 5000/2 with the quantity mass^2 and unit gram^2.

The correct equation should have been a dimensionless quantity 1 and no unit 1, which you can see when evaluation the expression with the private predicate eval_/2:

?- units:eval_({A == 5000*gram / (2 * gram)}, E).
E = {A==5000/2} * 1[1].

Most logic programming instance of the CLP(X) schema dont support
global variables. Only CLP(B), the sat solver by Markus Triska has something.
But CLP(FD) and CLP(BNR) dont have such a feature, which would be

a first step towards a tighter integration of units. So Sympy can do:

from sympy.physics.units import *

(5000 * gram) / (2 * gram)
# 2500

5 * gram + 1000 * gram
# 1005*gram

https://stackoverflow.com/q/48661433/28555451

Whereas SWI-Prolog can only do:

?- qeval({A == 5000*gram / (2 * gram)}).
A = 2500.

?- qeval({A == 5*gram + 1000*gram}).
ERROR: Domain error: `q{q:1,u:1,v:_7394}' expected, found `5+1000 * kind(isq:mass)[si:gram]'

Unless you would juggle with the variable A, hand another variable A’ to
clp(BNR) and solve a unit less problem, and then reconstruct A with unit
in pretty printed form as a Prolog term 1005*gram.

Or did I make a parenthesis error again?

BTW: You could even beat Sympy if you could handle:

(5 * kilogram) / (2 * gram)
# (5 * kilogram) / (2 * gram)

5 * gram + 1 * kilogram
# 5*gram + kilogram

Not sure why Sympy supposedly can’t do these examples…

Yeah, since clpBNR does not have an assignment operator, we need to specify the units for both side of the equation.
The best we can do for now is to infer the unit before hand:

?- units:eval_(5*gram + 1000*gram, E), qeval({A*E.u == E}).
E = 5+1000 * isq:mass[si:gram],
A = 1005.

The units:eval_/2 predicate is the internal units library expression parser predicate.
Unfortunately, we have to evaluate the expression 5*gram + 1000*gram first, because swi-prolog dict functional notation is eager and will expand before qeval has done it’s job.

It can :slight_smile:

?- units:eval_(5*gram + 1*kilogram, E), qeval({A*E.u == E}).
E = 5+10^3 * isq:mass[si:gram],
A = 1005.

And by the way, you don’t need clpBNR for the examples you gave:

?- qeval(X is 5*gram + 1*kilogram).
X = 1005 * isq:mass[si:gram].

Crazy, it turns out you can just do this:

?- qeval({Z is X*gram + Y*kilogram}), X = 1, Z.v = 1001.
Z = 1001 * isq:mass[si:gram],
X = Y, Y = 1.