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