1/*
    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(bv_q,
   41	[
   42	    allvars/2,
   43	    backsubst/3,
   44	    backsubst_delta/4,
   45	    basis_add/2,
   46	    dec_step/2,
   47	    deref/2,
   48	    deref_var/2,
   49	    detach_bounds/1,
   50	    detach_bounds_vlv/5,
   51	    determine_active_dec/1,
   52	    determine_active_inc/1,
   53	    dump_var/6,
   54	    dump_nz/5,
   55	    export_binding/1,
   56	    get_or_add_class/2,
   57	    inc_step/2,
   58	    intro_at/3,
   59	    iterate_dec/2,
   60	    lb/3,
   61	    pivot_a/4,
   62	    pivot/5,
   63	    rcbl_status/6,
   64	    reconsider/1,
   65	    same_class/2,
   66	    solve/1,
   67	    solve_ord_x/3,
   68	    ub/3,
   69	    unconstrained/4,
   70	    var_intern/2,
   71	    var_intern/3,
   72	    var_with_def_assign/2,
   73	    var_with_def_intern/4,
   74	    maximize/1,
   75	    minimize/1,
   76	    sup/2,
   77	    sup/4,
   78	    inf/2,
   79	    inf/4,
   80	    'solve_<'/1,
   81	    'solve_=<'/1,
   82	    'solve_=\\='/1,
   83	    log_deref/4
   84	]).   85:- use_module(store_q,
   86	[
   87	    add_linear_11/3,
   88	    add_linear_f1/4,
   89	    add_linear_ff/5,
   90	    delete_factor/4,
   91	    indep/2,
   92	    isolate/3,
   93	    nf2sum/3,
   94	    nf_rhs_x/4,
   95	    nf_substitute/4,
   96	    normalize_scalar/2,
   97	    mult_hom/3,
   98	    mult_linear_factor/3
   99	]).  100:- use_module('../clpqr/class',
  101	[
  102	    class_allvars/2,
  103	    class_basis/2,
  104	    class_basis_add/3,
  105	    class_basis_drop/2,
  106	    class_basis_pivot/3,
  107	    class_new/5
  108	]).  109:- use_module(ineq_q,
  110	[
  111	    ineq/4
  112	]).  113:- use_module(nf_q,
  114	[
  115	    {}/1,
  116	    split/3,
  117	    wait_linear/3
  118	]).  119:- use_module(bb_q,
  120	[
  121	    vertex_value/2
  122	]).  123:- use_module(library(ordsets),
  124	[
  125	    ord_add_element/3
  126	]).  127
  128% For the rhs maint. the following events are important:
  129%
  130%	-) introduction of an indep var at active bound B
  131%	-) narrowing of active bound
  132%	-) swap active bound
  133%	-) pivot
  134%
  135
  136% a variables bound (L/U) can have the states:
  137%
  138%	-) t_none	no bounds
  139%	-) t_l		inactive lower bound
  140%	-) t_u		inactive upper bound
  141%	-) t_L		active lower bound
  142%	-) t_U		active upper bound
  143%	-) t_lu		inactive lower and upper bound
  144%	-) t_Lu		active lower bound and inactive upper bound
  145%	-) t_lU		inactive lower bound and active upper bound
  146
  147% ----------------------------------- deref -----------------------------------
  148%
  149
  150% deref(Lin,Lind)
  151%
  152% Makes a linear equation of the form [v(I,[])|H] into a solvable linear
  153% equation.
  154% If the variables are new, they are initialized with the linear equation X=X.
  155
  156deref(Lin,Lind) :-
  157	split(Lin,H,I),
  158	normalize_scalar(I,Nonvar),
  159	length(H,Len),
  160	log_deref(Len,H,[],Restd),
  161	add_linear_11(Nonvar,Restd,Lind).
  162
  163% log_deref(Len,[Vs|VsTail],VsTail,Res)
  164%
  165% Logarithmically converts a linear equation in normal form ([v(_,_)|_]) into a
  166% linear equation in solver form ([I,R,K*X|_]). Res contains the result, Len is
  167% the length of the part to convert and [Vs|VsTail] is a difference list
  168% containing the equation in normal form.
  169
  170log_deref(0,Vs,Vs,Lin) :-
  171	!,
  172	Lin = [0,0].
  173log_deref(1,[v(K,[X^1])|Vs],Vs,Lin) :-
  174	!,
  175	deref_var(X,Lx),
  176	mult_linear_factor(Lx,K,Lin).
  177log_deref(2,[v(Kx,[X^1]),v(Ky,[Y^1])|Vs],Vs,Lin) :-
  178	!,
  179	deref_var(X,Lx),
  180	deref_var(Y,Ly),
  181	add_linear_ff(Lx,Kx,Ly,Ky,Lin).
  182log_deref(N,V0,V2,Lin) :-
  183	P is N >> 1,
  184	Q is N - P,
  185	log_deref(P,V0,V1,Lp),
  186	log_deref(Q,V1,V2,Lq),
  187	add_linear_11(Lp,Lq,Lin).
  188
  189% deref_var(X,Lin)
  190%
  191% Returns the equation of variable X. If X is a new variable, a new equation
  192% X = X is made.
  193
  194deref_var(X,Lin) :-
  195	(   get_attr(X,clpqr_itf,Att)
  196	->  (   \+ arg(1,Att,clpq)
  197	    ->  throw(error(permission_error('mix CLP(Q) variables with',
  198		'CLP(R) variables:',X),context(_)))
  199	    ;   arg(4,Att,lin(Lin))
  200	    ->  true
  201	    ;   setarg(2,Att,type(t_none)),
  202		setarg(3,Att,strictness(0)),
  203		Lin = [0,0,l(X*1,Ord)],
  204		setarg(4,Att,lin(Lin)),
  205		setarg(5,Att,order(Ord))
  206	    )
  207	;   Lin = [0,0,l(X*1,Ord)],
  208	    put_attr(X,clpqr_itf,t(clpq,type(t_none),strictness(0),
  209		lin(Lin),order(Ord),n,n,n,n,n,n))
  210	).
  211
  212% TODO
  213%
  214%
  215
  216var_with_def_assign(Var,Lin) :-
  217	Lin = [I,_|Hom],
  218	(   Hom = []
  219	->  % X=k
  220	    Var = I
  221	;   Hom = [l(V*K,_)|Cs]
  222	->  (   Cs = [],
  223		K =:= 1,
  224		I =:= 0
  225	    ->	% X=Y
  226		Var = V
  227	    ;	% general case
  228		var_with_def_intern(t_none,Var,Lin,0)
  229	    )
  230	).
  231
  232% var_with_def_intern(Type,Var,Lin,Strictness)
  233%
  234% Makes Lin the linear equation of new variable Var, makes all variables of
  235% Lin, and Var of the same class and bounds Var by type(Type) and
  236% strictness(Strictness)
  237
  238var_with_def_intern(Type,Var,Lin,Strict) :-
  239	put_attr(Var,clpqr_itf,t(clpq,type(Type),strictness(Strict),lin(Lin),
  240	    order(_),n,n,n,n,n,n)),	% check uses
  241	Lin = [_,_|Hom],
  242	get_or_add_class(Var,Class),
  243	same_class(Hom,Class).
  244
  245% TODO
  246%
  247%
  248
  249var_intern(Type,Var,Strict) :-
  250	put_attr(Var,clpqr_itf,t(clpq,type(Type),strictness(Strict),
  251	    lin([0,0,l(Var*1,Ord)]),order(Ord),n,n,n,n,n,n)),
  252	get_or_add_class(Var,_Class).
  253
  254% TODO
  255%
  256%
  257
  258var_intern(Var,Class) :-	% for ordered/1 but otherwise free vars
  259	get_attr(Var,clpqr_itf,Att),
  260	arg(2,Att,type(_)),
  261	arg(4,Att,lin(_)),
  262	!,
  263	get_or_add_class(Var,Class).
  264var_intern(Var,Class) :-
  265	put_attr(Var,clpqr_itf,t(clpq,type(t_none),strictness(0),
  266	    lin([0,0,l(Var*1,Ord)]),order(Ord),n,n,n,n,n,n)),
  267	get_or_add_class(Var,Class).
  268
  269% -----------------------------------------------------------------------------
  270
  271% export_binding(Lst)
  272%
  273% Binds variables X to Y where Lst contains elements of the form [X-Y].
  274
  275export_binding([]).
  276export_binding([X-Y|Gs]) :-
  277	Y = X,
  278	export_binding(Gs).
  279
  280% 'solve_='(Nf)
  281%
  282% Solves linear equation Nf = 0 where Nf is in normal form.
  283
  284'solve_='(Nf) :-
  285	deref(Nf,Nfd),	% dereferences and turns Nf into solvable form Nfd
  286	solve(Nfd).
  287
  288% 'solve_=\\='(Nf)
  289%
  290% Solves linear inequality Nf =\= 0 where Nf is in normal form.
  291
  292'solve_=\\='(Nf) :-
  293	deref(Nf,Lind),	% dereferences and turns Nf into solvable form Lind
  294	Lind = [Inhom,_|Hom],
  295	(   Hom = []
  296	->  Inhom =\= 0
  297	;   % make new variable Nz = Lind
  298	    var_with_def_intern(t_none,Nz,Lind,0),
  299	    % make Nz nonzero
  300	    get_attr(Nz,clpqr_itf,Att),
  301	    setarg(8,Att,nonzero)
  302	).
  303
  304% 'solve_<'(Nf)
  305%
  306% Solves linear inequality Nf < 0 where Nf is in normal form.
  307
  308'solve_<'(Nf) :-
  309	split(Nf,H,I),
  310	ineq(H,I,Nf,strict).
  311
  312% 'solve_=<'(Nf)
  313%
  314% Solves linear inequality Nf =< 0 where Nf is in normal form.
  315
  316'solve_=<'(Nf) :-
  317	split(Nf,H,I),
  318	ineq(H,I,Nf,nonstrict).
  319
  320maximize(Term) :-
  321	minimize(-Term).
  322
  323%
  324% This is NOT coded as minimize(Expr) :- inf(Expr,Expr).
  325%
  326% because the new version of inf/2 only visits
  327% the vertex where the infimum is assumed and returns
  328% to the 'current' vertex via backtracking.
  329% The rationale behind this construction is to eliminate
  330% all garbage in the solver data structures produced by
  331% the pivots on the way to the extremal point caused by
  332% {inf,sup}/{2,4}.
  333%
  334% If we are after the infimum/supremum for minimizing/maximizing,
  335% this strategy may have adverse effects on performance because
  336% the simplex algorithm is forced to re-discover the
  337% extremal vertex through the equation {Inf =:= Expr}.
  338%
  339% Thus the extra code for {minimize,maximize}/1.
  340%
  341% In case someone comes up with an example where
  342%
  343%   inf(Expr,Expr)
  344%
  345% outperforms the provided formulation for minimize - so be it.
  346% Both forms are available to the user.
  347%
  348minimize(Term) :-
  349	wait_linear(Term,Nf,minimize_lin(Nf)).
  350
  351% minimize_lin(Lin)
  352%
  353% Minimizes the linear expression Lin. It does so by making a new
  354% variable Dep and minimizes its value.
  355
  356minimize_lin(Lin) :-
  357	deref(Lin,Lind),
  358	var_with_def_intern(t_none,Dep,Lind,0),
  359	determine_active_dec(Lind),
  360	iterate_dec(Dep,Inf),
  361	{ Dep =:= Inf }.
  362
  363sup(Expression,Sup) :-
  364	sup(Expression,Sup,[],[]).
  365
  366sup(Expression,Sup,Vector,Vertex) :-
  367	inf(-Expression,-Sup,Vector,Vertex).
  368
  369inf(Expression,Inf) :-
  370	inf(Expression,Inf,[],[]).
  371
  372inf(Expression,Inf,Vector,Vertex) :-
  373	% wait until Expression becomes linear, Nf contains linear Expression
  374	% in normal form
  375	wait_linear(Expression,Nf,inf_lin(Nf,Inf,Vector,Vertex)).
  376
  377inf_lin(Lin,_,Vector,_) :-
  378	deref(Lin,Lind),
  379	var_with_def_intern(t_none,Dep,Lind,0),	% make new variable Dep = Lind
  380	determine_active_dec(Lind),	% minimizes Lind
  381	iterate_dec(Dep,Inf),
  382	vertex_value(Vector,Values),
  383	nb_setval(inf,[Inf|Values]),
  384	fail.
  385inf_lin(_,Infimum,_,Vertex) :-
  386	nb_current(inf,L),
  387	nb_delete(inf),
  388	assign([Infimum|Vertex],L).
  389
  390% assign(L1,L2)
  391%
  392% The elements of L1 are pairwise assigned to the elements of L2
  393% by means of asserting {X =:= Y} where X is an element of L1 and Y
  394% is the corresponding element of L2.
  395
  396assign([],[]).
  397assign([X|Xs],[Y|Ys]) :-
  398	{X =:= Y},		  % more defensive/expressive than X=Y
  399	assign(Xs,Ys).
  400
  401% --------------------------------- optimization ------------------------------
  402%
  403% The _sn(S) =< 0 row might be temporarily infeasible.
  404% We use reconsider/1 to fix this.
  405%
  406%   s(S) e [_,0] = d +xi ... -xj, Rhs > 0 so we want to decrease s(S)
  407%
  408%   positive xi would have to be moved towards their lower bound,
  409%   negative xj would have to be moved towards their upper bound,
  410%
  411%   the row s(S) does not limit the lower bound of xi
  412%   the row s(S) does not limit the upper bound of xj
  413%
  414%   a) if some other row R is limiting xk, we pivot(R,xk),
  415%      s(S) will decrease and get more feasible until (b)
  416%   b) if there is no limiting row for some xi: we pivot(s(S),xi)
  417%					    xj: we pivot(s(S),xj)
  418%      which cures the infeasibility in one step
  419%
  420
  421
  422% iterate_dec(OptVar,Opt)
  423%
  424% Decreases the bound on the variables of the linear equation of OptVar as much
  425% as possible and returns the resulting optimal bound in Opt. Fails if for some
  426% variable, a status of unlimited is found.
  427
  428iterate_dec(OptVar,Opt) :-
  429	get_attr(OptVar,clpqr_itf,Att),
  430	arg(4,Att,lin([I,R|H])),
  431	dec_step(H,Status),
  432	(   Status = applied
  433	->  iterate_dec(OptVar,Opt)
  434	;   Status = optimum,
  435	    Opt is R + I
  436	).
  437
  438% iterate_inc(OptVar,Opt)
  439%
  440% Increases the bound on the variables of the linear equation of OptVar as much
  441% as possible and returns the resulting optimal bound in Opt. Fails if for some
  442% variable, a status of unlimited is found.
  443
  444iterate_inc(OptVar,Opt) :-
  445	get_attr(OptVar,clpqr_itf,Att),
  446	arg(4,Att,lin([I,R|H])),
  447	inc_step(H,Status),
  448	(   Status = applied
  449	->  iterate_inc(OptVar,Opt)
  450	;   Status = optimum,
  451	    Opt is R + I
  452	).
  453
  454%
  455% Status = {optimum,unlimited(Indep,DepT),applied}
  456% If Status = optimum, the tables have not been changed at all.
  457% Searches left to right, does not try to find the 'best' pivot
  458% Therefore we might discover unboundedness only after a few pivots
  459%
  460
  461
  462dec_step_cont([],optimum,Cont,Cont).
  463dec_step_cont([l(V*K,OrdV)|Vs],Status,ContIn,ContOut) :-
  464	get_attr(V,clpqr_itf,Att),
  465	arg(2,Att,type(W)),
  466	arg(6,Att,class(Class)),
  467	(   dec_step_2_cont(W,l(V*K,OrdV),Class,Status,ContIn,ContOut)
  468	->  true
  469	;   dec_step_cont(Vs,Status,ContIn,ContOut)
  470	).
  471
  472inc_step_cont([],optimum,Cont,Cont).
  473inc_step_cont([l(V*K,OrdV)|Vs],Status,ContIn,ContOut) :-
  474	get_attr(V,clpqr_itf,Att),
  475	arg(2,Att,type(W)),
  476	arg(6,Att,class(Class)),
  477	(   inc_step_2_cont(W,l(V*K,OrdV),Class,Status,ContIn,ContOut)
  478	->  true
  479	;   inc_step_cont(Vs,Status,ContIn,ContOut)
  480	).
  481
  482dec_step_2_cont(t_U(U),l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  483	K > 0,
  484	(   lb(Class,OrdV,Vub-Vb-_)
  485	->  % found a lower bound
  486	    Status = applied,
  487	    pivot_a(Vub,V,Vb,t_u(U)),
  488	    replace_in_cont(ContIn,Vub,V,ContOut)
  489	;   Status = unlimited(V,t_u(U)),
  490	    ContIn = ContOut
  491	).
  492dec_step_2_cont(t_lU(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  493	K > 0,
  494	Init is L - U,
  495	class_basis(Class,Deps),
  496	lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  497	pivot_b(Vub,V,Vb,t_lu(L,U)),
  498	replace_in_cont(ContIn,Vub,V,ContOut).
  499dec_step_2_cont(t_L(L),l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  500	K < 0,
  501	(   ub(Class,OrdV,Vub-Vb-_)
  502	->  Status = applied,
  503	    pivot_a(Vub,V,Vb,t_l(L)),
  504	    replace_in_cont(ContIn,Vub,V,ContOut)
  505	;   Status = unlimited(V,t_l(L)),
  506	    ContIn = ContOut
  507	).
  508dec_step_2_cont(t_Lu(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  509	K < 0,
  510	Init is U - L,
  511	class_basis(Class,Deps),
  512	ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  513	pivot_b(Vub,V,Vb,t_lu(L,U)),
  514	replace_in_cont(ContIn,Vub,V,ContOut).
  515dec_step_2_cont(t_none,l(V*_,_),_,unlimited(V,t_none),Cont,Cont).
  516
  517
  518
  519inc_step_2_cont(t_U(U),l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  520	K < 0,
  521	(   lb(Class,OrdV,Vub-Vb-_)
  522	->  Status = applied,
  523	    pivot_a(Vub,V,Vb,t_u(U)),
  524	    replace_in_cont(ContIn,Vub,V,ContOut)
  525	;   Status = unlimited(V,t_u(U)),
  526	    ContIn = ContOut
  527	).
  528inc_step_2_cont(t_lU(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  529	K < 0,
  530	Init is L - U,
  531	class_basis(Class,Deps),
  532	lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  533	pivot_b(Vub,V,Vb,t_lu(L,U)),
  534	replace_in_cont(ContIn,Vub,V,ContOut).
  535inc_step_2_cont(t_L(L),l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  536	K > 0,
  537	(   ub(Class,OrdV,Vub-Vb-_)
  538	->  Status = applied,
  539	    pivot_a(Vub,V,Vb,t_l(L)),
  540	    replace_in_cont(ContIn,Vub,V,ContOut)
  541	;   Status = unlimited(V,t_l(L)),
  542	    ContIn = ContOut
  543	).
  544inc_step_2_cont(t_Lu(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  545	K > 0,
  546	Init is U - L,
  547	class_basis(Class,Deps),
  548	ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  549	pivot_b(Vub,V,Vb,t_lu(L,U)),
  550	replace_in_cont(ContIn,Vub,V,ContOut).
  551inc_step_2_cont(t_none,l(V*_,_),_,unlimited(V,t_none),Cont,Cont).
  552
  553replace_in_cont([],_,_,[]).
  554replace_in_cont([H1|T1],X,Y,[H2|T2]) :-
  555	(   H1 == X
  556	->  H2 = Y,
  557	    T1 = T2
  558	;   H2 = H1,
  559	    replace_in_cont(T1,X,Y,T2)
  560	).
  561
  562dec_step([],optimum).
  563dec_step([l(V*K,OrdV)|Vs],Status) :-
  564	get_attr(V,clpqr_itf,Att),
  565	arg(2,Att,type(W)),
  566	arg(6,Att,class(Class)),
  567	(   dec_step_2(W,l(V*K,OrdV),Class,Status)
  568	->  true
  569	;   dec_step(Vs,Status)
  570	).
  571
  572dec_step_2(t_U(U),l(V*K,OrdV),Class,Status) :-
  573	K > 0,
  574	(   lb(Class,OrdV,Vub-Vb-_)
  575	->  % found a lower bound
  576	    Status = applied,
  577	    pivot_a(Vub,V,Vb,t_u(U))
  578	;   Status = unlimited(V,t_u(U))
  579	).
  580dec_step_2(t_lU(L,U),l(V*K,OrdV),Class,applied) :-
  581	K > 0,
  582	Init is L - U,
  583	class_basis(Class,Deps),
  584	lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  585	pivot_b(Vub,V,Vb,t_lu(L,U)).
  586dec_step_2(t_L(L),l(V*K,OrdV),Class,Status) :-
  587	K < 0,
  588	(   ub(Class,OrdV,Vub-Vb-_)
  589	->  Status = applied,
  590	    pivot_a(Vub,V,Vb,t_l(L))
  591	;   Status = unlimited(V,t_l(L))
  592	).
  593dec_step_2(t_Lu(L,U),l(V*K,OrdV),Class,applied) :-
  594	K < 0,
  595	Init is U - L,
  596	class_basis(Class,Deps),
  597	ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  598	pivot_b(Vub,V,Vb,t_lu(L,U)).
  599dec_step_2(t_none,l(V*_,_),_,unlimited(V,t_none)).
  600
  601inc_step([],optimum).	% if status has not been set yet: no changes
  602inc_step([l(V*K,OrdV)|Vs],Status) :-
  603	get_attr(V,clpqr_itf,Att),
  604	arg(2,Att,type(W)),
  605	arg(6,Att,class(Class)),
  606	(   inc_step_2(W,l(V*K,OrdV),Class,Status)
  607	->  true
  608	;   inc_step(Vs,Status)
  609	).
  610
  611inc_step_2(t_U(U),l(V*K,OrdV),Class,Status) :-
  612	K < 0,
  613	(   lb(Class,OrdV,Vub-Vb-_)
  614	->  Status = applied,
  615	    pivot_a(Vub,V,Vb,t_u(U))
  616	;   Status = unlimited(V,t_u(U))
  617	).
  618inc_step_2(t_lU(L,U),l(V*K,OrdV),Class,applied) :-
  619	K < 0,
  620	Init is L - U,
  621	class_basis(Class,Deps),
  622	lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  623	pivot_b(Vub,V,Vb,t_lu(L,U)).
  624inc_step_2(t_L(L),l(V*K,OrdV),Class,Status) :-
  625	K > 0,
  626	(   ub(Class,OrdV,Vub-Vb-_)
  627	->  Status = applied,
  628	    pivot_a(Vub,V,Vb,t_l(L))
  629	;   Status = unlimited(V,t_l(L))
  630	).
  631inc_step_2(t_Lu(L,U),l(V*K,OrdV),Class,applied) :-
  632	K > 0,
  633	Init is U - L,
  634	class_basis(Class,Deps),
  635	ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  636	pivot_b(Vub,V,Vb,t_lu(L,U)).
  637inc_step_2(t_none,l(V*_,_),_,unlimited(V,t_none)).
  638
  639% ------------------------- find the most constraining row --------------------
  640%
  641% The code for the lower and the upper bound are dual versions of each other.
  642% The only difference is in the orientation of the comparisons.
  643% Indeps are ruled out by their types.
  644% If there is no bound, this fails.
  645%
  646% *** The actual lb and ub on an indep variable X are [lu]b + b(X), where b(X)
  647% is the value of the active bound.
  648%
  649% Nota bene: We must NOT consider infeasible rows as candidates to
  650%	     leave the basis!
  651%
  652% ub(Class,OrdX,Ub)
  653%
  654% See lb/3: this is similar
  655
  656ub(Class,OrdX,Ub) :-
  657	class_basis(Class,Deps),
  658	ub_first(Deps,OrdX,Ub).
  659
  660% ub_first(Deps,X,Dep-W-Ub)
  661%
  662% Finds the tightest upperbound for variable X from the linear equations of
  663% basis variables Deps, and puts the resulting bound in Ub. Dep is the basis
  664% variable that generates the bound, and W is bound of that variable that has
  665% to be activated to achieve this.
  666
  667ub_first([Dep|Deps],OrdX,Tightest) :-
  668	(   get_attr(Dep,clpqr_itf,Att),
  669	    arg(2,Att,type(Type)),
  670	    arg(4,Att,lin(Lin)),
  671	    ub_inner(Type,OrdX,Lin,W,Ub),
  672	    Ub >= 0
  673	->  ub(Deps,OrdX,Dep-W-Ub,Tightest)
  674	;   ub_first(Deps,OrdX,Tightest)
  675	).
  676
  677% ub(Deps,OrdX,TightestIn,TightestOut)
  678%
  679% See lb/4: this is similar
  680
  681ub([],_,T0,T0).
  682ub([Dep|Deps],OrdX,T0,T1) :-
  683	(   get_attr(Dep,clpqr_itf,Att),
  684	    arg(2,Att,type(Type)),
  685	    arg(4,Att,lin(Lin)),
  686	    ub_inner(Type,OrdX,Lin,W,Ub),
  687	    T0 = _-Ubb,
  688	    Ub < Ubb,
  689	    Ub >= 0
  690	->  ub(Deps,OrdX,Dep-W-Ub,T1)	% tighter bound, use new bound
  691	;   ub(Deps,OrdX,T0,T1)	% no tighter bound, keep current one
  692	).
  693
  694% ub_inner(Type,OrdX,Lin,W,Ub)
  695%
  696% See lb_inner/5: this is similar
  697
  698ub_inner(t_l(L),OrdX,Lin,t_L(L),Ub) :-
  699	nf_rhs_x(Lin,OrdX,Rhs,K),
  700	K < 0,
  701	Ub is (L - Rhs) rdiv K.
  702ub_inner(t_u(U),OrdX,Lin,t_U(U),Ub) :-
  703	nf_rhs_x(Lin,OrdX,Rhs,K),
  704	K > 0,
  705	Ub is (U - Rhs) rdiv K.
  706ub_inner(t_lu(L,U),OrdX,Lin,W,Ub) :-
  707	nf_rhs_x(Lin,OrdX,Rhs,K),
  708	(   K < 0 % use lowerbound
  709	->  W = t_Lu(L,U),
  710	    Ub = (L - Rhs) rdiv K
  711	;   K > 0 % use upperbound
  712	->  W = t_lU(L,U),
  713	    Ub = (U - Rhs) rdiv K
  714	).
  715
  716% lb(Class,OrdX,Lb)
  717%
  718% Returns in Lb how much we can lower the upperbound of X without violating
  719% a bound of the basisvariables.
  720% Lb has the form Dep-W-Lb with Dep the variable whose bound is violated when
  721% lowering the bound for X more, W the actual bound that has to be activated
  722% and Lb the amount that the upperbound can be lowered.
  723% X has ordering OrdX and class Class.
  724
  725lb(Class,OrdX,Lb) :-
  726	class_basis(Class,Deps),
  727	lb_first(Deps,OrdX,Lb).
  728
  729% lb_first(Deps,OrdX,Tightest)
  730%
  731% Returns in Tightest how much we can lower the upperbound of X without
  732% violating a bound of Deps.
  733% Tightest has the form Dep-W-Lb with Dep the variable whose bound is violated
  734% when lowering the bound for X more, W the actual bound that has to be
  735% activated and Lb the amount that the upperbound can be lowered. X has
  736% ordering attribute OrdX.
  737
  738lb_first([Dep|Deps],OrdX,Tightest) :-
  739	(   get_attr(Dep,clpqr_itf,Att),
  740	    arg(2,Att,type(Type)),
  741	    arg(4,Att,lin(Lin)),
  742	    lb_inner(Type,OrdX,Lin,W,Lb),
  743	    Lb =< 0 % Lb > 0 means a violated bound
  744	->  lb(Deps,OrdX,Dep-W-Lb,Tightest)
  745	;   lb_first(Deps,OrdX,Tightest)
  746	).
  747
  748% lb(Deps,OrdX,TightestIn,TightestOut)
  749%
  750% See lb_first/3: this one does the same thing, but is used for the steps after
  751% the first one and remembers the tightest bound so far.
  752
  753lb([],_,T0,T0).
  754lb([Dep|Deps],OrdX,T0,T1) :-
  755	(   get_attr(Dep,clpqr_itf,Att),
  756	    arg(2,Att,type(Type)),
  757	    arg(4,Att,lin(Lin)),
  758	    lb_inner(Type,OrdX,Lin,W,Lb),
  759	    T0 = _-Lbb,
  760	    Lb > Lbb,	% choose the least lowering, others might violate
  761			% bounds
  762	    Lb =< 0	% violation of a bound (without lowering)
  763	->  lb(Deps,OrdX,Dep-W-Lb,T1)
  764	;   lb(Deps,OrdX,T0,T1)
  765	).
  766
  767% lb_inner(Type,X,Lin,W,Lb)
  768%
  769% Returns in Lb how much lower we can make X without violating a bound
  770% by using the linear equation Lin of basis variable B which has type
  771% Type and which has to activate a bound (type W) to do so.
  772%
  773% E.g. when B has a lowerbound L, then L should always be smaller than I + R.
  774% So a lowerbound of X (which has scalar K in Lin), could be at most
  775% (L-(I+R))/K lower than its upperbound (if K is positive).
  776% Also note that Lb should always be smaller than 0, otherwise the row is
  777% not feasible.
  778% X has ordering attribute OrdX.
  779
  780lb_inner(t_l(L),OrdX,Lin,t_L(L),Lb) :-
  781	nf_rhs_x(Lin,OrdX,Rhs,K), % if linear equation Lin contains the term
  782				  % X*K, Rhs is the right hand side of that
  783				  % equation
  784	K > 0,
  785	Lb is (L - Rhs) rdiv K.
  786lb_inner(t_u(U),OrdX,Lin,t_U(U),Lb) :-
  787	nf_rhs_x(Lin,OrdX,Rhs,K),
  788	K < 0, % K < 0
  789	Lb is (U - Rhs) rdiv K.
  790lb_inner(t_lu(L,U),OrdX,Lin,W,Lb) :-
  791	nf_rhs_x(Lin,OrdX,Rhs,K),
  792	(   K < 0
  793	->  W = t_lU(L,U),
  794	    Lb is (U - Rhs) rdiv K
  795	;   K > 0
  796	->  W = t_Lu(L,U),
  797	    Lb is (L - Rhs) rdiv K
  798	).
  799
  800% ---------------------------------- equations --------------------------------
  801%
  802% backsubstitution will not make the system infeasible, if the bounds on the
  803% indep vars are obeyed, but some implied values might pop up in rows where X
  804% occurs
  805%	-) special case X=Y during bs -> get rid of dependend var(s), alias
  806%
  807
  808solve(Lin) :-
  809	Lin = [I,_|H],
  810	solve(H,Lin,I,Bindings,[]),
  811	export_binding(Bindings).
  812
  813% solve(Hom,Lin,I,Bind,BindT)
  814%
  815% Solves a linear equation Lin = [I,_|H] = 0 and exports the generated bindings
  816
  817solve([],_,I,Bind0,Bind0) :-
  818	!,
  819	I =:= 0.
  820solve(H,Lin,_,Bind0,BindT) :-
  821	sd(H,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT),
  822	get_attr(Selected,clpqr_itf,Att),
  823	arg(5,Att,order(Ord)),
  824	isolate(Ord,Lin,Lin1),	% Lin = 0 => Selected = Lin1
  825	(   Category = 1 % classless variable, no bounds
  826	->  setarg(4,Att,lin(Lin1)),
  827	    Lin1 = [Inhom,_|Hom],
  828	    bs_collect_binding(Hom,Selected,Inhom,Bind0,BindT),
  829	    eq_classes(NV,NVT,ClassesUniq)
  830	;   Category = 2 % class variable, no bounds
  831	->  arg(6,Att,class(NewC)),
  832	    class_allvars(NewC,Deps),
  833	    (   ClassesUniq = [_] % rank increasing
  834	    ->	bs_collect_bindings(Deps,Ord,Lin1,Bind0,BindT)
  835	    ;   Bind0 = BindT,
  836		bs(Deps,Ord,Lin1)
  837	    ),
  838	    eq_classes(NV,NVT,ClassesUniq)
  839	;   Category = 3 % classless variable, all variables in Lin and
  840			 % Selected are bounded
  841	->  arg(2,Att,type(Type)),
  842	    setarg(4,Att,lin(Lin1)),
  843	    deactivate_bound(Type,Selected),
  844	    eq_classes(NV,NVT,ClassesUniq),
  845	    basis_add(Selected,Basis),
  846	    undet_active(Lin1),	% we can't tell which bound will likely be a
  847				% problem at this point
  848	    Lin1 = [Inhom,_|Hom],
  849	    bs_collect_binding(Hom,Selected,Inhom,Bind0,Bind1),	% only if
  850								% Hom = []
  851	    rcbl(Basis,Bind1,BindT) % reconsider entire basis
  852	;   Category = 4 % class variable, all variables in Lin and Selected
  853			 % are bounded
  854	->  arg(2,Att,type(Type)),
  855	    arg(6,Att,class(NewC)),
  856	    class_allvars(NewC,Deps),
  857	    (   ClassesUniq = [_] % rank increasing
  858	    ->	bs_collect_bindings(Deps,Ord,Lin1,Bind0,Bind1)
  859	    ;   Bind0 = Bind1,
  860		bs(Deps,Ord,Lin1)
  861	    ),
  862	    deactivate_bound(Type,Selected),
  863	    basis_add(Selected,Basis),
  864	    % eq_classes( NV, NVT, ClassesUniq),
  865	    %  4 -> var(NV)
  866	    equate(ClassesUniq,_),
  867	    undet_active(Lin1),
  868	    rcbl(Basis,Bind1,BindT)
  869	).
  870
  871%
  872% Much like solve, but we solve for a particular variable of type t_none
  873%
  874
  875% solve_x(H,Lin,I,X,[Bind|BindT],BindT)
  876%
  877%
  878
  879solve_x(Lin,X) :-
  880	Lin = [I,_|H],
  881	solve_x(H,Lin,I,X,Bindings,[]),
  882	export_binding(Bindings).
  883
  884solve_x([],_,I,_,Bind0,Bind0) :-
  885	!,
  886	I =:= 0.
  887solve_x(H,Lin,_,X,Bind0,BindT) :-
  888	sd(H,[],ClassesUniq,9-9-0,_,NV,NVT),
  889	get_attr(X,clpqr_itf,Att),
  890	arg(5,Att,order(OrdX)),
  891	isolate(OrdX,Lin,Lin1),
  892	(   arg(6,Att,class(NewC))
  893	->  class_allvars(NewC,Deps),
  894	    (   ClassesUniq = [_] % rank increasing
  895	    ->	bs_collect_bindings(Deps,OrdX,Lin1,Bind0,BindT)
  896	    ;   Bind0 = BindT,
  897		bs(Deps,OrdX,Lin1)
  898	    ),
  899	    eq_classes(NV,NVT,ClassesUniq)
  900	;   setarg(4,Att,lin(Lin1)),
  901	    Lin1 = [Inhom,_|Hom],
  902	    bs_collect_binding(Hom,X,Inhom,Bind0,BindT),
  903	    eq_classes(NV,NVT,ClassesUniq)
  904	).
  905
  906% solve_ord_x(Lin,OrdX,ClassX)
  907%
  908% Does the same thing as solve_x/2, but has the ordering of X and its class as
  909% input. This also means that X has a class which is not sure in solve_x/2.
  910
  911solve_ord_x(Lin,OrdX,ClassX) :-
  912	Lin = [I,_|H],
  913	solve_ord_x(H,Lin,I,OrdX,ClassX,Bindings,[]),
  914	export_binding(Bindings).
  915
  916solve_ord_x([],_,I,_,_,Bind0,Bind0) :-
  917	I =:= 0.
  918solve_ord_x([_|_],Lin,_,OrdX,ClassX,Bind0,BindT) :-
  919	isolate(OrdX,Lin,Lin1),
  920	Lin1 = [_,_|H1],
  921	sd(H1,[],ClassesUniq1,9-9-0,_,NV,NVT), % do sd on Lin without X, then
  922					       % add class of X
  923	ord_add_element(ClassesUniq1,ClassX,ClassesUniq),
  924	class_allvars(ClassX,Deps),
  925	(   ClassesUniq = [_] % rank increasing
  926	->  bs_collect_bindings(Deps,OrdX,Lin1,Bind0,BindT)
  927	;   Bind0 = BindT,
  928	    bs(Deps,OrdX,Lin1)
  929	),
  930	eq_classes(NV,NVT,ClassesUniq).
  931
  932% sd(H,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT)
  933
  934% sd(Hom,ClassesIn,ClassesOut,PreferenceIn,PreferenceOut,[NV|NVTail],NVTail)
  935%
  936% ClassesOut is a sorted list of the different classes that are either in
  937% ClassesIn or that are the classes of the variables in Hom. Variables that do
  938% not belong to a class yet, are put in the difference list NV.
  939
  940sd([],Class0,Class0,Preference0,Preference0,NV0,NV0).
  941sd([l(X*K,_)|Xs],Class0,ClassN,Preference0,PreferenceN,NV0,NVt) :-
  942	get_attr(X,clpqr_itf,Att),
  943	(   arg(6,Att,class(Xc)) % old: has class
  944	->  NV0 = NV1,
  945	    ord_add_element(Class0,Xc,Class1),
  946	    (   arg(2,Att,type(t_none))
  947	    ->  preference(Preference0,2-X-K,Preference1)
  948		    % has class, no bounds => category 2
  949	    ;   preference(Preference0,4-X-K,Preference1)
  950		    % has class, is bounded => category 4
  951	    )
  952	;   % new: has no class
  953	    Class1 = Class0,
  954	    NV0 = [X|NV1], % X has no class yet, add to list of new variables
  955	    (   arg(2,Att,type(t_none))
  956	    ->  preference(Preference0,1-X-K,Preference1)
  957		    % no class, no bounds => category 1
  958	    ;   preference(Preference0,3-X-K,Preference1)
  959		    % no class, is bounded => category 3
  960	    )
  961	),
  962	sd(Xs,Class1,ClassN,Preference1,PreferenceN,NV1,NVt).
  963
  964%
  965% A is best sofar, B is current
  966% smallest prefered
  967preference(A,B,Pref) :-
  968	A = Px-_-_,
  969	B = Py-_-_,
  970	(   Px < Py
  971	->  Pref = A
  972	;   Pref = B
  973	).
  974
  975% eq_classes(NV,NVTail,Cs)
  976%
  977% Attaches all classless variables NV to a new class and equates all other
  978% classes with this class. The equate operation only happens after attach_class
  979% because the unification of classes can bind the tail of the AllVars attribute
  980% to a nonvar and then the attach_class operation wouldn't work.
  981
  982eq_classes(NV,_,Cs) :-
  983	var(NV),
  984	!,
  985	equate(Cs,_).
  986eq_classes(NV,NVT,Cs) :-
  987	class_new(Su,clpq,NV,NVT,[]), % make a new class Su with NV as the variables
  988	attach_class(NV,Su), % attach the variables NV to Su
  989	equate(Cs,Su).
  990
  991equate([],_).
  992equate([X|Xs],X) :- equate(Xs,X).
  993
  994%
  995% assert: none of the Vars has a class attribute yet
  996%
  997attach_class(Xs,_) :-
  998	var(Xs), % Tail
  999	!.
 1000attach_class([X|Xs],Class) :-
 1001	get_attr(X,clpqr_itf,Att),
 1002	setarg(6,Att,class(Class)),
 1003	attach_class(Xs,Class).
 1004
 1005% unconstrained(Lin,Uc,Kuc,Rest)
 1006%
 1007% Finds an unconstrained variable Uc (type(t_none)) in Lin with scalar Kuc and
 1008% removes it from Lin to return Rest.
 1009
 1010unconstrained(Lin,Uc,Kuc,Rest) :-
 1011	Lin = [_,_|H],
 1012	sd(H,[],_,9-9-0,Category-Uc-_,_,_),
 1013	Category =< 2,
 1014	get_attr(Uc,clpqr_itf,Att),
 1015	arg(5,Att,order(OrdUc)),
 1016	delete_factor(OrdUc,Lin,Rest,Kuc).
 1017
 1018%
 1019% point the vars in Lin into the same equivalence class
 1020% maybe join some global data
 1021%
 1022same_class([],_).
 1023same_class([l(X*_,_)|Xs],Class) :-
 1024	get_or_add_class(X,Class),
 1025	same_class(Xs,Class).
 1026
 1027% get_or_add_class(X,Class)
 1028%
 1029% Returns in Class the class of X if X has one, or a new class where X now
 1030% belongs to if X didn't have one.
 1031
 1032get_or_add_class(X,Class) :-
 1033	get_attr(X,clpqr_itf,Att),
 1034	arg(1,Att,CLP),
 1035	(   arg(6,Att,class(ClassX))
 1036	->  ClassX = Class
 1037	;   setarg(6,Att,class(Class)),
 1038	    class_new(Class,CLP,[X|Tail],Tail,[])
 1039	).
 1040
 1041% allvars(X,Allvars)
 1042%
 1043% Allvars is a list of all variables in the class to which X belongs.
 1044
 1045allvars(X,Allvars) :-
 1046	get_attr(X,clpqr_itf,Att),
 1047	arg(6,Att,class(C)),
 1048	class_allvars(C,Allvars).
 1049
 1050% deactivate_bound(Type,Variable)
 1051%
 1052% The Type of the variable is changed to reflect the deactivation of its
 1053% bounds.
 1054% t_L(_) becomes t_l(_), t_lU(_,_) becomes t_lu(_,_) and so on.
 1055
 1056deactivate_bound(t_l(_),_).
 1057deactivate_bound(t_u(_),_).
 1058deactivate_bound(t_lu(_,_),_).
 1059deactivate_bound(t_L(L),X) :-
 1060	get_attr(X,clpqr_itf,Att),
 1061	setarg(2,Att,type(t_l(L))).
 1062deactivate_bound(t_Lu(L,U),X) :-
 1063	get_attr(X,clpqr_itf,Att),
 1064	setarg(2,Att,type(t_lu(L,U))).
 1065deactivate_bound(t_U(U),X) :-
 1066	get_attr(X,clpqr_itf,Att),
 1067	setarg(2,Att,type(t_u(U))).
 1068deactivate_bound(t_lU(L,U),X) :-
 1069	get_attr(X,clpqr_itf,Att),
 1070	setarg(2,Att,type(t_lu(L,U))).
 1071
 1072% intro_at(X,Value,Type)
 1073%
 1074% Variable X gets new type Type which reflects the activation of a bound with
 1075% value Value. In the linear equations of all the variables belonging to the
 1076% same class as X, X is substituted by [0,Value,X] to reflect the new active
 1077% bound.
 1078
 1079intro_at(X,Value,Type) :-
 1080	get_attr(X,clpqr_itf,Att),
 1081	arg(5,Att,order(Ord)),
 1082	arg(6,Att,class(Class)),
 1083	setarg(2,Att,type(Type)),
 1084	(   Value =:= 0
 1085	->  true
 1086	;   backsubst_delta(Class,Ord,X,Value)
 1087	).
 1088
 1089% undet_active(Lin)
 1090%
 1091% For each variable in the homogene part of Lin, a bound is activated
 1092% if an inactive bound exists. (t_l(L) becomes t_L(L) and so on)
 1093
 1094undet_active([_,_|H]) :-
 1095	undet_active_h(H).
 1096
 1097% undet_active_h(Hom)
 1098%
 1099% For each variable in homogene part Hom, a bound is activated if an
 1100% inactive bound exists (t_l(L) becomes t_L(L) and so on)
 1101
 1102undet_active_h([]).
 1103undet_active_h([l(X*_,_)|Xs]) :-
 1104	get_attr(X,clpqr_itf,Att),
 1105	arg(2,Att,type(Type)),
 1106	undet_active(Type,X),
 1107	undet_active_h(Xs).
 1108
 1109% undet_active(Type,Var)
 1110%
 1111% An inactive bound of Var is activated if such exists
 1112% t_lu(L,U) is arbitrarily chosen to become t_Lu(L,U)
 1113
 1114undet_active(t_none,_).	% type_activity
 1115undet_active(t_L(_),_).
 1116undet_active(t_Lu(_,_),_).
 1117undet_active(t_U(_),_).
 1118undet_active(t_lU(_,_),_).
 1119undet_active(t_l(L),X) :- intro_at(X,L,t_L(L)).
 1120undet_active(t_u(U),X) :- intro_at(X,U,t_U(U)).
 1121undet_active(t_lu(L,U),X) :- intro_at(X,L,t_Lu(L,U)).
 1122
 1123% determine_active_dec(Lin)
 1124%
 1125% Activates inactive bounds on the variables of Lin if such bounds exist.
 1126% If the type of a variable is t_none, this fails. This version is aimed
 1127% to make the R component of Lin as small as possible in order not to violate
 1128% an upperbound (see reconsider/1)
 1129
 1130determine_active_dec([_,_|H]) :-
 1131	determine_active(H,-1).
 1132
 1133% determine_active_inc(Lin)
 1134%
 1135% Activates inactive bounds on the variables of Lin if such bounds exist.
 1136% If the type of a variable is t_none, this fails. This version is aimed
 1137% to make the R component of Lin as large as possible in order not to violate
 1138% a lowerbound (see reconsider/1)
 1139
 1140determine_active_inc([_,_|H]) :-
 1141	determine_active(H,1).
 1142
 1143% determine_active(Hom,S)
 1144%
 1145% For each variable in Hom, activates its bound if it is not yet activated.
 1146% For the case of t_lu(_,_) the lower or upper bound is activated depending on
 1147% K and S:
 1148% If sign of K*S is negative, then lowerbound, otherwise upperbound.
 1149
 1150determine_active([],_).
 1151determine_active([l(X*K,_)|Xs],S) :-
 1152	get_attr(X,clpqr_itf,Att),
 1153	arg(2,Att,type(Type)),
 1154	determine_active(Type,X,K,S),
 1155	determine_active(Xs,S).
 1156
 1157determine_active(t_L(_),_,_,_).
 1158determine_active(t_Lu(_,_),_,_,_).
 1159determine_active(t_U(_),_,_,_).
 1160determine_active(t_lU(_,_),_,_,_).
 1161determine_active(t_l(L),X,_,_) :- intro_at(X,L,t_L(L)).
 1162determine_active(t_u(U),X,_,_) :- intro_at(X,U,t_U(U)).
 1163determine_active(t_lu(L,U),X,K,S) :-
 1164	KS is K*S,
 1165	(   KS < 0
 1166	->  intro_at(X,L,t_Lu(L,U))
 1167	;   KS > 0
 1168	->  intro_at(X,U,t_lU(L,U))
 1169	).
 1170
 1171%
 1172% Careful when an indep turns into t_none !!!
 1173%
 1174
 1175detach_bounds(V) :-
 1176	get_attr(V,clpqr_itf,Att),
 1177	arg(2,Att,type(Type)),
 1178	arg(4,Att,lin(Lin)),
 1179	arg(5,Att,order(OrdV)),
 1180	arg(6,Att,class(Class)),
 1181	setarg(2,Att,type(t_none)),
 1182	setarg(3,Att,strictness(0)),
 1183	(   indep(Lin,OrdV)
 1184	->  (   ub(Class,OrdV,Vub-Vb-_)
 1185	    ->	% exchange against thightest
 1186		class_basis_drop(Class,Vub),
 1187		pivot(Vub,Class,OrdV,Vb,Type)
 1188	    ;   lb(Class,OrdV,Vlb-Vb-_)
 1189	    ->  class_basis_drop(Class,Vlb),
 1190		pivot(Vlb,Class,OrdV,Vb,Type)
 1191	    ;   true
 1192	    )
 1193	;   class_basis_drop(Class,V)
 1194	).
 1195
 1196detach_bounds_vlv(OrdV,Lin,Class,Var,NewLin) :-
 1197	(   indep(Lin,OrdV)
 1198	->  Lin = [_,R|_],
 1199	    (   ub(Class,OrdV,Vub-Vb-_)
 1200	    ->  % in verify_lin, class might contain two occurrences of Var,
 1201		% but it doesn't matter which one we delete
 1202		class_basis_drop(Class,Var),
 1203		pivot_vlv(Vub,Class,OrdV,Vb,R,NewLin)
 1204	    ;   lb(Class,OrdV,Vlb-Vb-_)
 1205	    ->  class_basis_drop(Class,Var),
 1206		pivot_vlv(Vlb,Class,OrdV,Vb,R,NewLin)
 1207	    ;   NewLin = Lin
 1208	    )
 1209	;   NewLin = Lin,
 1210	    class_basis_drop(Class,Var)
 1211	).
 1212
 1213% ----------------------------- manipulate the basis --------------------------
 1214
 1215% basis_drop(X)
 1216%
 1217% Removes X from the basis of the class to which X belongs.
 1218
 1219basis_drop(X) :-
 1220	get_attr(X,clpqr_itf,Att),
 1221	arg(6,Att,class(Cv)),
 1222	class_basis_drop(Cv,X).
 1223
 1224% basis(X,Basis)
 1225%
 1226% Basis is the basis of the class to which X belongs.
 1227
 1228basis(X,Basis) :-
 1229	get_attr(X,clpqr_itf,Att),
 1230	arg(6,Att,class(Cv)),
 1231	class_basis(Cv,Basis).
 1232
 1233% basis_add(X,NewBasis)
 1234%
 1235% NewBasis is the result of adding X to the basis of the class to which X
 1236% belongs.
 1237
 1238basis_add(X,NewBasis) :-
 1239	get_attr(X,clpqr_itf,Att),
 1240	arg(6,Att,class(Cv)),
 1241	class_basis_add(Cv,X,NewBasis).
 1242
 1243% basis_pivot(Leave,Enter)
 1244%
 1245% Removes Leave from the basis of the class to which it belongs, and adds
 1246% Enter to that basis.
 1247
 1248basis_pivot(Leave,Enter) :-
 1249	get_attr(Leave,clpqr_itf,Att),
 1250	arg(6,Att,class(Cv)),
 1251	class_basis_pivot(Cv,Enter,Leave).
 1252
 1253% ----------------------------------- pivot -----------------------------------
 1254
 1255% pivot(Dep,Indep)
 1256%
 1257% The linear equation of variable Dep, is transformed into one of variable
 1258% Indep, containing Dep. Then, all occurrences of Indep in linear equations are
 1259% substituted by this new definition
 1260
 1261%
 1262% Pivot ignoring rhs and active states
 1263%
 1264
 1265pivot(Dep,Indep) :-
 1266	get_attr(Dep,clpqr_itf,AttD),
 1267	arg(4,AttD,lin(H)),
 1268	arg(5,AttD,order(OrdDep)),
 1269	get_attr(Indep,clpqr_itf,AttI),
 1270	arg(5,AttI,order(Ord)),
 1271	arg(5,AttI,class(Class)),
 1272	delete_factor(Ord,H,H0,Coeff),
 1273	K is -1 rdiv Coeff,
 1274	add_linear_ff(H0,K,[0,0,l(Dep* -1,OrdDep)],K,Lin),
 1275	backsubst(Class,Ord,Lin).
 1276
 1277% pivot_a(Dep,Indep,IndepT,DepT)
 1278%
 1279% Removes Dep from the basis, puts Indep in, and pivots the equation of
 1280% Dep to become one of Indep. The type of Dep becomes DepT (which means
 1281% it gets deactivated), the type of Indep becomes IndepT (which means it
 1282% gets activated)
 1283
 1284
 1285pivot_a(Dep,Indep,Vb,Wd) :-
 1286	basis_pivot(Dep,Indep),
 1287	get_attr(Indep,clpqr_itf,Att),
 1288	arg(2,Att,type(Type)),
 1289	arg(5,Att,order(Ord)),
 1290	arg(6,Att,class(Class)),
 1291	pivot(Dep,Class,Ord,Vb,Type),
 1292	get_attr(Indep,clpqr_itf,Att2), %changed?
 1293	setarg(2,Att2,type(Wd)).
 1294
 1295pivot_b(Vub,V,Vb,Wd) :-
 1296	(   Vub == V
 1297	->  get_attr(V,clpqr_itf,Att),
 1298	    arg(5,Att,order(Ord)),
 1299	    arg(6,Att,class(Class)),
 1300	    setarg(2,Att,type(Vb)),
 1301	    pivot_b_delta(Vb,Delta), % nonzero(Delta)
 1302	    backsubst_delta(Class,Ord,V,Delta)
 1303	;   pivot_a(Vub,V,Vb,Wd)
 1304	).
 1305
 1306pivot_b_delta(t_Lu(L,U),Delta) :- Delta is L-U.
 1307pivot_b_delta(t_lU(L,U),Delta) :- Delta is U-L.
 1308
 1309% select_active_bound(Type,Bound)
 1310%
 1311% Returns the bound that is active in Type (if such exists, 0 otherwise)
 1312
 1313select_active_bound(t_L(L),L).
 1314select_active_bound(t_Lu(L,_),L).
 1315select_active_bound(t_U(U),U).
 1316select_active_bound(t_lU(_,U),U).
 1317select_active_bound(t_none,0).
 1318%
 1319% for project.pl
 1320%
 1321select_active_bound(t_l(_),0).
 1322select_active_bound(t_u(_),0).
 1323select_active_bound(t_lu(_,_),0).
 1324
 1325
 1326% pivot(Dep,Class,IndepOrd,DepAct,IndAct)
 1327%
 1328% See pivot/2.
 1329% In addition, variable Indep with ordering IndepOrd has an active bound IndAct
 1330
 1331%
 1332%
 1333% Pivot taking care of rhs and active states
 1334%
 1335pivot(Dep,Class,IndepOrd,DepAct,IndAct) :-
 1336	get_attr(Dep,clpqr_itf,Att),
 1337	arg(4,Att,lin(H)),
 1338	arg(5,Att,order(DepOrd)),
 1339	setarg(2,Att,type(DepAct)),
 1340	select_active_bound(DepAct,AbvD), % New current value for Dep
 1341	select_active_bound(IndAct,AbvI), % Old current value of Indep
 1342	delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ...
 1343	AbvDm is -AbvD,
 1344	AbvIm is -AbvI,
 1345	add_linear_f1([0,AbvIm],Coeff,H0,H1),
 1346	K is -1 rdiv Coeff,
 1347	add_linear_ff(H1,K,[0,AbvDm,l(Dep* -1,DepOrd)],K,H2),
 1348	    % Indep = -1/Coeff*... + 1/Coeff*Dep
 1349	add_linear_11(H2,[0,AbvIm],Lin),
 1350	backsubst(Class,IndepOrd,Lin).
 1351
 1352% Rewrite Dep = ... + Coeff*Indep + ...
 1353% into Indep = ... + -1/Coeff*Dep + ...
 1354%
 1355% For backsubstitution, old current value of Indep must be removed from RHS
 1356% New current value of Dep must be added to RHS
 1357% For solving: old current value of Indep should be out of RHS
 1358
 1359pivot_vlv(Dep,Class,IndepOrd,DepAct,AbvI,Lin) :-
 1360	get_attr(Dep,clpqr_itf,Att),
 1361	arg(4,Att,lin(H)),
 1362	arg(5,Att,order(DepOrd)),
 1363	setarg(2,Att,type(DepAct)),
 1364	select_active_bound(DepAct,AbvD), % New current value for Dep
 1365	delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ...
 1366	AbvDm is -AbvD,
 1367	AbvIm is -AbvI,
 1368	add_linear_f1([0,AbvIm],Coeff,H0,H1),
 1369	K is -1 rdiv Coeff,
 1370	add_linear_ff(H1,K,[0,AbvDm,l(Dep* -1,DepOrd)],K,Lin),
 1371	    % Indep = -1/Coeff*... + 1/Coeff*Dep
 1372	add_linear_11(Lin,[0,AbvIm],SubstLin),
 1373	backsubst(Class,IndepOrd,SubstLin).
 1374
 1375% backsubst_delta(Class,OrdX,X,Delta)
 1376%
 1377% X with ordering attribute OrdX, is substituted in all linear equations of
 1378% variables in the class Class, by linear equation [0,Delta,l(X*1,OrdX)]. This
 1379% reflects the activation of a bound.
 1380
 1381backsubst_delta(Class,OrdX,X,Delta) :-
 1382	backsubst(Class,OrdX,[0,Delta,l(X*1,OrdX)]).
 1383
 1384% backsubst(Class,OrdX,Lin)
 1385%
 1386% X with ordering OrdX is substituted in all linear equations of variables in
 1387% the class Class, by linear equation Lin
 1388
 1389backsubst(Class,OrdX,Lin) :-
 1390	class_allvars(Class,Allvars),
 1391	bs(Allvars,OrdX,Lin).
 1392
 1393% bs(Vars,OrdV,Lin)
 1394%
 1395% In all linear equations of the variables Vars, variable V with ordering
 1396% attribute OrdV is substituted by linear equation Lin.
 1397%
 1398% valid if nothing will go ground
 1399%
 1400
 1401bs(Xs,_,_) :-
 1402	var(Xs),
 1403	!.
 1404bs([X|Xs],OrdV,Lin) :-
 1405	(   get_attr(X,clpqr_itf,Att),
 1406	    arg(4,Att,lin(LinX)),
 1407	    nf_substitute(OrdV,Lin,LinX,LinX1) % does not change attributes
 1408	->  setarg(4,Att,lin(LinX1)),
 1409	    bs(Xs,OrdV,Lin)
 1410	;   bs(Xs,OrdV,Lin)
 1411	).
 1412
 1413%
 1414% rank increasing backsubstitution
 1415%
 1416
 1417% bs_collect_bindings(Deps,SelectedOrd,Lin,Bind,BindT)
 1418%
 1419% Collects bindings (of the form [X-I] where X = I is the binding) by
 1420% substituting Selected in all linear equations of the variables Deps (which
 1421% are of the same class), by Lin. Selected has ordering attribute SelectedOrd.
 1422%
 1423% E.g. when V = 2X + 3Y + 4, X = 3V + 2Z and Y = 4X + 3
 1424% we can substitute V in the linear equation of X: X = 6X + 9Y + 2Z + 12
 1425% we can't substitute V in the linear equation of Y of course.
 1426
 1427bs_collect_bindings(Xs,_,_,Bind0,BindT) :-
 1428	var(Xs),
 1429	!,
 1430	Bind0 = BindT.
 1431bs_collect_bindings([X|Xs],OrdV,Lin,Bind0,BindT) :-
 1432	(   get_attr(X,clpqr_itf,Att),
 1433	    arg(4,Att,lin(LinX)),
 1434	    nf_substitute(OrdV,Lin,LinX,LinX1) % does not change attributes
 1435	->  setarg(4,Att,lin(LinX1)),
 1436	    LinX1 = [Inhom,_|Hom],
 1437	    bs_collect_binding(Hom,X,Inhom,Bind0,Bind1),
 1438	    bs_collect_bindings(Xs,OrdV,Lin,Bind1,BindT)
 1439	;   bs_collect_bindings(Xs,OrdV,Lin,Bind0,BindT)
 1440	).
 1441
 1442% bs_collect_binding(Hom,Selected,Inhom,Bind,BindT)
 1443%
 1444% Collects binding following from Selected = Hom + Inhom.
 1445% If Hom = [], returns the binding Selected-Inhom (=0)
 1446%
 1447bs_collect_binding([],X,Inhom) --> [X-Inhom].
 1448bs_collect_binding([_|_],_,_) --> [].
 1449
 1450%
 1451% reconsider the basis
 1452%
 1453
 1454% rcbl(Basis,Bind,BindT)
 1455%
 1456%
 1457
 1458rcbl([],Bind0,Bind0).
 1459rcbl([X|Continuation],Bind0,BindT) :-
 1460	(   rcb_cont(X,Status,Violated,Continuation,NewContinuation) % have a culprit
 1461	->  rcbl_status(Status,X,NewContinuation,Bind0,BindT,Violated)
 1462	;   rcbl(Continuation,Bind0,BindT)
 1463	).
 1464
 1465rcb_cont(X,Status,Violated,ContIn,ContOut) :-
 1466	get_attr(X,clpqr_itf,Att),
 1467	arg(2,Att,type(Type)),
 1468	arg(4,Att,lin([I,R|H])),
 1469	(   Type = t_l(L) % case 1: lowerbound: R + I should always be larger
 1470			  % than the lowerbound
 1471	->  R + I =< L,
 1472	    Violated = l(L),
 1473	    inc_step_cont(H,Status,ContIn,ContOut)
 1474	;   Type = t_u(U) % case 2: upperbound: R + I should always be smaller
 1475			  % than the upperbound
 1476	->  R + I >= U,
 1477	    Violated = u(U),
 1478	    dec_step_cont(H,Status,ContIn,ContOut)
 1479	;   Type = t_lu(L,U) % case 3: check both
 1480	->  At is R + I,
 1481	    (   At =< L
 1482	    ->	Violated = l(L),
 1483		inc_step_cont(H,Status,ContIn,ContOut)
 1484	    ;   At >= U
 1485	    ->	Violated = u(U),
 1486		dec_step_cont(H,Status,ContIn,ContOut)
 1487	    )
 1488	). % other types imply nonbasic variable or unbounded variable
 1489
 1490
 1491
 1492%
 1493% reconsider one element of the basis
 1494% later: lift the binds
 1495%
 1496reconsider(X) :-
 1497	rcb(X,Status,Violated),
 1498	!,
 1499	rcbl_status(Status,X,[],Binds,[],Violated),
 1500	export_binding(Binds).
 1501reconsider(_).
 1502
 1503%
 1504% Find a basis variable out of its bound or at its bound
 1505% Try to move it into whithin its bound
 1506%   a) impossible -> fail
 1507%   b) optimum at the bound -> implied value
 1508%   c) else look at the remaining basis variables
 1509%
 1510%
 1511% Idea: consider a variable V with linear equation Lin.
 1512% When a bound on a variable X of Lin gets activated, its value, multiplied
 1513% with the scalar of X, is added to the R component of Lin.
 1514% When we consider the lowerbound of V, it must be smaller than R + I, since R
 1515% contains at best the lowerbounds of the variables in Lin (but could contain
 1516% upperbounds, which are of course larger). So checking this can show the
 1517% violation of a bound of V. A similar case works for the upperbound.
 1518
 1519rcb(X,Status,Violated) :-
 1520	get_attr(X,clpqr_itf,Att),
 1521	arg(2,Att,type(Type)),
 1522	arg(4,Att,lin([I,R|H])),
 1523	(   Type = t_l(L) % case 1: lowerbound: R + I should always be larger
 1524			  % than the lowerbound
 1525	->  R + I =< L,
 1526	    Violated = l(L),
 1527	    inc_step(H,Status)
 1528	;   Type = t_u(U) % case 2: upperbound: R + I should always be smaller
 1529			  % than the upperbound
 1530	->  R + I >= U,
 1531	    Violated = u(U),
 1532	    dec_step(H,Status)
 1533	;   Type = t_lu(L,U) % case 3: check both
 1534	->  At is R + I,
 1535	    (   At =< L
 1536	    ->	Violated = l(L),
 1537		inc_step(H,Status)
 1538	    ;   At >= U
 1539	    ->	Violated = u(U),
 1540		dec_step(H,Status)
 1541	    )
 1542	). % other types imply nonbasic variable or unbounded variable
 1543
 1544% rcbl_status(Status,X,Continuation,[Bind|BindT],BindT,Violated)
 1545%
 1546%
 1547
 1548rcbl_status(optimum,X,Cont,B0,Bt,Violated) :- rcbl_opt(Violated,X,Cont,B0,Bt).
 1549rcbl_status(applied,X,Cont,B0,Bt,Violated) :- rcbl_app(Violated,X,Cont,B0,Bt).
 1550rcbl_status(unlimited(Indep,DepT),X,Cont,B0,Bt,Violated) :-
 1551	rcbl_unl(Violated,X,Cont,B0,Bt,Indep,DepT).
 1552
 1553%
 1554% Might reach optimum immediately without changing the basis,
 1555% but in general we must assume that there were pivots.
 1556% If the optimum meets the bound, we backsubstitute the implied
 1557% value, solve will call us again to check for further implied
 1558% values or unsatisfiability in the rank increased system.
 1559%
 1560rcbl_opt(l(L),X,Continuation,B0,B1) :-
 1561	get_attr(X,clpqr_itf,Att),
 1562	arg(2,Att,type(Type)),
 1563	arg(3,Att,strictness(Strict)),
 1564	arg(4,Att,lin(Lin)),
 1565	Lin = [I,R|_],
 1566	Opt is R + I,
 1567	(   L < Opt
 1568	->  narrow_u(Type,X,Opt), % { X =< Opt }
 1569	    rcbl(Continuation,B0,B1)
 1570	;   L =:= Opt,
 1571	    Strict /\ 2 =:= 0, % meets lower
 1572	    Mop is -Opt,
 1573	    normalize_scalar(Mop,MopN),
 1574	    add_linear_11(MopN,Lin,Lin1),
 1575	    Lin1 = [Inhom,_|Hom],
 1576	    (   Hom = []
 1577	    ->  rcbl(Continuation,B0,B1) % would not callback
 1578	    ;   solve(Hom,Lin1,Inhom,B0,B1)
 1579	    )
 1580	).
 1581rcbl_opt(u(U),X,Continuation,B0,B1) :-
 1582	get_attr(X,clpqr_itf,Att),
 1583	arg(2,Att,type(Type)),
 1584	arg(3,Att,strictness(Strict)),
 1585	arg(4,Att,lin(Lin)),
 1586	Lin = [I,R|_],
 1587	Opt is R + I,
 1588	(   U > Opt
 1589	->  narrow_l(Type,X,Opt), % { X >= Opt }
 1590	    rcbl(Continuation,B0,B1)
 1591	;   U =:= Opt,
 1592	    Strict /\ 1 =:= 0, % meets upper
 1593	    Mop is -Opt,
 1594	    normalize_scalar(Mop,MopN),
 1595	    add_linear_11(MopN,Lin,Lin1),
 1596	    Lin1 = [Inhom,_|Hom],
 1597	    (   Hom = []
 1598	    ->  rcbl(Continuation,B0,B1) % would not callback
 1599	    ;   solve(Hom,Lin1,Inhom,B0,B1)
 1600	    )
 1601	).
 1602
 1603%
 1604% Basis has already changed when this is called
 1605%
 1606rcbl_app(l(L),X,Continuation,B0,B1) :-
 1607	get_attr(X,clpqr_itf,Att),
 1608	arg(4,Att,lin([I,R|H])),
 1609	(   R + I > L % within bound now
 1610	->  rcbl(Continuation,B0,B1)
 1611	;   inc_step(H,Status),
 1612	    rcbl_status(Status,X,Continuation,B0,B1,l(L))
 1613	).
 1614rcbl_app(u(U),X,Continuation,B0,B1) :-
 1615	get_attr(X,clpqr_itf,Att),
 1616	arg(4,Att,lin([I,R|H])),
 1617	(   R + I < U % within bound now
 1618	->  rcbl(Continuation,B0,B1)
 1619	;   dec_step(H,Status),
 1620	    rcbl_status(Status,X,Continuation,B0,B1,u(U))
 1621	).
 1622%
 1623% This is never called for a t_lu culprit
 1624%
 1625rcbl_unl(l(L),X,Continuation,B0,B1,Indep,DepT) :-
 1626	pivot_a(X,Indep,t_L(L),DepT), % changes the basis
 1627	rcbl(Continuation,B0,B1).
 1628rcbl_unl(u(U),X,Continuation,B0,B1,Indep,DepT) :-
 1629	pivot_a(X,Indep,t_U(U),DepT), % changes the basis
 1630	rcbl(Continuation,B0,B1).
 1631
 1632% narrow_u(Type,X,U)
 1633%
 1634% Narrows down the upperbound of X (type Type) to U.
 1635% Fails if Type is not t_u(_) or t_lu(_)
 1636
 1637narrow_u(t_u(_),X,U) :-
 1638	get_attr(X,clpqr_itf,Att),
 1639	setarg(2,Att,type(t_u(U))).
 1640narrow_u(t_lu(L,_),X,U) :-
 1641	get_attr(X,clpqr_itf,Att),
 1642	setarg(2,Att,type(t_lu(L,U))).
 1643
 1644% narrow_l(Type,X,L)
 1645%
 1646% Narrows down the lowerbound of X (type Type) to L.
 1647% Fails if Type is not t_l(_) or t_lu(_)
 1648
 1649narrow_l( t_l(_),    X, L) :-
 1650	get_attr(X,clpqr_itf,Att),
 1651	setarg(2,Att,type(t_l(L))).
 1652
 1653narrow_l( t_lu(_,U), X, L) :-
 1654	get_attr(X,clpqr_itf,Att),
 1655	setarg(2,Att,type(t_lu(L,U))).
 1656
 1657% ----------------------------------- dump ------------------------------------
 1658
 1659% dump_var(Type,Var,I,H,Dump,DumpTail)
 1660%
 1661% Returns in Dump a representation of the linear constraint on variable
 1662% Var which has linear equation H + I and has type Type.
 1663
 1664dump_var(t_none,V,I,H) -->
 1665	!,
 1666	(   {
 1667		H = [l(W*K,_)],
 1668		V == W,
 1669		I =:= 0,
 1670		K =:= 1
 1671	    }
 1672	->  % indep var
 1673	    []
 1674	;   {nf2sum(H,I,Sum)},
 1675	    [V = Sum]
 1676	).
 1677dump_var(t_L(L),V,I,H) -->
 1678	!,
 1679	dump_var(t_l(L),V,I,H).
 1680% case lowerbound: V >= L or V > L
 1681% say V >= L, and V = K*V1 + ... + I, then K*V1 + ... + I >= L
 1682% and K*V1 + ... >= L-I and V1 + .../K = (L-I)/K
 1683dump_var(t_l(L),V,I,H) -->
 1684	!,
 1685	{
 1686	    H = [l(_*K,_)|_], % avoid 1 >= 0
 1687	    get_attr(V,clpqr_itf,Att),
 1688	    arg(3,Att,strictness(Strict)),
 1689	    Sm is Strict /\ 2,
 1690	    Kr is 1 rdiv K,
 1691	    Li is Kr*(L - I),
 1692	    mult_hom(H,Kr,H1),
 1693	    nf2sum(H1,0,Sum),
 1694	    (   K > 0 % K > 0
 1695	    ->	dump_strict(Sm,Sum >= Li,Sum > Li,Result)
 1696	    ;   dump_strict(Sm,Sum =< Li,Sum < Li,Result)
 1697	    )
 1698	},
 1699	[Result].
 1700dump_var(t_U(U),V,I,H) -->
 1701	!,
 1702	dump_var(t_u(U),V,I,H).
 1703dump_var(t_u(U),V,I,H) -->
 1704	!,
 1705	{
 1706	    H = [l(_*K,_)|_], % avoid 0 =< 1
 1707	    get_attr(V,clpqr_itf,Att),
 1708	    arg(3,Att,strictness(Strict)),
 1709	    Sm is Strict /\ 1,
 1710	    Kr is 1 rdiv K,
 1711	    Ui is Kr*(U-I),
 1712	    mult_hom(H,Kr,H1),
 1713	    nf2sum(H1,0.0,Sum),
 1714	    (   K > 0
 1715	    ->	dump_strict(Sm,Sum =< Ui,Sum < Ui,Result)
 1716	    ;   dump_strict(Sm,Sum >= Ui,Sum > Ui,Result)
 1717	    )
 1718	},
 1719	[Result].
 1720dump_var(t_Lu(L,U),V,I,H) -->
 1721	!,
 1722	dump_var(t_l(L),V,I,H),
 1723	dump_var(t_u(U),V,I,H).
 1724dump_var(t_lU(L,U),V,I,H) -->
 1725	!,
 1726	dump_var(t_l(L),V,I,H),
 1727	dump_var(t_u(U),V,I,H).
 1728dump_var(t_lu(L,U),V,I,H) -->
 1729	!,
 1730	dump_var(t_l(L),V,I,H),
 1731	dump_var(t_U(U),V,I,H).
 1732dump_var(T,V,I,H) --> % should not happen
 1733	[V:T:I+H].
 1734
 1735% dump_strict(FilteredStrictness,Nonstrict,Strict,Res)
 1736%
 1737% Unifies Res with either Nonstrict or Strict depending on FilteredStrictness.
 1738% FilteredStrictness is the component of strictness related to the bound: 0
 1739% means nonstrict, 1 means strict upperbound, 2 means strict lowerbound,
 1740% 3 is filtered out to either 1 or 2.
 1741
 1742dump_strict(0,Result,_,Result).
 1743dump_strict(1,_,Result,Result).
 1744dump_strict(2,_,Result,Result).
 1745
 1746% dump_nz(V,H,I,Dump,DumpTail)
 1747%
 1748% Returns in Dump a representation of the nonzero constraint of variable V
 1749% which has linear
 1750% equation H + I.
 1751
 1752dump_nz(_,H,I) -->
 1753	{
 1754	    H = [l(_*K,_)|_],
 1755	    Kr is 1 rdiv K,
 1756	    I1 is -Kr*I,
 1757	    mult_hom(H,Kr,H1),
 1758	    nf2sum(H1,0,Sum)
 1759	},
 1760	[Sum =\= I1].
 1761
 1762		 /*******************************
 1763		 *	       SANDBOX		*
 1764		 *******************************/
 1765:- multifile
 1766	sandbox:safe_primitive/1. 1767
 1768sandbox:safe_primitive(bv_q:inf(_,_)).
 1769sandbox:safe_primitive(bv_q:inf(_,_,_,_)).
 1770sandbox:safe_primitive(bv_q:sup(_,_)).
 1771sandbox:safe_primitive(bv_q:sup(_,_,_,_)).
 1772sandbox:safe_primitive(bv_q:maximize(_)).
 1773sandbox:safe_primitive(bv_q:minimize(_))