% fougau.pl =================================================================
% constraint handling rules for linear arithmetic
% fouriers algorithm for inequalities and gaussian elemination for equalities
% thom fruehwirth 950405-06, 980312
% 961107 Christian Holzbaur, SICStus mods
% CHOOSE one of the propagation rules below and comment out the others!
% completeness and termination depends on propagation rule used
:- use_module( library(chr)).
:- use_module( library(lists), [append/3]).
:- ensure_loaded('math-utilities').
handler fougau.
% auxiliary constraint to delay a goal G until it is ground
constraints check/1.
check(G) <=> ground(G) | G.
constraints {}/1.
{ C,Cs } <=> { C }, { Cs }.
{A =< B} <=> ground(A),ground(B) | A=<B.
{A >= B} <=> ground(A),ground(B) | A>=B.
{A < B} <=> ground(A),ground(B) | A<B.
{A > B} <=> ground(A),ground(B) | A>B.
{A =\= B} <=> ground(A),ground(B) | A=\=B.
{A =< B} <=> normalize(B,A,P,C), eq(P,C,'>'('=')).
{A >= B} <=> normalize(A,B,P,C), eq(P,C,'>'('=')).
{A < B} <=> normalize(B,A,P,C), eq(P,C,'>'('>')).
{A > B} <=> normalize(A,B,P,C), eq(P,C,'>'('>')).
%{A < B} <=> normalize(B,A+1,P,C), eq(P,C,(>=)). % adopt to integer
%{A > B} <=> normalize(A,B+1,P,C), eq(P,C,(>=)). % adopt to integer
{A =\= B} <=> normalize(A+X,B,P,C), eq(P,C,(=:=)), check(X=\=0).
{A =:= B} <=> ground(A),ground(B) | X is A-B, zero(X). % handle imprecision of reals
{A =:= B} <=> var(A), ground(B) | A is B.
{B =:= A} <=> var(A), ground(B) | A is B.
{A =:= B} <=> unconstrained(A),var(B) | A=B.
{B =:= A} <=> unconstrained(A),var(B) | A=B.
{A =:= B} <=> normalize(A,B,P,C), eq(P,C,(=:=)).
constraints eq/3.
% eq(P,C,R)
% P is a polynomial (list of monomials variable*coefficient),
% C is a numeric constant and R is the relation between P and C
% simplify single equation
zero @ eq([],C1,(=:=)) <=> zero(C1).
zero @ eq([],C1,'>'('=')) <=> C1>=0.
zero @ eq([],C1,'>'('>')) <=> C1>0.
unify @ eq([X*C2],C1,(=:=)) <=> nonground(X),nonzero(C2) | is_div(C1,C2,X).
%, integer(X) % if integers only
simplify @ eq(P0,C1,R) <=> delete(X*C2,P0,P),ground(X) | % R any relation
%, integer(X), % if integers only
is_mul(X,C2,XC2),
C is XC2+C1,
eq(P,C,R).
/*
% must use if you unify variables of equations with each other
unified @ eq(P0,C1,R1) <=>
append(P1,[X*C2|P2],P0),var(X),delete(Y*C3,P2,P3),X==Y
|
C23 is C2+C3,
append(P1,[X*C23|P3],P4),
sort1(P4,P5),
eq(P5,C1,R1).
*/
%(1) remove redundant inequation
% -1 (change in number of constraints)
red_poly @ eq([X*C1X|P1],C1,'>'(R1)) \ eq([X*C2X|P2],C2,'>'(R2)) <=>
C is C2X/C1X, % explicit because of call_explicit bug
C>0, % same sign
C1C is C1*C,
C1C=<C2, % remove right one
stronger(C1X,C1C,R1,C2X,C2,R2), % remove right one if C1C=:= C2
same_poly(P1,P2,C)
|
true.
%(2) equate opposite inequations
% -1
opp_poly @ eq([X*C1X|P1],C1,'>'(R1)), eq([X*C2X|P2],C2,'>'(R2)) <=>
C is C2X/C1X,
C<0, % different sign
C1C is C1*C,
C1C>=C2, % applicable?
same_poly(P1,P2,C)
|
Z is C1C-C2, zero(Z), % must have identical constants
(R1)=('='), (R2)=('='), % fail if one of R's is ('>')
eq([X*C1X|P1],C1,(=:=)).
%(3) usual equation replacement (like math-gauss.chr)
% 0
/*
elimin_eager @ eq([X*C2X|PX],C1X,(=:=)) \ eq(P0,C1,R) <=> % R any relation
extract(X*C2,P0,P)
|
is_div(C2,C2X,CX),
mult_const(eq0(C1X,PX),CX,P2),
add_eq0(eq0(C1,P),P2,eq0(C3,P3)),
sort1(P3,P4),
eq(P4,C3,R).
*/
elimin_lazy @ eq([X*C2X|PX],C1X,(=:=)) \ eq([X*C2|P],C1,R) <=>
is_div(C2,C2X,CX),
mult_const(eq0(C1X,PX),CX,P2),
add_eq0(eq0(C1,P),P2,eq0(C3,P3)),
sort1(P3,P4),
eq(P4,C3,R).
% choose one of the propagation rules below and comment out the others!
% completeness and termination depends on propagation rule used
%(4) propagate, transitive closure of inequations, various versions
% +1
/*
% complete, but too costly, propagate_lazy is as good, can loop
propagate_eager @ eq([X*C2X|PX],C1X,'>'(R1)), eq(P0,C1,'>'(R2)) ==>
extract(X*C2,P0,P),
is_div(C2,C2X,CX),
CX>0 % different sign?
|
combine_ineqs(R1,R2,R3),
mult_const(eq0(C1X,PX),CX,P2),
add_eq0(eq0(C1,P),P2,eq0(C3,P3)),
sort1(P3,P4),
eq(P4,C3,'>'(R3)).
% complete, may loop
propagate_lazy @ eq([X*C2X|PX],C1X,'>'(R1)), eq([X*C2|P],C1,'>'(R2)) ==>
is_div(C2,C2X,CX),
CX>0 % different sign?
|
combine_ineqs(R1,R2,R3),
mult_const(eq0(C1X,PX),CX,P2),
add_eq0(eq0(C1,P),P2,eq0(C3,P3)),
sort1(P3,P4),
eq(P4,C3,'>'(R3)).
*/
/*
% incomplete, number of variables does not increase, may loop
propagate_pair @ eq([X*C2X|PX],C1X,'>'(R1)), eq(P0,C1,'>'(R2)) ==>
not(PX=[_,_,_|_]), % single variable or pair of variables only
extract(X*C2,P0,P),
is_div(C2,C2X,CX),
CX>0 % different sign?
|
combine_ineqs(R1,R2,R3),
mult_const(eq0(C1X,PX),CX,P2),
add_eq0(eq0(C1,P),P2,eq0(C3,P3)),
sort1(P3,P4),
eq(P4,C3,'>'(R3)).
*/
% incomplete, is interval reasoning, number of variables decreases, loop free
propagate_single @ eq([X*C2X],C1X,'>'(R1)), eq(P0,C1,'>'(R2)) ==> % single variable only
(P0=[V*C2|P],V==X ; P0=[VC,V*C2|PP],V==X,P=[VC|PP]), % only first or second variable
is_div(C2,C2X,CX),
CX>0 % different sign?
|
combine_ineqs(R1,R2,R3),
is_mul(C1X,CX,C1XCX),
C3 is C1+C1XCX,
eq(P,C3,'>'(R3)).
/*
% incomplete, ignore inequations until they are sufficiently simplified
%propagate_never @ eq([X*C2X|PX],C1X,'>'(R1)), eq([X*C2|P],C1,'>'(R2)) ==>
% fail | true.
*/
% handle nonlinear equations ------------------------------------------------
operator(450,xfx,eqnonlin).
constraints (eqnonlin)/2.
non_linear @ X eqnonlin A <=> ground(A) | A1 is A, {X=:=A1}.
non_linear @ X eqnonlin A*B <=> ground(A) | A1 is A, {X=:=A1*B}.
non_linear @ X eqnonlin B*A <=> ground(A) | A1 is A, {X=:=A1*B}.
% labeling, useful really only for integers ---------------------------------
%label_with eq([XC],C1,'>'('=')) if true.
%eq([XC],C1,'>'('=')) :- eq([XC],C1,(=:=)) ; eq([XC],C1,'>'('>')).
% auxiliary predicates --------------------------------------------------------
% combine two inequalities
combine_ineqs(('='),('='),('=')):- !.
combine_ineqs(_,_,('>')).
same_poly([],[],C).
same_poly([X*C1|P1],[X*C2|P2],C) ?-
%X==Y,
C4 is C*C1-C2, zero(C4),
same_poly(P1,P2,C).
stronger(C1X,C1C,R1,C2X,C2,R2):-
C1C=:=C2 ->
\+ (R1=('='),R2=('>')),
C1A is abs(C1X)+1/abs(C1X), C2A is abs(C2X)+1/abs(C2X),
C1A=<C2A
;
true.
/* end of file math-fougau.chr ----------------------------------------------*/