Formulating quadratic equations in clpfd

Question

CLPFD-systems are not primarily targeted to handle quadratic equations efficiently, nevertheless, are there better ways to formulate problems like the following?

It seems the problem boils down to equations like the following. SWI with library(clpfd) gave:

?- time( ((L+7)^2#=L^2+189, L in 0..10000000) ).
% 252,169,718 inferences, 87208.554 CPU in 87445.038 seconds (100% CPU, 2892 Lips)
ERROR: Out of local stack

But now the latest version in SWI gives

?- time( ((L+7)^2#=L^2+189, L in 0..10000000) ).
% 3,805,545,940 inferences, 868.635 CPU in 870.311 seconds (100% CPU, 4381063 Lips)
L = 10.

and in SICStus 4.3beta7 I get:

| ?- statistics(runtime,_).
yes
| ?- (L+7)*(L+7)#=L*L+189, L in 0..10000000.
L = 10 ? ;
no
| ?- statistics(runtime,[_,T_ms]).
T_ms = 2550 ? ;
no

Daube

I also noticed - solving some Project Euler problem - that plain arithmetic can be much faster than CLP(FD). Seems that Markus has got an interesting problem to solve... BTW the time to find the first solution to (L+7)*(L+7)#=L*L+189, L in 0..10000, label([L]) is much smaller than the time required when the ^2 operator is used... – Farina

@CapelliC: * in place of ^2 leads to less consistency. Also, the cost is proportional to the size of the domain for L. – Daube

Is actually solving the equation out of the question? What is your use case exactly? I am asking because it feels maybe a bit wasteful to not try to simplify before doing anything else. – Dulcia

@Boris: Solving the equation is always an option but at what price? – Daube

It's difficult to know the "price" when the use case is not known. If those are indeed just quadratic equations, it's fairly trivial to simplify them (or solve them) programmatically. – Dulcia

@Boris one might not use this in the quadratic case since there are precise algebraic solutions, but (I think) Ulrich is offering this by way of example. I played with it a bit in GNU and also saw the huge inconsistency between using X^2 versus X*X (in GNU, the X^2 was far better behaved. I think in general, for algebraic equations, one might need to apply some maths to prune the search tree a bit. Otherwise, the Prolog FD implementation(s) will struggle. – Amice

@repeat That was a year ago! :p I can't say I really remember what I did. I'm thinking I probably ran Ulrich's equation, but I don't recall how I measured it. I just tried it again and couldn't get the function to run more than once in the X * X case without resulting in a segfault (gprolog V 1.4.2). – Amice

Answer 1

To solve this quickly, a constraint solver that has already a primitive X #= Y^2 constraint, might also implement the following rules:

Rule #1:

X #= (Y+C)^2 --> H #= Y^2, X #= H + 2*C*Y + C^2
% where C is an integer and X,Y variables

Rule #2:

X #= Z^2, Y #= Z^2 --> X #= Z^2, X #= Y.
% where X,Y,Z are variables

The above rules will reduce the equation to a linear equation, which is anyway directly solved by a decent CLP(FD) system. Here you see the difference between a system with and without rule #2:

Without Rule #2:

?- (L+7)*(L+7) #= L*L+189, stored.
_C #= 140+_B-14*L.
_B #>= 0.
_C #>= 0.
_B #= L*L.
_C #= L*L.
Yes

With Rule #2:

?- (L+7)*(L+7) #= L*L+189, stored.
L = 10
?- statistics(time, S), (L+7)*(L+7) #= L*L+189, statistics(time, T), U is T-S.
U = 3

But rule #2 looks a little ad hoc to me. Not yet sure whether one should keep it.

Bye

Christa

Answer 2

clpBNR fares better:

% *** clpBNR v0.9.13alpha ***.
Welcome to SWI-Prolog (threaded, 64 bits, version 8.4.3)

?- time(( { (L+7)**2 == L**2+189 }, L::integer(0, 10_000_000), solve(L) )).
% 1,269,854 inferences, 0.163 CPU in 0.164 seconds (100% CPU, 7770826 Lips)
L = 10 ;
% 271,125,405 inferences, 40.874 CPU in 40.957 seconds (100% CPU, 6633143 Lips)
false.

Comparing clpfd and clpBNR using a range that succeeds in a reasonable timescale:

?- use_module(library(clpfd)).
?- time( ((L+7)^2#=L^2+189, L in 0..90_000) ).
% 13,609,198 inferences, 11.482 CPU in 11.497 seconds (100% CPU, 1185274 Lips)
L = 10.

?- time(( { (L+7)**2 == L**2+189 }, L::integer(0, 90_000), solve(L) )).
% 540,550 inferences, 0.107 CPU in 0.107 seconds (100% CPU, 5058454 Lips)
L = 10 ;
% 2,105,638 inferences, 0.320 CPU in 0.321 seconds (100% CPU, 6578814 Lips)
false.

Juratory

Better than what? Could you enter the same query for comparison on the same machine? – Daube