1%
    2% CLP(BNR) == Constraints On Boolean, Integer, and Real Intervals
    3%
    4/*	The MIT License (MIT)
    5 *
    6 *	Copyright (c) 2019-2025 Rick Workman
    7 *
    8 *	Permission is hereby granted, free of charge, to any person obtaining a copy
    9 *	of this software and associated documentation files (the "Software"), to deal
   10 *	in the Software without restriction, including without limitation the rights
   11 *	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
   12 *	copies of the Software, and to permit persons to whom the Software is
   13 *	furnished to do so, subject to the following conditions:
   14 *
   15 *	The above copyright notice and this permission notice shall be included in all
   16 *	copies or substantial portions of the Software.
   17 *
   18 *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   19 *	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   20 *	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   21 *	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   22 *	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
   23 *	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   24 *	SOFTWARE.
   25 */
   26
   27:- module(clpBNR,          % SWI module declaration
   28	[
   29	op(700, xfx, ::),
   30	(::)/2,                % declare interval
   31	{}/1,                  % define constraint
   32	add_constraint/1,      % more primitive define single constraint, bypass simplify 
   33	interval/1,            % filter for clpBNR constrained var
   34	interval_degree/2,     % number of constraints on clpBNR constrained var
   35	interval_goals/2,      % list of goals to build clpBNR constrained var
   36	list/1,                % O(1) list filter (also for compatibility)
   37	domain/2, range/2,     % get type and bounds (domain)
   38	delta/2,               % width (span) of an interval or numeric (also arithmetic function)
   39	midpoint/2,            % midpoint of an interval (or numeric) (also arithmetic function)
   40	median/2,              % median of an interval (or numeric) (also arithmetic function)
   41	lower_bound/1,         % narrow interval to point equal to lower bound
   42	upper_bound/1,         % narrow interval to point equal to upper bound
   43
   44	% additional constraint operators
   45	op(200, fy, ~),        % boolean 'not'
   46	op(500, yfx, and),     % boolean 'and'
   47	op(500, yfx, or),      % boolean 'or'
   48	op(500, yfx, nand),    % boolean 'nand'
   49	op(500, yfx, nor),     % boolean 'nor'
   50	% op(500, yfx, xor),   % boolean 'xor', but operator already defined (P=400) for arithmetic
   51	op(700, xfx, <>),      % integer not equal
   52	op(700, xfx, <=),      % included (one way narrowing)
   53
   54	% utilities
   55	print_interval/1, print_interval/2,      % pretty print interval with optional stream
   56	small/1, small/2,      % defines small interval width based on precision value
   57	solve/1, solve/2,      % solve (list of) intervals using split to find point solutions
   58	splitsolve/1, splitsolve/2,   % solve (list of) intervals using split
   59	absolve/1, absolve/2,  % absolve (list of) intervals, narrows by nibbling bounds
   60	enumerate/1,           % "enumerate" integers
   61	global_minimum/2,      % find interval containing global minimum(s) for an expression
   62	global_minimum/3,      % global_minimum/2 with definable precision
   63	global_maximum/2,      % find interval containing global maximum(s) for an expression
   64	global_maximum/3,      % global_maximum/2 with definable precision
   65	global_minimize/2,     % global_minimum/2 plus narrow vars to found minimizers
   66	global_minimize/3,     % global_minimum/3 plus narrow vars to found minimizers
   67	global_maximize/2,     % global_maximum/2 plus narrow vars to found maximizers
   68	global_maximize/3,     % global_maximum/3 plus narrow vars to found maximizers
   69	nb_setbounds/2,        % non-backtracking set bounds (use with branch and bound)
   70	partial_derivative/3,  % differentiate Exp wrt. X and simplify
   71	clpStatistics/0,       % reset
   72	clpStatistic/1,        % get selected
   73	clpStatistics/1,       % get all defined in a list
   74	watch/2,               % enable monitoring of changes for interval or (nested) list of intervals
   75	trace_clpBNR/1         % enable/disable tracing of clpBNR narrowing ops
   76	]).   77
   78/*		missing(?) functionality from original CLP(BNR): utility accumulate/2.		*/
   79
   80/* supported interval relations:
   81
   82+	-	*	/                         %% arithmetic
   83** ^                                  %% includes real exponent, odd/even integer
   84abs                                   %% absolute value
   85sqrt                                  %% positive square root
   86min	max                               %% binary min/max
   87==	is	<>	=<	>=	<	>             %% comparison (`is` synonym for `==`)
   88<=	                                  %% included (one way narrowing)
   89and	or	nand	nor	xor	->	,         %% boolean (`,` synonym for `and`)
   90-	~                                 %% unary negate and not
   91exp	log                               %% exp/ln
   92sin	asin	cos	acos	tan	atan      %% trig functions
   93integer                               %% must be an integral value, e.g., 1 and 1.0 are both integral
   94sig                                   %% signum of real, (-1,0,+1)
   95
   96*/

  117version("0.12.3").
  118
  119% support various optimizations via goal expansion
  120:- discontiguous clpBNR:goal_expansion/2.  121
  122% debug feature control and messaging
  123:- if(exists_source(swish(lib/swish_debug))).  124	:- create_prolog_flag(clpBNR_swish, true, [access(read_only)]).  125	:- use_module(swish(lib/swish_debug)).  126	:- use_module(library(http/html_write)).  127:- else.  128	:- use_module(library(debug)).  129:- endif.  130
  131:- use_module(library(prolog_versions)).  % SWIP dependency enforcement
  132
  133:- require_prolog_version('9.1.22',       % properly exported arithmetic functions
  134	          [ rational   % require rational number support, implies bounded=false
  135	          ]).  136
  137%
  138% Optimize arithmetic, but not debug.
  139%
  140set_optimize_flags_ :-      % called at start of load
  141	set_prolog_flag(optimise,true),              % scoped to file/module
  142	current_prolog_flag(optimise_debug,ODflag),  % save; restore in initialization
  143	nb_linkval('$optimise_debug_save',ODflag),
  144	set_prolog_flag(optimise_debug,false).       % so debug/3, debugging/1 don't get "optimized"
  145
  146restore_optimize_flags_ :-  % called at module initialization (after load)
  147	nb_getval('$optimise_debug_save',ODflag), nb_delete('$optimise_debug_save'),
  148	set_prolog_flag(optimise_debug,ODflag).
  149
  150:- set_optimize_flags_.  151
  152% local debug and trace message support
  153debug_clpBNR_(FString,Args) :- debug(clpBNR,FString,Args).
  154
  155% sandboxing for SWISH
  156:- multifile(sandbox:safe_prolog_flag/1).  157:- multifile(sandbox:safe_global_variable/1).  158:- multifile(sandbox:safe_primitive/1).  159:- multifile(sandbox:safe_meta/2).  160
  161current_node_(Node) :-  % look back to find current Op being executed for debug messages
  162	prolog_current_frame(Frame),  % this is a little grungy, but necessary to get intervals
  163	prolog_frame_attribute(Frame,parent_goal,doNode_(Args,Op,_,_,_,_,_)),
  164	map_constraint_op_(Op,Args,Node),
  165	!.
  166
  167sandbox:safe_primitive(clpBNR:current_node_(_Node)). 
  168
  169%
  170% statistics
  171%
  172
  173% assign,increment/read global counter (assumed to be ground value so use _linkval)
  174g_assign(G,V)  :- nb_linkval(G,V).
  175g_inc(G)       :- nb_getval(G,N), N1 is N+1, nb_linkval(G,N1).
  176g_incb(G,I)    :- nb_getval(G,N), N1 is N+I, b_setval(G,N1).    % undone on backtrack
  177g_read(G,V)    :- nb_getval(G,V).
  178
  179sandbox:safe_global_variable('$clpBNR:thread_init_done').
  180sandbox:safe_global_variable('$clpBNR:userTime').
  181sandbox:safe_global_variable('$clpBNR:inferences').
  182sandbox:safe_global_variable('$clpBNR:gc_time').
  183
  184%  
  185% Global var statistics are per thread and therefore must be "lazily" initialized
  186% Also ensures that thread copies of flags are properly set.
  187% This exception handler will be invoked the first time '$clpBNR:thread_init_done' is read
  188% Public predicates ::, {}, clpStatistics/1 and range read this global before proceeding. 
  189%
  190user:exception(undefined_global_variable,'$clpBNR:thread_init_done',retry) :- !,
  191	set_prolog_flags,     % initialize/check environment flags  
  192	clpStatistics,        % defines remaining globals 
  193	g_assign('$clpBNR:thread_init_done',1).
  194
  195%
  196% Define custom clpBNR flags when module loaded
  197%
  198:- create_prolog_flag(clpBNR_iteration_limit,3000,[type(integer),keep(true)]).  199:- create_prolog_flag(clpBNR_default_precision,6,[type(integer),keep(true)]).  200:- create_prolog_flag(clpBNR_verbose,false,[type(boolean),keep(true)]).  201
  202sandbox:safe_prolog_flag(clpBNR_iteration_limit,_).
  203sandbox:safe_prolog_flag(clpBNR_default_precision,_).
  204sandbox:safe_prolog_flag(clpBNR_verbose,_).
  205%
  206% Set public flags (see module/thread initialization)
  207%
  208set_prolog_flags :-
  209	set_prolog_flag(prefer_rationals, true),           % enable rational arithmetic
  210	(current_prolog_flag(max_rational_size,_)
  211	 -> true                                           % already defined, leave as is
  212	 ;  set_prolog_flag(max_rational_size, 16)         % rational size in bytes before ..
  213	),
  214	set_prolog_flag(max_rational_size_action, float),  % conversion to float
  215
  216	set_prolog_flag(float_overflow,infinity),          % enable IEEE continuation values
  217	set_prolog_flag(float_zero_div,infinity),
  218	set_prolog_flag(float_undefined,nan),
  219	set_prolog_flag(write_attributes,portray).         % thread-local, init to 'portray'

 clpStatistics is det

Resets clpBNR statistics - always succeeds.

clpBNR collects a number of "operational measurements" on a per-thread basis and combines them with some system statistics for subsequent querying. clpBNR measurements include:

narrowingOps	number of interval primitives called
narrowingFails	number of interval primitive failures
node_count	number of nodes in clpBNR constraint network
max_iterations	maximum number of iterations before throttling occurs (max/limit)

System statistics included in clpStatistics:

userTime	from statistics:cputime
gcTime	from statistics:garbage_collection.Time
globalStack	from statistics:globalused/statistics:global
trailStack	from statistics:trailused/statistics:trail
localStack	from statistics:localused/statistics:local
inferences	from statistics:inferences

/

 clpStatistic(?S) is nondet

Succeeds if S unifies with a clpStatistic value; otherwise fails. On backtracking all values that unify with S will be generated. Examples:

?- clpStatistics, X::real, {X**4-4*X**3+4*X**2-4*X+3==0}, clpStatistic(narrowingOps(Ops)).
Ops = 2245,
X::real(-1.509169756145379, 4.18727500493995).

?- clpStatistics, X::real, {X**4-4*X**3+4*X**2-4*X+3==0}, clpStatistic(S).
S = userTime(0.02277600000000035),
X::real(-1.509169756145379, 4.18727500493995) ;
S = gcTime(0.0),
X::real(-1.509169756145379, 4.18727500493995) ;
S = globalStack(43696/524256),
X::real(-1.509169756145379, 4.18727500493995) ;
S = trailStack(664/133096),
X::real(-1.509169756145379, 4.18727500493995) ;
S = localStack(1864/118648),
X::real(-1.509169756145379, 4.18727500493995) ;
S = inferences(86215),
X::real(-1.509169756145379, 4.18727500493995) ;
S = narrowingOps(2245),
X::real(-1.509169756145379, 4.18727500493995) ;
S = narrowingFails(0),
X::real(-1.509169756145379, 4.18727500493995) ;
S = node_count(9),
X::real(-1.509169756145379, 4.18727500493995) ;
S = max_iterations(2245/3000),
X::real(-1.509169756145379, 4.18727500493995).

*/

 clpStatistics(?Ss) is semidet

Succeeds if Ss unifies with a list of clpStatistic's values; otherwise fails. Example:

?- clpStatistics, X::real, {X**4-4*X**3+4*X**2-4*X+3==0}, clpStatistics(Ss).
Ss = [userTime(0.023398999999999504), gcTime(0.001), globalStack(19216/131040), trailStack(1296/133096), localStack(2184/118648), inferences(82961), narrowingOps(2245), narrowingFails(0), node_count(9), max_iterations(2245/3000)],
X::real(-1.509169756145379, 4.18727500493995).

/

  285:- discontiguous clpBNR:clpStatistics/0, clpBNR:clpStatistic/1.  286
  287clpStatistics :-
  288	% garbage_collect,  % ? do gc before time snapshots
  289	statistics(cputime,T), g_assign('$clpBNR:userTime',T),   % thread based
  290	statistics(inferences,I), g_assign('$clpBNR:inferences',I),
  291	statistics(garbage_collection,[_,_,G,_]), g_assign('$clpBNR:gc_time',G),
  292	fail.  % backtrack to reset other statistics.
  293
  294clpStatistic(_) :- g_read('$clpBNR:thread_init_done',0).  % ensures per-thread initialization then fails
  295
  296clpStatistic(userTime(T)) :- statistics(cputime,T1), g_read('$clpBNR:userTime',T0), T is T1-T0.
  297
  298clpStatistic(gcTime(G)) :- statistics(garbage_collection,[_,_,G1,_]), g_read('$clpBNR:gc_time',G0), G is (G1-G0)/1000.0.
  299
  300clpStatistic(globalStack(U/T)) :- statistics(globalused,U), statistics(global,T).
  301
  302clpStatistic(trailStack(U/T)) :- statistics(trailused,U), statistics(trail,T).
  303
  304clpStatistic(localStack(U/T)) :- statistics(localused,U), statistics(local,T).
  305
  306clpStatistic(inferences(I)) :- statistics(inferences,I1), g_read('$clpBNR:inferences',I0), I is I1-I0.

 list(?X:list) is semidet

Succeeds if X is a list; otherwise fails. Note: not equivalent to is_list/1 but executes in O(1) time. This filter is provided for historical compatability. /

  314list(X) :- compound(X) ->  X=[_|_] ; X==[].
  315
  316%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  317%
  318%  SWI-Prolog implementation of IA
  319%
  320%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  321%
  322% Intervals are constrained (attributed) variables.
  323%
  324% Current interval bounds updates via setarg(Val) which is backtrackable

 interval(?X:interval) is semidet

Succeeds if X is an interval, i.e., a variable with a clpBNR attribute; otherwise fails. /

  330interval(Int) :- get_attr(Int, clpBNR, _).

 interval_degree(?X:numeric, ?N:integer) is semidet

Succeeds if X is numeric and N = number of clpBNR constraints on X; otherwise fails. If X is a number, N = 0. Examples:

?- {X==Y+1}, interval_degree(X,N).
N = 1,
X::real(-1.0Inf, 1.0Inf),
Y::real(-1.0Inf, 1.0Inf).

?- interval_degree(42,N).
N = 0.

/

  346interval_degree(X, N) :-
  347	number(X)
  348	-> N = 0
  349	;  interval_object(X, _, _, Nodelist), % fail if not a number or interval
  350	   len_nodelist(Nodelist,0,N).         % current number of elements in (indefinite) Nodelist
  351
  352len_nodelist(T,N,N) :- var(T), !.          % end of indefinite list
  353len_nodelist([_|T],Nin,N) :- 
  354	Nout is Nin+1,
  355	len_nodelist(T,Nout,N).
  356
  357% internal abstraction
  358interval_object(Int, Type, Val, Nodelist) :-
  359	get_attr(Int, clpBNR, interval(Type, Val, Nodelist, _)).
  360
  361% define a new or replace existing clpBNR attribute
  362define_interval(Int, Type, Val, Nodelist, Flags) :-
  363	put_attr(Int, clpBNR, interval(Type, Val, Nodelist, Flags)).
  364
  365% update constraints in Nodelist
  366update_interval_nodelist_(Int,NodeList) :-
  367	get_attr(Int, clpBNR, Def),
  368	setarg(3,Def,NodeList).
  369
  370% flags (list)  abstraction
  371get_interval_flags_(Int, Flags) :-
  372	get_attr(Int, clpBNR, Attr),
  373	arg(4,Attr,Flags).     % interval(_, _, _, Flags)).
  374
  375set_interval_flags_(Int, Flags) :-  % flags assumed to be ground so no copy required
  376	get_attr(Int, clpBNR, Attr),
  377	setarg(4,Attr,Flags).  % interval(_, _, _, Flags)).
  378
  379%
  380% Interval value constants
  381%
  382universal_interval((-1.0Inf,1.0Inf)).
  383
  384empty_interval((1.0Inf,-1.0Inf)).
  385
  386new_universal_interval(boolean,Int) :- !,
  387	define_interval(Int, integer, (0,1), _Nodelist, []).
  388new_universal_interval(Type,Int) :-
  389	universal_interval(UI),
  390	define_interval(Int, Type, UI, _Nodelist, []).
  391
  392% Finite intervals 
  393finite_interval(real,    (-1.0e+16,1.0e+16)).
  394finite_interval(integer, I) :- 
  395	(current_prolog_flag(min_tagged_integer,L), current_prolog_flag(max_tagged_integer,H)
  396	 -> I = (L,H)  % SWIP: use tagged limits for finite default
  397	 ;  I = (-10000000000,10000000000)  % otherwise just big ints
  398	).
  399finite_interval(boolean, (0,1)).
  400
  401finite_int((L,H)) :- \+(L = -1.0Inf ; H = 1.0Inf).    % a test for a finite interval
  402
  403zero_int((L,H)) :- L==H, L=:=0.                       % a test for point interval zero
  404% used by primitive op so expand inline
  405goal_expansion(zero_int(I), (I=(L,H), L==H, L=:=0)).  % SWIP inline optimization
  406
  407% legal integer bounds values
  408integerBnd(1.0Inf).
  409integerBnd(-1.0Inf).
  410integerBnd(B) :- integer(B).
  411
  412% precise bounds values - Note: assumes external commit so cuts not req'd in first two clauses
  413preciseBnd(1.0Inf).
  414preciseBnd(-1.0Inf).
  415preciseBnd(1.5NaN) :- !, fail.
  416preciseBnd(B) :- preciseNumber(B). 
  417
  418preciseNumber(B) :-
  419	rational(B) -> true 
  420	; number(B), 0 is cmpr(B,rationalize(B)). % rational(B)=:=rationalize(B), fails if float not precise
  421
  422:- include(clpBNR/ia_primitives).  % interval arithmetic relations via evalNode/4.

 nb_setbounds(?X:interval, +Bs:number_list) is semidet

Succeeds if X is an interval and can be narrowed to the bounds Bs = [L,U]; otherwise fails. On backtracking, this value is not undone.

Caution: this predicate is non-logical and intended for specialized use case, e.g., some branch-and-bound algorithms (narrow to current solution, then backtrack to next solution). /

  431%
  432%  non-backtracking set bounds for use with branch and bound
  433%
  434nb_setbounds(Int, [L,U]) :- 
  435	number(L), number(U),
  436	get_attr(Int, clpBNR, Def),
  437	arg(2, Def, Val),             % WAM code
  438	^(Val,(L,U),NewVal),          % new range is intersection (from ia_primitives)
  439	nb_setarg(2, Def, NewVal).
  440
  441%
  442% get current value
  443%
  444getValue(Int, Val) :- 
  445	number(Int)
  446	 -> Val=(Int,Int)                                   % numeric point value
  447	 ;  get_attr(Int, clpBNR, interval(_, Val, _, _)).  % interval, optimized for SWIP
  448
  449%
  450% set interval value (assumes bounds check has already been done)
  451% Note: putValue_ cannot modify the new value since the narrowing op is still in the Agenda
  452%	and may not be run again. Insert `integral` op for integer bounds adjustment.
  453%
  454putValue_(New, Int, NodeList) :- 
  455	get_attr(Int, clpBNR, Def)             % still an interval ?
  456	 -> (debugging(clpBNR) -> check_monitor_(Int, New, Def) ; true),
  457	    Def = interval(Type,_,Nodes,_),    % get Type and Nodes for queueing
  458	    New = (L,H),
  459	    ( 0 is cmpr(L,H)                   % point interval ->
  460	     -> setarg(3,Def,_NL),             % clear node list (so nothing done in unify)
  461	        pointValue_(L,H,Int),          % unify, avoiding -0.0 and preferring rational (invokes hook)
  462	        NodeList = Nodes               % still add Nodes to Agenda
  463	     ;  setarg(2,Def,New),             % update value in interval (backtrackable)
  464	        % if integer has non-integral bounds, schedule `integral` op to fix it
  465	        (Type == integer
  466	         -> ( integerBnd(L), integerBnd(H) -> NodeList = Nodes ; NodeList = [node(integral,_,0,$(Int))|_] )
  467	         ;  NodeList = Nodes
  468	        )
  469	    )
  470	 ; true.  % no longer an interval, NodeList is empty
  471
  472pointValue_(-0.0,_,0.0) :-!.
  473pointValue_(L,H,Int) :- (rational(L) -> Int = L ; Int = H).

 range(?X, ?Bs:number_list) is semidet

Succeeds if X is numeric and Bs unifies with a list containing the lower and upper bound of X; otherwise fails. If X is a logic variable range(X,[2,3]) is equivalent to X::real(2,3). If X is a number the lower and upper bounds are the same. Examples:

?- X::integer(1,10), range(X,Bs).
Bs = [1, 10],
X::integer(1, 10).

?- range(42,Bs).
Bs = [42, 42].

?- range(X,[2,3]).
X::real(2, 3).

/

  491%
  492%  range(Int, Bounds) for compatibility 
  493%
  494range(Int, [L,H]) :- getValue(Int, (IL,IH)), !,  % existing interval or number =>  number(IL)
  495	(var(L) -> L=IL ; non_empty(L,IL)),          % range check (L=<IL,IH=<H), no narrowing
  496	(var(H) -> H=IH ; non_empty(IH,H)).
  497range(Int, [L,H]) :- var(Int),  % for other var(Int), declare it to a real interval with specified bounds
  498	Int::real(L,H).

 domain(?X:interval, ?Dom) is semidet

Succeeds if X is an interval and Dom unifies with the domain of X; otherwise fails. Dom is a compound term with a functor specifying the type (real or integer) and the arguments specifying the bounds. If X has a domain of integer(0,1), Dom will be boolean. Examples:

?- range(X,[2,3]), domain(X,Dom).
Dom = real(2, 3),
X::real(2, 3).

?- X::integer(0,1),domain(X,Dom).
Dom = boolean,
X::boolean.

?- domain(X,Dom).
false.

Note: unlike range/2, domain/2 will not change X. /

  518domain(Int, Dom) :-
  519	interval_object(Int, Type, Val, _),
  520	interval_domain_(Type, Val, Dom).
  521
  522interval_domain_(integer,(0,1),boolean) :- !.  % integer(0,1) is synonymous with boolean
  523interval_domain_(T,(L,H),Dom) :- Dom=..[T,L,H].
  524
  525:- use_module(library(arithmetic), [arithmetic_function/1]).   % extended arithmetic functions

 delta(?X:numeric, ?W:number) is semidet

Succeeds if X is numeric and W unifies with the width of X (upper bound-lowerbound); otherwise fails. Examples:

?- X:: real(1r2,5r3),delta(X,D).
D = 7r6,
X::real(0.5, 1.6666666666666667).

?- delta(42,W).
W = 0.

delta is also available as an arithmetic function:

?- X::real(1r2,pi), W is delta(X).
W = 2.6415926535897936,
X::real(0.5, 3.1415926535897936).

/

  546%
  547%  delta(Int, Wid) width/span of an interval or numeric value, can be infinite
  548%
  549:- arithmetic_function(delta/1).  550
  551delta(Int, Wid) :-
  552	getValue(Int,(L,H)),
  553	Wid is roundtoward(H-L,to_positive).

 midpoint(?X:numeric, ?M:number) is semidet

Succeeds if X is numeric and M unifies with the midpoint of X; otherwise fails. Examples:

?- X:: real(1r2,5r3), midpoint(X,M).
M = 13r12,
X::real(0.5, 1.6666666666666667).

?- midpoint(42,M).
M = 42.

midpoint is also available as an arithmetic function:

?- X::real(1r2,pi), M is midpoint(X).
M = 1.8207963267948968,
X::real(0.5, 3.1415926535897936).

/

  574%
  575%  midpoint(Int, Wid) midpoint of an interval or numeric value
  576% based on:
  577%	Frédéric Goualard. How do you compute the midpoint of an interval?. 
  578%	ACM Transactions on Mathematical Software, Association for Computing Machinery, 2014,
  579%	40 (2), 10.1145/2493882. hal-00576641v1
  580% Exception, single infinite bound treated as largest finite FP value
  581%
  582:- arithmetic_function(midpoint/1).  583
  584midpoint(Int, Mid) :-
  585	getValue(Int,(L,H)),
  586	midpoint_(L,H,Mid).
  587
  588midpoint_(L,H,M)       :- L =:= -H, !, M=0.              % symmetric including (-inf,inf)
  589midpoint_(-1.0Inf,H,M) :- !, M is nexttoward(-1.0Inf,0)/2 + H/2.
  590midpoint_(L,1.0Inf,M)  :- !, M is L/2 + nexttoward(1.0Inf,0)/2.
  591midpoint_(L,H,M)       :- M1 is L/2 + H/2, M=M1.        % general case

 median(?X:numeric, ?M:float) is semidet

Succeeds if X is numeric and M unifies with the median of X; otherwise fails. The median is 0 if the domain of X contains 0; otherwise it is the floating point value which divides the interval into two sub-domains each containing approximately equal numbers of floating point values. Examples:

?- X:: real(1r2,5r3), median(X,M).
M = 0.9128709291752769,
X::real(0.5, 1.6666666666666667).

?- median(42,M).
M = 42.0.

median is also available as an arithmetic function:

?- X::real(1r2,pi), M is median(X).
M = 1.2533141373155003,
X::real(0.5, 3.1415926535897936).

/

  612%
  613% median(Int,Med) from CLP(RI)
  614% Med = 0 if Int contains 0, else a number which divides Int into equal
  615% numbers of FP values. Med is always a float
  616%
  617:- arithmetic_function(median/1).  618
  619median(Int, Med) :-
  620	getValue(Int,(L,H)),
  621	median_bound_(lo,L,FL),
  622	median_bound_(hi,H,FH),
  623	median_(FL,FH,Med), !.
  624	
  625median_bound_(lo,B,FB) :- B=:=0, FB is nexttoward(B,1.0).
  626median_bound_(lo,-1.0Inf,FB) :- FB is nexttoward(-1.0Inf,0.0).
  627median_bound_(lo,B,FB) :- FB is roundtoward(float(B), to_negative).
  628
  629median_bound_(hi,B,FB) :- B=:=0, !, FB is nexttoward(B,-1.0).
  630median_bound_(hi,1.0Inf,FB) :- FB is nexttoward(1.0Inf,0.0).
  631median_bound_(hi,B,FB) :- FB is roundtoward(float(B), to_positive).
  632
  633median_(B,B,B).                          % point interval
  634median_(L,H,0.0) :- L < 0.0, H > 0.0.    % contains 0: handles (-inf,inf)
  635median_(L,H,M)   :- M is copysign(sqrt(abs(L))*sqrt(abs(H)),L).      % L and H have same sign
  636
  637%
  638%  lower_bound and upper_bound
  639%

 lower_bound(?X:numeric) is semidet

Succeeds if X is numeric and unifies with the lower bound of its domain. Examples:

?- X::integer(1,10),lower_bound(X).
X = 1.

?- X = 42, lower_bound(X).
X = 42.

Note that lower_bound will unify X with a number on success, but it may fail if this value is inconsistent with current constraints. /

  653lower_bound(Int) :-
  654	getValue(Int,(L,_H)),
  655	Int=L.

 upper_bound(?X:numeric) is semidet

Succeeds if X is numeric and unifies with the upper bound of its domain. Examples:

?- X::integer(1,10),upper_bound(X).
X = 10.

?- X = 42, upper_bound(X).
X = 42.

Note that upper_bound will unify X with a number on success, but it may fail if this value is inconsistent with current constraints. /

  670upper_bound(Int) :-
  671	getValue(Int,(_L,H)),
  672	Int=H.

 ::(-X:numeric_List, ?Dom) is semidet

Succeeds if variable X has domain Dom; otherwise fails. If Dom, or either bound of Dom, is a variable (or missing), it will be unified with the default value depending on its type. Default domains are real(-1.0e+16, 1.0e+16) and integer(-72057594037927936, 72057594037927935). Examples:

?- X::real(-pi/2,pi/2).
X::real(-1.5707963267948968, 1.5707963267948968).

?- X::real, Y::integer.
X::real(-1.0e+16, 1.0e+16),
Y::integer(-72057594037927936, 72057594037927935).

?- Y::integer(1,_), Y::Dom.
Dom = integer(1, 72057594037927935),
Y::integer(1, 72057594037927935).

?- B::boolean.
B::boolean.

?- 42::Dom.
false.

Note that bounds can be defined using arithmetic expressions.

Alternatively, the first argument may be a list of variables:

?- [B1,B2,B3]::boolean.
B1::boolean,
B2::boolean,
B3::boolean.

?- length(Vs,3), Vs::real(-1,1).
Vs = [_A, _B, _C],
_A::real(-1, 1),
_B::real(-1, 1),
_C::real(-1, 1).

Numeric strings can be used to define real domains with explicit precision:

?- X::"0.5".
X::real(0.44999999999999996, 0.5500000000000002).

?- X::"0.5000", domain(X,Dom).
Dom = real(0.49994999999999995, 0.5000500000000001),
X:: 0.500... .

?- X::"100".
X::real(99.94999999999999, 100.05000000000001).

?- X::"0.5e3".
false.

Note that floating point exponent format in numeric strings is not supported . /

  728Rs::Dom :- list(Rs),!,                    % list of vars
  729	intervals_(Rs,Dom).
  730
  731R::Dom :-                                 % single var
  732	g_read('$clpBNR:thread_init_done',_), % ensures per-thread initialization 
  733	(var(Dom)                             % domain undefined
  734	 -> (var(R) -> int_decl_(real,_,R) ; true),  % default or domain query (if interval(R) else fail)
  735	    domain(R,Dom)                     % extract domain
  736	 ;  (string(Dom)
  737	     -> Type = real,                  % Dom = numeric string
  738	        Val = Dom
  739	     ;  Dom=..[Type|Bounds],          % Dom = Type(L,H)
  740	        Val=..[','|Bounds]
  741	    ),
  742	    int_decl_(Type,Val,R)
  743	).
  744
  745intervals_([],_Dom).
  746intervals_([Int|Ints],Dom) :-
  747	copy_term(Dom,CDom),              % need a fresh copy of domain for each interval
  748	Int::CDom, !,
  749	intervals_(Ints,Dom).
  750
  751int_decl_(boolean,_,R) :- !,          % boolean is integer; 0=false, 1=true, ignore any bounds.
  752	int_decl_(integer,(0,1),R).
  753int_decl_(Type,(','),R) :- !,         % no bounds, fill with vars
  754	int_decl_(Type,(_,_),R).
  755int_decl_(Type,Val,R) :- string(Val), !,
  756	number_string(Num,Val),           % rounded numeric string
  757	(sub_string(Val,_,1,FracLen,'.')
  758	 -> sub_string(Val,_,FracLen,0,FracString),  % contains '.'
  759	    number_string(N,FracString), integer(N)  % Frac must be integer (no exponent)
  760	 ;  integer(Num), FracLen = 1                % no '.', Num must be an integer
  761	),
  762	% from number of fractional digits calculate rounding value as a float 
  763	Rnd is 0.5/10**FracLen,           % should be precise
  764	int_decl_(Type,(Num-Rnd,Num+Rnd),R).
  765int_decl_(Type,Val,R) :- interval_object(R,CType,CVal,_NL), !,  % already interval
  766	(Type = CType, Val = CVal         % query, no change
  767	 -> true
  768	 ;	Val = (L,H),                  % changing type,bounds?
  769	    lower_bound_val_(Type,L,IL),
  770	    upper_bound_val_(Type,H,IH),
  771	    applyType_(Type, R, T/T, Agenda),           % coerce reals to integers (or no-op).
  772	    ^(CVal,(IL,IH),New),          % low level functional primitive
  773	    updateValue_(CVal, New, R, 1, Agenda, NewAgenda),  % update value (Agenda empty if no value change)
  774	    stable_(NewAgenda)            % then execute Agenda
  775	).
  776int_decl_(Type,(L,H),R) :- var(R), !, % new interval (R can't be existing interval)
  777	lower_bound_val_(Type,L,IL),
  778	upper_bound_val_(Type,H,IH),
  779	C is cmpr(IL,IH),  % compare bounds
  780	(C == 0
  781	 -> (rational(IL) -> R=IL ; R = IH)  % point range, can unify now
  782	 ;  C == -1,                         % necessary condition: IL < IH
  783	    define_interval(R, Type, (IL,IH), _Nodelist, [])
  784	).
  785int_decl_(Type,(L,H),R) :- constant_type_(Type,R), % R already a point value, check range
  786	lower_bound_val_(Type,L,IL), non_empty(IL,R),  % IL=<R,
  787	upper_bound_val_(Type,H,IH), non_empty(R,IH).  % R=<IH.
  788
  789lower_bound_val_(Type,L,L) :- var(L), !,   % unspecified bound, make it finite
  790	finite_interval(Type,(L,_)).
  791lower_bound_val_(real,L,IL) :-             % real: evaluate and round outward (if float)
  792	((L == pi ; L == e)
  793	 -> IL is roundtoward(L,to_negative)
  794	 ;  Lv is L,
  795	    (preciseBnd(Lv) -> IL=Lv ; IL is nexttoward(Lv,-1.0Inf)) 
  796	).
  797lower_bound_val_(integer,L,IL) :-          % integer: make integer
  798	IL is ceiling(L).
  799lower_bound_val_(boolean,L,IL) :-          % boolean: narrow to L=0
  800	IL is max(0,ceiling(L)).
  801
  802upper_bound_val_(Type,H,H) :- var(H), !,   % unspecified bound, make it finite
  803	finite_interval(Type,(_,H)).
  804upper_bound_val_(real,H,IH) :-             % real: evaluate and round outward (if float)
  805	((H == pi ; H == e)
  806	 -> IH is roundtoward(H,to_positive)
  807	 ;  Hv is H,
  808	    (preciseBnd(Hv) -> IH=Hv ; IH is nexttoward(Hv,1.0Inf)) 
  809	).
  810upper_bound_val_(integer,H,IH) :-          % integer: make integer
  811	IH is floor(H).
  812upper_bound_val_(boolean,H,IH) :-          % boolean: narrow to H=1
  813	IH is min(1,floor(H)).
  814
  815constant_type_(real,C) :- number(C).
  816constant_type_(integer,C) :- integer(C), !.
  817constant_type_(integer,C) :- float(C), float_class(C,infinite).
  818
  819applyType_(NewType, Int, Agenda, NewAgenda) :-      % narrow Int to Type
  820	get_attr(Int,clpBNR,interval(Type,Val,Nodelist,Flags)),
  821	(NewType @< Type    % standard order of atoms:  boolean @< integer @< real
  822	 -> (debugging(clpBNR) -> check_monitor_(Int, NewType, interval(Type,Val,Nodelist,Flags)) ; true),
  823	    Val = (L,H),
  824	    lower_bound_val_(NewType,L,IL),
  825	    upper_bound_val_(NewType,H,IH),
  826	    (IL == IH
  827	     -> Int=IL  % narrowed to point
  828	     ;  define_interval(Int,integer,(IL,IH),Nodelist,Flags),  % set Type (only case allowed)
  829	        linkNodeList_(Nodelist, Agenda, NewAgenda)
  830	    )
  831	 ; NewAgenda = Agenda
  832	).
  833
  834%
  835% this goal gets triggered whenever an interval is unified, valid for a numeric value or another interval
  836%
  837attr_unify_hook(IntDef, Num) :-         % unify an interval with a numeric
  838	number(Num),
  839	IntDef = interval(Type,(L,H),Nodelist,_Flags),
  840	constant_type_(Type,Num),                % numeric consistent with type
  841	% L=<Num, Num=<H, assumes L < H
  842	cmpr(L,Num) + cmpr(Num,H) < 0, !,        % and in range (not NaN)
  843	(debugging(clpBNR) -> monitor_unify_(IntDef, Num) ; true),
  844	(var(Nodelist)
  845	 -> true                                 % nothing to do
  846	 ;  linkNodeList_(Nodelist, T/T, Agenda),
  847	    stable_(Agenda)                      % broadcast change
  848	).
  849
  850attr_unify_hook(interval(Type1,V1,Nodelist1,Flags1), Int) :-   % unifying two intervals
  851	get_attr(Int, clpBNR, interval(Type2,V2,Nodelist2,Flags2)),  %%SWIP attribute def.
  852	^(V1,V2,R),                     % unified range is intersection (from ia_primitives),
  853	mergeValues_(Type1, Type2, NewType, R, NewR), !,
  854	mergeNodes_(Nodelist2,Nodelist1,Newlist),  % unified node list is a merge of two lists
  855	mergeFlags_(Flags1,Flags2,Flags),
  856	(debugging(clpBNR) -> monitor_unify_(interval(Type1,V1,_,Flags), Int) ; true),
  857	% update new type, value and constraint list, undone on backtracking
  858	define_interval(Int,NewType,NewR,Newlist,Flags),
  859	(var(Newlist)
  860	 -> true                                 % nothing to do
  861	 ;  linkNodeList_(Newlist, T/T, Agenda),
  862	    stable_(Agenda)                      % broadcast change
  863	).
  864
  865attr_unify_hook(interval(Type,Val,_Nodelist,_Flags), V) :-   % new value out of range
  866	g_inc('$clpBNR:evalNodeFail'),  % count of primitive call failures
  867	debugging(clpBNR),             % fail immediately unless debugging
  868	debug_clpBNR_('Failed to unify ~w(~w) with ~w',[Type,Val,V]),
  869	fail.
  870
  871% awkward monitor case because original interval gone
  872monitor_unify_(IntDef, Update) :-  % debugging, manufacture temporary interval
  873	IntDef = interval(Type,Val,Nodelist,Flags),
  874	define_interval(Temp, Type, Val, Nodelist, Flags),
  875	check_monitor_(Temp, Update, IntDef).
  876
  877% if types identical, result type and bounds unchanged,
  878% else one is integer so result type integer, and integer bounds applied
  879mergeValues_(T, T, T, R, R) :- !.
  880mergeValues_(_, _, integer, (L,H), (IL,IH)) :-
  881	lower_bound_val_(integer,L,IL),     % type check bounds
  882	upper_bound_val_(integer,H,IH),
  883	non_empty(IL,IH).                   % still non-empty
  884
  885% optimize for one or both lists (dominant case)
  886mergeFlags_([],Flags2,Flags2) :- !.
  887mergeFlags_(Flags1,[],Flags1) :- !.
  888mergeFlags_([F1|Flags1],Flags2,Flags) :-   % discard if F in Flags2 
  889	functor(F1,N,1),                       % ignore value
  890	functor(F2,N,1),
  891	memberchk(F2,Flags2), !,
  892	mergeFlags_(Flags1,Flags2,Flags).
  893mergeFlags_([F1|Flags1],Flags2,[F1|Flags]) :-  % keep F, not in Flags2 
  894	mergeFlags_(Flags1,Flags2,Flags).
  895
  896% merge two node lists removing duplicates based on operation and arguments
  897mergeNodes_([N],NodeList,NodeList) :- var(N),!.         % end of list
  898mergeNodes_([node(Op,_,_,Ops)|Ns],NodeList,NewList) :-  % if same Op and Ops, discard
  899	matching_node_(NodeList,Op,Ops), !,
  900	mergeNodes_(Ns,NodeList,NewList).
  901mergeNodes_([N|Ns],NodeList,[N|NewList]) :-             % not a duplicate, retain
  902	mergeNodes_(Ns,NodeList,NewList).
  903
  904matching_node_([node(Op,_,_,NOps)|_Ns],Op,Ops) :-
  905	NOps==Ops, !.  % identical args
  906matching_node_([N|Ns],Op,Ops) :-
  907	nonvar(N),     % not end of list
  908	matching_node_(Ns,Op,Ops).

 {+Constraints} is semidet

Succeeds if Constraints is a sequence of one or more boolean expressions (typically equalities and inequalities) defining a set of valid and consistent constraints; otherwise fails. The ',' binary operator is used to syntactically separate individual constraints and has semantics and (similar to the use of ',' in clause bodies). Arithmetic expressions are expressed using the same set of operators used in functional Prolog arithmetic (listed below) with addition of boolean operators to support that sub-domain of reals.

Table of supported interval relations:

+ - * /	arithmetic
** ^	includes real exponent, odd/even integer
abs	absolute value
sqrt	positive square root
min max	binary min/max
== is <> =\= =< >= < >	comparison (is and =\= synonyms for == and <>)
<=	included (one way narrowing)
and or nand nor xor -> ,	boolean (`,` synonym for and)
- ~	unary negate and not (boolean)
exp log	exp/ln
sin asin cos acos tan atan	trig functions
integer	must be an integer value
sig	signum of real (-1,0,+1)

clpBNR defines the following additional operators for use in constraint expressions:

op(200, fy, ~)	boolean 'not'
op(500, yfx, and)	boolean 'and'
op(500, yfx, or)	boolean 'or'
op(500, yfx, nand)	boolean 'nand'
op(500, yfx, nor)	boolean 'nor'

Note that the comparison operators <>, =\=, '<' and '>' are unsound (due to incompleteness) over the real domain but sound over the integer domain. Strict inequality (<> and =\=) is disallowed for type real (will be converted to type integer) but < and > are supported for reals since they may be useful for things like branch and bound searches (with caution). The boolean functions are restricted to type 'boolean', constraining their argument types accordingly. Some examples:

?- {X == Y+1, Y >= 1}.
X::real(2, 1.0Inf),
Y::real(1, 1.0Inf).

?- {X == cos(X)}.
X:: 0.73908513321516... .

?- X::real, {X**4-4*X**3+4*X**2-4*X+3==0}.
X::real(-1.509169756145379, 4.18727500493995).

?- {A or B, C and D}.
C = D, D = 1,
A::boolean,
B::boolean.

Note that any variable in a constraint expression with no domain will be assigned the most general value consistent with the operator types, e.g., real(-1.0Inf,1.0Inf), boolean, etc.

Numeric strings can be used in constraint expressions for domains with explicit precision:

?- {X == pi+"0.5"}.
X::real(3.591592653589793, 3.691592653589794).

/

  965{Cons} :-
  966	g_read('$clpBNR:thread_init_done',_),     % ensures per-thread initialization
  967	term_variables(Cons, CVars),              % predeclare variables to maintain ordering
  968	declare_vars_(CVars),                     % convert any plain vars to intervals
  969	addConstraints_(Cons,T/T,Agenda),         % add constraints
  970	stable_(Agenda).                          % then execute Agenda

 add_constraint(+Constraint) is semidet

Succeeds if Constraint is a supported constraint defined by {}/1. (Conjunction of contraints is supported using `,` operator).

add_constraint/1 is a primitive version of {}/1, without expression simplification, reducing overhead in some scenarios. */

  979add_constraint(Con) :-
  980	g_read('$clpBNR:thread_init_done',_),     % ensures per-thread initialization
  981	term_variables(Con, CVars),               % predeclare variables to maintain ordering
  982	declare_vars_(CVars),
  983	constraint_(Con),    % a constraint is a boolean expression that evaluates to true
  984	buildConstraint_(Con,T/T,Agenda),
  985	stable_(Agenda).
  986
  987% **Obsolete**, use public `add_constraint`
  988% low overhead version for internal use (also used by arithmetic_types pack)
  989constrain_(C) :- 
  990	add_constraint(C).
  991
  992declare_vars_([]).
  993declare_vars_([CV|CVars]) :-
  994	% if not already an interval, constrain to real(-inf,inf)
  995	(interval(CV) -> true ; new_universal_interval(real,CV)),
  996	declare_vars_(CVars).
  997
  998%%:- include(clpBNR/ia_simplify).  % simplifies constraints to a hopefully more efficient equivalent
  999:- use_module(clpBNR/simplify_constraint). 1000
 1001addConstraints_([],Agenda,Agenda) :- !.
 1002addConstraints_([C|Cs],Agenda,NewAgenda) :-
 1003	nonvar(C),
 1004	addConstraints_(C,Agenda,NextAgenda), !,
 1005	addConstraints_(Cs,NextAgenda,NewAgenda).
 1006addConstraints_((C,Cs),Agenda,NewAgenda) :-  % Note: comma as operator
 1007	nonvar(C),
 1008	addConstraints_(C,Agenda,NextAgenda), !,
 1009	addConstraints_(Cs,NextAgenda,NewAgenda).
 1010addConstraints_(C,Agenda,NewAgenda) :-
 1011	constraint_(C),    % a constraint is a boolean expression that evaluates to true
 1012	(simplify(C,CS)-> true ; CS = C),  % possible rewrite
 1013	buildConstraint_(CS, Agenda, NewAgenda).
 1014
 1015buildConstraint_(C,Agenda,NewAgenda) :-
 1016	debug_clpBNR_('Add ~p',{C}),
 1017	% need to catch errors from ground expression evaluation
 1018	catch(build_(C, 1, boolean, Agenda, NewAgenda),_Err,fail), !.
 1019buildConstraint_(C,_Agenda,_NewAgenda) :-
 1020	debug_clpBNR_('{} failure due to bad or inconsistent constraint: ~p',{C}),
 1021	fail.
 1022
 1023%
 1024% build a node from an expression
 1025%
 1026build_(Var, Var, VarType, Agenda, NewAgenda) :-         % already an interval or new variable
 1027	var(Var), !,
 1028	(interval_object(Var,Type,_,_)                      % already an interval?
 1029	 -> (Type \== VarType                               % type narrowing?
 1030	     -> applyType_(VarType, Var, Agenda, NewAgenda) % coerces exisiting intervals to required type
 1031	     ;  NewAgenda = Agenda                          % nothing to do   
 1032	    )
 1033	 ;  new_universal_interval(VarType,Var),            % implicit interval creation
 1034	    NewAgenda = Agenda                              % nothing to schedule
 1035	).
 1036build_(::(L,H), Int, VarType, Agenda, Agenda) :-        % hidden :: feature: interval bounds literal (without fuzzing)
 1037	number(L), number(H), !,
 1038	C is cmpr(L,H),  % compare bounds
 1039	(C == 0
 1040	 -> (rational(L) -> Int=L ; Int=H)                  % point value, if either bound rational (precise), use it
 1041	 ;	C == -1,                                        % necessary condition: L < H
 1042	    (VarType = real -> true ; true),                % if undefined Type, use 'real' (efficient ignore/1)
 1043	    define_interval(Int, VarType, (L,H), _Nodelist, [])  % create clpBNR interval
 1044	).
 1045build_(NumStr, Int, _VarType, Agenda, Agenda) :-
 1046	string(NumStr), !,                                  % not necessary but optimal
 1047	int_decl_(real,NumStr,Int).                         % rounded numeric string
 1048build_(Num, Int, VarType, Agenda, Agenda) :-            % atomic value representing a numeric
 1049	atomic(Num), !,                                     % includes inf, nan, pi, e
 1050	% imprecise numbers will be fuzzed
 1051	(preciseBnd(Num) -> Int = Num ;  int_decl_(VarType,(Num,Num),Int)).
 1052build_(Exp, Num, _, Agenda, Agenda) :-                  % pre-compile any ground precise Exp
 1053	ground(Exp),
 1054	safe_(Exp),                                         % safe to evaluate using is/2
 1055	Num is Exp,
 1056	preciseBnd(Num),                                    % precise result
 1057	!.
 1058build_(Exp, Z, _, Agenda, NewAgenda) :-                 % deconstruct to primitives
 1059	Exp =.. [F|Args],
 1060	fmap_(F,Op,[Z|Args],ArgsR,Types), !,                % supported arithmetic op
 1061	build_args_(ArgsR,Objs,Types,Agenda,ObjAgenda),
 1062	newNode_(Op,Objs,ObjAgenda,NewAgenda).
 1063build_(Exp, Z, _, Agenda, NewAgenda) :-                 % user defined
 1064	Exp =.. [Prim|Args],
 1065	chk_primitive_(Prim),
 1066	build_args_([Z|Args],Objs,_Types,Agenda,ObjAgenda),
 1067	newNode_(user(Prim),Objs,ObjAgenda,NewAgenda).
 1068
 1069build_args_([],[],_,Agenda,Agenda).
 1070build_args_([Arg|ArgsR],[Obj|Objs],[Type|Types],Agenda,NewAgenda) :-
 1071	(var(Type) -> Type=real ; true),                    % var if user defined
 1072	build_(Arg,Obj,Type,Agenda,NxtAgenda),
 1073	build_args_(ArgsR,Objs,Types,NxtAgenda,NewAgenda).
 1074
 1075chk_primitive_(Prim) :-  % wraps safe usage of unsafe current_predicate/2
 1076	UsrHead =..[Prim,'$op',_,_,_],
 1077	current_predicate(_,clpBNR:UsrHead).
 1078
 1079sandbox:safe_primitive(clpBNR:chk_primitive_(_Prim)).
 1080
 1081% to invoke user defined primitive
 1082call_user_primitive(Prim, P, InArgs, OutArgs) :-  % wraps unsafe meta call/N
 1083	call(clpBNR:Prim, '$op', InArgs, OutArgs, P).
 1084
 1085% really unsafe, but in a pengine a user can't define any predicates in another module, so this is safe
 1086sandbox:safe_meta(clpBNR:call_user_primitive(_Prim, _P, _InArgs, _OutArgs), []).
 1087
 1088% only called when argument(s) ground
 1089safe_(_ xor _) :- !,                              % clpBNR xor incompatible with `is` xor
 1090	fail.
 1091safe_(integer(_)) :- !,                           % clpBNR integer incompatible with `is` integer
 1092	fail.
 1093safe_(asin(_)) :- !,                              % asin multi-valued, can't use `is`
 1094	fail.
 1095safe_(acos(_)) :- !,                              % acos multi-valued, can't use `is`
 1096	fail.
 1097safe_(atan(_)) :- !,                              % atan multi-valued, can't use `is`
 1098	fail.
 1099safe_(_ ** E) :-                                  % ** multi-valued for rational even denominator
 1100	rational(E,_N,D),
 1101	0 is D mod 2, !,
 1102	fail.
 1103safe_(F) :- 
 1104	current_arithmetic_function(F),               % evaluable by is/2
 1105	F =.. [_Op|Args],
 1106	safe_args_(Args).
 1107
 1108safe_args_([]).
 1109safe_args_([A|Args]) :-
 1110	(atomic(A) -> \+ string(A) ; safe_(A)),       % numeric strings aren't safe
 1111	safe_args_(Args).
 1112
 1113%  a constraint must evaluate to a boolean 
 1114constraint_(C) :- nonvar(C), C =..[Op|_], fmap_(Op,_,_,_,[boolean|_]), !.
 1115
 1116%  map constraint operator to primitive/arity/types
 1117fmap_(+,    add,   ZXY,     ZXY,     [real,real,real]).
 1118fmap_(-,    add,   [Z,X,Y], [X,Z,Y], [real,real,real]).     % note subtract before minus 
 1119fmap_(*,    mul,   ZXY,     ZXY,     [real,real,real]).
 1120fmap_(/,    mul,   [Z,X,Y], [X,Z,Y], [real,real,real]).
 1121fmap_(**,   pow,   ZXY,     ZXY,     [real,real,real]).
 1122fmap_(^,    pow,   ZXY,     ZXY,     [real,real,real]).
 1123fmap_(min,  min,   ZXY,     ZXY,     [real,real,real]).
 1124fmap_(max,  max,   ZXY,     ZXY,     [real,real,real]).
 1125fmap_(==,   eq,    ZXY,     ZXY,     [boolean,real,real]).  % strict equality
 1126fmap_(=:=,  eq,    ZXY,     ZXY,     [boolean,real,real]).  % Prolog compatible, strict equality
 1127fmap_(is,   eq,    ZXY,     ZXY,     [boolean,real,real]).
 1128fmap_(<>,   ne,    ZXY,     ZXY,     [boolean,integer,integer]).
 1129fmap_(=\=,  ne,    ZXY,     ZXY,     [boolean,integer,integer]).  % Prolog compatible
 1130fmap_(=<,   le,    ZXY,     ZXY,     [boolean,real,real]).
 1131fmap_(>=,   le,    [Z,X,Y], [Z,Y,X], [boolean,real,real]).
 1132fmap_(<,    lt,    ZXY,     ZXY,     [boolean,real,real]).
 1133fmap_(>,    lt,    [Z,X,Y], [Z,Y,X], [boolean,real,real]).
 1134fmap_(<=,   in,    ZXY,     ZXY,     [boolean,real,real]).  % inclusion/subinterval
 1135
 1136fmap_(and,  and,   ZXY,     ZXY,     [boolean,boolean,boolean]).
 1137fmap_(',',  and,   ZXY,     ZXY,     [boolean,boolean,boolean]).  % for usability
 1138fmap_(or,   or,    ZXY,     ZXY,     [boolean,boolean,boolean]).
 1139fmap_(nand, nand,  ZXY,     ZXY,     [boolean,boolean,boolean]).
 1140fmap_(nor,  nor,   ZXY,     ZXY,     [boolean,boolean,boolean]).
 1141fmap_(xor,  xor,   ZXY,     ZXY,     [boolean,boolean,boolean]).
 1142fmap_(->,   imB,   ZXY,     ZXY,     [boolean,boolean,boolean]).
 1143
 1144fmap_(sqrt, sqrt,  ZX,      ZX,      [real,real]).          % pos. root version vs. **(1/2)
 1145fmap_(-,    minus, ZX,      ZX,      [real,real]).
 1146fmap_(~,    not,   ZX,      ZX,      [boolean,boolean]).
 1147fmap_(integer,int, ZX,      ZX,      [boolean,real]).
 1148fmap_(sgn,  sgn,   ZX,      ZX,      [integer,real]).       % signum, Z::integer(-1,1)
 1149fmap_(exp,  exp,   ZX,      ZX,      [real,real]).
 1150fmap_(log,  exp,   [Z,X],   [X,Z],   [real,real]).
 1151fmap_(abs,  abs,   ZX,      ZX,      [real,real]).
 1152fmap_(sin,  sin,   ZX,      ZX,      [real,real]).
 1153fmap_(asin, sin,   [Z,X],   [X,Z],   [real,real]).
 1154fmap_(cos,  cos,   ZX,      ZX,      [real,real]).
 1155fmap_(acos, cos,   [Z,X],   [X,Z],   [real,real]).
 1156fmap_(tan,  tan,   ZX,      ZX,      [real,real]).
 1157fmap_(atan, tan,   [Z,X],   [X,Z],   [real,real]).
 1158
 1159% reverse map node info to a facsimile of the original constraint
 1160map_constraint_op_(integral,$(V),integral(V)) :- !.
 1161map_constraint_op_(user(Func),Args,C) :- !,
 1162	remap_(Func,Args,C).
 1163map_constraint_op_(Op,Args,C) :-
 1164	fmap_(COp,Op,_,_,_),
 1165	remap_(COp,Args,C),
 1166	!.
 1167
 1168remap_(Op,$(Z,X,Y),C) :- constraint_(Op), Z==1, !,  % simplification for binary constraints
 1169	C=..[Op,X,Y]. 
 1170remap_(Op,$(Z,X),C) :- constraint_(Op), Z==1, !,    % simplification for unary constraints (~)
 1171	C=..[Op,X].
 1172remap_(Op,$(Z,X,Y),Z==C) :- !,
 1173	C=..[Op,X,Y].
 1174remap_(Op,$(Z,X),Z==C) :-
 1175	C=..[Op,X].
 1176
 1177%
 1178% Node constructor
 1179%
 1180newNode_(eq, [Z,X,Y], Agenda, Agenda) :- Z==1, !,   % optimize for X==Y, no node required
 1181	(number(X), number(Y) -> 0 is cmpr(X,Y) ; X=Y). % heavy lifting done by attr_unify_hook.
 1182newNode_(mul, [Z,X,Y], Agenda, Agenda) :-

 1183	(number(X), X =:= 0, interval_object(Y, _, Yv, _), finite_int(Yv), Z = X
 1184	;number(Y), Y =:= 0, interval_object(X, _, Xv, _), finite_int(Xv), Z = Y
 1185	), !

 1185.
 1186newNode_(add, [Z,X,Y], Agenda, NewAgenda) :-

 1187	(number(X), X =:= 0, newNode_(eq, [1,Z,Y], Agenda, NewAgenda)
 1188	;number(Y), Y =:= 0, newNode_(eq, [1,Z,X], Agenda, NewAgenda)
 1189	), !

 1189.
 1190newNode_(Op, Objs, Agenda, NewAgenda) :-
 1191	Args =.. [$|Objs],  % store arguments as $/N where N=1..3
 1192	NewNode = node(Op, _P, 0, Args),  % L=0
 1193	addNode_(Objs,NewNode),
 1194	% increment count of added nodes, will be decremented on backtracking/failure
 1195	g_incb('$clpBNR:node_count',1),
 1196	linkNode_(Agenda, NewNode, NewAgenda).
 1197
 1198addNode_([],_Node).
 1199addNode_([Arg|Args],Node) :-
 1200	(number(Arg) 
 1201	 -> true                                          % if Arg a number, nothing to do
 1202	 ;  interval_object(Arg, _Type, _Val, Nodelist),  % else add Node to Arg's Nodelist
 1203	    newmember(Nodelist, Node)
 1204	),
 1205	addNode_(Args,Node).
 1206
 1207sandbox:safe_global_variable('$clpBNR:node_count').
 1208
 1209clpStatistics :-
 1210	g_assign('$clpBNR:node_count',0),  % reset/initialize node count to 0
 1211	fail.  % backtrack to reset other statistics.
 1212
 1213clpStatistic(node_count(C)) :-
 1214	g_read('$clpBNR:node_count',C).
 1215
 1216% extend list with N
 1217newmember([X|Xs],N) :- 
 1218	(nonvar(X)
 1219	 -> newmember(Xs,N) % not end of (indefinite) list.
 1220     ;  X = N           % end of list, insert N and new tail Xs
 1221    ).
 1222
 1223%
 1224% Process Agenda to narrow intervals (fixed point iteration)
 1225%
 1226stable_([]/[]) :- !.  % small optimization when agenda is empty
 1227stable_(Agenda) :-
 1228	current_prolog_flag(clpBNR_iteration_limit,Ops),  % budget for current operation
 1229	stableLoop_(Agenda,Ops),
 1230	!.  % achieved stable state with empty Agenda -> commit.
 1231
 1232stableLoop_([]/[], OpsLeft) :- !,           % terminate successfully when agenda comes to an end
 1233	g_read('$clpBNR:iterations',Cur),        % maintain "low" water mark (can be negative)
 1234	(OpsLeft<Cur -> g_assign('$clpBNR:iterations',OpsLeft) ; true),
 1235	(OpsLeft<0 -> E is -OpsLeft, debug_clpBNR_('Iteration throttle limit exceeded by ~w ops.',E) ; true).
 1236stableLoop_([Node|Agenda]/T, OpsLeft) :-
 1237	Node = node(Op,P,_,Args),  % if node on queue ignore link bit (was: Node = node(Op,P,1,Args))
 1238	doNode_(Args, Op, P, OpsLeft, DoOver, Agenda/T, NxtAgenda),  % undoable on backtrack
 1239	nb_setarg(3,Node,0),                    % reset linked bit
 1240	% if doNode_ requested DoOver and above Ops threshold, put Node back at the end of the queue 
 1241	(atom(DoOver), OpsLeft > 0 -> linkNode_(NxtAgenda,Node,NewAgenda) ; NewAgenda = NxtAgenda),
 1242	RemainingOps is OpsLeft-1,              % decrement OpsLeft (can go negative)
 1243	stableLoop_(NewAgenda,RemainingOps).
 1244
 1245% support for max_iterations statistic
 1246sandbox:safe_global_variable('$clpBNR:iterations').
 1247
 1248clpStatistics :-
 1249	current_prolog_flag(clpBNR_iteration_limit,L), 
 1250	g_assign('$clpBNR:iterations',L),  % set "low" water mark to upper limit
 1251	fail.  % backtrack to reset other statistics.
 1252
 1253clpStatistic(max_iterations(O/L)) :-
 1254	g_read('$clpBNR:iterations',Ops),
 1255	current_prolog_flag(clpBNR_iteration_limit,L),
 1256	O is L-Ops.  % convert "low" water mark to high water mark
 1257
 1258%
 1259% doNode_/7 : Evaluate a node and add new nodes to end of queue. `evalNode` primitives can
 1260%	 fail, resulting in eventual failure of `stable_`, i.e., inconsistent constraint set.
 1261%
 1262% Note: primitives operate on interval values (sets) only, unaware of common arguments,
 1263% so additional consistency checks required on update.
 1264%
 1265doNode_($(ZArg,XArg,YArg), Op, P, OpsLeft, DoOver, Agenda, NewAgenda) :-  % Arity 3 Op
 1266	(var(P)                                          % check persistent bit
 1267	 -> getValue(ZArg,ZVal),
 1268	    getValue(XArg,XVal),
 1269	    getValue(YArg,YVal),
 1270	    evalNode(Op, P, $(ZVal,XVal,YVal), $(NZVal,NXVal,NYVal)),  % can fail
 1271	    % enforce consistency for common arguments by intersecting and redoing as required.
 1272	    (var(ZArg)                % if Z == X and/or Y
 1273	     -> (ZArg==XArg -> consistent_value_(NZVal,NXVal,NZ1,DoOver) ; NZ1 = NZVal),
 1274	        (ZArg==YArg -> consistent_value_(NZ1,  NYVal,NZ2,DoOver) ; NZ2 = NZ1),
 1275	        updateValue_(ZVal, NZ2, ZArg, OpsLeft, Agenda, AgendaZ)
 1276	     ;  AgendaZ = Agenda
 1277	    ),
 1278	    (var(XArg), XArg==YArg    % if X==Y
 1279	     -> consistent_value_(NXVal,NYVal,NVal,DoOver),
 1280	        updateValue_(XVal, NVal,  XArg, OpsLeft, AgendaZ, NewAgenda)  % only one update needed
 1281	     ;  updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, AgendaZX),
 1282	        updateValue_(YVal, NYVal, YArg, OpsLeft, AgendaZX, NewAgenda)
 1283	    )
 1284	 ; % P = p, trim persistent nodes from all arguments
 1285	    trim_op_(ZArg), trim_op_(XArg), trim_op_(YArg),  
 1286	    g_incb('$clpBNR:node_count',-1),  % decrement node count
 1287	    NewAgenda = Agenda
 1288	).
 1289
 1290doNode_($(ZArg,XArg), Op, P, OpsLeft, DoOver, Agenda, NewAgenda) :-       % Arity 2 Op
 1291	(var(P)                                          % check persistent bit
 1292	 -> getValue(ZArg,ZVal),
 1293	    getValue(XArg,XVal),
 1294	    evalNode(Op, P, $(ZVal,XVal), $(NZVal,NXVal)),      % can fail
 1295	    % enforce consistency for common arguments by intersecting and redoing as required.
 1296	    (var(ZArg), ZArg==XArg    % if Z==X
 1297	     -> consistent_value_(NZVal,NXVal,NVal,DoOver),     % consistent value, DoOver if needed
 1298	        updateValue_(ZVal, NVal,  ZArg, OpsLeft, Agenda, NewAgenda)  % only one update needed
 1299	     ;  updateValue_(ZVal, NZVal, ZArg, OpsLeft, Agenda, AgendaZ),
 1300	        updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, NewAgenda)
 1301	    )
 1302	 ; % P = p, trim persistent nodes from all arguments
 1303	    trim_op_(ZArg), trim_op_(XArg),
 1304	    g_incb('$clpBNR:node_count',-1),  % decrement node count
 1305	    NewAgenda = Agenda
 1306	).
 1307 
 1308doNode_($(Arg), Op, P, _OpsLeft, _, Agenda, NewAgenda) :-                 % Arity 1 Op
 1309	(var(P)                                          % check persistent bit
 1310	 -> getValue(Arg,Val),
 1311	    evalNode(Op, P, $(Val), $(NVal)),                   % can fail causing stable_ to fail => backtracking
 1312	    updateValue_(Val, NVal, Arg, 1, Agenda,NewAgenda)   % always update value regardless of OpsLeft limiter	
 1313	 ;  % P = p, trim persistent nodes from argument
 1314	    trim_op_(Arg),
 1315	    g_incb('$clpBNR:node_count',-1),  % decrement node count
 1316	    NewAgenda = Agenda
 1317	).
 1318
 1319consistent_value_(Val,Val,Val,_) :- !.                       % same value
 1320consistent_value_(Val1,Val2,Val,true) :- ^(Val1,Val2,Val).   % different values, intersect
 1321
 1322%
 1323% remove any persistent nodes from Arg
 1324%	called whenever a persistent node is encountered in FP iteration
 1325%
 1326trim_op_(Arg) :-
 1327	number(Arg)
 1328	 -> true                         % if a number, nothing to trim
 1329	 ;  get_attr(Arg, clpBNR, Def),  % an interval
 1330	    arg(3,Def,NList),            % Def = interval(_, _, NList, _),
 1331	    trim_persistent_(NList,TrimList),
 1332	    % if trimmed list empty, set to a new unshared var to avoid cycles(?) on backtracking
 1333	    (var(TrimList) -> setarg(3,Def,_) ; setarg(3,Def,TrimList)).  % write trimmed node list
 1334
 1335trim_persistent_(T,T) :- var(T), !.    % end of indefinite node list
 1336trim_persistent_([node(_,P,_,_)|Ns],TNs) :- nonvar(P), !, trim_persistent_(Ns,TNs).
 1337trim_persistent_([N|Ns],[N|TNs]) :- trim_persistent_(Ns,TNs).
 1338
 1339%
 1340% Any changes in interval values should come through here.
 1341% Note: This captures all updated state for undoing on backtracking
 1342%
 1343updateValue_(Old, New, Int, OpsLeft, Agenda, NewAgenda) :-  % set interval value to New
 1344	Old \== New,                                 % if value changes
 1345	(OpsLeft>0 ; propagate_if_(Old, New)), !,    % if OpsLeft >0 or narrowing sufficent
 1346	putValue_(New, Int, Nodelist),               % update value (may fail)
 1347	linkNodeList_(Nodelist, Agenda, NewAgenda).  % then propagate change
 1348updateValue_(_, _, _, _, Agenda, Agenda).        % otherwise just continue with Agenda
 1349
 1350% propgate if old range large enough and sufficient narrowing (> 10%)
 1351propagate_if_((OL,OH), (NL,NH)) :- 
 1352	OH-OL > 1E-32,           % large enough
 1353	(NH-NL)/(OH-OL) < 0.9.   % sufficient narrowing
 1354
 1355linkNodeList_(X,      List, List) :- var(X), !.  % end of indefinite list
 1356linkNodeList_([X|Xs], List, NewList) :-
 1357	(arg(3,X,Linked), Linked == 1                % test linked flag (only SWIP VM codes)
 1358	 -> linkNodeList_(Xs, List, NewList)         % add rest of the node list
 1359	 ;  linkNode_(List, X, NextList),            % not linked add it to list
 1360	    linkNodeList_(Xs, NextList, NewList)     % add rest of the node list
 1361	).
 1362
 1363linkNode_(List, X, NextList) :-                      % add to list
 1364	arg(1,X,Op), arg(4,X,Args),  % X = node(Op,_,_,Args)
 1365	(promote_(Args,Op)  % promote some primitives (put at start of work queue)
 1366	 -> linkNodeHead_(List, X, NextList)
 1367	  ; linkNodeTail_(List, X, NextList)
 1368	),
 1369	setarg(3,X,1).                                   % set linked bit
 1370
 1371% promotion heuristic which seems generally effective
 1372promote_($(_),_).                                       % promote low arity primitives
 1373promote_($(_,_),_).
 1374promote_($(Z,X,Y),add) :- \+ (var(Z), var(X), var(Y)).  % conditionally promote 
 1375promote_($(Z,X,Y),mul) :- \+ (var(Z), var(X), var(Y)).
 1376promote_($(Z,X,Y),pow) :- \+ (var(Z), var(X), var(Y)).
 1377
 1378linkNodeHead_(List/Tail, X, [X|List]/Tail).
 1379linkNodeTail_(List/[X|Tail], X, List/Tail).
 1380
 1381:- include(clpBNR/ia_utilities).  % print,solve, etc.

 watch(+X:interval_List, +Action:atom) is semidet

Succeeds if X is an interval and Action is an atom; otherwise fails. If successful, and Action is not none, a watchpoint is placed on X. Watchpoints are only "actioned" when the debug topic clpBNR is enabled. If Action = log a debug message is printed when the interval doman narrows. If Action = trace the debugger is invoked. If Action = none the watchpoint is removed. /

 1388watch(Int,Action) :-
 1389	atom(Action), 
 1390	current_module(clpBNR),  % this just marks watch/2 as unsafe regarding body
 1391	get_interval_flags_(Int,Flags), !,
 1392	remove_(Flags,watch(_),Flags1),
 1393	(Action = none -> Flags2=Flags1 ; Flags2=[watch(Action)|Flags1]),
 1394	set_interval_flags_(Int,Flags2).
 1395watch(Ints,Action) :-
 1396	list(Ints),
 1397	watch_list_(Ints,Action).
 1398
 1399remove_([],_,[]).
 1400remove_([X|Xs],X,Xs) :- !.
 1401remove_([X|Xs],X,[X|Ys]) :-
 1402	remove_(Xs,X,Ys).
 1403
 1404watch_list_([],_Action).
 1405watch_list_([Int|Ints],Action) :-
 1406	watch(Int,Action),
 1407	watch_list_(Ints,Action).
 1408
 1409% check if watch enabled on this interval
 1410check_monitor_(Int, Update, interval(_Type,_Val,_Nodelist,Flags)) :-
 1411	(memberchk(watch(Action), Flags)
 1412	 -> once(monitor_action_(Action, Update, Int))
 1413	 ; true
 1414	).
 1415
 1416%
 1417% for monitoring changes - all actions defined here
 1418%
 1419monitor_action_(trace, Update, Int) :-  !, % log changes to console and enter debugger
 1420	monitor_action_(log, Update, Int),
 1421	trace.
 1422monitor_action_(log, Update, Int) :-  var(Update), !,  % log interval unify
 1423	debug_clpBNR_('Unify ~p with ~p',[Int,Update]).
 1424monitor_action_(log, Update, Int) :-  number(Update), !,    % log interval unify with number
 1425	domain(Int,Dom),
 1426	debug_clpBNR_('Unify _?{~p} with ~p',[Dom,Update]).
 1427monitor_action_(log, integer, Int) :-  !,  % change type to integer
 1428	debug_clpBNR_('Set type of  ~p to ~p',[Int,integer]).
 1429monitor_action_(log, Val, Int) :-  !,  % narrow range
 1430	debug_clpBNR_('Set value of ~p to (~p)',[Int,Val]).
 1431monitor_action_(_, _, _).  % default to noop (as in 'none')
 1432
 1433sandbox:safe_primitive(clpBNR:watch(_Int,Action)) :- % watch(X,trace) is unsafe.
 1434	Action \= trace.
 1435% only safe because watch(X,trace) is unsafe.
 1436sandbox:safe_primitive(clpBNR:monitor_action_(_Action, _Update, _Int)).