1/*  $Id$
    2
    3    Part of CLP(Q) (Constraint Logic Programming over Rationals)
    4
    5    Author:        Leslie De Koninck
    6    E-mail:        Leslie.DeKoninck@cs.kuleuven.be
    7    WWW:           http://www.swi-prolog.org
    8		   http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09
    9    Copyright (C): 2006, K.U. Leuven and
   10		   1992-1995, Austrian Research Institute for
   11		              Artificial Intelligence (OFAI),
   12			      Vienna, Austria
   13
   14    This software is based on CLP(Q,R) by Christian Holzbaur for SICStus
   15    Prolog and distributed under the license details below with permission from
   16    all mentioned authors.
   17
   18    This program is free software; you can redistribute it and/or
   19    modify it under the terms of the GNU General Public License
   20    as published by the Free Software Foundation; either version 2
   21    of the License, or (at your option) any later version.
   22
   23    This program is distributed in the hope that it will be useful,
   24    but WITHOUT ANY WARRANTY; without even the implied warranty of
   25    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   26    GNU General Public License for more details.
   27
   28    You should have received a copy of the GNU Lesser General Public
   29    License along with this library; if not, write to the Free Software
   30    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
   31
   32    As a special exception, if you link this library with other files,
   33    compiled with a Free Software compiler, to produce an executable, this
   34    library does not by itself cause the resulting executable to be covered
   35    by the GNU General Public License. This exception does not however
   36    invalidate any other reasons why the executable file might be covered by
   37    the GNU General Public License.
   38*/
   39
   40:- module(store_q,
   41	[
   42	    add_linear_11/3,
   43	    add_linear_f1/4,
   44	    add_linear_ff/5,
   45	    normalize_scalar/2,
   46	    delete_factor/4,
   47	    mult_linear_factor/3,
   48	    nf_rhs_x/4,
   49	    indep/2,
   50	    isolate/3,
   51	    nf_substitute/4,
   52	    mult_hom/3,
   53	    nf2sum/3,
   54	    nf_coeff_of/3,
   55	    renormalize/2	
   56	]).   57
   58% normalize_scalar(S,[N,Z])
   59%
   60% Transforms a scalar S into a linear expression [S,0]
   61
   62normalize_scalar(S,[S,0]).
   63
   64% renormalize(List,Lin)
   65%
   66% Renormalizes the not normalized linear expression in List into
   67% a normalized one. It does so to take care of unifications.
   68% (e.g. when a variable X is bound to a constant, the constant is added to
   69% the constant part of the linear expression; when a variable X is bound to
   70% another variable Y, the scalars of both are added)
   71
   72renormalize([I,R|Hom],Lin) :-
   73	length(Hom,Len),
   74	renormalize_log(Len,Hom,[],Lin0),
   75	add_linear_11([I,R],Lin0,Lin).
   76
   77% renormalize_log(Len,Hom,HomTail,Lin)
   78%
   79% Logarithmically renormalizes the homogene part of a not normalized
   80% linear expression. See also renormalize/2.
   81
   82renormalize_log(1,[Term|Xs],Xs,Lin) :-
   83	!,
   84	Term = l(X*_,_),
   85	renormalize_log_one(X,Term,Lin).
   86renormalize_log(2,[A,B|Xs],Xs,Lin) :-
   87	!,
   88	A = l(X*_,_),
   89	B = l(Y*_,_),
   90	renormalize_log_one(X,A,LinA),
   91	renormalize_log_one(Y,B,LinB),
   92	add_linear_11(LinA,LinB,Lin).
   93renormalize_log(N,L0,L2,Lin) :-
   94	P is N>>1,
   95	Q is N-P,
   96	renormalize_log(P,L0,L1,Lp),
   97	renormalize_log(Q,L1,L2,Lq),
   98	add_linear_11(Lp,Lq,Lin).
   99
  100% renormalize_log_one(X,Term,Res)
  101%
  102% Renormalizes a term in X: if X is a nonvar, the term becomes a scalar.
  103
  104renormalize_log_one(X,Term,Res) :-
  105	var(X),
  106	Term = l(X*K,_),
  107	get_attr(X,clpqr_itf,Att),
  108	arg(5,Att,order(OrdX)), % Order might have changed
  109	Res = [0,0,l(X*K,OrdX)].
  110renormalize_log_one(X,Term,Res) :-
  111	nonvar(X),
  112	Term = l(X*K,_),
  113	Xk is X*K,
  114	normalize_scalar(Xk,Res).
  115
  116% ----------------------------- sparse vector stuff ---------------------------- %
  117
  118% add_linear_ff(LinA,Ka,LinB,Kb,LinC)
  119%
  120% Linear expression LinC is the result of the addition of the 2 linear expressions
  121% LinA and LinB, each one multiplied by a scalar (Ka for LinA and Kb for LinB).
  122
  123add_linear_ff(LinA,Ka,LinB,Kb,LinC) :-
  124	LinA = [Ia,Ra|Ha],
  125	LinB = [Ib,Rb|Hb],
  126	LinC = [Ic,Rc|Hc],
  127	Ic is Ia*Ka+Ib*Kb,
  128	Rc is Ra*Ka+Rb*Kb,
  129 	add_linear_ffh(Ha,Ka,Hb,Kb,Hc).
  130
  131% add_linear_ffh(Ha,Ka,Hb,Kb,Hc)
  132%
  133% Homogene part Hc is the result of the addition of the 2 homogene parts Ha and Hb,
  134% each one multiplied by a scalar (Ka for Ha and Kb for Hb)
  135
  136add_linear_ffh([],_,Ys,Kb,Zs) :- mult_hom(Ys,Kb,Zs).
  137add_linear_ffh([l(X*Kx,OrdX)|Xs],Ka,Ys,Kb,Zs) :-
  138	add_linear_ffh(Ys,X,Kx,OrdX,Xs,Zs,Ka,Kb).
  139
  140% add_linear_ffh(Ys,X,Kx,OrdX,Xs,Zs,Ka,Kb)
  141%
  142% Homogene part Zs is the result of the addition of the 2 homogene parts Ys and
  143% [l(X*Kx,OrdX)|Xs], each one multiplied by a scalar (Ka for [l(X*Kx,OrdX)|Xs] and Kb for Ys)
  144
  145add_linear_ffh([],X,Kx,OrdX,Xs,Zs,Ka,_) :- mult_hom([l(X*Kx,OrdX)|Xs],Ka,Zs).
  146add_linear_ffh([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs,Ka,Kb) :-
  147	compare(Rel,OrdX,OrdY),
  148	(   Rel = (=)
  149	->  Kz is Kx*Ka+Ky*Kb,
  150	    (   Kz =:= 0
  151	    ->  add_linear_ffh(Xs,Ka,Ys,Kb,Zs)
  152	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  153		add_linear_ffh(Xs,Ka,Ys,Kb,Ztail)
  154	    )
  155	;   Rel = (<)
  156	->  Zs = [l(X*Kz,OrdX)|Ztail],
  157	    Kz is Kx*Ka,
  158	    add_linear_ffh(Xs,Y,Ky,OrdY,Ys,Ztail,Kb,Ka)
  159	;   Rel = (>)
  160	->  Zs = [l(Y*Kz,OrdY)|Ztail],
  161	    Kz is Ky*Kb,
  162	    add_linear_ffh(Ys,X,Kx,OrdX,Xs,Ztail,Ka,Kb)
  163     	).
  164
  165% add_linear_f1(LinA,Ka,LinB,LinC)
  166%
  167% special case of add_linear_ff with Kb = 1
  168
  169add_linear_f1(LinA,Ka,LinB,LinC) :-
  170	LinA = [Ia,Ra|Ha],
  171	LinB = [Ib,Rb|Hb],
  172	LinC = [Ic,Rc|Hc],
  173	Ic is Ia*Ka+Ib,
  174	Rc is Ra*Ka+Rb,
  175	add_linear_f1h(Ha,Ka,Hb,Hc).
  176
  177% add_linear_f1h(Ha,Ka,Hb,Hc)
  178%
  179% special case of add_linear_ffh/5 with Kb = 1
  180
  181add_linear_f1h([],_,Ys,Ys).
  182add_linear_f1h([l(X*Kx,OrdX)|Xs],Ka,Ys,Zs) :-
  183	add_linear_f1h(Ys,X,Kx,OrdX,Xs,Zs,Ka).
  184
  185% add_linear_f1h(Ys,X,Kx,OrdX,Xs,Zs,Ka)
  186%
  187% special case of add_linear_ffh/8 with Kb = 1
  188
  189add_linear_f1h([],X,Kx,OrdX,Xs,Zs,Ka) :- mult_hom([l(X*Kx,OrdX)|Xs],Ka,Zs).
  190add_linear_f1h([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs,Ka) :-
  191	compare(Rel,OrdX,OrdY),
  192	(   Rel = (=)
  193	->  Kz is Kx*Ka+Ky,
  194	    (   Kz =:= 0
  195	    ->  add_linear_f1h(Xs,Ka,Ys,Zs)
  196	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  197		add_linear_f1h(Xs,Ka,Ys,Ztail)
  198	    )
  199	;   Rel = (<)
  200	->  Zs = [l(X*Kz,OrdX)|Ztail],
  201	    Kz is Kx*Ka,
  202	    add_linear_f1h(Xs,Ka,[l(Y*Ky,OrdY)|Ys],Ztail)
  203 	;   Rel = (>)
  204	->  Zs = [l(Y*Ky,OrdY)|Ztail],
  205	    add_linear_f1h(Ys,X,Kx,OrdX,Xs,Ztail,Ka)
  206	).
  207
  208% add_linear_11(LinA,LinB,LinC)
  209%
  210% special case of add_linear_ff with Ka = 1 and Kb = 1
  211
  212add_linear_11(LinA,LinB,LinC) :-
  213	LinA = [Ia,Ra|Ha],
  214	LinB = [Ib,Rb|Hb],
  215	LinC = [Ic,Rc|Hc],
  216	Ic is Ia+Ib,
  217	Rc is Ra+Rb,
  218	add_linear_11h(Ha,Hb,Hc).
  219
  220% add_linear_11h(Ha,Hb,Hc)
  221%
  222% special case of add_linear_ffh/5 with Ka = 1 and Kb = 1
  223
  224add_linear_11h([],Ys,Ys).
  225add_linear_11h([l(X*Kx,OrdX)|Xs],Ys,Zs) :-
  226	add_linear_11h(Ys,X,Kx,OrdX,Xs,Zs).
  227
  228% add_linear_11h(Ys,X,Kx,OrdX,Xs,Zs)
  229%
  230% special case of add_linear_ffh/8 with Ka = 1 and Kb = 1
  231
  232add_linear_11h([],X,Kx,OrdX,Xs,[l(X*Kx,OrdX)|Xs]).
  233add_linear_11h([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs) :-
  234	compare(Rel,OrdX,OrdY),
  235	(   Rel = (=)
  236	->  Kz is Kx+Ky,
  237	    (   Kz =:= 0
  238	    ->  add_linear_11h(Xs,Ys,Zs)
  239	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  240		add_linear_11h(Xs,Ys,Ztail)
  241	    )
  242	;   Rel = (<)
  243	->  Zs = [l(X*Kx,OrdX)|Ztail],
  244	    add_linear_11h(Xs,Y,Ky,OrdY,Ys,Ztail)
  245	;   Rel = (>)
  246	->  Zs = [l(Y*Ky,OrdY)|Ztail],
  247	    add_linear_11h(Ys,X,Kx,OrdX,Xs,Ztail)
  248	).
  249
  250% mult_linear_factor(Lin,K,Res)
  251%
  252% Linear expression Res is the result of multiplication of linear
  253% expression Lin by scalar K
  254
  255mult_linear_factor(Lin,K,Mult) :-
  256	K =:= 1,
  257	!,
  258	Mult = Lin.
  259mult_linear_factor(Lin,K,Res) :-
  260	Lin = [I,R|Hom],
  261	Res = [Ik,Rk|Mult],
  262	Ik is I*K,
  263	Rk is R*K,
  264	mult_hom(Hom,K,Mult).
  265
  266% mult_hom(Hom,K,Res)
  267%
  268% Homogene part Res is the result of multiplication of homogene part
  269% Hom by scalar K
  270
  271mult_hom([],_,[]).
  272mult_hom([l(A*Fa,OrdA)|As],F,[l(A*Fan,OrdA)|Afs]) :-
  273	Fan is F*Fa,
  274	mult_hom(As,F,Afs).
  275
  276% nf_substitute(Ord,Def,Lin,Res)
  277%
  278% Linear expression Res is the result of substitution of Var in
  279% linear expression Lin, by its definition in the form of linear
  280% expression Def
  281
  282nf_substitute(OrdV,LinV,LinX,LinX1) :-
  283	delete_factor(OrdV,LinX,LinW,K),
  284	add_linear_f1(LinV,K,LinW,LinX1).
  285
  286% delete_factor(Ord,Lin,Res,Coeff)
  287%
  288% Linear expression Res is the result of the deletion of the term
  289% Var*Coeff where Var has ordering Ord from linear expression Lin
  290
  291delete_factor(OrdV,Lin,Res,Coeff) :-
  292	Lin = [I,R|Hom],
  293	Res = [I,R|Hdel],
  294	delete_factor_hom(OrdV,Hom,Hdel,Coeff).
  295
  296% delete_factor_hom(Ord,Hom,Res,Coeff)
  297%
  298% Homogene part Res is the result of the deletion of the term
  299% Var*Coeff from homogene part Hom
  300
  301delete_factor_hom(VOrd,[Car|Cdr],RCdr,RKoeff) :-
  302	Car = l(_*Koeff,Ord),
  303	compare(Rel,VOrd,Ord),
  304	(   Rel= (=)
  305	->  RCdr = Cdr,
  306	    RKoeff=Koeff
  307	;   Rel= (>)
  308	->  RCdr = [Car|RCdr1],
  309	    delete_factor_hom(VOrd,Cdr,RCdr1,RKoeff)
  310	).
  311
  312
  313% nf_coeff_of(Lin,OrdX,Coeff)
  314%
  315% Linear expression Lin contains the term l(X*Coeff,OrdX)
  316
  317nf_coeff_of([_,_|Hom],VOrd,Coeff) :-
  318	nf_coeff_hom(Hom,VOrd,Coeff).
  319
  320% nf_coeff_hom(Lin,OrdX,Coeff)
  321%
  322% Linear expression Lin contains the term l(X*Coeff,OrdX) where the
  323% order attribute of X = OrdX
  324
  325nf_coeff_hom([l(_*K,OVar)|Vs],OVid,Coeff) :-
  326	compare(Rel,OVid,OVar),
  327	(   Rel = (=)
  328	->  Coeff = K
  329	;   Rel = (>)
  330	->  nf_coeff_hom(Vs,OVid,Coeff)
  331	).
  332
  333% nf_rhs_x(Lin,OrdX,Rhs,K)
  334%
  335% Rhs = R + I where Lin = [I,R|Hom] and l(X*K,OrdX) is a term of Hom
  336
  337nf_rhs_x(Lin,OrdX,Rhs,K) :-
  338	Lin = [I,R|Tail],
  339	nf_coeff_hom(Tail,OrdX,K),
  340	Rhs is R+I.	% late because X may not occur in H
  341
  342% isolate(OrdN,Lin,Lin1)
  343%
  344% Linear expression Lin1 is the result of the transformation of linear expression
  345% Lin = 0 which contains the term l(New*K,OrdN) into an equivalent expression Lin1 = New.
  346
  347isolate(OrdN,Lin,Lin1) :-
  348	delete_factor(OrdN,Lin,Lin0,Coeff),
  349	K is -1 rdiv Coeff,
  350	mult_linear_factor(Lin0,K,Lin1).
  351
  352% indep(Lin,OrdX)
  353%
  354% succeeds if Lin = [0,_|[l(X*1,OrdX)]]
  355
  356indep(Lin,OrdX) :-
  357	Lin = [I,_|[l(_*K,OrdY)]],
  358	OrdX == OrdY,
  359	K =:= 1,
  360	I =:= 0.
  361
  362% nf2sum(Lin,Sofar,Term)
  363%
  364% Transforms a linear expression into a sum
  365% (e.g. the expression [5,_,[l(X*2,OrdX),l(Y*-1,OrdY)]] gets transformed into 5 + 2*X - Y)
  366
  367nf2sum([],I,I).
  368nf2sum([X|Xs],I,Sum) :-
  369	(   I =:= 0
  370	->  X = l(Var*K,_),
  371 	    (   K =:= 1
  372	    ->  hom2sum(Xs,Var,Sum)
  373	    ;   K =:= -1
  374	    ->  hom2sum(Xs,-Var,Sum)
  375	    ;	hom2sum(Xs,K*Var,Sum)
  376	    )
  377	;   hom2sum([X|Xs],I,Sum)
  378 	).
  379
  380% hom2sum(Hom,Sofar,Term)
  381%
  382% Transforms a linear expression into a sum
  383% this predicate handles all but the first term
  384% (the first term does not need a concatenation symbol + or -)
  385% see also nf2sum/3
  386
  387hom2sum([],Term,Term).
  388hom2sum([l(Var*K,_)|Cs],Sofar,Term) :-
  389	(   K =:= 1
  390	->  Next = Sofar + Var
  391	;   K =:= -1
  392	->  Next = Sofar - Var
  393	;   K < 0
  394	->  Ka is -K,
  395	    Next = Sofar - Ka*Var
  396	;   Next = Sofar + K*Var
  397	),
  398	hom2sum(Cs,Next,Term)