% 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 -----------------------------------------------
