1/*	The MIT License (MIT)
    2 *
    3 *	Copyright (c) 2019-2024 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% non-empty interval test
   54%
   55goal_expansion(non_empty(L,H), 1 =\= cmpr(L,H)).  % SWIP inline optimization for non_empty/2
   56
   57non_empty((L,H)) :- non_empty(L,H).
   58non_empty(L,H) :- 1 =\= cmpr(L,H).
   59
   60%
   61% interval category: non-negative (includes zero), non-positive, or split (contains 0)
   62%
   63intCase(C, (L,H)) :-
   64	( L>=0 -> C=p            % non-negative
   65	; H=<0 -> C=n            % non-positive
   66	;         C=s            % split
   67	).
   68
   69%
   70% Z := X ^ Y  (interval intersection)
   71%
   72^((Xl,Xh), (Yl,Yh), (Zl,Zh)) :-
   73	Zl is maxr(Xl, Yl), Zh is minr(Xh, Yh),  % eliminates NaN's, prefers rationals
   74	non_empty(Zl,Zh).
   75
   76%
   77% Z := X \/ Y (interval union)
   78%
   79union_((Xl,Xh),(Yl,Yh),(Zl,Zh)) :-
   80	Zl is minr(Xl,Yl), Zh is maxr(Xh,Yh).  % eliminates NaN's, prefers rationals
   81
   82%
   83% to_float/2 converts an non-float bounds to floats
   84%  useful when when float arguments will be produced on narrowing, e.g., elementary functions
   85%  avoids incorrect rounding on argument conversion (monotonic decreasing segments)
   86%
   87to_float((L,H),(FL,FH)) :-
   88	( float(L) -> FL = L ; FL is roundtoward(float(L), to_negative) ), 
   89	( float(H) -> FH = H ; FH is roundtoward(float(H), to_positive) ). 
   90
   91%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   92%
   93% Relational Operations (narrowing operations)
   94%
   95:- discontiguous clpBNR:narrowing_op/4.   96:- style_check(-singleton).  % for narrowing_op
   97%
   98%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   99
  100% integral (contract to integer bounds) - optimized version of (int,P,$(1,X),_)
  101
  102narrowing_op(integral, P, $(X), $(NewX)) :-
  103	integral_(X,NewX).
  104	
  105integral_((Xl,Xh),(NXl,NXh)) :-
  106	NXl is ceiling(Xl), NXh is floor(Xh),
  107	non_empty(NXl,NXh).
  108
  109%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  110
  111% Z == integer(X) % (Z boolean) 
  112%  true if X is an integer point, false if it doesn't contain an integer
  113%  else if true, narrows X to integer bounds
  114
  115narrowing_op(int, P, $(Z,X), $(NewZ,NewX)) :-
  116	int_(Z,X,P,NewZ,NewX).
  117	
  118int_(Z,X,p,NewZ,X) :-                       % persistent cases
  119	X = (Xl,Xh),
  120	(0 is cmpr(Xl,Xh)                   % point value?
  121	 -> (integer(Xl) -> B = 1 ; B = 0)  % yes, true if integer else false
  122	 ;  ceiling(Xl) > floor(Xh),        % no, and X doesn't contain an integer
  123	    B = 0 
  124	),
  125	!,  % succeed, persistent
  126	^(Z,(B,B),NewZ).
  127	
  128int_(Z,X,_P,NewZ,NewX) :-                   % possible narrowing case
  129	(Z = (1,1) 
  130	 -> integral_(X,NewX),              % Z true, narrow X to integers
  131	    NewZ = Z
  132	 ;  NewX = X,                       % Z false or indeterminate, no change to X
  133	    ^(Z,(0,1),NewZ)                 % Z boolean
  134	).
  135
  136%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  137
  138% Z==(X==Y)  % (Z boolean)
  139
  140narrowing_op(eq, P, $(Z,X,Y), Out) :-
  141	(^(X,Y,New)                       % intersect?
  142	 -> (Z = (1,1)
  143	     -> Out = $(Z,New,New), 
  144	        New = (L,H),
  145	        (0 is cmpr(L,H) -> P = p ; true)   % persistent if a point
  146	     ; % intersect and Z false or indeterminate, all things possible
  147	        ( X = (Xl,Xh), Y = (Yl,Yh),
  148	         0 is cmpr(Xl,Xh) \/ cmpr(Yl,Yh) \/ cmpr(Xl,Yl)  % X and Y same point
  149	         -> ^(Z,(1,1),NewZ),      % Z true
  150	            P = p                 % persistant
  151	         ;  ^(Z,(0,1),NewZ)       % ensure Z boolean
  152	        ),
  153	        Out = $(NewZ,X,Y)
  154	    )
  155	 ;  % X and Y do not intersect  
  156	    ^(Z, (0,0), NewZ),            % fail primitive if Z cannot be false
  157	    Out = $(NewZ,X,Y), P = p      % persistent
  158	 ).
  159
  160%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  161
  162% Z==(X<>Y)  % (Z boolean, X,Y integer)
  163
  164narrowing_op(ne, P, $(Z,X,Y), Out) :-
  165	(^(X,Y,(L,H))
  166	 -> (Z = (1,1)
  167	     -> ne_int_(X,Y,NewX),        % narrow X if Y is a point and a bound of X (always succeeds but may not narrow
  168	        ne_int_(Y,NewX,NewY),     % narrow Y if NewX is a point and a bound of Y (always succeeds but may not narrow)	
  169	        Out = $(Z,NewX,NewY)
  170	     ;  (X = Y, L = H -> B = (0,0), P = p ;  B = (0,1)),
  171	        ^(Z,B,NewZ),
  172	        Out = $(NewZ,X,Y)
  173	     )
  174	 ;  % X and Y do not intersect
  175	    ^(Z, (1,1), NewZ),            % true and 
  176	    P = p,                        % persistent
  177	    Out = $(NewZ,X,Y)
  178	).
  179
  180% narrow X and Y given Y <> X,  where Y and X are integer intervals (enforced elsewhere)
  181% integers so standard comparison works.
  182ne_int_((X,H), (X,X), (NewL,H)) :- !,  % X is a point,  and low bound of Y
  183	NewL is X+1, NewL=<H.
  184ne_int_((L,X), (X,X), (L,NewH)) :- !,  % X is a point,  and high bound of Y
  185	NewH is X-1, L=<NewH.
  186ne_int_(Y, _X, Y).                      % no narrowing
  187
  188%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  189
  190% Z==(X=<Y)  % (Z boolean)
  191
  192narrowing_op(le, P, $(Z, X, Y), Out) :-
  193	X = (Xl,Xh), Y = (Yl,Yh),
  194	(non_empty(Xh,Yl)                                  % Xh =< Yl
  195	 -> ^(Z,(1,1),NewZ), P = p, Out = $(NewZ,X,Y)      % persistent X =< Y
  196	 ;  (-1 is cmpr(Yh,Xl)                             % Yh < Xl
  197	     -> ^(Z,(0,0),NewZ), P = p, Out = $(NewZ,X,Y)  % persistent Y < X
  198	     ;  (Z = (1,1)                                 % if {X =< Y}, can narrow X and Y
  199	         -> NXh is minr(Xh,Yh),                    % NewX := (Xl,NXh) 
  200	            NYl is maxr(Xl,Yl),                    % NewY := (NYl,Yh)
  201	            Out = $(Z,(Xl,NXh),(NYl,Yh))
  202	         ;  ^(Z,(0,1),NewZ),                       % Z false or indeterminate
  203	            Out = $(NewZ,X,Y)
  204	        )
  205	    )
  206	).
  207
  208%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  209
  210% Z==(X<Y)  % (Z boolean) generally unsound on reals but possibly useful for minima/maxima
  211
  212narrowing_op(lt, P, $(Z, X, Y), Out) :-
  213	X = (Xl,Xh), Y = (Yl,Yh),
  214	(-1 is cmpr(Xh,Yl)                                 % Xh < Yl
  215	 -> ^(Z,(1,1),NewZ), P = p, Out = $(NewZ,X,Y)      % persistent X < Y
  216	 ;  (1 =\= cmpr(Yh,Xl)                             % Yh =< Xl
  217	     -> ^(Z,(0,0),NewZ), P = p, Out = $(NewZ,X,Y)  % persistent Y =< X	
  218	     ;  (Z = (1,1)                                 % if {X < Y}, can narrow X and Y
  219	         -> next_lt_(Yh,-1,YhD),                   % YhD is next downward value from Yh
  220	            NXh is minr(Xh,YhD),                   % NewX := (Xl,NXh)
  221	            next_lt_(Xl,1,XlU),                    % NXlU is next upward value from NXl
  222	            NYl is maxr(XlU,Yl),                   % NewY := (NYl,Yh)
  223	            (NXh < NYl -> P=p ; true),             % persistent due to increment?
  224	            Out = $(Z,(Xl,NXh),(NYl,Yh))
  225	         ;  ^(Z,(0,1),NewZ),                       % Z false or indeterminate
  226	            Out = $(NewZ,X,Y)
  227	        )
  228	    )
  229	).
  230
  231% Note: can't narrow an infinite bound, minimize change to bound
  232% Repeat, not sound on reals (uses nexttoward, missing values between floats)
  233next_lt_( 1.0Inf, _,  1.0Inf) :- !.
  234next_lt_(-1.0Inf, _, -1.0Inf) :- !.
  235next_lt_(V, -1, NV) :- NV is maxr(V-1,nexttoward(V,-1.0Inf)).  % integers will get sorted out with `integral`
  236next_lt_(V,  1, NV) :- NV is minr(V+1,nexttoward(V, 1.0Inf)).
  237
  238%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  239
  240% Z==(X<=Y)  % inclusion, constrains X to be subinterval of Y (Z boolean)
  241
  242narrowing_op(in, P, $(Z, X, Y), $(NewZ, NewX, Y)):-
  243	% Only two cases: either X and Y intersect or they don't.
  244	(^(X,Y,X1)               % X1 is intersection of X and Y
  245	 -> (Z == (1,1) 
  246	     -> NewX = X1        % Z true, X must be subinterval of Y
  247	     ;  NewX = X,        % Z is false or indeterminate, X unchanged
  248	        ^(Z,(0,1),NewZ)  % Z boolean
  249	    )
  250	 ;  ^(Z,(0,0),NewZ),     % X and Y don't intersect, Z must be false
  251	 	NewX = X,            % no change in X
  252	    P=p                  % persistent, since X and Y can never intersect
  253	).
  254
  255%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  256
  257% Z==X+Y
  258
  259narrowing_op(add, _, $((Zl,Zh), (Xl,Xh), (Yl,Yh)), $((NZl,NZh), (NXl,NXh), (NYl,NYh))) :-
  260	NZl is maxr(Zl, roundtoward( Xl+Yl,  to_negative)),       % NewZ := Z ^ (X+Y),
  261	NZh is minr(Zh, roundtoward( Xh+Yh,  to_positive)),
  262	non_empty(NZl,NZh),
  263	% Note: subtraction done by adding minus values so rounding mode consistent
  264	% during any numeric type conversion.
  265	NXl is maxr(Xl, roundtoward(NZl+(-Yh),  to_negative)),    % NewX := X ^ (NZ-Y),
  266	NXh is minr(Xh, roundtoward(NZh+(-Yl),  to_positive)),
  267	non_empty(NXl,NXh),
  268	NYl is maxr(Yl, roundtoward(NZl+(-NXh), to_negative)),    % NewY := Y ^ (NZ-NX).
  269	NYh is minr(Yh, roundtoward(NZh+(-NXl), to_positive)),
  270	non_empty(NYl,NYh).
  271
  272%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  273
  274% Z==X*Y
  275
  276narrowing_op(mul, _, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  277	intCase(Cx,X),
  278	intCase(Cy,Y),
  279	multCase(Cx,Cy,X,Y,Z,NewZ),                       % NewZ := Z ^ (X*Y)
  280	intCase(CNz,NewZ),
  281	odivCase(CNz,Cx,NewZ,X,Y,Yp),                     % Yp := Y ^ (Z/X),
  282	intCase(Cyp,Yp),
  283	odivCase(CNz,Cyp,NewZ,Yp,X,NewX),                 % NewX := X ^ (Z/Yp),
  284	% if X narrowed it may be necessary to recalculate Y due to non-optimal ordering.
  285	(Y = Yp, \+(X = NewX)                   % if Y didn't narrow and X did
  286	 -> intCase(CNx,NewX),
  287	    odivCase(CNz,CNx,NewZ,NewX,Y,NewY)            % re calculate: NewY := Y ^ NewZ/NewX
  288	 ;  NewY = Yp
  289	).
  290
  291%
  292% * cases ("Interval Arithmetic: from Principles to Implementation", Fig. 3)
  293%
  294% Note: for mixed sign, mixed mode, one of the rounding directions will be incorrect.
  295%  For now assume(?) this will at worst cancel outward rounding
  296%multCase(z,_, X, _, _, X) :- !.  % X==0
  297%multCase(_,z, _, Y, _, Y).       % Y==0
  298
  299multCase(p,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- 
  300	NZl is maxr(Zl,roundtoward(Xl*Yl,to_negative)),
  301	NZh is minr(Zh,roundtoward(Xh*Yh,to_positive)),
  302	non_empty(NZl,NZh).
  303multCase(p,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- 
  304	NZl is maxr(Zl,roundtoward(Xh*Yl,to_negative)),
  305	NZh is minr(Zh,roundtoward(Xl*Yh,to_positive)),
  306	non_empty(NZl,NZh).
  307multCase(n,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- 
  308	NZl is maxr(Zl,roundtoward(Xl*Yh,to_negative)),
  309	NZh is minr(Zh,roundtoward(Xh*Yl,to_positive)),
  310	non_empty(NZl,NZh).
  311multCase(n,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
  312	NZl is maxr(Zl,roundtoward(-Xh * -Yh,to_negative)),  % negate for any conversion rounding 
  313	NZh is minr(Zh,roundtoward(-Xl * -Yl,to_positive)),
  314	non_empty(NZl,NZh).
  315
  316multCase(p,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
  317	NZl is maxr(Zl,roundtoward(Xh*Yl,to_negative)),
  318	NZh is minr(Zh,roundtoward(Xh*Yh,to_positive)),
  319	non_empty(NZl,NZh).
  320multCase(n,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
  321	NZl is maxr(Zl,roundtoward(Xl*Yh,to_negative)),
  322	NZh is minr(Zh,roundtoward(-Xl * -Yl,to_positive)),  % negate for any conversion rounding
  323	non_empty(NZl,NZh).
  324multCase(s,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
  325	NZl is maxr(Zl,roundtoward(Xl*Yh,to_negative)),
  326	NZh is minr(Zh,roundtoward(Xh*Yh,to_positive)),
  327	non_empty(NZl,NZh).
  328multCase(s,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
  329	NZl is maxr(Zl,roundtoward(Xh*Yl,to_negative)),
  330	NZh is minr(Zh,roundtoward(-Xl * -Yl,to_positive)),  % negate for any conversion rounding
  331	non_empty(NZl,NZh).
  332
  333multCase(s,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
  334	NZl is maxr(Zl,minr(roundtoward(Xl*Yh,to_negative),roundtoward(Xh*Yl,to_negative))),
  335	NZh is minr(Zh,maxr(roundtoward(-Xl* -Yl,to_positive),roundtoward(Xh*Yh,to_positive))),
  336	non_empty(NZl,NZh).
  337
  338%
  339% / cases ("Interval Arithmetic: from Principles to Implementation", Fig. 4)
  340% Tricky handling the "zero" cases - returns universal interval when no narowing possible:
  341%	1. denominator is zero or contains zero (z and s)
  342%	2. denominator is bounded by zero (p and n) and numerator contains zero (see numeric guards)
  343% Numeric guards check for 0/0, e.g., p,p -> Xh+Yh>0
  344
  345odivCase(p,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  346	sign(Yl)+sign(Xl) > 0  % X > 0 or Y > 0
  347	 -> NZl is maxr(Zl,roundtoward(Xl/Yh,to_negative)),
  348	    NZh is minr(Zh,roundtoward(Xh/max(0.0,Yl),to_positive)),
  349	    non_empty(NZl,NZh)
  350	 ;  NZl = Zl, NZh = Zh.
  351odivCase(p,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  352	sign(Xl)-sign(Yh) > 0  % X > 0 or Y < 0
  353	 -> NZl is maxr(Zl,roundtoward(-Xh / -min(-0.0,Yh),to_negative)),  % use min to preserve sign
  354	    NZh is minr(Zh,roundtoward(-Xl / -Yl          ,to_positive)),
  355	    non_empty(NZl,NZh)
  356	 ;  NZl = Zl, NZh = Zh.
  357odivCase(p,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  358	Xl > 0     % X > 0
  359	 -> ZIh is roundtoward(Xl/Yh,to_negative),  % (ZIh,inf)
  360	    ZIl is roundtoward(Xl/Yl,to_positive),  % (-inf,ZIl)
  361	    ( 1 is cmpr(Zl,ZIl) -> NZl is maxr(Zl,ZIh) ; NZl = Zl),  % Zl > ZIl ?
  362	    (-1 is cmpr(Zh,ZIh) -> NZh is minr(Zh,ZIl) ; NZh = Zh),  % Zh > ZIl ?
  363	    non_empty(NZl,NZh)
  364	 ;  NZl = Zl, NZh = Zh.
  365
  366odivCase(n,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
  367	sign(Yl)-sign(Xh) > 0  % X < 0 or Y > 0
  368	 -> NZl is maxr(Zl,roundtoward(-Xl / -max(0.0,Yl),to_negative)),
  369	    NZh is minr(Zh,roundtoward(-Xh / -Yh,to_positive)),
  370	    non_empty(NZl,NZh)
  371	 ;  NZl = Zl, NZh = Zh.
  372odivCase(n,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  373	sign(Yh)+sign(Xh) < 0  % X < 0 or Y < 0
  374	 -> NZl is maxr(Zl,roundtoward(Xh/Yl,to_negative)),
  375	    NZh is minr(Zh,roundtoward(Xl/min(-0.0,Yh),to_positive)),  % use min to preserve sign
  376	    non_empty(NZl,NZh)
  377	 ;  NZl = Zl, NZh = Zh.
  378odivCase(n,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  379	Xh < 0     % X < 0
  380	 -> ZIh is roundtoward( Xh /  Yl,to_negative),  % (ZIh,inf)
  381	    ZIl is roundtoward(-Xh / -Yh,to_positive),  % (-inf,ZIl)
  382	    ( 1 is cmpr(Zl,ZIh) -> NZl is maxr(Zl,ZIl) ; NZl = Zl),  % Zl > ZIh ?
  383	    (-1 is cmpr(Zh,ZIl) -> NZh is minr(Zh,ZIh) ; NZh = Zh),  % Zh > ZIl ?
  384	    non_empty(NZl,NZh)
  385	 ;  NZl = Zl, NZh = Zh.
  386
  387odivCase(s,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  388	Yl > 0     % Y > 0
  389	 -> NZl is maxr(Zl,roundtoward(-Xl / -Yl,to_negative)),
  390	    NZh is minr(Zh,roundtoward( Xh /  Yl,to_positive)),
  391	    non_empty(NZl,NZh)
  392	 ;  NZl = Zl, NZh = Zh.
  393odivCase(s,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
  394	Yh < 0     % Y < 0
  395	 -> NZl is maxr(Zl,roundtoward(Xh/Yh,to_negative)),
  396	    NZh is minr(Zh,roundtoward(Xl/Yh,to_positive)),
  397	    non_empty(NZl,NZh)
  398	 ;  NZl = Zl, NZh = Zh.
  399
  400odivCase(s,s,       _,       _,       Z,         Z).  % both contain 0, no narrowing
  401
  402%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  403
  404% Z==min(X,Y)    Z==max(X,Y)
  405
  406narrowing_op(min, _, $((Zl,Zh),(Xl,Xh),(Yl,Yh)), New) :-
  407	Z1l is maxr(Zl,minr(Xl,Yl)),   % Z1 := Z ^ min(X,Y),
  408	Z1h is minr(Zh,minr(Xh,Yh)),
  409	(   Yh < Xl -> Case = 1        % narrow Y
  410	; ( Xh < Yl -> Case = 2        % narrow X
  411	;   Case = 3                   % narrow both                  
  412	  )
  413	),
  414	minimax(Case, (Zl,1.0Inf), (Z1l,Z1h),(Xl,Xh),(Yl,Yh), New).
  415
  416narrowing_op(max, _, $((Zl,Zh),(Xl,Xh),(Yl,Yh)), New) :-
  417	Z1l is maxr(Zl,maxr(Xl,Yl)),   % Z1 := Z ^ max(X,Y),
  418	Z1h is minr(Zh,maxr(Xh,Yh)),
  419	(   Xh < Yl -> Case = 1        % narrow Y
  420	; ( Yh < Xl -> Case = 2        % narrow X
  421	;   Case = 3                   % narrow both                  
  422	  )
  423	),
  424	minimax(Case, (-1.0Inf,Zh), (Z1l,Z1h),(Xl,Xh),(Yl,Yh), New).
  425	
  426minimax(1, Zpartial, Z, X, Y, $(Z, X, NewY)) :-     % Case 1, narrow Y
  427	^(Y,Zpartial,NewY).                    % NewY := Y ^ Zpartial.
  428minimax(2, Zpartial, Z, X, Y, $(Z, NewX, Y)) :-     % Case 2, narrow X
  429	^(X,Zpartial,NewX).                    % NewX := X ^ Zpartial,
  430minimax(3, Zpartial, Z, X, Y, $(Z, NewX, NewY)) :-  % Case 3, overlap
  431	^(X,Zpartial,NewX),                    % NewX := X ^ Zpartial,
  432	^(Y,Zpartial,NewY).                    % NewY := Y ^ Zpartial.
  433
  434%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  435
  436% Z==abs(X)
  437
  438narrowing_op(abs, _, $(Z,X), $(NewZ, NewX)) :-
  439	intCase(Cx,X),
  440	absCase(Cx,X,Z,NewZ), !,
  441	non_empty(NewZ),
  442	intersect_regions_(X,NewZ,NewX),
  443	non_empty(NewX).
  444
  445%
  446% intersect and join Z with positive and negative regions of X
  447%
  448intersect_regions_(Z,X,NewZ) :-
  449	(^(Z,X,ZPos) -> true ; empty_interval(ZPos)),  % positive Z region
  450	(-(X,Z,ZNeg) -> true ; empty_interval(ZNeg)),  % negative Z region
  451	union_(ZPos,ZNeg,NewZ).
  452
  453% NZ := Z ^ -X (unary minus)
  454-((Xl,Xh), (Zl,Zh), (NZl,NZh)) :-
  455	NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl),
  456	non_empty(NZl,NZh).
  457
  458%
  459% abs(X) cases
  460%
  461absCase(p, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,Xl),  NZh is minr(Zh,Xh).
  462absCase(n, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl).
  463absCase(s, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,0),   NZh is minr(Zh,maxr(-Xl,Xh)).
  464
  465%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  466
  467% Z== -X
  468
  469narrowing_op(minus, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
  470	NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl),    % NewZ := Z ^ -X, no empty check required
  471	NXl is maxr(Xl,-NZh), NXh is minr(Xh,-NZl).  % NewX := X ^ -Z, no empty check required
  472
  473%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  474
  475% Z== sqrt(X)  (Z is positive root of positive X)
  476
  477narrowing_op(sqrt, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
  478	Xh>=0, Xl1 is maxr(0,Xl),  % narrow X to positive range
  479	NZl is maxr(maxr(0,Zl),roundtoward(sqrt(Xl1),to_negative)),
  480	NZh is minr(Zh,roundtoward(sqrt(Xh),to_positive)),
  481	non_empty(NZl,NZh),
  482	NXh is minr(Xh,roundtoward(NZh*NZh,to_positive)),
  483	NXl is maxr(Xl,-NXh),
  484	non_empty(NXl,NXh).
  485
  486%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  487
  488% Z== exp(X)
  489
  490narrowing_op(exp, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
  491	NZl is maxr(Zl,roundtoward(exp(Xl),to_negative)),
  492	NZh is minr(Zh,roundtoward(exp(Xh),to_positive)),
  493	non_empty(NZl,NZh),
  494	NXl is maxr(Xl,roundtoward(log(NZl),to_negative)),
  495	NXh is minr(Xh,roundtoward(log(NZh),to_positive)),
  496	non_empty(NXl,NXh).
  497
  498%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  499
  500% Z== X**Y
  501
  502narrowing_op(pow, P, $(Z,X,Y), $(NewZ,NewX,NewY)) :-
  503	powr_(P, Z,X,Y, NewZ,NewX,NewY).
  504
  505powr_(p, Z,X,Y, NewZ,NewX,NewY) :-                  % (persistant) special cases
  506	Z = (Zl,Zh), X = (Xl,Xh), Y = (Yl,Yh), 
  507	( equals_one(Xl), equals_one(Xh) ->             % 1**Y == Z
  508		!, ^(Z,(1,1),NewZ), NewX = X, NewY = Y      %  ==> Z := (1,1)
  509	; equals_one(Zl), equals_one(Zh) ->             % X**Y == 1 , X doesn't contain 1 or -1
  510		((Xl =< -1,-1 =< Xh ; Xl =< 1,1 =< Xh)
  511		 -> fail                                    % fail to use general pt_powrCase
  512		 ;  !, NewZ = Z, NewX = X, ^(Y,(0,0),NewY)  %  ==> Y := (0,0)
  513		)
  514	; equals_zero(Yl), equals_zero(Yh) ->           % X**0 == Z
  515		!, ^(Z,(1,1),NewZ), NewX = X, NewY = Y      %  ==> Z := (1,1)
  516	).
  517
  518powr_(_, Z,X,(R,R), NewZ,NewX,NewY) :- rational(R,N,D), !,   % R = rational point exponent
  519	N1 is N mod 2, D1 is D rem 2,  % use `mod` for negative N
  520	pt_powrCase(N1,D1,Z,X,R,NewZ), non_empty(NewZ),
  521	Ri is 1/R,
  522	pt_powrCase(D1,N1,X,NewZ,Ri,NewX), non_empty(NewX),
  523	NewY = (R,R).
  524
  525powr_(_, Z,X,Y, NewZ,NewX,NewY) :-  % interval Y
  526	intCase(C,X),
  527	powr_intY_(C, Z,X,Y, NewZ,NewX,NewY).
  528
  529equals_zero(0).
  530equals_zero(0.0).
  531equals_zero(-0.0).
  532
  533equals_one(1).
  534equals_one(1.0).
  535
  536% Assumption : even numerator and denominator can't occur (non-normalized rational).
  537pt_powrCase(N1,D1,Z,X,R,NewZ) :-  R < 0, !,      % R negative
  538	Rp is -R,
  539	universal_interval(UI),
  540	intCase(Cz,Z),
  541	odivCase(p,Cz,(1,1),Z,UI,Zi),                % Zi := UI ^ 1/Z
  542	pt_powrCase(N1,D1,Zi,X,Rp,NZi),
  543	intCase(CNzi,NZi),
  544	odivCase(p,CNzi,(1,1),NZi,Z,NewZ).           % NewZ := Z ^ 1/NZi
  545
  546pt_powrCase(1,0,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :-  % R positive, numerator odd, denominator even
  547	% reflect about X axis to include positive and negative roots
  548	Xh >=0,   % some part of X must be positive
  549	Zl1 is maxr(0, roundtoward(Xl**R,to_negative)),  % neg Xl ==> nan
  550	Zh1 is roundtoward(Xh**R,to_positive),
  551	Zpl is maxr(Zl,Zl1),  Zph is minr(Zh,Zh1),       % positive case
  552	Znl is maxr(Zl,-Zh1), Znh is minr(Zh,-Zl1),      % negative case
  553	( 1 is cmpr(Znl,Znh) -> NZl is Zpl, NZh is Zph   % Znl > Znh -> negative case is empty
  554	; 1 is cmpr(Zpl,Zph) -> NZl is Znl, NZh is Znh   % Zpl > Zph -> positive case is empty
  555	; NZl is minr(Zpl,Znl), NZh is maxr(Zph,Znh)     % else union of positive and negative cases
  556	).
  557
  558pt_powrCase(0,1,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :-  % R positive, numerator even, denominator odd
  559	% reflect about Z axis
  560	( Xh < 0 -> Xl1 is -Xh, Xh1 is -Xl           % negative X, negate interval
  561	;(Xl > 0 -> Xl1 = Xl, Xh1 = Xh               % positive X, As is
  562	;           Xl1 = 0, Xh1 is maxr(-Xl,Xh)     % split
  563	)),
  564	% NZl can't be negative
  565	NZl is maxr(Zl, roundtoward(Xl1**R,to_negative)),  % Xl1 known to be positive
  566	NZh is minr(Zh, roundtoward(Xh1**R,to_positive)).
  567
  568pt_powrCase(1,1,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :-  % R positive, numerator odd, denominator odd
  569	% continuous monotonic
  570	NZl is maxr(Zl, roundtoward(Xl**R,to_negative)),
  571	NZh is minr(Zh, roundtoward(Xh**R,to_positive)).
  572
  573% Y is an interval
  574powr_intY_(p, Z,X,Y, NewZ,NewX,NewY) :-              % positive X, interval Y 
  575	powr_prim_(Z,X,Y,NewZ), non_empty(NewZ),         % NewZ := X**Y
  576	universal_interval(UI),
  577	intCase(Cy,Y),
  578	odivCase(p,Cy,(1,1),Y,UI,Yi),                    % Yi := UI ^ 1/Y
  579	powr_prim_(X,NewZ,Yi,NewX), non_empty(NewX),     % NewX := NewZ**(1/Y)
  580	ln(NewZ,Ynum), intCase(Cyn,Ynum),
  581	ln(NewX,Yden), intCase(Cyd,Yden),
  582	odivCase(Cyn,Cyd,Ynum,Yden,Y,NewY).              % NewY := Y ^ log(NewZ)/log(NewX)
  583
  584powr_intY_(n, Z,X,Y, NewZ,NewX,NewY) :-              % negative X, interval Y
  585	% In this case Y may contain an rational which is legal for negative X, e.g.,
  586	% to calculate an odd power or root. However it still may be able to bound Z.
  587	% This case flips X to positive, calculates upper bound Zh1, and then intersects
  588	% Z with (-Zh1,Zh1). 
  589	X = (Xl,Xh), X1l is -Xh, X1h is -Xl,
  590	powr_intY_(p, Z,(X1l,X1h),Y, (_,Z1h),(X2l,X2h),NewY),  % use code for positive case
  591	X3l is -X2h, X3h is -X2l, ^(X,(X3l,X3h),NewX),    % flip back for NewX
  592	Z1l is -Z1h, ^(Z,(Z1l,Z1h),NewZ).                 % maximal range = (-Z1h,Z1h)  
  593
  594powr_intY_(s, Z,X,Y, NewZ,NewX,NewY):-                % split X, interval Y
  595	% union of positive and negative regions
  596	X = (Xl,Xh), Xfl is -Xl,  % flip negative to positive for calculation
  597	powr_intY_(p, Z,(0,Xfl),Y, (_,Znh),NXn,NYn),      % Znl always 0,
  598	powr_intY_(p, Z,(0,Xh), Y, (_,Zph),NXp,NYp),      % Zpl always 0
  599	% union and intersect with original
  600	NZl is -Znh, NZh is maxr(Znh,Zph), ^(Z,(NZl,NZh),NewZ),
  601	union_(NXn,NXp,NX), ^(X,NX,NewX),
  602	union_(NYn,NYp,NY), ^(Y,NY,NewY).
  603	
  604powr_prim_((Zl,Zh), (Xl,Xh), (Yl,Yh), (NZl,NZh)) :-  % X assumed positive
  605	Zll is roundtoward(Xl**Yl,to_negative),
  606	Zlh is roundtoward(Xl**Yh,to_positive),
  607	Zhl is roundtoward(Xh**Yl,to_negative),
  608	Zhh is roundtoward(Xh**Yh,to_positive),
  609	NZl is maxr(Zl,minr(Zll, minr(Zlh, minr(Zhl, Zhh)))), % intersect with min of all
  610	NZh is minr(Zh,maxr(Zll, maxr(Zlh, maxr(Zhl, Zhh)))). % intersect with max of all
  611
  612ln((Xl,Xh), (Zl,Zh)) :-
  613	Zl is roundtoward(log(Xl),to_negative),
  614	Zh is roundtoward(log(Xh),to_positive).
  615
  616%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  617
  618% Z== sin(X)
  619
  620narrowing_op(sin, _, $(Z,X), $(NewZ, NewX)) :-
  621	sin_(X,S,Z1), ^(Z,Z1,NewZ),
  622	asin_(NewZ,S,X1), ^(X,X1,NewX).
  623
  624sin_((X,X),S,NZ) :-  % point infinity special case
  625	abs(X) =:= 1.0Inf, !,
  626	S = (-10,10), NZ = (-1.0,1.0). 
  627sin_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
  628	Sl is round(Xl/pi),                    % determine monotonic sector of width pi
  629	Sh is round(Xh/pi),                    % sector 0 is (-pi/2,pi/2)
  630	to_float((Xl,Xh),(FXl,FXh)),           % convert now to float
  631	sinCase(FXl,FXh,Sl,Sh,NZl,NZh).
  632
  633sinCase(Xl,Xh,S,S,NZl,NZh) :- !,           % same sector
  634	( 0 is S mod 2                         % monotonically increasing ?
  635	 -> NZl is roundtoward(sin(Xl), to_negative),
  636	    NZh is roundtoward(sin(Xh), to_positive)
  637	 ;  NZl is roundtoward(sin(Xh), to_negative),
  638	    NZh is roundtoward(sin(Xl), to_positive)
  639	).
  640sinCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !,  % adjacent sectors 
  641	( 0 is Sl mod 2                        % monotonically increasing ?
  642	 -> NZl is minr(roundtoward(sin(Xl), to_negative),  % Z contains +1
  643	                roundtoward(sin(Xh), to_negative)),
  644	    NZh =  1.0          
  645	 ; 	NZl = -1.0,                                      % Z contains -1
  646	    NZh is maxr(roundtoward(sin(Xl), to_positive),
  647	                roundtoward(sin(Xh), to_positive))
  648	).
  649sinCase(_,_,_,_,-1.0,1.0).                  % span multiple sectors
  650
  651asin_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
  652	asinCase(Xl,Xh,Sl,Sh,NZl,NZh).
  653
  654% Need to maintain adequate width applying offset, otherwise loss of precision
  655%	leads to unintentional rounding. This is done by fuzzing offset. 
  656asinCase(Xl,Xh,S,S,NZl,NZh) :- !,           % same sector
  657	Zl is roundtoward(asin(Xl),to_negative),
  658	Zh is roundtoward(asin(Xh),to_positive),
  659	Offs is S*pi,
  660	bounded_number(Offl,Offh,Offs),         % fuzz offset
  661	( 0 is S mod 2                          % monotonically increasing ?
  662	 -> NZl is roundtoward(Offl+Zl,to_negative),
  663	    NZh is roundtoward(Offh+Zh,to_positive)
  664	 ;  NZl is roundtoward(Offl-Zh,to_negative),
  665	    NZh is roundtoward(Offh-Zl,to_positive)
  666	).
  667asinCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !,  % adjacent sector
  668	Offl is nexttoward(Sl*pi,-inf), Offh is nexttoward(Sh*pi,inf),  % fuzz offset
  669	( 0 is Sl mod 2                         % monotonically increasing ?
  670	 -> Z is roundtoward(asin(Xl),to_negative),
  671	    NZl is roundtoward(Offl+Z,to_negative),
  672	    NZh is roundtoward(Offh-Z,to_positive)
  673	 ;  Z is roundtoward(asin(Xh),to_negative),
  674	    NZl is roundtoward(Offl-Z,to_negative),
  675	    NZh is roundtoward(Offh+Z,to_positive)
  676	).
  677asinCase(_,_,_,_,-1.0Inf,1.0Inf).           % span multiple sectors
  678
  679%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  680
  681% Z== cos(X)
  682
  683narrowing_op(cos, _, $(Z,X), $(NewZ, NewX)) :-
  684	cos_(X,S,Z1), ^(Z,Z1,NewZ),
  685	acos_(NewZ,S,X1), ^(X,X1,NewX).
  686
  687cos_((X,X),S,NZ) :-  % point infinity special case
  688	abs(X) =:= 1.0Inf, !,
  689	S = (-10,10), NZ = (-1.0,1.0). 
  690cos_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
  691	Sl is floor(Xl/pi),
  692	Sh is floor(Xh/pi),  % sector 0 is (0,pi)
  693	to_float((Xl,Xh),(FXl,FXh)),             % convert now to float
  694	cosCase(FXl,FXh,Sl,Sh,NZl,NZh).
  695
  696cosCase(Xl,Xh,S,S,NZl,NZh) :- !,             % same sector
  697	( 1 is S mod 2                           % monotonically increasing ?
  698	 -> NZl is roundtoward(cos(Xl), to_negative),
  699	    NZh is roundtoward(cos(Xh), to_positive)
  700	 ;  NZl is roundtoward(cos(Xh), to_negative),
  701	    NZh is roundtoward(cos(Xl), to_positive)
  702	).
  703cosCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !,  % adjacent sector
  704	( 1 is Sl mod 2                          % monotonically increasing ?
  705	 -> NZl is minr(roundtoward(cos(Xl), to_negative),  % Z contains +1
  706	                roundtoward(cos(Xh), to_negative)),
  707	    NZh =  1.0          
  708	 ; 	NZl = -1.0,                                      % Z contains -1
  709	    NZh is maxr(roundtoward(cos(Xl), to_positive),
  710	                roundtoward(cos(Xh), to_positive))
  711	).
  712cosCase(_,_,_,_,-1.0,1.0).                  % span multiple sectors
  713
  714acos_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
  715	acosCase(Xl,Xh,Sl,Sh,NZl,NZh).
  716
  717% Need to maintain adequate width applying offset, otherwise loss of precision
  718%	leads to unintentional rounding. This is done by fuzzing offset. 
  719acosCase(Xl,Xh,S,S,NZl,NZh) :- !,           % same sector
  720	Zl is roundtoward(acos(Xh),to_negative),
  721	Zh is roundtoward(acos(Xl),to_positive),
  722	( 1 is S mod 2                          % monotonically increasing ?
  723	 -> Offs is (S+1)*pi,
  724	 	bounded_number(Offl,Offh,Offs),     % fuzz offset
  725	    NZl is roundtoward(Offl-Zh,to_negative),
  726        NZh is roundtoward(Offh-Zl,to_positive)               
  727	 ;  Offs is S*pi,
  728	 	bounded_number(Offl,Offh,Offs),     % fuzz offset
  729	    NZl is roundtoward(Offl+Zl,to_negative),
  730	    NZh is roundtoward(Offh+Zh,to_positive)
  731  	).
  732acosCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !,  % adjacent sector
  733	( 1 is Sl mod 2                         % monotonically increasing ?
  734	 -> Z is roundtoward(acos(Xl),to_positive),  % acos is monotone decreasing
  735	    Offl is nexttoward((Sl+1)*pi,-inf), NZl is roundtoward(Offl-Z,to_negative),
  736	    Offh is nexttoward(Sh*pi,inf) , NZh is roundtoward(Offh+Z,to_positive)
  737	 ;  Z is roundtoward(acos(Xh),to_positive),  % acos is monotone decreasing
  738	    Offl is nexttoward(Sl*pi,-inf) , NZl is roundtoward(Offl+Z,to_negative),
  739	    Offh is nexttoward((Sh+1)*pi,inf), NZh is roundtoward(Offh-Z,to_positive)	
  740  	).
  741acosCase(_,_,_,_,-1.0Inf,1.0Inf).           % span multiple sectors
  742
  743%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  744
  745% Z== tan(X)
  746
  747narrowing_op(tan, _, $(Z,X), $(NewZ, NewX)) :-
  748	X = (Xl,Xh),
  749	Sl is round(Xl/pi),  % sector 0 is (-pi/2,pi/2)
  750	Sh is round(Xh/pi),
  751	tanCase(Z,X,Sl,Sh,NewZ,NewX).
  752	
  753% Need to maintain adequate width applying offset, otherwise loss of precision
  754%	leads to unintentional rounding. This is done by fuzzing offset. 
  755tanCase(Z,(X,X),_,_,Z,(X,X)) :-    % point infinity special case
  756	abs(X) =:= 1.0Inf, !.
  757tanCase((Zl,Zh),(Xl,Xh),S,S,(NZl,NZh),(NXl,NXh)) :-  !,   % same sector
  758	Z1l is roundtoward(tan(Xl),to_negative),
  759	Z1h is roundtoward(tan(Xh),to_positive),
  760	( non_empty(Z1l,Z1h)
  761	 -> ^((Zl,Zh),(Z1l,Z1h),(NZl,NZh)),                   % good range
  762	 	Offs is S*pi,
  763	 	bounded_number(Offl,Offh,Offs),                   % fuzz offset
  764	    X1l is Offl+roundtoward(atan(NZl),to_negative),
  765	    X1h is Offh+roundtoward(atan(NZh),to_positive),
  766	    ^((Xl,Xh),(X1l,X1h),(NXl,NXh))
  767	 ;  % same sector range check failed, assume rounding crossed singularity
  768	    ((Z1l =< Zh ; Z1h >= Zl) -> true),                % disjoint ranges must intersect Z
  769	    NZl = Zl, NZh = Zh, NXl = Xl, NXh = Xh            % nothing changes
  770	).
  771tanCase(Z,(Xl,Xh),Sl,Sh,(NZl,NZh),(NXl,NXh)) :-  Sh-Sl =:= 1, !, % spans one singularity
  772	Sg is pi*(Sl+1r2),               % closest to singularity
  773	bounded_number(XSl,XSh,Sg),      % use bounded number to straddle
  774	(tanCase(Z,(Xl,XSl),Sl,Sl,(NZ1l,NZ1h),(NXl,NX1h))
  775	 -> (tanCase(Z,(XSh,Xh),Sh,Sh,_NZ2,(_,NXh))
  776	     -> NZl = -1.0Inf, NZh = 1.0Inf                   % both sectors, includes singularity
  777	     ;  NXh = NX1h, NZl = NZ1l, NZh = NZ1h            % only lo sector
  778	    )
  779	 ;  tanCase(Z,(XSh,Xh),Sh,Sh,(NZl,NZh),(NXl,NXh))     % only hi sector
  780	).
  781tanCase(Z,X,Sl,Sh,Z,X).                                   % spans multiple singularities
  782
  783%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  784
  785% Z== ~X boolean negation (Z and X boolean)
  786
  787narrowing_op(not, P, $(Z,X), $(NewZ, NewX)) :-
  788	notB_(X,Z,Z1,P), ^(Z,Z1,NewZ),
  789	notB_(NewZ,X,X1,P), ^(X,X1,NewX).
  790
  791notB_((0,0), _, (1,1), p) :- !.
  792notB_((1,1), _, (0,0), p) :- !.
  793notB_(    _, Z, NewZ,  _) :- ^(Z,(0,1),NewZ).
  794
  795%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  796
  797% Z==X and Y  boolean 'and'
  798
  799narrowing_op(and, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  800	andB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
  801
  802andB_rel_(Z,(1,1),(1,1), NewZ,(1,1),(1,1), p) :- !, ^(Z,(1,1),NewZ).  % all cases persistent (points) except last
  803andB_rel_(Z,(0,0),Y,     NewZ,(0,0),Y, p)     :- !, ^(Z,(0,0),NewZ), ^(Y,(0,1),NewY).
  804andB_rel_(Z,X,(0,0),     NewZ,X,(0,0), p)     :- !, ^(Z,(0,0),NewZ), ^(X,(0,1),NewX).
  805andB_rel_((1,1),X,Y,     (1,1),NewX,NewY, p)  :- !, ^(X,(1,1),NewX), ^(Y,(1,1),NewY).
  806andB_rel_((0,0),X,(1,1), (0,0),NewX,(1,1), p) :- !, ^(X,(0,0),NewX).
  807andB_rel_((0,0),(1,1),Y, (0,0),(1,1),NewY, p) :- !, ^(Y,(0,0),NewY).
  808andB_rel_(Z,X,Y,         NewZ,NewX,NewY, P)   :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
  809
  810%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  811
  812% Z==X and Y  boolean 'nand'
  813
  814narrowing_op(nand, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  815	nandB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
  816
  817nandB_rel_(Z,(1,1),(1,1), NewZ,(1,1),(1,1), p) :- !, ^(Z,(0,0),NewZ).  % all cases persistent except last
  818nandB_rel_(Z,(0,0),Y,     NewZ,(0,0),Y, p)     :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY).
  819nandB_rel_(Z,X,(0,0),     NewZ,X,(0,0), p)     :- !, ^(Z,(1,1),NewZ), ^(X,(0,1),NewX).
  820nandB_rel_((0,0),X,Y,     (0,0),NewX,NewY, p)  :- !, ^(X,(1,1),NewX), ^(Y,(1,1),NewY).
  821nandB_rel_((1,1),X,(1,1), (1,1),NewX,(1,1), p) :- !, ^(X,(0,0),NewX).
  822nandB_rel_((1,1),(1,1),Y, (1,1),(1,1),NewY, p) :- !, ^(Y,(0,0),NewY).
  823nandB_rel_(Z,X,Y,         NewZ,NewX,NewY, P)   :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
  824
  825%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  826
  827% Z==X or Y  boolean 'or'
  828
  829narrowing_op(or, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  830	orB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
  831
  832orB_rel_(Z,(0,0),(0,0), NewZ,(0,0),(0,0), p) :- !, ^(Z,(0,0),NewZ).  % all cases persistent (points) except last
  833orB_rel_(Z,(1,1),Y,     NewZ,(1,1),Y, p)     :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY).
  834orB_rel_(Z,X,(1,1),     NewZ,X,(1,1), p)     :- !, ^(Z,(1,1),NewZ), ^(X,(0,1),NewX).
  835orB_rel_((0,0),X,Y,     (0,0),NewX,NewY, p)  :- !, ^(X,(0,0),NewX), ^(Y,(0,0),NewY).
  836orB_rel_((1,1),X,(0,0), (1,1),NewX,(0,0), p) :- !, ^(X,(1,1),NewX).
  837orB_rel_((1,1),(0,0),Y, (1,1),(0,0),NewY, p) :- !, ^(Y,(1,1),NewY).
  838orB_rel_(Z,X,Y,         NewZ,NewX,NewY, P)   :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
  839
  840%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  841
  842% Z==X nor Y  boolean 'nor'
  843
  844narrowing_op(nor, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  845	norB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
  846
  847norB_rel_(Z,(0,0),(0,0), NewZ,(0,0),(0,0), p) :- !, ^(Z,(1,1),NewZ).  % all cases persistent (points) except last
  848norB_rel_(Z,(1,1),Y,     NewZ,(1,1),Y, p)     :- !, ^(Z,(0,0),NewZ), ^(Y,(0,1),NewY).
  849norB_rel_(Z,X,(1,1),     NewZ,X,(1,1), p)     :- !, ^(Z,(0,0),NewZ), ^(X,(0,1),NewX).
  850norB_rel_((1,1),X,Y,     (1,1),NewX,NewY, p)  :- !, ^(X,(0,0),NewX), ^(Y,(0,0),NewY).
  851norB_rel_((0,0),X,(0,0), (0,0),NewX,(0,0), p) :- !, ^(X,(1,1),NewX).
  852norB_rel_((0,0),(0,0),Y, (0,0),(0,0),NewY, p) :- !, ^(Y,(1,1),NewY).
  853norB_rel_(Z,X,Y,         NewZ,NewX,NewY, P)   :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
  854
  855%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  856
  857% Z==X xor Y  boolean 'xor'
  858
  859narrowing_op(xor, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  860	xorB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
  861
  862xorB_rel_(Z,(B,B),(B,B), NewZ,(B,B),(B,B), p) :- !, ^(Z,(0,0),NewZ).  % all cases persistent (points) except last
  863xorB_rel_(Z,(X,X),(Y,Y), NewZ,(X,X),(Y,Y), p) :- !, ^(Z,(1,1),NewZ).
  864xorB_rel_((B,B),X,(B,B), (B,B),NewX,(B,B), p) :- !, ^(X,(0,0),NewX).
  865xorB_rel_((Z,Z),X,(Y,Y), (Z,Z),NewX,(Y,Y), p) :- !, ^(X,(1,1),NewX).
  866xorB_rel_((B,B),(B,B),Y, (B,B),(B,B),NewY, p) :- !, ^(Y,(0,0),NewY).
  867xorB_rel_((Z,Z),(X,X),Y, (Z,Z),(X,X),NewY, p) :- !, ^(Y,(1,1),NewY).
  868xorB_rel_(Z,X,Y,         NewZ,NewX,NewY, P)   :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
  869
  870%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  871
  872% Z==X -> Y  boolean 'implies'
  873
  874narrowing_op(imB, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
  875	imB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
  876
  877imB_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
  878imB_rel_(Z,(0,0),Y, NewZ,(0,0),Y,    p) :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY).        % X=0, Z=1, persistant
  879imB_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
  880imB_rel_(Z,X,(1,1), NewZ,X,(1,1),    P) :- !, ^(Z,(1,1),NewZ).                      % Y=1, Z=1, persistant
  881imB_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
  882imB_rel_(Z,X,Y,     NewZ,NewX,NewY,  P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY).  % still indeterminate
  883
  884%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  885
  886:- style_check(+singleton).  % for narrowing_op
  887
  888% User defined IA primitives
  889
  890narrowing_op(user(Prim), P, InArgs, OutArgs) :-
  891	call_user_primitive(Prim, P, InArgs, OutArgs).  % from main module body