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