1/*	The MIT License (MIT)
    2 *
    3 *	Copyright (c) 2019-2026 Rick Workman
    4 *
    5 *	Permission is hereby granted, free of charge, to any person obtaining a copy
    6 *	of this software and associated documentation files (the "Software"), to deal
    7 *	in the Software without restriction, including without limitation the rights
    8 *	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
    9 *	copies of the Software, and to permit persons to whom the Software is
   10 *	furnished to do so, subject to the following conditions:
   11 *
   12 *	The above copyright notice and this permission notice shall be included in all
   13 *	copies or substantial portions of the Software.
   14 *
   15 *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   16 *	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   17 *	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   18 *	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   19 *	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
   20 *	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   21 *	SOFTWARE.
   22 */
   23
   24%
   25% API function for executing primitives
   26%
   27% evalNode(+primitive_name, ?persistence_flag, +$(inputs..), -$(outputs..))
   28%
   29evalNode(Op, P, Is, R) :-
   30	g_inc('$clpBNR:evalNode'),  % count of primitive calls
   31	narrowing_op(Op, P, Is, R),
   32	!.
   33evalNode(_Op, _, _, _):-
   34	g_inc('$clpBNR:evalNodeFail'),  % count of primitive call failures
   35	debugging(clpBNR),              % fail immediately unless debug=true
   36	current_node_(C),
   37	debug_clpBNR_("** fail ** ~w.",C),
   38	fail.
   39
   40sandbox:safe_global_variable('$clpBNR:evalNode').
   41sandbox:safe_global_variable('$clpBNR:evalNodeFail').
   42
   43clpStatistics :-
   44	g_assign('$clpBNR:evalNode',0),
   45	g_assign('$clpBNR:evalNodeFail',0),
   46	fail.  % backtrack to reset other statistics.
   47
   48clpStatistic(narrowingOps(C)) :- g_read('$clpBNR:evalNode',C).
   49
   50clpStatistic(narrowingFails(C)) :- g_read('$clpBNR:evalNodeFail',C).
   51
   52%
   53% Ops histogram data collection
   54%
   55clpStatistics :- 
   56	% initialize Ops count dictionary
   57	findall(Op-0,(clause(clpBNR:narrowing_op(Op,_,_,_),_),atom(Op)),Data),
   58	dict_create(Dict,'$clpBNR:ops_dict',Data),
   59	g_assign('$clpBNR:ops_dict',Dict),
   60	fail.  % backtrack to reset other statistics.
   61
   62clpStatistic(opCounts(Cs)) :-   % retrieve counts from dictionary
   63	g_read('$clpBNR:ops_dict',Dict),
   64	dict_pairs(Dict,'$clpBNR:ops_dict',Raw),
   65	exclude_zeros(Raw,Cs).       % exclude Ops with no calls
   66
   67exclude_zeros([],[]).
   68exclude_zeros([_-0|Raw],Cs)   :- !, exclude_zeros(Raw,Cs).
   69exclude_zeros([C|Raw],[C|Cs]) :- exclude_zeros(Raw,Cs).
   70
   71% use library(prolog_wrap) to conditionally gather Op counts statistics	
   72counts_clpBNR(Bool)  :-                  % query or already in defined state
   73	( current_predicate_wrapper(clpBNR:evalNode(_Op, _P, _Is, _R), 
   74	                            'clpBNR:evalNode', _Wrapped, _Body)
   75	 -> Bool = true ; Bool = false
   76	),
   77	!.
   78counts_clpBNR(true)  :-                  % turn on wrapper
   79	wrap_predicate(clpBNR:evalNode(Op, _P, _Is, _R),
   80	                   'clpBNR:evalNode', 
   81	                   Wrapped, evalNode_wrap_(Wrapped, Op)).
   82counts_clpBNR(false) :-                  % turn off wrapper
   83	unwrap_predicate(clpBNR:evalNode/4,  %(Op, P, Is, R),
   84	                   'clpBNR:evalNode').
   85
   86evalNode_wrap_(Wrapped, Op) :-
   87	 % increment count for Op entry in Dict 
   88	g_read('$clpBNR:ops_dict',Dict),
   89	get_dict(Op,Dict,N), N1 is N+1, nb_link_dict(Op,Dict,N1),
   90	Wrapped.                             % execute wrapped evalNode
   91
   92sandbox:safe_global_variable('$clpBNR:ops_dict').
   93
   94%
   95% non-empty interval test
   96%
   97goal_expansion(non_empty(L,H), 1 =\= cmpr(L,H)).  % SWIP inline optimization for non_empty/2
   98
   99non_empty((L,H)) :- non_empty(L,H).
  100non_empty(L,H) :- 1 =\= cmpr(L,H).
  101
  102%
  103% interval category: non-negative (includes zero), non-positive, or split (contains 0)
  104%
  105intCase(C, (L,H)) :-
  106	( L>=0 -> C=p            % non-negative
  107	; H=<0 -> C=n            % non-positive
  108	;         C=s            % split
  109	).
  110
  111%
  112% Z := X ^ Y  (interval intersection)
  113%
  114^((Xl,Xh), (Yl,Yh), (Zl,Zh)) :-
  115	Zl is maxr(Xl, Yl), Zh is minr(Xh, Yh),  % eliminates NaN's, prefers rationals
  116	non_empty(Zl,Zh).
  117
  118%
  119% Z := X \/ Y (interval union)
  120%
  121union_((Xl,Xh),(Yl,Yh),(Zl,Zh)) :-
  122	Zl is minr(Xl,Yl), Zh is maxr(Xh,Yh).  % eliminates NaN's, prefers rationals
  123
  124%
  125% (Xl,Xh) contains (Yl,Yh)
  126%
  127int_contains((Xl,Xh),(Yl,Yh)) :-
  128	Xl is minr(Xl,Yl), Xh is maxr(Xh,Yh).
  129
  130%
  131% to_float/2 converts an non-float bounds to floats
  132%  useful when when float arguments will be produced on narrowing, e.g., elementary functions
  133%  avoids incorrect rounding on argument conversion (monotonic decreasing segments)
  134%
  135to_float((L,H),(FL,FH)) :-
  136	( float(L) -> FL = L ; FL is roundtoward(float(L), to_negative) ), 
  137	( float(H) -> FH = H ; FH is roundtoward(float(H), to_positive) ). 
  138
  139%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  140%
  141% Relational Operations (narrowing operations)
  142%
  143:- discontiguous clpBNR:narrowing_op/4.  144:- style_check(-singleton).  % for narrowing_op
  145%
  146%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  147
  148% integral (contract to integer bounds) - optimized version of (int,P,$(1,X),_)
  149
  150narrowing_op(integral, P, $(X), $(NewX)) :-
  151	integral_(X,NewX).
  152	
  153integral_((Xl,Xh),(NXl,NXh)) :-
  154	NXl is ceiling(Xl), NXh is floor(Xh),
  155	non_empty(NXl,NXh).
  156
  157%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  158
  159% Z == integer(X) % (Z boolean) 
  160%  true if X is an integral point, false if it doesn't contain an integer
  161%  else if true, narrows X to integer bounds
  162
  163narrowing_op(int, P, $(Z,X), $(NewZ,NewX)) :-
  164	int_(Z,X,P,NewZ,NewX).
  165	
  166int_(Z,X,p,NewZ,X) :-                       % persistent cases
  167	X = (Xl,Xh),
  168	(0 is cmpr(Xl,Xh)                   % point value?
  169	 -> (0 is cmpr(Xl,integer(Xl)) -> B = 1 ; B = 0)  % yes, true if integer else false
  170	 ;  ceiling(Xl) > floor(Xh),        % no, and X doesn't contain an integer
  171	    B = 0 
  172	),
  173	!,  % succeed, persistent
  174	^(Z,(B,B),NewZ).
  175	
  176int_(Z,X,P,NewZ,NewX) :-                % possible narrowing case
  177	(Z = (1,1) 
  178	 -> integral_(X,NewX),              % Z true, narrow X to integers
  179	    (NewX = (N,N) -> P=p ; true),   % persistent if a point
  180	    NewZ = Z
  181	 ;  NewX = X,                       % Z false or indeterminate, no change to X
  182	    ^(Z,(0,1),NewZ)                 % Z boolean
  183	).
  184
  185%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  189narrowing_op(sgn, P, $(Z,X), $((NZl,NZh),NewX)) :-
  190	X = (Xl,Xh),
  191	SXl is cmpr(Xl,0), SXh is cmpr(Xh,0),
  192	^(Z,(SXl,SXh),(NZl,NZh)),
  193	(NZl >= 0 -> ^(X, (0,1.0Inf), X1)   ; X1=X),
  194	(NZh =< 0 -> ^(X1,(-1.0Inf,0),NewX) ; NewX = X1),
  195	(NZl == NZh -> P=p ; true).  % sgn is a point
  196
  197%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  198
  199% Z==(X==Y)  % (Z boolean)
  200
  201narrowing_op(eq, P, $(Z,X,Y), Out) :-
  202	(^(X,Y,New)                       % intersect?
  203	 -> (Z = (1,1)
  204	     -> Out = $(Z,New,New), 
  205	        New = (L,H),
  206	        (0 is cmpr(L,H) -> P = p ; true)   % persistent if a point
  207	     ; % intersect and Z false or indeterminate, all things possible
  208	        ( X = (Xl,Xh), Y = (Yl,Yh),
  209	         0 is cmpr(Xl,Xh) \/ cmpr(Yl,Yh) \/ cmpr(Xl,Yl)  % X and Y same point
  210	         -> ^(Z,(1,1),NewZ),      % Z true
  211	            P = p                 % persistant
  212	         ;  ^(Z,(0,1),NewZ)       % ensure Z boolean
  213	        ),
  214	        Out = $(NewZ,X,Y)
  215	    )
  216	 ;  % X and Y do not intersect  
  217	    ^(Z, (0,0), NewZ),            % fail primitive if Z cannot be false
  218	    Out = $(NewZ,X,Y), P = p      % persistent
  219	 ).
  220
  221%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  222
  223% Z==(X<>Y)  % (Z boolean, X,Y integer)
  224
  225narrowing_op(ne, P, $(Z,X,Y), Out) :-
  226	(^(X,Y,(L,H))
  227	 -> (Z = (1,1)
  228	     -> ne_int_(X,Y,NewX),        % narrow X if Y is a point and a bound of X (always succeeds but may not narrow
  229	        ne_int_(Y,NewX,NewY),     % narrow Y if NewX is a point and a bound of Y (always succeeds but may not narrow)	
  230	        Out = $(Z,NewX,NewY)
  231	     ;  (X = Y, L = H -> B = (0,0), P = p ;  B = (0,1)),
  232	        ^(Z,B,NewZ),
  233	        Out = $(NewZ,X,Y)
  234	     )
  235	 ;  % X and Y do not intersect
  236	    ^(Z, (1,1), NewZ),            % true and 
  237	    P = p,                        % persistent
  238	    Out = $(NewZ,X,Y)
  239	).
  240
  241% narrow X and Y given Y <> X,  where Y and X are integer intervals (enforced elsewhere)
  242% integers so standard comparison works.
  243ne_int_((X,H), (X,X), (NewL,H)) :- !,  % X is a point,  and low bound of Y
  244	NewL is X+1, NewL=<H.
  245ne_int_((L,X), (X,X), (L,NewH)) :- !,  % X is a point,  and high bound of Y
  246	NewH is X-1, L=<NewH.
  247ne_int_(Y, _X, Y).                      % no narrowing
  248
  249%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  250
  251% Z==(X=<Y)  % (Z boolean)
  252
  253narrowing_op(le, P, $(Z, X, Y), Out) :-
  254	X = (Xl,Xh), Y = (Yl,Yh),
  255	(non_empty(Xh,Yl)                                  % Xh =< Yl
  256	 -> ^(Z,(1,1),NewZ), P = p, Out = $(NewZ,X,Y)      % persistent X =< Y
  257	 ;  (-1 is cmpr(Yh,Xl)                             % Yh < Xl
  258	     -> ^(Z,(0,0),NewZ), P = p, Out = $(NewZ,X,Y)  % persistent Y < X
  259	     ;  le_int_(Z,Xl,Xh,Yl,Yh,P,Out)               % possibly non-persistent narrowing cases
  260	    )
  261	).
  262	
  263le_int_((1,1),Xl,Xh,Yl,Yh,P,$((1,1),(Xl,NXh),(NYl,Yh))) :- !,  % common case, narrow X and Y
  264	NXh is minr(Xh,Yh),                               % NewX := (Xl,NXh) 
  265	NYl is maxr(Xl,Yl),                               % NewY := (NYl,Yh)
  266	(non_empty(NXh,NYl) -> P=p ; true).               % now persistant?
  267le_int_(Z,Xl,Xh,Yl,Yh,P,$(NewZ,NewX,NewY)) :- !,
  268	(Z = (0,0)
  269	 -> le_int_((1,1),Yl,Yh,Xl,Xh,P,$(_,NewY,NewX)),  % false so Y =< X (Y < X is unsafe)
  270	 	NewZ = Z
  271	 ;  NewX = (Xl,Xh), NewY = (Yl,Yh),               % indeterminate, X and Y unchanged
  272	   ^(Z,(0,1),NewZ)                                % Z::boolean
  273	).
  274
  275%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  276
  277% Z==(X<Y)  % (Z boolean) generally unsound on reals but possibly useful for minima/maxima
  278
  279narrowing_op(lt, P, $(Z, X, Y), Out) :-
  280	X = (Xl,Xh), Y = (Yl,Yh),
  281	(-1 is cmpr(Xh,Yl)                                 % Xh < Yl
  282	 -> ^(Z,(1,1),NewZ), P = p, Out = $(NewZ,X,Y)      % persistent X < Y
  283	 ;  (1 =\= cmpr(Yh,Xl)                             % Yh =< Xl
  284	     -> ^(Z,(0,0),NewZ), P = p, Out = $(NewZ,X,Y)  % persistent Y =< X	
  285	     ;  lt_int_(Z,Xl,Xh,Yl,Yh,P,Out)               % possibly non-persistent narrowing cases
  286	    )
  287	).
  288
  289lt_int_((1,1),Xl,Xh,Yl,Yh,P,$((1,1),(Xl,NXh),(NYl,Yh))) :- !,  % common case, narrow X and Y
  290	next_lt_(Yh,-1,YhD),                              % YhD is next downward value from Yh
  291	NXh is minr(Xh,YhD),                              % NewX := (Xl,NXh) 
  292	next_lt_(Xl,1,XlU),                               % NXlU is next upward value from NXl
  293	NYl is maxr(XlU,Yl),                              % NewY := (NYl,Yh)
  294	(non_empty(NXh,NYl) -> P=p ; true).               % now persistent?
  295lt_int_(Z,Xl,Xh,Yl,Yh,P,$(NewZ,NewX,NewY)) :- !,
  296	(Z = (0,0)
  297	 -> lt_int_((1,1),Yl,Yh,Xl,Xh,P,$(_,NewY,NewX)),  % false so Y < X
  298	    NewZ = Z
  299	 ;  NewX = (Xl,Xh), NewY = (Yl,Yh),               % indeterminate, X and Y unchanged
  300	   ^(Z,(0,1),NewZ)                                % Z::boolean
  301	).
  302
  303% Note: can't narrow an infinite bound, minimize change to bound
  304% Repeat, not sound on reals (uses nexttoward, missing values between floats)
  305next_lt_( 1.0Inf, _,  1.0Inf) :- !.
  306next_lt_(-1.0Inf, _, -1.0Inf) :- !.
  307next_lt_(V, -1, NV) :- NV is maxr(V-1,nexttoward(V,-1.0Inf)).  % integers will get sorted out with `integral`
  308next_lt_(V,  1, NV) :- NV is minr(V+1,nexttoward(V, 1.0Inf)).
  309
  310%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  311
  312% Z==(X<=Y)  % inclusion, constrains X to be subinterval of Y (Z boolean)
  313
  314narrowing_op(in, P, $(Z, X, Y), $(NewZ, NewX, Y)):-
  315	% Only two cases: either X and Y intersect or they don't.
  316	(^(X,Y,X1)               % X1 is intersection of X and Y
  317	 -> (Z = (1,1) 
  318	     -> NewX = X1        % Z true, X must be subinterval of Y
  319	     ;  NewX = X,        % Z is false or indeterminate, X unchanged
  320	        ^(Z,(0,1),NewZ)  % Z boolean
  321	    )
  322	 ;  ^(Z,(0,0),NewZ),     % X and Y don't intersect, Z must be false
  323	    NewX = X,            % no change in X
  324	    P=p                  % persistent, since X and Y can never intersect
  325	).
  326
  327%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  328
  329% Z==X+Y
  330
  331narrowing_op(add, _, $((Zl,Zh), (Xl,Xh), (Yl,Yh)), $((NZl,NZh), (NXl,NXh), (NYl,NYh))) :-
  332	NZl is maxr(Zl, roundtoward( Xl+Yl,  to_negative)),       % NewZ := Z ^ (X+Y),
  333	NZh is minr(Zh, roundtoward( Xh+Yh,  to_positive)),
  334	non_empty(NZl,NZh),
  335	% Note: subtraction done by adding minus values so rounding mode consistent
  336	% during any numeric type conversion.
  337	NXl is maxr(Xl, roundtoward(NZl+(-Yh),  to_negative)),    % NewX := X ^ (NZ-Y),
  338	NXh is minr(Xh, roundtoward(NZh+(-Yl),  to_positive)),
  339	non_empty(NXl,NXh),
  340	NYl is maxr(Yl, roundtoward(NZl+(-NXh), to_negative)),    % NewY := Y ^ (NZ-NX).
  341	NYh is minr(Yh, roundtoward(NZh+(-NXl), to_positive)),
  342	non_empty(NYl,NYh).
  343
  344%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  345
  346% Z==X*Y
  347
  348narrowing_op(mul, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  349	intCase(Cx,X),
  350	intCase(Cy,Y),
  351	multCase(Cx,Cy,X,Y,Z1), ^(Z,Z1,NewZ),           % NewZ := Z ^ (X*Y)
  352	% Uses totality condition from 
  353	%   https://www.academia.edu/52359860/A_unified_framework_for_interval_constraints_and_interval_arithmetic
  354	% to avoid unnecessary evaluation of NewX and NewY
  355	(int_contains(Z,Z1)

  356	 -> NewX = X, NewY = Y
  357	 ; 
  358	    intCase(CNz,NewZ),
  359	    odivCase(CNz,Cx,NewZ,X,Y,Yp),               % Yp := Y ^ (Z/X),
  360	    intCase(Cyp,Yp),
  361	    odivCase(CNz,Cyp,NewZ,Yp,X,NewX),           % NewX := X ^ (Z/Yp),
  362	    % if X narrowed it may be necessary to recalculate Y due to non-optimal ordering.
  363	    (Y == Yp, X \== NewX                        % if Y didn't narrow and X did
  364	     -> intCase(CNx,NewX),
  365	        odivCase(CNz,CNx,NewZ,NewX,Y,NewY)      % re calculate: NewY := Y ^ NewZ/NewX
  366	     ;  NewY = Yp,
  367	        ( ((zero_int(NewX), finite_int(NewY)) ; (zero_int(NewY), finite_int(NewX)))
  368	         -> P=p  % if X or Y = 0, and other finite, persistent 
  369	         ;  true
  370	        )
  371	    )
  372	)

  372.
  373/*  general case, but sub-par performance since some re-calculations unnecessary
  374narrowing_op(mul, _, $(Z,X,Y), New) :-
  375	...
  376	Nxt = $(NewZ,NewX,Yp),
  377	($(Z,X,Y) == Nxt                   
  378	 -> New = Nxt                       % no change, exit
  379	 ;  narrowing_op(mul, _, Nxt, New)  % something narrowed, repeat until no change
  380	). 
  381*/
  382
  383%
  384% * cases ("Interval Arithmetic: from Principles to Implementation", Fig. 3)
  385%
  386% Note: Use negation to ensure proper rounding control during any conversions.
  387% This is quite tricky because negation inverts direction of rounding - take care.
  388%multCase(z,_, X, _, _, X) :- !.  % X==0
  389%multCase(_,z, _, Y, _, Y).       % Y==0
  390multCase(p,p, (Xl,Xh), (Yl,Yh), (NZl,NZh)):- 
  391	NZl is roundtoward(Xl*Yl,to_negative),
  392	NZh is roundtoward(Xh*Yh,to_positive).
  393multCase(p,n, (Xl,Xh), (Yl,Yh), (NZl,NZh)):- 
  394	NZl is -roundtoward(Xh * -Yl,to_positive),
  395	NZh is -roundtoward(Xl * -Yh,to_negative).
  396multCase(n,p, (Xl,Xh), (Yl,Yh), (NZl,NZh)):- 
  397	NZl is -roundtoward(-Xl * Yh,to_positive),
  398	NZh is -roundtoward(-Xh * Yl,to_negative).
  399multCase(n,n, (Xl,Xh), (Yl,Yh), (NZl,NZh)):-
  400	NZl is roundtoward(-Xh * -Yh,to_negative),  % negate for any conversion rounding 
  401	NZh is roundtoward(-Xl * -Yl,to_positive).
  402
  403multCase(p,s, (Xl,Xh), (Yl,Yh), (NZl,NZh)):-
  404	NZl is -roundtoward(Xh * -Yl,to_positive),
  405	NZh is roundtoward(Xh *  Yh,to_positive).
  406multCase(n,s, (Xl,Xh), (Yl,Yh), (NZl,NZh)):-
  407	NZl is -roundtoward(-Xl *  Yh,to_positive),
  408	NZh is roundtoward(-Xl * -Yl,to_positive).
  409multCase(s,p, (Xl,Xh), (Yl,Yh), (NZl,NZh)):-
  410	NZl is -roundtoward(-Xl * Yh,to_positive),
  411	NZh is roundtoward( Xh * Yh,to_positive).
  412multCase(s,n, (Xl,Xh), (Yl,Yh), (NZl,NZh)):-
  413	NZl is -roundtoward( Xh * -Yl,to_positive),
  414	NZh is roundtoward(-Xl * -Yl,to_positive).
  415
  416multCase(s,s, (Xl,Xh), (Yl,Yh), (NZl,NZh)):-
  417	NZl is minr(-roundtoward( Xl * -Yh,to_positive),-roundtoward(-Xh * Yl,to_positive)),
  418	NZh is maxr( roundtoward(-Xl * -Yl,to_positive), roundtoward( Xh * Yh,to_positive)).
  419
  420%
  421% / cases ("Interval Arithmetic: from Principles to Implementation", Fig. 4)
  422% Tricky handling the "zero" cases - returns universal interval when no narowing possible:
  423%	1. denominator is zero or contains zero (z and s)
  424%	2. denominator is bounded by zero (p and n) and numerator contains zero (see numeric guards)
  425% Numeric guards check for 0/0, e.g., p,p -> Xh+Yh>0
  426
  427odivCase(p,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  428	sign(Yl)+sign(Xl) > 0  % X > 0 or Y > 0
  429	 -> NZl is maxr(Zl,roundtoward(Xl/Yh,to_negative)),
  430	    NZh is minr(Zh,roundtoward(Xh/max(0.0,Yl),to_positive)),
  431	    non_empty(NZl,NZh)
  432	 ;  NZl = Zl, NZh = Zh.
  433odivCase(p,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  434	sign(Xl)-sign(Yh) > 0  % X > 0 or Y < 0
  435	 -> NZl is maxr(Zl,roundtoward(-Xh / -min(-0.0,Yh),to_negative)),  % use min to preserve sign
  436	    NZh is minr(Zh,roundtoward(-Xl / -Yl          ,to_positive)),
  437	    non_empty(NZl,NZh)
  438	 ;  NZl = Zl, NZh = Zh.
  439odivCase(p,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  440	Xl > 0     % X > 0
  441	 -> ZIh is roundtoward(Xl/Yh,to_negative),  % (ZIh,inf)
  442	    ZIl is roundtoward(Xl/Yl,to_positive),  % (-inf,ZIl)
  443	    ( 1 is cmpr(Zl,ZIl) -> NZl is maxr(Zl,ZIh) ; NZl = Zl),  % Zl > ZIl ?
  444	    (-1 is cmpr(Zh,ZIh) -> NZh is minr(Zh,ZIl) ; NZh = Zh),  % Zh > ZIl ?
  445	    non_empty(NZl,NZh)
  446	 ;  NZl = Zl, NZh = Zh.
  447
  448odivCase(n,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
  449	sign(Yl)-sign(Xh) > 0  % X < 0 or Y > 0
  450	 -> NZl is maxr(Zl,roundtoward(-Xl / -max(0.0,Yl),to_negative)),  % if Yx=0, preserve sign
  451	    NZh is minr(Zh,roundtoward(-Xh / -max(0.0,Yh),to_positive)),
  452	    non_empty(NZl,NZh)
  453	 ;  NZl = Zl, NZh = Zh.
  454odivCase(n,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  455	sign(Yh)+sign(Xh) < 0  % X < 0 or Y < 0
  456	 -> NZl is maxr(Zl,roundtoward(Xh/Yl,to_negative)),
  457	    NZh is minr(Zh,roundtoward(Xl/min(-0.0,Yh),to_positive)),  % use min to preserve sign
  458	    non_empty(NZl,NZh)
  459	 ;  NZl = Zl, NZh = Zh.
  460odivCase(n,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  461	Xh < 0     % X < 0
  462	 -> ZIh is roundtoward( Xh /  Yl,to_negative),  % (ZIh,inf)
  463	    ZIl is roundtoward(-Xh / -Yh,to_positive),  % (-inf,ZIl)
  464	    ( 1 is cmpr(Zl,ZIh) -> NZl is maxr(Zl,ZIl) ; NZl = Zl),  % Zl > ZIh ?
  465	    (-1 is cmpr(Zh,ZIl) -> NZh is minr(Zh,ZIh) ; NZh = Zh),  % Zh > ZIl ?
  466	    non_empty(NZl,NZh)
  467	 ;  NZl = Zl, NZh = Zh.
  468
  469odivCase(s,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  470	Yl > 0     % Y > 0
  471	 -> NZl is maxr(Zl,roundtoward(-Xl / -Yl,to_negative)),
  472	    NZh is minr(Zh,roundtoward( Xh /  Yl,to_positive)),
  473	    non_empty(NZl,NZh)
  474	 ;  NZl = Zl, NZh = Zh.
  475odivCase(s,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  476	Yh < 0     % Y < 0
  477	 -> NZl is maxr(Zl,roundtoward(Xh/Yh,to_negative)),
  478	    NZh is minr(Zh,roundtoward(Xl/Yh,to_positive)),
  479	    non_empty(NZl,NZh)
  480	 ;  NZl = Zl, NZh = Zh.
  481
  482odivCase(s,s,       _,       _,       Z,         Z).  % both contain 0, no narrowing
  483
  484%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  485
  486% Z==min(X,Y)    Z==max(X,Y)
  487
  488narrowing_op(min, _, $((Zl,Zh),(Xl,Xh),(Yl,Yh)), New) :-
  489	Z1l is maxr(Zl,minr(Xl,Yl)),  % Z1 := Z ^ min(X,Y),
  490	Z1h is minr(Zh,minr(Xh,Yh)),
  491	minimax((Zl,1.0Inf), (Z1l,Z1h),(Xl,Xh),(Yl,Yh), New).
  492
  493narrowing_op(max, _, $((Zl,Zh),(Xl,Xh),(Yl,Yh)), New) :-
  494	Z1l is maxr(Zl,maxr(Xl,Yl)),  % Z1 := Z ^ max(X,Y),
  495	Z1h is minr(Zh,maxr(Xh,Yh)),
  496	minimax((-1.0Inf,Zh), (Z1l,Z1h),(Xl,Xh),(Yl,Yh), New).
  497
  498minimax(_, Z1,X,Y, NewXYZ) :-                  % Case 1, X not in Z1
  499	\+(^(Z1,X,_)), !,                          % _ := Z1 \^ X,
  500	NewXYZ = $(New, X, New),
  501	^(Y,Z1,New).                               % New := Y ^ Z1.
  502
  503minimax(_, Z1,X,Y, NewXYZ) :-                  % Case 2, Y not in Z1
  504	\+(^(Z1,Y,_)), !,                          % _ := Z1 \^ Y,
  505	NewXYZ = $(New, New, Y),
  506	^(X,Z1,New).                               % New := X ^ Z1.
  507
  508minimax(Zpartial, Z1,X,Y, $(Z1, NewX, NewY)) :- % Case 3, overlap
  509	^(X,Zpartial,NewX),                        % NewX := X ^ Zpartial,
  510	^(Y,Zpartial,NewY).                        % NewY := Y ^ Zpartial.
  511
  512%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  513
  514% Z==abs(X)
  515
  516narrowing_op(abs, _, $(Z,X), $(NewZ, NewX)) :-
  517	intCase(Cx,X),
  518	absCase(Cx,X,Z,NewZ), !,
  519	non_empty(NewZ),
  520	intersect_regions_(X,NewZ,NewX),
  521	non_empty(NewX).
  522
  523%
  524% intersect and join Z with positive and negative regions of X
  525%
  526intersect_regions_(Z,X,NewZ) :-
  527	(^(Z,X,ZPos) -> true ; empty_interval(ZPos)),  % positive Z region
  528	(-(X,Z,ZNeg) -> true ; empty_interval(ZNeg)),  % negative Z region
  529	union_(ZPos,ZNeg,NewZ).
  530
  531% NZ := Z ^ -X (unary minus)
  532-((Xl,Xh), (Zl,Zh), (NZl,NZh)) :-
  533	NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl),
  534	non_empty(NZl,NZh).
  535
  536%
  537% abs(X) cases
  538%
  539absCase(p, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,Xl),  NZh is minr(Zh,Xh).
  540absCase(n, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl).
  541absCase(s, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,0),   NZh is minr(Zh,maxr(-Xl,Xh)).
  542
  543%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  544
  545% Z== -X
  546
  547narrowing_op(minus, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
  548	NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl),    % NewZ := Z ^ -X, no empty check required
  549	NXl is maxr(Xl,-NZh), NXh is minr(Xh,-NZl).  % NewX := X ^ -Z, no empty check required
  550
  551%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  552
  553% Z== sqrt(X)  (Z is positive root of positive X)
  554
  555narrowing_op(sqrt, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
  556	Xh>=0, Xl1 is maxr(0,Xl),  % narrow X to positive range
  557	NZl is maxr(maxr(0,Zl),roundtoward(sqrt(Xl1),to_negative)),
  558	NZh is minr(Zh,roundtoward(sqrt(Xh),to_positive)),
  559	non_empty(NZl,NZh),
  560	NXh is minr(Xh,roundtoward(NZh*NZh,to_positive)),
  561	NXl is maxr(Xl,-NXh),
  562	non_empty(NXl,NXh).
  563
  564%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  565
  566% Z== exp(X)
  567
  568narrowing_op(exp, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
  569	NZl is maxr(Zl,roundtoward(exp(Xl),to_negative)),
  570	NZh is minr(Zh,roundtoward(exp(Xh),to_positive)),
  571	non_empty(NZl,NZh),
  572	NXl is maxr(Xl,roundtoward(log(NZl),to_negative)),
  573	NXh is minr(Xh,roundtoward(log(NZh),to_positive)),
  574	non_empty(NXl,NXh).
  575
  576%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  577
  578% Z== X**Y
  579
  580narrowing_op(pow, P, $(Z,X,Y), $(NewZ,NewX,NewY)) :-
  581	powr_(P, Z,X,Y, NewZ,NewX,NewY).
  582
  583powr_(p, Z,X,Y, NewZ,NewX,NewY) :-                  % (persistant) special cases
  584	Z = (Zl,Zh), X = (Xl,Xh), Y = (Yl,Yh), 
  585	( equals_one(Xl), equals_one(Xh) ->             % 1**Y == Z
  586		!, ^(Z,(1,1),NewZ), NewX = X, NewY = Y      %  ==> Z := (1,1)
  587	; equals_one(Zl), equals_one(Zh) ->             % X**Y == 1 , X doesn't contain 1 or -1
  588		((Xl =< -1,-1 =< Xh ; Xl =< 1,1 =< Xh)
  589		 -> fail                                    % fail to use general pt_powrCase
  590		 ;  !, NewZ = Z, NewX = X, ^(Y,(0,0),NewY)  %  ==> Y := (0,0)
  591		)
  592	; equals_zero(Yl), equals_zero(Yh) ->           % X**0 == Z
  593		!, ^(Z,(1,1),NewZ), NewX = X, NewY = Y      %  ==> Z := (1,1)
  594	).
  595
  596powr_(_, Z,X,(R,R), NewZ,NewX,NewY) :- rational(R,N,D), !,   % R = rational point exponent
  597	N1 is N mod 2, D1 is D rem 2,  % use `mod` for negative N
  598	pt_powrCase(N1,D1,Z,X,R,NewZ), non_empty(NewZ),
  599	Ri is 1/R,
  600	pt_powrCase(D1,N1,X,NewZ,Ri,NewX), non_empty(NewX),
  601	NewY = (R,R).
  602
  603powr_(_, Z,X,Y, NewZ,NewX,NewY) :-  % interval Y
  604	intCase(C,X),
  605	powr_intY_(C, Z,X,Y, NewZ,NewX,NewY).
  606
  607equals_zero(0).
  608equals_zero(0.0).
  609equals_zero(-0.0).
  610
  611equals_one(1).
  612equals_one(1.0).
  613
  614% Assumption : even numerator and denominator can't occur (non-normalized rational).
  615pt_powrCase(N1,D1,Z,X,R,NewZ) :-  R < 0, !,      % R negative
  616	Rp is -R,
  617	universal_interval(UI),
  618	intCase(Cz,Z),
  619	odivCase(p,Cz,(1,1),Z,UI,Zi),                % Zi := UI ^ 1/Z
  620	pt_powrCase(N1,D1,Zi,X,Rp,NZi),
  621	intCase(CNzi,NZi),
  622	odivCase(p,CNzi,(1,1),NZi,Z,NewZ).           % NewZ := Z ^ 1/NZi
  623
  624pt_powrCase(1,0,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :-  % R positive, numerator odd, denominator even
  625	% reflect about X axis to include positive and negative roots
  626	Xh >=0,   % some part of X must be positive
  627	Zl1 is maxr(0, roundtoward(Xl**R,to_negative)),  % neg Xl ==> nan
  628	Zh1 is roundtoward(Xh**R,to_positive),
  629	Zpl is maxr(Zl,Zl1),  Zph is minr(Zh,Zh1),       % positive case
  630	Znl is maxr(Zl,-Zh1), Znh is minr(Zh,-Zl1),      % negative case
  631	( 1 is cmpr(Znl,Znh) -> NZl is Zpl, NZh is Zph   % Znl > Znh -> negative case is empty
  632	; 1 is cmpr(Zpl,Zph) -> NZl is Znl, NZh is Znh   % Zpl > Zph -> positive case is empty
  633	; NZl is minr(Zpl,Znl), NZh is maxr(Zph,Znh)     % else union of positive and negative cases
  634	).
  635
  636pt_powrCase(0,1,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :-  % R positive, numerator even, denominator odd
  637	% reflect about Z axis
  638	( Xh < 0 -> Xl1 is -Xh, Xh1 is -Xl           % negative X, negate interval
  639	;(Xl > 0 -> Xl1 = Xl, Xh1 = Xh               % positive X, As is
  640	;           Xl1 = 0, Xh1 is maxr(-Xl,Xh)     % split
  641	)),
  642	% NZl can't be negative
  643	NZl is maxr(Zl, roundtoward(Xl1**R,to_negative)),  % Xl1 known to be positive
  644	NZh is minr(Zh, roundtoward(Xh1**R,to_positive)).
  645
  646pt_powrCase(1,1,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :-  % R positive, numerator odd, denominator odd
  647	% continuous monotonic
  648	NZl is maxr(Zl, roundtoward(Xl**R,to_negative)),
  649	NZh is minr(Zh, roundtoward(Xh**R,to_positive)).
  650
  651% Y is an interval
  652powr_intY_(p, Z,X,Y, NewZ,NewX,NewY) :-              % positive X, interval Y 
  653	powr_prim_(Z,X,Y,NewZ),                          % NewZ := X**Y
  654	universal_interval(UI),
  655	intCase(Cy,Y),
  656	odivCase(p,Cy,(1,1),Y,UI,Yi),                    % Yi := UI ^ 1/Y 
  657	powr_prim_(X,NewZ,Yi,NewX),                      % NewX := NewZ**(1/Y) (no narrow if 0 in Y)
  658	ln(NewZ,Ynum), intCase(Cyn,Ynum),
  659	ln(NewX,Yden), intCase(Cyd,Yden),
  660	odivCase(Cyn,Cyd,Ynum,Yden,Y,NewY).              % NewY := Y ^ log(NewZ)/log(NewX)
  661
  662powr_intY_(n, Z,X,Y, NewZ,NewX,NewY) :-              % negative X, interval Y
  663	% In this case Y may contain an rational which is legal for negative X, e.g.,
  664	% to calculate an odd power or root. However it still may be able to bound Z.
  665	% This case flips X to positive, calculates upper bound Zh1, and then intersects
  666	% Z with (-Zh1,Zh1). 
  667	X = (Xl,Xh), X1l is -Xh, X1h is -Xl,
  668	powr_intY_(p, Z,(X1l,X1h),Y, (_,Z1h),(X2l,X2h),NewY),  % use code for positive case
  669	X3l is -X2h, X3h is -X2l, ^(X,(X3l,X3h),NewX),    % flip back for NewX
  670	Z1l is -Z1h, ^(Z,(Z1l,Z1h),NewZ).                 % maximal range = (-Z1h,Z1h)  
  671
  672powr_intY_(s, Z,X,Y, NewZ,NewX,NewY):-                % split X, interval Y
  673	% union of positive and negative regions
  674	X = (Xl,Xh), Xfl is -Xl,  % flip negative to positive for calculation
  675	powr_intY_(p, Z,(0,Xfl),Y, (_,Znh),NXn,NYn),      % Znl always 0,
  676	powr_intY_(p, Z,(0,Xh), Y, (_,Zph),NXp,NYp),      % Zpl always 0
  677	% union and intersect with original
  678	NZl is -Znh, NZh is maxr(Znh,Zph), ^(Z,(NZl,NZh),NewZ),
  679	union_(NXn,NXp,NX), ^(X,NX,NewX),
  680	union_(NYn,NYp,NY), ^(Y,NY,NewY).
  681	
  682powr_prim_((Zl,Zh), (Xl,Xh), (Yl,Yh), NewZ) :-        % X assumed positive
  683	Zll is roundtoward(Xl**Yl,to_negative),
  684	Zlh is roundtoward(Xl**Yh,to_positive),
  685	Zhl is roundtoward(Xh**Yl,to_negative),
  686	Zhh is roundtoward(Xh**Yh,to_positive),
  687	NZl is maxr(Zl,minr(Zll, minr(Zlh, minr(Zhl, Zhh)))), % intersect with min of all
  688	NZh is minr(Zh,maxr(Zll, maxr(Zlh, maxr(Zhl, Zhh)))), % intersect with max of all
  689	(non_empty(NZl,NZh) -> NewZ = (NZl,NZh) ; NewZ = (Zl,Zh)).  % on empty, no narrowing 
  690
  691ln((Xl,Xh), (Zl,Zh)) :-
  692	Zl is roundtoward(log(Xl),to_negative),
  693	Zh is roundtoward(log(Xh),to_positive).
  694
  695%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  696
  697% Z== sin(X)
  698
  699narrowing_op(sin, _, $(Z,X), $(NewZ, NewX)) :-
  700	sin_(X,S,Z,NewZ),
  701	asin_(NewZ,S,X,NewX).
  702
  703sin_((Xl,Xh),(Sl,Sh),Z,NZ) :-
  704	(Xl == Xh, abs(Xl) =:= 1.0Inf
  705	 -> Sl = -10, Sh = 10,                 %  point infinity special case
  706	    ^(Z,(-1.0,1.0),NZ)
  707	 ;  Sl is floor(roundtoward(Xl/pi+0.5,to_positive)), % determine monotonic sector of width pi
  708	    Sh is floor(roundtoward(Xh/pi+0.5,to_negative)), % sector 0 is (-pi/2,pi/2) - see asin(x)
  709	    ( Sh-Sl =< 2                       % more than 2 sectors -> no narrowing
  710	     -> to_float((Xl,Xh),FX),          % convert now to float
  711	        sinCase_(FX,Sl,Sh,Z,(1,-1),NZ),
  712	        non_empty(NZ)
  713	     ; ^(Z,(-1.0,1.0),NZ) 
  714	    )
  715	).
  716
  717sinCase_(X,Sl,Sh,Z,ZIn,NZ) :-
  718	(Sl == Sh
  719	 -> sinSector_(Sl,X,Z,ZIn,NZ)          % single/last sector 
  720	 ;  XSh is pi*(Sl+0.5),                % spanning sectors =<2
  721	    X = (Xl,Xh),
  722	    sinSector_(Sl,(Xl,XSh),Z,ZIn,ZOut),
  723	    SNl is Sl+1,
  724	    sinCase_((XSh,Xh),SNl,Sh,Z,ZOut,NZ)
  725	).
  726
  727sinSector_(S,(Xl,Xh),(Zl,Zh),(ZlIn,ZhIn),(ZlOut,ZhOut)) :-
  728	( 0 is S mod 2                         % monotonically increasing ?
  729	 -> ZSl is roundtoward(sin(Xl), to_negative),
  730	    ZSh is roundtoward(sin(Xh), to_positive)
  731	 ;  ZSl is roundtoward(sin(Xh), to_negative),
  732	    ZSh is roundtoward(sin(Xl), to_positive)
  733	),
  734	( ^((Zl,Zh),(ZSl,ZSh),(Z1l,Z1h)) % union with previous Z*In
  735	 -> ZlOut is minr(ZlIn,Z1l), ZhOut is maxr(ZhIn,Z1h)
  736	 ;  ZlOut = ZlIn,            ZhOut = ZhIn 
  737	).
  738
  739asin_(Z,S,X,NX) :- 
  740	( Z = (-1.0,1.0)
  741	 -> NX = X                             % no narrowing of X possible based on Z
  742	 ;  S = (Sl,Sh),
  743	    (Sl == Sh
  744	     -> asinSector_(Sl,Z,XS),          % same sector
  745	        ^(X,XS,NX)
  746	     ;  ( Sh-Sl =< 2                   % X spans sectors
  747	         -> asinCaseLo_(Sl,Z,X,XSl),
  748	            asinCaseHi_(Sh,Z,X,XSh),
  749	            ^(X,(XSl,XSh),NX)
  750	         ;  NX = X
  751	        )
  752	    ) 
  753	).
  754
  755asinCaseLo_(S,Z,X,NXl) :-
  756	asinSector_(S,Z,XS),
  757	( ^(X,XS,(NXl,_))
  758	 -> true
  759	 ;  X = (_,Xh), Xh >= pi*(S+0.5),      % next sector in range 
  760	    S1 is S+1,
  761	    asinCaseLo_(S1,Z,X,NXl)
  762	).
  763
  764asinCaseHi_(S,Z,X,NXh) :-
  765	asinSector_(S,Z,XS),
  766	( ^(X,XS,(_,NXh))
  767	 -> true
  768	 ;  X = (Xl,_), Xl =< pi*(S-0.5),      % previous sector in range 
  769	    S1 is S-1, 
  770	    asinCaseHi_(S1,Z,X,NXh)
  771	).
  772
  773asinSector_(S,(Zl,Zh),(XSl,XSh)) :-
  774	Offs is S*pi,
  775    bounded_number(Offl,Offh,Offs),     % fuzz offset
  776%%	Offl is nexttoward(S*pi,-inf), Offh is nexttoward(S*pi,inf),  % fuzz offset
  777%%	Offl is roundtoward(S*pi,to_negative), Offh is roundtoward(S*pi,to_positive),  % fuzz offset
  778	(0 is S mod 2
  779	 -> XSl is roundtoward(Offl+asin(Zl),to_negative),  % monotonically increasing
  780	    XSh is roundtoward(Offh+asin(Zh),to_positive)
  781	 ;  XSl is roundtoward(Offl-asin(Zh),to_negative),  % monotonically decreasing
  782	    XSh is roundtoward(Offh-asin(Zl),to_positive)
  783	).
  784
  785%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  786
  787% Z== cos(X)
  788
  789narrowing_op(cos, _, $(Z,X), $(NewZ, NewX)) :-
  790	cos_(X,S,Z,NewZ),
  791	acos_(NewZ,S,X,NewX).
  792
  793cos_((Xl,Xh),(Sl,Sh),Z,NZ) :-
  794	(Xl == Xh, abs(Xl) =:= 1.0Inf
  795	 -> Sl = -10, Sh = 10,                 %  point infinity special case
  796	    ^(Z,(-1.0,1.0),NZ)
  797	 ;  Sl is floor(roundtoward(Xl/pi,to_positive)),  % determine monotonic sector of width pi
  798	    Sh is floor(roundtoward(Xh/pi,to_negative)),  % sector 0 is (0,pi) - see acos(x)
  799	    ( Sh-Sl =< 2                           % more than 2 sectors -> no narrowing
  800	     -> to_float((Xl,Xh),FX),       % convert now to float
  801	        cosCase_(FX,Sl,Sh,Z,(1,-1),NZ),
  802	        non_empty(NZ)
  803	     ; ^(Z,(-1.0,1.0),NZ)
  804	    )
  805	).
  806
  807cosCase_(X,Sl,Sh,Z,ZIn,NZ) :- 
  808	( Sl == Sh
  809	 -> cosSector_(Sl,X,Z,ZIn,NZ)          % single/last sector 
  810	 ;  SNl is Sl+1,                       % spanning sectors =<2
  811	    XSh is pi*(Sl+1),
  812	    X = (Xl,Xh),
  813	    cosSector_(Sl,(Xl,XSh),Z,ZIn,ZOut),
  814	    cosCase_((XSh,Xh),SNl,Sh,Z,ZOut,NZ)
  815	).
  816
  817cosSector_(S,(Xl,Xh),(Zl,Zh),(ZlIn,ZhIn),(ZlOut,ZhOut)) :-
  818	( 1 is S mod 2                         % monotonically increasing ?
  819	 -> ZSl is roundtoward(cos(Xl), to_negative),
  820	    ZSh is roundtoward(cos(Xh), to_positive)
  821	 ;  ZSl is roundtoward(cos(Xh), to_negative),
  822	    ZSh is roundtoward(cos(Xl), to_positive)
  823	),
  824	( ^((Zl,Zh),(ZSl,ZSh),(Z1l,Z1h)) % union with previous Z*In
  825	 -> ZlOut is minr(ZlIn,Z1l), ZhOut is maxr(ZhIn,Z1h)
  826	 ;  ZlOut = ZlIn,            ZhOut = ZhIn 
  827	).
  828
  829acos_(Z,S,X,NX) :- 
  830	( Z = (-1.0,1.0)
  831	 -> NX = X                             % no narrowing of X possible based on Z
  832	 ;  S = (Sl,Sh),
  833	    (Sl == Sh
  834	     -> acosSector_(Sl,Z,XS),          % same sector
  835	        ^(X,XS,NX)
  836	     ;  ( Sh-Sl =< 2                   % X spans sectors
  837	         -> acosCaseLo_(Sl,Z,X,XSl),
  838	            acosCaseHi_(Sh,Z,X,XSh),
  839	            ^(X,(XSl,XSh),NX)
  840	         ;  NX = X
  841	        )
  842	    ) 
  843	).
  844
  845acosCaseLo_(S,Z,X,NXl) :-
  846	acosSector_(S,Z,XS),
  847	( ^(X,XS,(NXl,_))
  848	 -> true
  849	 ;  S1 is S+1,                         % next sector in range 
  850	    X = (_,Xh), Xh >= pi*S1,
  851	    acosCaseLo_(S1,Z,X,NXl)
  852	).
  853
  854acosCaseHi_(S,Z,X,NXh) :-
  855	acosSector_(S,Z,XS),
  856	( ^(X,XS,(_,NXh))
  857	 -> true
  858	 ;  X = (Xl,_), Xl =< pi*S,            % previous sector in range 
  859	    S1 is S-1,
  860	    acosCaseHi_(S1,Z,X,NXh)
  861	).
  862
  863acosSector_(S,(Zl,Zh),(XSl,XSh)) :-
  864	X1l is roundtoward(acos(Zh),to_negative),
  865	X1h is roundtoward(acos(Zl),to_positive),
  866	(1 is S mod 2
  867	 -> Offs is (S+1)*pi,
  868	    bounded_number(Offl,Offh,Offs),     % fuzz offset
  869	    XSl is roundtoward(Offl-X1h,to_negative),  % monotonically increasing
  870	    XSh is roundtoward(Offh-X1l,to_positive)
  871	 ;  Offs is S*pi,
  872	    bounded_number(Offl,Offh,Offs),     % fuzz offset
  873	    XSl is roundtoward(Offl+X1l,to_negative),  % monotonically decreasing
  874	    XSh is roundtoward(Offh+X1h,to_positive)
  875	).
  876
  877%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  878
  879% Z== tan(X)
  880
  881narrowing_op(tan, _, $(Z,X), $(NewZ, NewX)) :-
  882	X = (Xl,Xh),
  883	Sl is round(Xl/pi),  % sector 0 is (-pi/2,pi/2)
  884	Sh is round(Xh/pi),
  885	tanCase(Z,X,Sl,Sh,NewZ,NewX).
  886	
  887% Need to maintain adequate width applying offset, otherwise loss of precision
  888%	leads to unintentional rounding. This is done by fuzzing offset. 
  889tanCase(Z,(X,X),_,_,Z,(X,X)) :-    % point infinity special case
  890	abs(X) =:= 1.0Inf, !.
  891tanCase((Zl,Zh),(Xl,Xh),S,S,(NZl,NZh),NX) :-  !,          % same sector
  892	Z1l is roundtoward(tan(Xl),to_negative),
  893	Z1h is roundtoward(tan(Xh),to_positive),
  894	( non_empty(Z1l,Z1h)
  895	 -> ^((Zl,Zh),(Z1l,Z1h),(NZl,NZh)),                   % good range
  896	    atan_((NZl,NZh),S,X1),
  897	    ^((Xl,Xh),X1,NX)
  898	 ;  % same sector range check failed, assume rounding crossed singularity
  899	    ((Z1l =< Zh ; Z1h >= Zl) -> true),                % disjoint ranges must intersect Z
  900	    NZl = Zl, NZh = Zh, NX = (Xl,Xh)                  % nothing changes
  901	).
  902tanCase(Z,(Xl,Xh),Sl,Sh,(NZl,NZh),(NXl,NXh)) :-  Sh-Sl =:= 1,  % spans one singularity
  903	Sg is pi*(Sl+0.5),               % closest to singularity (integer Sl and 0.5 is precise float)
  904	bounded_number(XSl,XSh,Sg),      % use bounded number to straddle
  905	(tanCase(Z,(Xl,XSl),Sl,Sl,(NZ1l,NZ1h),(NXl,NX1h))
  906	 -> (tanCase(Z,(XSh,Xh),Sh,Sh,_NZ2,(_,NXh))
  907	     -> fail                                          % both sectors, includes singularity
  908	     ;  NXh = NX1h, NZl = NZ1l, NZh = NZ1h            % only lo sector
  909	    )
  910	 ;  tanCase(Z,(XSh,Xh),Sh,Sh,(NZl,NZh),(NXl,NXh))     % only hi sector
  911	),
  912	!.
  913tanCase(Z,X,Sl,Sh,Z,NX) :-  Sh-Sl =:= 2,

  914	% X is odd multiple of (-pi/2,pi/2)	?
  915	X = (Xl,Xh),
  916	LB is rational(roundtoward(Xl/pi + 0.5, to_positive)), integer(LB),	
  917	UB is rational(roundtoward(Xh/pi + 0.5, to_negative)), integer(UB),
  918	!,
  919	S is (Sl+Sh)/2,
  920	atan_(Z,S,X1), ^(X,X1,NX)

  920.
  921tanCase(Z,X,Sl,Sh,Z,X).                                   % spans multiple singularities
  922
  923atan_((Zl,Zh),S,(Xl,Xh)) :-
  924	Offs is S*pi,
  925	bounded_number(Offl,Offh,Offs),                       % fuzz offset
  926	Xl is Offl+roundtoward(atan(Zl),to_negative),
  927	Xh is Offh+roundtoward(atan(Zh),to_positive).
  928
  929%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  930
  931% Z== ~X boolean negation (Z and X boolean)
  932
  933narrowing_op(not, P, $(Z,X), $(NewZ, NewX)) :-
  934	notB_(X,Z,Z1,P), ^(Z,Z1,NewZ),
  935	notB_(NewZ,X,X1,P), ^(X,X1,NewX).
  936
  937notB_((0,0), _, (1,1), p) :- !.
  938notB_((1,1), _, (0,0), p) :- !.
  939notB_(    _, Z, NewZ,  _) :- ^(Z,(0,1),NewZ).
  940
  941%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  942
  943% Z==X and Y  boolean 'and'
  944
  945narrowing_op(and, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  946	andB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
  947
  948andB_rel_(Z,(1,1),(1,1), NewZ,(1,1),(1,1), p) :- !, ^(Z,(1,1),NewZ).  % all cases persistent (points) except last
  949andB_rel_(Z,(0,0),Y,     NewZ,(0,0),Y, p)     :- !, ^(Z,(0,0),NewZ), ^(Y,(0,1),NewY).
  950andB_rel_(Z,X,(0,0),     NewZ,X,(0,0), p)     :- !, ^(Z,(0,0),NewZ), ^(X,(0,1),NewX).
  951andB_rel_((1,1),X,Y,     (1,1),NewX,NewY, p)  :- !, ^(X,(1,1),NewX), ^(Y,(1,1),NewY).
  952andB_rel_((0,0),X,(1,1), (0,0),NewX,(1,1), p) :- !, ^(X,(0,0),NewX).
  953andB_rel_((0,0),(1,1),Y, (0,0),(1,1),NewY, p) :- !, ^(Y,(0,0),NewY).
  954andB_rel_(Z,X,Y,         NewZ,NewX,NewY, P)   :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
  955
  956%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  957
  958% Z==X and Y  boolean 'nand'
  959
  960narrowing_op(nand, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  961	nandB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
  962
  963nandB_rel_(Z,(1,1),(1,1), NewZ,(1,1),(1,1), p) :- !, ^(Z,(0,0),NewZ).  % all cases persistent except last
  964nandB_rel_(Z,(0,0),Y,     NewZ,(0,0),Y, p)     :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY).
  965nandB_rel_(Z,X,(0,0),     NewZ,X,(0,0), p)     :- !, ^(Z,(1,1),NewZ), ^(X,(0,1),NewX).
  966nandB_rel_((0,0),X,Y,     (0,0),NewX,NewY, p)  :- !, ^(X,(1,1),NewX), ^(Y,(1,1),NewY).
  967nandB_rel_((1,1),X,(1,1), (1,1),NewX,(1,1), p) :- !, ^(X,(0,0),NewX).
  968nandB_rel_((1,1),(1,1),Y, (1,1),(1,1),NewY, p) :- !, ^(Y,(0,0),NewY).
  969nandB_rel_(Z,X,Y,         NewZ,NewX,NewY, P)   :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
  970
  971%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  972
  973% Z==X or Y  boolean 'or'
  974
  975narrowing_op(or, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  976	orB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
  977
  978orB_rel_(Z,(0,0),(0,0), NewZ,(0,0),(0,0), p) :- !, ^(Z,(0,0),NewZ).  % all cases persistent (points) except last
  979orB_rel_(Z,(1,1),Y,     NewZ,(1,1),Y, p)     :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY).
  980orB_rel_(Z,X,(1,1),     NewZ,X,(1,1), p)     :- !, ^(Z,(1,1),NewZ), ^(X,(0,1),NewX).
  981orB_rel_((0,0),X,Y,     (0,0),NewX,NewY, p)  :- !, ^(X,(0,0),NewX), ^(Y,(0,0),NewY).
  982orB_rel_((1,1),X,(0,0), (1,1),NewX,(0,0), p) :- !, ^(X,(1,1),NewX).
  983orB_rel_((1,1),(0,0),Y, (1,1),(0,0),NewY, p) :- !, ^(Y,(1,1),NewY).
  984orB_rel_(Z,X,Y,         NewZ,NewX,NewY, P)   :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
  985
  986%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  987
  988% Z==X nor Y  boolean 'nor'
  989
  990narrowing_op(nor, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  991	norB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
  992
  993norB_rel_(Z,(0,0),(0,0), NewZ,(0,0),(0,0), p) :- !, ^(Z,(1,1),NewZ).  % all cases persistent (points) except last
  994norB_rel_(Z,(1,1),Y,     NewZ,(1,1),Y, p)     :- !, ^(Z,(0,0),NewZ), ^(Y,(0,1),NewY).
  995norB_rel_(Z,X,(1,1),     NewZ,X,(1,1), p)     :- !, ^(Z,(0,0),NewZ), ^(X,(0,1),NewX).
  996norB_rel_((1,1),X,Y,     (1,1),NewX,NewY, p)  :- !, ^(X,(0,0),NewX), ^(Y,(0,0),NewY).
  997norB_rel_((0,0),X,(0,0), (0,0),NewX,(0,0), p) :- !, ^(X,(1,1),NewX).
  998norB_rel_((0,0),(0,0),Y, (0,0),(0,0),NewY, p) :- !, ^(Y,(1,1),NewY).
  999norB_rel_(Z,X,Y,         NewZ,NewX,NewY, P)   :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
 1000
 1001%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1002
 1003% Z==X xor Y  boolean 'xor'
 1004
 1005narrowing_op(xor, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
 1006	xorB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
 1007
 1008xorB_rel_(Z,(B,B),(B,B), NewZ,(B,B),(B,B), p) :- !, ^(Z,(0,0),NewZ).  % all cases persistent (points) except last
 1009xorB_rel_(Z,(X,X),(Y,Y), NewZ,(X,X),(Y,Y), p) :- !, ^(Z,(1,1),NewZ).
 1010xorB_rel_((B,B),X,(B,B), (B,B),NewX,(B,B), p) :- !, ^(X,(0,0),NewX).
 1011xorB_rel_((Z,Z),X,(Y,Y), (Z,Z),NewX,(Y,Y), p) :- !, ^(X,(1,1),NewX).
 1012xorB_rel_((B,B),(B,B),Y, (B,B),(B,B),NewY, p) :- !, ^(Y,(0,0),NewY).
 1013xorB_rel_((Z,Z),(X,X),Y, (Z,Z),(X,X),NewY, p) :- !, ^(Y,(1,1),NewY).
 1014xorB_rel_(Z,X,Y,         NewZ,NewX,NewY, P)   :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
 1015
 1016%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1017
 1018% Z==X -> Y  boolean 'implies'
 1019
 1020narrowing_op(imB, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
 1021	imB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
 1022
 1023imB_rel_(Z,(1,1),Y, New,(1,1),New,   P) :- !, ^(Z,Y,New), (New=(B,B) -> P=p ; true).   % X=1, Z=Y, persistant if point
 1024imB_rel_(Z,(0,0),Y, NewZ,(0,0),Y,    p) :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY).        % X=0, Z=1, persistant
 1025imB_rel_(Z,X,(0,0), NewZ,NewX,(0,0), P) :- !, notB_(X,Z,NewZ,P), notB_(NewZ,X,NewX,P). % Y=0, Z=~X, persistant if point
 1026imB_rel_(Z,X,(1,1), NewZ,X,(1,1),    P) :- !, ^(Z,(1,1),NewZ).                      % Y=1, Z=1, persistant
 1027imB_rel_((0,0),X,Y, (0,0),NewX,NewY, p) :- !, ^(X,(1,1),NewX), ^(Y,(0,0),NewY).     % Z=0,X=1,Y=0, persistant
 1028imB_rel_(Z,X,Y,     NewZ,NewX,NewY,  P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
 1029
 1030%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1031
 1032:- style_check(+singleton).  % for narrowing_op
 1033
 1034% User defined IA primitives
 1035
 1036narrowing_op(user(Prim), P, InArgs, OutArgs) :-
 1037	call_user_primitive(Prim, P, InArgs, OutArgs).  % from main module body