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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  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,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).