301 lines
7.0 KiB
Perl
301 lines
7.0 KiB
Perl
|
% math-utilities.pl ===========================================================
|
||
|
|
% auxiliary predicates for math*.pl constraint solvers
|
||
|
|
% thom fruehwirth 1991-92, revised 930518,931223,940304
|
||
|
|
% 961030 christian holzbaur, SICStus adaption
|
||
|
|
|
||
|
|
:- use_module( library('chr/getval')).
|
||
|
|
:- use_module( library('chr/matching')).
|
||
|
|
:- use_module( library('chr/ordering'), [globalize/1,var_compare/3]).
|
||
|
|
|
||
|
|
|
||
|
|
% SETTINGS --------------------------------------------------------------------
|
||
|
|
|
||
|
|
% for use in is/2: precision, slack variables, simulated infimum, etc.
|
||
|
|
|
||
|
|
% Code works with flag prefer_rationals on or off
|
||
|
|
% and with float_precision single or double
|
||
|
|
|
||
|
|
% adapt precision for zero/1 test
|
||
|
|
:- ( current_module(eclipse) ->
|
||
|
|
get_flag(float_precision,G)
|
||
|
|
;
|
||
|
|
G = double
|
||
|
|
),
|
||
|
|
(G==single -> setval(precision,1.0e-06),setval(mprecision,-1.0e-06)
|
||
|
|
;
|
||
|
|
G==double -> setval(precision,1.0e-12),setval(mprecision,-1.0e-12)
|
||
|
|
).
|
||
|
|
|
||
|
|
slack(X,X). % :- X>=0.
|
||
|
|
|
||
|
|
inf( 3.40282e38).
|
||
|
|
minf( -3.40282e38).
|
||
|
|
sup( 1.0e-45).
|
||
|
|
msup( -1.0e-45).
|
||
|
|
|
||
|
|
:- multifile portray/1.
|
||
|
|
|
||
|
|
portray( X) :- math_portray( X, Xp), print( Xp).
|
||
|
|
|
||
|
|
|
||
|
|
% PRETTY PRINT ---------------------------------------------------------------
|
||
|
|
|
||
|
|
% for math-gauss.pl and math-elim.pl
|
||
|
|
math_portray(equals(P,C),P1=:=0):- zero(C),!,
|
||
|
|
make_poly(P,P1).
|
||
|
|
math_portray(equals(P,C),P1=:=C1):-!,
|
||
|
|
MC is (-C),
|
||
|
|
avoid_float(MC,C1),
|
||
|
|
make_poly(P,P1).
|
||
|
|
% for math-fougau.pl
|
||
|
|
math_portray(eq(P,C,(=:=)),P1=:=C1):-!,
|
||
|
|
MC is (-C),
|
||
|
|
avoid_float(MC,C1),
|
||
|
|
make_poly(P,P1).
|
||
|
|
math_portray(eq(P,C,'>'('=')),P1>=C1):-!,
|
||
|
|
MC is (-C),
|
||
|
|
avoid_float(MC,C1),
|
||
|
|
make_poly(P,P1).
|
||
|
|
math_portray(eq(P,C,'>'('>')),P1>C1):-!,
|
||
|
|
MC is (-C),
|
||
|
|
avoid_float(MC,C1),
|
||
|
|
make_poly(P,P1).
|
||
|
|
% for all three math*pl solvers
|
||
|
|
math_portray(eqnonlin(X,(E)),X=:=E):-!.
|
||
|
|
|
||
|
|
|
||
|
|
make_poly([],0).
|
||
|
|
make_poly([X*C],-CX):- C<0,!,
|
||
|
|
C1 is (-C),
|
||
|
|
avoid_float(C1,C2),
|
||
|
|
make_mono(C2,X,CX).
|
||
|
|
make_poly([X*C],CX):-!,
|
||
|
|
avoid_float(C,C1),
|
||
|
|
make_mono(C1,X,CX).
|
||
|
|
make_poly([X*C|P],P1-CX):- C<0,!,
|
||
|
|
C1 is (-C),
|
||
|
|
avoid_float(C1,C2),
|
||
|
|
make_mono(C2,X,CX),
|
||
|
|
make_poly(P,P1).
|
||
|
|
make_poly([X*C|P],P1+CX):-
|
||
|
|
avoid_float(C,C1),
|
||
|
|
make_mono(C1,X,CX),
|
||
|
|
make_poly(P,P1).
|
||
|
|
|
||
|
|
make_mono(C,X,CX):- nonvar(X),X=slack(Y),!,make_mono(C,Y,CX).
|
||
|
|
make_mono(C,X,CX1):- nonvar(X),number(X),!,CX is C*X,avoid_float(CX,CX1).
|
||
|
|
make_mono(1,X,X):-!.
|
||
|
|
% make_mono(1_1,X,X):-!.
|
||
|
|
make_mono(C,X,C*X).
|
||
|
|
|
||
|
|
|
||
|
|
% AUXILIARY PREDICATES -------------------------------------------------------
|
||
|
|
|
||
|
|
nonground( X) :- ground( X), !, fail.
|
||
|
|
nonground( _).
|
||
|
|
|
||
|
|
%
|
||
|
|
% sort X*K,slack(_)*K with globalized Xs
|
||
|
|
%
|
||
|
|
sort1(A,B):-
|
||
|
|
msort(A,C),
|
||
|
|
((C=[X*_|_],nonvar(X),X=slack(_))->A=B;B=C). % slacks unordered why?
|
||
|
|
|
||
|
|
msort( L, S) :-
|
||
|
|
length( L, Len),
|
||
|
|
msort( Len, L, [], S).
|
||
|
|
|
||
|
|
msort( 0, L, L, []) :- !.
|
||
|
|
msort( 1, [X|L], L, [X]) :- !.
|
||
|
|
msort( N, L0, L2, S) :-
|
||
|
|
P is N>>1,
|
||
|
|
Q is N-P,
|
||
|
|
msort( P, L0, L1, Sp),
|
||
|
|
msort( Q, L1, L2, Sq),
|
||
|
|
merge( Sp, Sq, S).
|
||
|
|
|
||
|
|
merge( [], B, B) :- !.
|
||
|
|
merge( A, [], A) :- !.
|
||
|
|
merge( [A|As], [B|Bs], Res) :-
|
||
|
|
cmp( R, A, B),
|
||
|
|
merge( R, A, As, B, Bs, Res).
|
||
|
|
|
||
|
|
merge( =, A, As, _, Bs, [A|Rest]) :- merge( As, Bs, Rest).
|
||
|
|
merge( <, A, As, B, Bs, [A|Rest]) :- merge( As, [B|Bs], Rest).
|
||
|
|
merge( >, A, As, B, Bs, [B|Rest]) :- merge( [A|As], Bs, Rest).
|
||
|
|
|
||
|
|
cmp( R, X, Y) :- var(X), var(Y), !, var_compare( R, X, Y).
|
||
|
|
cmp( R, X, _) :- var(X), !, R = (<).
|
||
|
|
cmp( R, _, Y) :- var(Y), !, R = (>).
|
||
|
|
cmp( R, X, Y) :-
|
||
|
|
functor( X, Fx, Ax),
|
||
|
|
functor( Y, Fy, Ay),
|
||
|
|
compare( Rr, Ax/Fx, Ay/Fy),
|
||
|
|
( Rr = (=),
|
||
|
|
Ax > 0 ->
|
||
|
|
cmp_args( 1,Ax, X, Y, R)
|
||
|
|
;
|
||
|
|
R = Rr
|
||
|
|
).
|
||
|
|
|
||
|
|
cmp_args( N,M, _, _, R) :- N>M, !, R = (=).
|
||
|
|
cmp_args( N,M, X, Y, R) :-
|
||
|
|
arg( N, X, Ax),
|
||
|
|
arg( N, Y, Ay),
|
||
|
|
cmp( Rr, Ax, Ay),
|
||
|
|
( Rr = (=) ->
|
||
|
|
N1 is N+1,
|
||
|
|
cmp_args( N1,M, X, Y, R)
|
||
|
|
;
|
||
|
|
R = Rr
|
||
|
|
).
|
||
|
|
|
||
|
|
|
||
|
|
rev([],L,L).
|
||
|
|
rev([X|L1],L2,L3):- rev(L1,[X|L2],L3).
|
||
|
|
|
||
|
|
extract(X*C2,P0,P) ?- delete(Y*C2,P0,P),X==Y,!.
|
||
|
|
|
||
|
|
delete( X, [X|L], L).
|
||
|
|
delete( Y, [X|Xs], [X|Xt]) :-
|
||
|
|
delete( Y, Xs, Xt).
|
||
|
|
|
||
|
|
zero( slack(S)) ?- !, zero( S).
|
||
|
|
zero(C):-
|
||
|
|
float(C) ->
|
||
|
|
getval(precision,P),
|
||
|
|
getval(mprecision,MP),
|
||
|
|
MP < C, % cope with imprecision
|
||
|
|
C < P
|
||
|
|
;
|
||
|
|
C=:=0.
|
||
|
|
|
||
|
|
nonzero(C):- zero(C), !, fail.
|
||
|
|
nonzero(_).
|
||
|
|
|
||
|
|
unwrap( slack(S), X) ?- !, X=S.
|
||
|
|
unwrap( X, X).
|
||
|
|
|
||
|
|
is_div( C1, C2, C3) :-
|
||
|
|
unwrap( C1, C11),
|
||
|
|
unwrap( C2, C21),
|
||
|
|
unwrap( C3, C31),
|
||
|
|
is_divu( C11, C21, C31).
|
||
|
|
|
||
|
|
is_divu(C1,C2,C3):- zero(C1),!,C3=0.
|
||
|
|
is_divu(C1,C2,C3):- X is -(C1/C2), % minus here to get sign needed in handlers
|
||
|
|
avoid_float(X,C3).
|
||
|
|
|
||
|
|
is_mul( C1, C2, C3) :-
|
||
|
|
unwrap( C1, C11),
|
||
|
|
unwrap( C2, C21),
|
||
|
|
unwrap( C3, C31),
|
||
|
|
is_mulu( C11, C21, C31).
|
||
|
|
|
||
|
|
is_mulu(C1,C2,C3):- zero(C1),!,C3=0.
|
||
|
|
is_mulu(C1,C2,C3):- zero(C2),!,C3=0.
|
||
|
|
is_mulu(C1,C2,C3):- X is C1*C2,
|
||
|
|
avoid_float(X,C3).
|
||
|
|
|
||
|
|
avoid_float(X,C3):-
|
||
|
|
float(X) -> Y is round(X),Z is X-Y,(zero(Z)-> C3 is integer(Y);C3=X) ; C3=X.
|
||
|
|
|
||
|
|
|
||
|
|
simplifyable(X*C,P,P1):- delete(X*C,P,P1),ground(X),!.
|
||
|
|
|
||
|
|
|
||
|
|
% HANDLING SLACK VARIABLES ----------------------------------------------------
|
||
|
|
|
||
|
|
all_slacks([]).
|
||
|
|
all_slacks([slack(Sl)*C|P]) ?- % check_slack(Sl),
|
||
|
|
all_slacks(P).
|
||
|
|
|
||
|
|
all_slacks([],_).
|
||
|
|
all_slacks([slack(Sl)*C|P],S) ?- % check_slack(Sl),
|
||
|
|
sign(C,S),
|
||
|
|
all_slacks(P,S).
|
||
|
|
|
||
|
|
check_slack( S) :- find_constraint( S, basic(_)#_), !.
|
||
|
|
check_slack( _) :- raise_exception( slack).
|
||
|
|
|
||
|
|
sign(C,0):- zero(C),!.
|
||
|
|
sign(C,S):- C>0 -> S=1 ; S=(-1).
|
||
|
|
|
||
|
|
all_zeroes([]).
|
||
|
|
all_zeroes([slack(0)*C|P]) :-
|
||
|
|
all_zeroes(P).
|
||
|
|
|
||
|
|
|
||
|
|
% COMPUTING WITH POLYNOMIALS -------------------------------------------------
|
||
|
|
|
||
|
|
% gets rounded constant C from is_div/3
|
||
|
|
mult_const(eq0(C1,P1),C,eq0(0 ,[])):- C=:=0,!.
|
||
|
|
mult_const(eq0(C1,P1),C,eq0(C1,P1)):- C=:=1,!.
|
||
|
|
mult_const(eq0(C1,P1),C2,eq0(C,P)):-
|
||
|
|
(zero(C1) -> C=0 ; C is C1*C2),
|
||
|
|
mult_const1(P1,C2,P).
|
||
|
|
mult_const1([],C,[]).
|
||
|
|
mult_const1([Xi*Ci|Poly],C,PolyR):-
|
||
|
|
(zero(Ci) -> PolyR=NPoly ; NCi is Ci*C,PolyR=[Xi*NCi|NPoly]),
|
||
|
|
mult_const1(Poly,C,NPoly).
|
||
|
|
|
||
|
|
% gets input from const_mult/3
|
||
|
|
add_eq0(eq0(C1,P1),eq0(C2,P2),eq0(C,P0)):-
|
||
|
|
Ci is C1+C2,
|
||
|
|
(zero(Ci) -> C=0 ; C=Ci),
|
||
|
|
add_eq1(P1,P2,P0).
|
||
|
|
% sort(P,P0).
|
||
|
|
add_eq1([],Poly,Poly):-!.
|
||
|
|
add_eq1(Poly,[],Poly):-!.
|
||
|
|
add_eq1([Xi1*Ci1|Poly1],Poly21,Poly):-
|
||
|
|
delete(Xi2*Ci2,Poly21,Poly2),Xi2==Xi1,
|
||
|
|
!,
|
||
|
|
Ci is Ci1+Ci2,
|
||
|
|
(zero(Ci) -> Poly=Poly3 ; Poly=[Xi1*Ci|Poly3]),
|
||
|
|
add_eq1(Poly1,Poly2,Poly3).
|
||
|
|
add_eq1([Xi1*Ci1|Poly1],Poly2,[Xi1*Ci1|Poly3]):-
|
||
|
|
add_eq1(Poly1,Poly2,Poly3).
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
normalize(A,B,P2,C1):-
|
||
|
|
normalize1(A-B,P),
|
||
|
|
P=eq0(C1,P1),rev(P1,[],P1R),globalize(P1R),
|
||
|
|
sort1(P1,P2).
|
||
|
|
|
||
|
|
normalize1(V,P) ?- var(V),!,
|
||
|
|
P=eq0(0,[V*1]).
|
||
|
|
normalize1(C,P) ?- ground(C),!,
|
||
|
|
C1 is C,P=eq0(C1,[]).
|
||
|
|
normalize1(slack(V),P) ?- !,
|
||
|
|
P=eq0(0,[slack(V)*1]).
|
||
|
|
normalize1((+E),P) ?-!,
|
||
|
|
normalize1(E,P).
|
||
|
|
normalize1((-E),P) ?-!,
|
||
|
|
normalize1(E,P1),
|
||
|
|
mult_const(P1,(-1),P).
|
||
|
|
normalize1(A*B,C) ?- ground(A),!,
|
||
|
|
normalize1(B,BN),
|
||
|
|
mult_const(BN,A,C).
|
||
|
|
normalize1(B*A,C) ?- ground(A),!,
|
||
|
|
normalize1(B,BN),
|
||
|
|
mult_const(BN,A,C).
|
||
|
|
normalize1(B/A,C) ?- ground(A),!,
|
||
|
|
normalize1(B,BN),
|
||
|
|
A1 is 1/A,
|
||
|
|
mult_const(BN,A1,C).
|
||
|
|
normalize1(A-B,C) ?- !,
|
||
|
|
normalize1(A,AN),
|
||
|
|
normalize1((-B),BN),
|
||
|
|
add_eq0(AN,BN,C).
|
||
|
|
normalize1(A+B,C) ?- !,
|
||
|
|
normalize1(A,AN),
|
||
|
|
normalize1(B,BN),
|
||
|
|
add_eq0(AN,BN,C).
|
||
|
|
normalize1(E,C) ?-
|
||
|
|
C=eq0(0,[CX*1]),
|
||
|
|
eqnonlin(CX,E). % add a nonlinear equation constraint
|
||
|
|
|
||
|
|
|
||
|
|
% end of file math-utilities.pl -----------------------------------------------
|