diff --git a/docs/CLP_BNR_Guide/CLP_BNR_Guide.html b/docs/CLP_BNR_Guide/CLP_BNR_Guide.html old mode 100755 new mode 100644 diff --git a/docs/CLP_BNR_Guide/CLP_BNR_Guide.myw b/docs/CLP_BNR_Guide/CLP_BNR_Guide.myw old mode 100755 new mode 100644 diff --git a/docs/CLP_BNR_Guide/lib/commonmark.min.js b/docs/CLP_BNR_Guide/lib/commonmark.min.js old mode 100755 new mode 100644 diff --git a/docs/CLP_BNR_Guide/lib/grit.js b/docs/CLP_BNR_Guide/lib/grit.js old mode 100755 new mode 100644 diff --git a/docs/Package/clpBNR-0.7.4.zip b/docs/Package/clpBNR-0.7.4.zip deleted file mode 100644 index 814c941..0000000 Binary files a/docs/Package/clpBNR-0.7.4.zip and /dev/null differ diff --git a/docs/Package/clpBNR-0.8.1.zip b/docs/Package/clpBNR-0.8.1.zip deleted file mode 100644 index 0bc8a8a..0000000 Binary files a/docs/Package/clpBNR-0.8.1.zip and /dev/null differ diff --git a/docs/Package/clpBNR-0.8.2.zip b/docs/Package/clpBNR-0.8.2.zip deleted file mode 100644 index 4501e28..0000000 Binary files a/docs/Package/clpBNR-0.8.2.zip and /dev/null differ diff --git a/docs/Package/clpBNR-0.8.zip b/docs/Package/clpBNR-0.8.zip deleted file mode 100644 index f1088de..0000000 Binary files a/docs/Package/clpBNR-0.8.zip and /dev/null differ diff --git a/docs/Package/clpBNR-0.9.0.zip b/docs/Package/clpBNR-0.9.0.zip deleted file mode 100644 index b6817ed..0000000 Binary files a/docs/Package/clpBNR-0.9.0.zip and /dev/null differ diff --git a/docs/Package/clpBNR-0.9.1.zip b/docs/Package/clpBNR-0.9.1.zip deleted file mode 100644 index 68384e0..0000000 Binary files a/docs/Package/clpBNR-0.9.1.zip and /dev/null differ diff --git a/docs/Package/clpBNR-0.9.2.zip b/docs/Package/clpBNR-0.9.2.zip deleted file mode 100644 index ed6e428..0000000 Binary files a/docs/Package/clpBNR-0.9.2.zip and /dev/null differ diff --git a/docs/Package/clpBNR-0.9.3.zip b/docs/Package/clpBNR-0.9.3.zip deleted file mode 100644 index e10a056..0000000 Binary files a/docs/Package/clpBNR-0.9.3.zip and /dev/null differ diff --git a/docs/Package/clpBNR-0.9.4.zip b/docs/Package/clpBNR-0.9.4.zip deleted file mode 100644 index e95567c..0000000 Binary files a/docs/Package/clpBNR-0.9.4.zip and /dev/null differ diff --git a/docs/Package/clpBNR-0.9.4/prolog/clpBNR.pl b/docs/Package/clpBNR-0.9.4/prolog/clpBNR.pl deleted file mode 100755 index f804a44..0000000 --- a/docs/Package/clpBNR-0.9.4/prolog/clpBNR.pl +++ /dev/null @@ -1,825 +0,0 @@ -% -% CLP(BNR) == Constraints On Boolean, Integer, and Real Intervals -% -/* The MIT License (MIT) - * - * Copyright (c) 2019,2020 Rick Workman - * - * Permission is hereby granted, free of charge, to any person obtaining a copy - * of this software and associated documentation files (the "Software"), to deal - * in the Software without restriction, including without limitation the rights - * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the Software is - * furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included in all - * copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - * SOFTWARE. - */ -:- module(clpBNR, % SWI module declaration - [ - op(700, xfx, ::), - op(150, xf, ...), % postfix op currently just for output - (::)/2, % declare interval - {}/1, % define constraint - interval/1, % filter for clpBNR constrained var - list/1, % for compatibility - domain/2, range/2, % for compatibility - delta/2, % width (span) of an interval or numeric (also arithmetic function) - midpoint/2, % midpoint of an interval (or numeric) (also arithmetic function) - median/2, % median of an interval (or numeric) (also arithmetic function) - lower_bound/1, % narrow interval to point equal to lower bound - upper_bound/1, % narrow interval to point equal to upper bound - - % additional constraint operators - op(200, fy, ~), % boolean 'not' - op(500, yfx, and), % boolean 'and' - op(500, yfx, or), % boolean 'or' - op(500, yfx, nand), % boolean 'nand' - op(500, yfx, nor), % boolean 'nor' - op(500, yfx, xor), % boolean 'xor' - op(700, xfx, <>), % integer not equal - op(700, xfx, <=), % set inclusion - op(700, xfx, =>), % set inclusion - - % utilities - print_interval/1, print_interval/2, % pretty print interval with optional stream - small/1, small/2, % defines small interval width based on precision value - solve/1, solve/2, % solve (list of) intervals using split to find point solutions - splitsolve/1, splitsolve/2, % solve (list of) intervals using split - absolve/1, absolve/2, % absolve (list of) intervals, narrows by nibbling bounds - enumerate/1, % "enumerate" integers - global_minimum/2, % find interval containing global minimums for an expression - global_minimum/3, % global_minimum/2 with definable precision - global_maximum/2, % find interval containing global minimums for an expression - global_maximum/3, % global_maximum/2 with definable precision - nb_setbounds/2, % non-backtracking set bounds (use with branch and bound) - partial_derivative/3, % differentiate Exp wrt. X and simplify - clpStatistics/0, % reset - clpStatistic/1, % get selected - clpStatistics/1, % get all defined in a list - watch/2 % enable monitoring of changes for interval or (nested) list of intervals - ]). - -/* missing(?) functionality: utility accumulate/2. */ - -/* supported interval relations: - -+ - * / %% arithmetic -** %% includes real exponent, odd/even integer -abs %% absolute value -sqrt %% positive square root -min max %% binary min/max -== is <> =< >= < > %% comparison -<= => %% inclusion -and or nand nor xor -> %% boolean -- ~ %% unary negate and not -exp log %% exp/ln -sin asin cos acos tan atan %% trig functions - -*/ - -:- use_module(library(arithmetic)). % for interval arithmetic functions -:- use_module(library(lists),[subtract/3,union/3]). % for flags - -:- style_check([-singleton, -discontiguous]). % :- discontiguous ... not reliable. - -% -% SWIP optimise control - set flag to true for compiled arithmetic -% -:- (current_prolog_flag(optimise,Opt), - nb_setval('clpBNR:temp',Opt), % save current value to restore on :- initialization/1 - set_prolog_flag(optimise,false) - ). -% -% Define debug_clpBNR_/2 before turning on optimizer removing debug calls. -% -debug_clpBNR_(FString,Args) :- debug(clpBNR,FString,Args). - -:- set_prolog_flag(optimise,true). - -current_node_(Node) :- % look back to find current Op being excuted for debug messages - prolog_current_frame(F), % this is a little grungy, but necessary to get intervals - prolog_frame_attribute(F,parent_goal,doNode_(Arg,Op,_,_,_,_)), - Arg =.. [_|Args], - Node=..[Op|Args]. - -% -% statistics -% - -% assign,increment/read global counter (assumed to be ground value so use _linkval) -g_assign(G,V) :- nb_linkval(G,V). -g_inc(G) :- nb_getval(G,N), N1 is N+1, nb_linkval(G,N1). -g_incb(G) :- nb_getval(G,N), N1 is N+1, b_setval(G,N1). % undone on backtrack -g_read(G,V) :- nb_getval(G,V). - -:- discontiguous - clpBNR:clpStatistics/0, clpBNR:clpStatistic/1, - sandbox:safe_global_variable/1. - -sandbox:safe_global_variable('clpBNR:userTime'). -sandbox:safe_global_variable('clpBNR:inferences'). -sandbox:safe_global_variable('clpBNR:gc_time'). - -clpStatistics :- - % garbage_collect, % ? do gc before time snapshots - statistics(cputime,T), g_assign('clpBNR:userTime',T), % thread based - statistics(inferences,I), g_assign('clpBNR:inferences',I), - statistics(garbage_collection,[_,_,G,_]), g_assign('clpBNR:gc_time',G), - fail. % backtrack to reset other statistics. - -clpStatistic(userTime(T)) :- statistics(cputime,T1), g_read('clpBNR:userTime',T0), T is T1-T0. - -clpStatistic(gcTime(G)) :- statistics(garbage_collection,[_,_,G1,_]), g_read('clpBNR:gc_time',G0), G is (G1-G0)/1000.0. - -clpStatistic(globalStack(U/T)) :- statistics(globalused,U), statistics(global,T). - -clpStatistic(trailStack(U/T)) :- statistics(trailused,U), statistics(trail,T). - -clpStatistic(localStack(U/T)) :- statistics(localused,U), statistics(local,T). - -clpStatistic(inferences(I)) :- statistics(inferences,I1), g_read('clpBNR:inferences',I0), I is I1-I0. - - -:- include(ia_primitives). % interval arithmetic relations via evalNode/4. - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% SWI-Prolog implementation of IA -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% Intervals are constrained (attributed) variables. -% -% Current interval bounds updates via setarg(Val) which is backtrackable -% -% interval(X) - filter -% -interval(X) :- get_attr(X, clpBNR, _). - -% internal abstraction -interval_object(Int, Type, Val, Nodelist) :- - get_attr(Int, clpBNR, interval(Type, Val, Nodelist, Flags)). - -% flags (list) abstraction -get_interval_flags_(Int,Flags) :- - get_attr(Int, clpBNR, interval(Type, Val, Nodelist, Flags)). - -set_interval_flags_(Int,Flags) :- % flags assumed to be ground so no copy required - get_attr(Int, clpBNR, interval(Type, Val, Nodelist, _)), - put_attr(Int, clpBNR, interval(Type, Val, Nodelist, Flags)). - -% -% Interval value constants -% -universal_interval((-1.0Inf,1.0Inf)). - -% Finite intervals - 64 bit IEEE reals, -finite_interval(real, (-1.0e+16,1.0e+16)). -finite_interval(integer, (L,H)) :- %% SWIP: - current_prolog_flag(bounded,false),!, % integers are unbounded, but use tagged limits for finite default - current_prolog_flag(min_tagged_integer,L), - current_prolog_flag(max_tagged_integer,H). -%finite_interval(boolean, (0,1)). - -% Empty (L>H) -%empty_interval([L,H]) :- universal_interval([H,L]). - -% -% non-backtracking set bounds for use with branch and bound -% -nb_setbounds(Int, [L,U]) :- - get_attr(Int, clpBNR, Def), - Def = interval(_, Val, _, _), - ^(Val,(L,U),NewVal), % new range is intersection (from ia_primitives) - nb_setarg(2, Def, NewVal). - -% -% get current value -% -getValue(Int, Val) :- - number(Int) - -> Val=(Int,Int) % numeric point value - ; get_attr(Int, clpBNR, interval(_, Val, _, _)). % interval, optimized for SWIP - -% -% set monitor action on an interval -% -watch(Int,Action) :- - atom(Action), - get_interval_flags_(Int,Flags), !, - lists:subtract(Flags,[watch(_)],Flags1), % remove any previous setting - (Action = none -> true ; set_interval_flags_(Int,[watch(Action)|Flags1])). -watch(Ints,Action) :- - is_list(Ints), - watch_list_(Ints,Action). - -watch_list_([],Action). -watch_list_([Int|Ints],Action) :- - watch(Int,Action), - watch_list_(Ints,Action). - -% check if watch enabled on this interval -check_monitor_(Int, Update, interval(Type,Val,Nodelist,Flags)) :- - (memberchk(watch(Action), Flags) - -> once(monitor_action_(Action, Update, Int)) % in ia_utilities - ; true - ). - -% -% set interval value (assumes bounds check has already been done) -% -putValue_(New, Int, NodeList) :- - get_attr(Int, clpBNR, Def), !, % still an interval - (debugging(clpBNR,true) -> check_monitor_(Int, New, Def) ; true), - Def = interval(Type,_,Nodes,_), % get Type and Nodes before setValue_ - setValue_(New,Int,Def), % set new value - queue_nodes_(Type,New,Int,Nodes,NodeList). % construct node list to schedule -putValue_((L,H), Num, NodeList) :- number(Num), % catch things like [0,-0.0] - L=:=Num, H=:=Num. - -setValue_((L,H),Int,Def) :- L=:=H, !, % narrowed to a point, unify with interval - setarg(3,Def,_NL), % clear node list (so nothing done in unify) - (rational(L) -> Int=L ; Int=H). % if either bound rational (precise), use it -setValue_(New,Int,Def) :- % update value in interval (backtrackable) - setarg(2,Def,New). - -queue_nodes_(real,_,_,Nodes,Nodes). % type real - just use Nodes -queue_nodes_(integer,(L,H),_,Nodes,Nodes) :- % type integer #1 with integral bounds - integral_(L), integral_(H), !. % run Nodes if bounds integers or infinities -queue_nodes_(integer,_,Int,_,[node(integral,_,0,$(Int))|_]). % type integer #2 - % schedule an "integral" node to re-adjust bounds (Nodes run on adjustment) - -integral_(1.0Inf). -integral_(-1.0Inf). -integral_(B) :- integer(B). - -% -% range(Int, Bounds) for compatability -% -% for interval(Int) and number(Int), check if value is (always) in specified range, unifying any vars with current value -range(Int, [L,H]) :- getValue(Int, (IL,IH)), number(IL), !, % existing interval - (var(L) -> L=IL ; L= H=IH ; IH= 0.0. % contains 0 (handles (-inf,inf) -median_(L,H,M) :- M is copysign(sqrt(abs(L))*sqrt(abs(H)),L). % L and H have same sign - -% -% lower_bound and upper_bound -% -lower_bound(Int) :- - getValue(Int,(L,H)), - Int=L. - -upper_bound(Int) :- - getValue(Int,(L,H)), - Int=H. - -% -% Interval declarations -% - -Ints::Dom :- is_list(Ints),!, - intervals_(Ints,Dom). - -R::Dom :- var(R), var(Dom), !, % declare R = real(L,H), Note: R can be interval - int_decl_(real,(_,_),R), - domain(R,Dom). - -R::Dom :- var(Dom), !, % domain query (if interval(R) else fail) - domain(R,Dom). % "domain" query, unify interval Type and Bounds - -R::Dom :- % interval(R) or number(R) and nonvar(Dom) - Dom=..[Type|Bounds], - (Bounds=[] -> Val=(_,_) ; Val=..[,|Bounds]), - int_decl_(Type,Val,R). - -int_decl_(boolean,_,R) :- !, % boolean is integer; 0=false, 1=true, ignore any bounds. - int_decl_(integer,(0,1),R). - -int_decl_(Type,(L,H),R) :- interval_object(R,IType,CVal,_NL), !, % already interval - lower_bound_val_(Type,L,IL), % changing type,bounds? - upper_bound_val_(Type,H,IH), - applyType_(Type, IType, R, T/T, Agenda), % coerce reals to integers (or no-op). - ^(CVal,(IL,IH),New), % low level functional primitive - updateValue_(CVal, New, R, 1, Agenda, NewAgenda), % update value (Agenda empty if no value change) - stable_(NewAgenda). % then execute Agenda - -int_decl_(Type,(L,H),R) :- var(R), !, % new interval (R can't be existing interval) - lower_bound_val_(Type,L,IL), - upper_bound_val_(Type,H,IH), - IL= - R=IL ; % point range, can unify now - put_attr(R, clpBNR, interval(Type, (IL,IH), _NL, [])) % attach clpBNR attribute - ). - -int_decl_(Type,(L,H),R) :- (Type=integer -> integer(R) ; number(R)), !, % R already a point value, check range - lower_bound_val_(Type,L,IL), IL= IL=Lv ; IL is nexttoward(Lv,-1.0Inf)). -lower_bound_val_(integer,L,IL) :- % integer: make integer, fail if inf - IL is ceiling(L), IL \= 1.0Inf. - -upper_bound_val_(Type,H,IH) :- var(H), !, % unspecified bound, make it finite - finite_interval(Type,(_,IH)). -upper_bound_val_(real,H,IH) :- % real: evaluate and round outward (if float) - Hv is H, Hv\= -1.0Inf, - (rational(Hv) -> IH=Hv ; IH is nexttoward(Hv,1.0Inf)). -upper_bound_val_(integer,H,IH) :- % integer: make integer, fail if -inf - IH is floor(H), IH \= -1.0Inf. - -applyType_(integer, real, Int, Agenda, NewAgenda) :- !, % narrow real to integer - get_attr(Int,clpBNR,interval(Type,Val,NodeList,Flags)), - (debugging(clpBNR,true) -> check_monitor_(Int, integer, interval(Type,Val,NodeList,Flags)) ; true), - Val = (L,H), - lower_bound_val_(integer,L,IL), - upper_bound_val_(integer,H,IH), - (IL=IH - -> Int=IL % narrowed to point - ; (put_attr(Int,clpBNR,interval(integer,(IL,IH),NodeList,Flags)), % set Type (only case allowed) - linkNodeList_(NodeList, Agenda, NewAgenda) - ) - ). -applyType_(Type,IType,Int, Agenda, Agenda). % anything else: no change - -% -% this goal gets triggered whenever an interval is unified, valid for a numeric value or another interval -% -attr_unify_hook(interval(Type,(L,H),Nodelist,Flags), V) :- % unify an interval with a numeric - (Type=integer -> integer(V) ; number(V)), % check that V is consistent with Type - L= monitor_unify_(interval(Type,(L,H),_,Flags), V) ; true), - linkNodeList_(Nodelist, T/T, Agenda), - stable_(Agenda). % broadcast change - -attr_unify_hook(interval(Type1,V1,Nodelist1,Flags1), Int) :- % unifying two intervals - get_attr(Int, clpBNR, interval(Type2,V2,Nodelist2,Flags2)), !, %%SWIP attribute def. - mergeType_(Type1, Type2, NewType), % unified Type=integer if either is an integer - ^(V1,V2,V), % unified range is intersection (from ia_primitives), - mergeNodes_(Nodelist2,Nodelist1,Newlist), % unified node list is a merge of two lists - mergeFlags_(Flags1,Flags2,Flags), - (debugging(clpBNR,true) -> monitor_unify_(interval(Type1,V1,_,Flags), Int) ; true), - % update new type, value and constraint list, undone on backtracking - put_attr(Int,clpBNR,interval(NewType,V,Newlist,Flags)), - linkNodeList_(Newlist, T/T, Agenda), - stable_(Agenda). % broadcast change - -attr_unify_hook(interval(Type,Val,Nodelist,Flags), V) :- % new value out of range - g_inc('clpBNR:evalNodeFail'), % count of primitive call failures - debugging(clpBNR, true), % fail immediately unless debug=true - debug_clpBNR_('Failed to unify ~w::~w with ~w',[Type,Val,V]), - fail. - -% awkward monitor case because original interval gone -monitor_unify_(IntDef, Update) :- % debbuging, manufacture temporary interval - put_attr(Temp,clpBNR,IntDef), - check_monitor_(Temp, Update, IntDef). - -% if both real, result type is real, else one is integer so result type integer -mergeType_(real, real, real) :- !. -mergeType_(_, _, integer). - -% optimize for one or both lists (dominant case) -mergeFlags_([],Flags2,Flags2) :- !. -mergeFlags_(Flags1,[],Flags1) :- !. -mergeFlags_(Flags1,Flags2,Flags) :- - lists:union(Flags1,Flags2,Flags). % ambiguous if both have same flag set to different values - -mergeNodes_([N],NodeList,NodeList) :- var(N),!. -mergeNodes_([N|Ns],NodeList,[N|NewList]) :- - N=node(Op,_,_,Ops), - notIn_(NodeList,Op,Ops), !, % test for equivalent node already in NodeList - mergeNodes_(Ns,NodeList,NewList). -mergeNodes_([N|Ns],NodeList,NewList) :- - mergeNodes_(Ns,NodeList,NewList). - -notIn_([node(NOp,_,_,NOps)|Ns],Op,Ops) :- % equivalent node(Op, _, _,Ops) ==> failure - NOp==Op, NOps==Ops, % avoid unification of ops and operands (Op may be more than just an atom) - !, fail. -notIn_([N|Ns],Op,Ops) :- % keep searching - nonvar(N), !, - notIn_(Ns,Op,Ops). -notIn_(_,_,_). % end of search - -% -% New Constraints use { ... } syntax. -% -{}. -{Cons} :- - addConstraints_(Cons,T/T,Agenda), % add constraints - stable_(Agenda). % then execute Agenda - -addConstraints_(C,Agenda,NewAgenda) :- - constraint_(C), !, % a constraint is a boolean expression that evaluates to true - simplify(C,CS), % optional - buildConstraint_(CS, Agenda, NewAgenda). -addConstraints_((C,Cs),Agenda,NewAgenda) :- % Note: comma as operator - nonvar(C), - addConstraints_(C,Agenda,NextAgenda), !, - addConstraints_(Cs,NextAgenda,NewAgenda). -addConstraints_([],Agenda,Agenda). -addConstraints_([C|Cs],Agenda,NewAgenda) :- - nonvar(C), - addConstraints_(C,Agenda,NextAgenda), !, - addConstraints_(Cs,NextAgenda,NewAgenda). - -% low overhead version for internal use -constrain_(C) :- - buildConstraint_(C,T/T,Agenda), - stable_(Agenda). - -buildConstraint_(C,Agenda,NewAgenda) :- - debug_clpBNR_('Add ~p',{C}), - build_(C, 1, boolean, Agenda, NewAgenda), !. - -:- include(ia_simplify). % simplifies constraints to a hopefully more efficient equivalent - -% -% build a node from an expression -% -build_(Int, Int, VarType, Agenda, NewAgenda) :- % existing interval object - interval_object(Int, Type, _, _), !, - applyType_(VarType, Type, Int, Agenda, NewAgenda). % coerces exsiting intervals to required type -build_(Var, Var, VarType, Agenda, Agenda) :- % implicit interval creation. - var(Var), !, - universal_interval(UI), - int_decl_(VarType,UI,Var). -build_(Num, Num, _, Agenda, Agenda) :- % rational numeric value is precise - rational(Num), !. -build_(Num, Int, VarType, Agenda, Agenda) :- % floating point constant, may not be precise - float(Num), !, - int_decl_(VarType,(Num,Num),Int). % may be fuzzed, so not a point -build_(pt(Num), Num, _, Agenda, Agenda) :- % point value, must be a number - number(Num), !. -build_(Exp, Num, VarType, Agenda, Agenda) :- % pre-compile ground Exp (including pi and e) - ground(Exp), - safe_(Exp), % safe to evaluate - L is roundtoward(Exp,to_negative), % rounding only affects float evaluation - (float(L) % not precise -> interval - -> (H is roundtoward(Exp,to_positive), int_decl_(VarType,(L,H),Num)) % if precise use L - ; Num=L % else use precise value - ), !. -build_(Exp, Z, _, Agenda, NewAgenda) :- % deconstruct to primitives - Exp =.. [F|Args], - fmap_(F,Op,[Z|Args],ArgsR,Types), !, % supported arithmetic op - build_args_(ArgsR,Objs,Types,Agenda,ObjAgenda), - newNode_(Op,Objs,ObjAgenda,NewAgenda). - -build_args_([],[],_,Agenda,Agenda). -build_args_([Arg|ArgsR],[Obj|Objs],[Type|Types],Agenda,NewAgenda) :- - build_(Arg,Obj,Type,Agenda,NxtAgenda), - build_args_(ArgsR,Objs,Types,NxtAgenda,NewAgenda). - -% only called when argument is ground -safe_(E) :- atomic(E), !. % all atomics, including [] -safe_([A|As]) :- !, - safe_(A), - safe_(As). -safe_(F) :- - current_arithmetic_function(F), % evaluable by is/2 - F =.. [Op|Args], - \+memberchk(Op,[**,sin,cos,tan,asin,acos,atan]), % unsafe operations due to rounding - safe_(Args). - -% a constraint must evaluate to a boolean -constraint_(C) :- nonvar(C), C =..[Op|_], fmap_(Op,_,_,_,[boolean|_]), !. - -% map constraint operator to primitive/arity/types -fmap_(+, add, ZXY, ZXY, [real,real,real]). -fmap_(-, add, [Z,X,Y], [X,Z,Y], [real,real,real]). % note subtract before minus -fmap_(*, mul, ZXY, ZXY, [real,real,real]). -fmap_(/, mul, [Z,X,Y], [X,Z,Y], [real,real,real]). -fmap_(**, pow, ZXY, ZXY, [real,real,real]). -fmap_(min, min, ZXY, ZXY, [real,real,real]). -fmap_(max, max, ZXY, ZXY, [real,real,real]). -fmap_(==, eq, ZXY, ZXY, [boolean,real,real]). % strict equality -fmap_(=:=, eq, ZXY, ZXY, [boolean,real,real]). % Prolog compatible, strict equality -fmap_(is, eq, ZXY, ZXY, [boolean,real,real]). -fmap_(<>, ne, ZXY, ZXY, [boolean,integer,integer]). -fmap_(=\=, ne, ZXY, ZXY, [boolean,integer,integer]). % Prolog compatible -fmap_(=<, le, ZXY, ZXY, [boolean,real,real]). -fmap_(>=, le, [Z,X,Y], [Z,Y,X], [boolean,real,real]). -fmap_(<, lt, ZXY, ZXY, [boolean,real,real]). -fmap_(>, lt, [Z,X,Y], [Z,Y,X], [boolean,real,real]). -fmap_(<=, in, ZXY, ZXY, [boolean,real,real]). % inclusion/subinterval -fmap_(=>, in, [Z,X,Y], [Z,Y,X], [boolean,real,real]). % inclusion/subinterval - -fmap_(and, and, ZXY, ZXY, [boolean,boolean,boolean]). -fmap_(or, or, ZXY, ZXY, [boolean,boolean,boolean]). -fmap_(nand, nand, ZXY, ZXY, [boolean,boolean,boolean]). -fmap_(nor, nor, ZXY, ZXY, [boolean,boolean,boolean]). -fmap_(xor, xor, ZXY, ZXY, [boolean,boolean,boolean]). -fmap_(->, imB, ZXY, ZXY, [boolean,boolean,boolean]). - -fmap_(sqrt,sqrt, ZX, ZX, [real,real]). % pos. root version vs. **(1/2) -fmap_(-, minus, ZX, ZX, [real,real]). -fmap_(~, not, ZX, ZX, [boolean,boolean]). -fmap_(exp, exp, ZX, ZX, [real,real]). -fmap_(log, exp, [Z,X], [X,Z], [real,real]). -fmap_(abs, abs, ZX, ZX, [real,real]). -fmap_(sin, sin, ZX, ZX, [real,real]). -fmap_(asin,sin, [Z,X], [X,Z], [real,real]). -fmap_(cos, cos, ZX, ZX, [real,real]). -fmap_(acos,cos, [Z,X], [X,Z], [real,real]). -fmap_(tan, tan, ZX, ZX, [real,real]). -fmap_(atan,tan, [Z,X], [X,Z], [real,real]). - -% reverse map from Op and Args (used by "verbose" top level output to reverse compile constraints) -remap_(Op,$(Z,X,Y),C) :- constraint_(Op), Z==1, !, % simplification for constraints - C=..[Op,X,Y]. -remap_(Op,$(Z,X,Y),Z==C) :- !, - C=..[Op,X,Y]. -remap_(Op,$(Z,X),Z==C) :- - C=..[Op,X]. - -% -% Node constructor -% -newNode_(Op, Objs, Agenda, NewAgenda) :- - Args =.. [$|Objs], % store arguments as $/N where N=1..3 - NewNode = node(Op, P, 0, Args), % L=0 - addNode_(Objs,NewNode), - % increment count of added nodes, will be decremented on backtracking/failure - g_incb('clpBNR:node_count'), - linkNode_(Agenda, NewNode, NewAgenda). - -addNode_([],_Node). -addNode_([Arg|Args],Node) :- - (interval_object(Arg, _Type, _Val, Nodelist) -> newmember(Nodelist, Node) ; true), - addNode_(Args,Node). - -clpStatistics :- - g_assign('clpBNR:node_count',0), % reset/initialize node count to 0 - fail. % backtrack to reset other statistics. - -clpStatistic(node_count(C)) :- - g_read('clpBNR:node_count',C). - -% extend list with X -newmember([X|Xs],N) :- - nonvar(X), !, % not end of (indefinite) list - newmember(Xs,N). -newmember([N|_],N). % end of list - -% -% Process Agenda to narrow intervals (fixed point iteration) -% -stable_(Agenda) :- - current_prolog_flag(clpBNR_iteration_limit,Ops), % budget for current operation - stableLoop_(Agenda,Ops), - !. % achieved stable state with empty Agenda -> commit. - -stableLoop_([]/[], OpsLeft) :- !, % terminate successfully when agenda comes to an end - g_read('clpBNR:iterations',Cur), % maintain "low" water mark (can be negative) - (OpsLeft g_assign('clpBNR:iterations',OpsLeft);true). -stableLoop_([Node|Agenda]/T, OpsLeft) :- - Node = node(Op,P,_,Args), % if node on queue ignore link bit (was: Node = node(Op,P,1,Args)) - doNode_(Args, Op, P, OpsLeft, Agenda/T, NewAgenda), % undoable on backtrack - nb_setarg(3,Node,0), % reset linked bit - RemainingOps is OpsLeft-1, % decrement OpsLeft (can go negative) - stableLoop_(NewAgenda,RemainingOps). - -% support for max_iterations statistic -sandbox:safe_global_variable('clpBNR:iterations'). - -clpStatistics :- - current_prolog_flag(clpBNR_iteration_limit,L), - g_assign('clpBNR:iterations',L), % set "low" water mark to upper limit - fail. % backtrack to reset other statistics. - -clpStatistic(max_iterations(O/L)) :- - g_read('clpBNR:iterations',Ops), - current_prolog_flag(clpBNR_iteration_limit,L), - O is L-Ops. % convert "low" water mark to high water mark - -% -% Execute a node on the queue -% Note: "special" ops like instantiate and integral not counted as narrowing op in clpStatistics -% -% Comment out the following to enable Op tracing: -goal_expansion(traceIntOp_(Op, Args, PrimArgs, New),true). - -doNode_($(ZArg,XArg,YArg), Op, P, OpsLeft, Agenda, NewAgenda) :- % Arity 3 Op - var(P), !, % check persistent bit - % nonground(Args,_), !, % not safe, unifications happens before attr_unify_hook - getValue(ZArg,ZVal), - getValue(XArg,XVal), - getValue(YArg,YVal), - evalNode(Op, P, $(ZVal,XVal,YVal), $(NZVal,NXVal,NYVal)), % can fail causing stable_ to fail => backtracking - traceIntOp_(Op, [ZArg,XArg,YArg], [ZVal,XVal,YVal], [NZVal,NXVal,NYVal]), % in ia_utilities - updateValue_(ZVal, NZVal, ZArg, OpsLeft, Agenda, AgendaZ), - updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, AgendaZX), - updateValue_(YVal, NYVal, YArg, OpsLeft, AgendaZX, NewAgenda). - -doNode_($(ZArg,XArg), Op, P, OpsLeft, Agenda, NewAgenda) :- % Arity 2 Op - var(P), !, % check persistent bit - % nonground(Args,_), !, % not safe, unifications happens before attr_unify_hook - getValue(ZArg,ZVal), - getValue(XArg,XVal), - evalNode(Op, P, $(ZVal,XVal), $(NZVal,NXVal)), % can fail causing stable_ to fail => backtracking - traceIntOp_(Op, [ZArg,XArg], [ZVal,XVal], [NZVal,NXVal]), % in ia_utilities - updateValue_(ZVal, NZVal, ZArg, OpsLeft, Agenda, AgendaZ), - updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, NewAgenda). - -doNode_($(Arg), Op, P, OpsLeft, Agenda, NewAgenda) :- % Arity 1 Op - var(P), !, % check persistent bit - % nonground(Args,_), !, % not safe, unifications happens before attr_unify_hook - getValue(Arg,Val), - evalNode(Op, P, $(Val), $(NVal)), % can fail causing stable_ to fail => backtracking - traceIntOp_(Op, [Arg], [Val], [NVal]), % in ia_utilities - updateValue_(Val, NVal, Arg, 1, Agenda,NewAgenda). % always update value regardless of OpsLeft limiter - -doNode_(Args, Op, p, _, Agenda, Agenda) :- % persistent bit "set", skip node and trim - trim_ops_(Args). - -% -% called whenever a persistent node is encountered in FP iteration -% remove it from any node arguments -% -trim_ops_($(Op1,Op2,Op3)) :- - trim_op_(Op1), - trim_op_(Op2), - trim_op_(Op3). -trim_ops_($(Op1,Op2)) :- - trim_op_(Op1), - trim_op_(Op2). -trim_ops_($(Op1)) :- - trim_op_(Op1). - -trim_op_(Arg) :- number(Arg), !. -trim_op_(Arg) :- - get_attr(Arg, clpBNR, Def), % an interval ? - Def = interval(_, _, NList, _), - trim_persistent_(NList,TrimList), - % if trimmed list empty, set to a new unshared var to avoid cycles(?) on backtracking - (var(TrimList) -> setarg(3,Def,_) ; setarg(3,Def,TrimList)). % write trimmed node list - -trim_persistent_(T,T) :- var(T), !. % end of indefinite list -trim_persistent_([node(_,P,_,_)|Ns],TNs) :- nonvar(P), !, trim_persistent_(Ns,TNs). -trim_persistent_([N|Ns],[N|TNs]) :- trim_persistent_(Ns,TNs). - -% -% Any changes in interval values should come through here. -% Note: This captures all updated state for undoing on backtracking -% -updateValue_(Old, Old, _, _, Agenda, Agenda) :- !. % no change in value (constant?) - -updateValue_(Old, New, Int, OpsLeft, Agenda, NewAgenda) :- % set interval value to New - % New = [NewL,NewH], (NewL > NewH -> trace ; true), - % NewL=0 or narrowing sufficent - putValue_(New, Int, Nodelist), % update value (may fail) - linkNodeList_(Nodelist, Agenda, NewAgenda). % then propagate change - -updateValue_(_, _, _, _, Agenda, Agenda). % otherwise just continue with Agenda - -propagate_if_(Ops, _, _) :- Ops>0, !. % check work limiter -propagate_if_(_, (OL,OH), (NL,NH)) :- (NH-NL)/(OH-OL) < 0.9. % any overflow in calculation will propagate - -linkNodeList_([X|Xs], List, NewList) :- - nonvar(X), !, % not end of list ... - (arg(3,X,1) % test linked flag - -> NextList = List % already linked - ; linkNode_(List, X, NextList) % not linked add it to list - ), - linkNodeList_(Xs, NextList, NewList). -linkNodeList_(_, List, List). % end of indefinite node list (don't make it definite) - -% Assumes persistant nodes are pre-trimmed (see putValue_/3) -linkNode_(List/[X|NextTail], X, List/NextTail) :- % add to list - setarg(3,X,1). % set linked bit - -:- include(ia_utilities). % print,solve, etc. - -% -% Get all defined statistics -% -clpStatistics(Ss) :- findall(S, clpStatistic(S), Ss). - -% end of reset chain succeeds. Need cut since predicate is "discontiguous". -clpStatistics :- !. - -clpBNR_version_("0.9.4"). - -:- initialization(( - % restore "optimise" flag - (nb_current('clpBNR:temp',Opt) - -> (nb_delete('clpBNR:temp'), set_prolog_flag(optimise,Opt)) - ; true - ), - - clpBNR_version_(Version), format("*** clpBNR v~walpha ***\n\n",[Version]), - (current_prolog_flag(bounded,true) - -> write("Error: clpBNR requires unbounded integers and rationals.\n") - ; true - ), - (current_prolog_flag(float_overflow,_) - -> true - ; write("Error: clpBNR requires support for IEEE arithmetic.\n") - ), - - % Set required arithmetic flags - set_prolog_flag(prefer_rationals, true), % enable rational arithmetic - set_prolog_flag(max_rational_size, 16), % rational size in bytes before .. - set_prolog_flag(max_rational_size_action, float), % conversion to float - - set_prolog_flag(float_overflow,infinity), % enable IEEE continuation values - set_prolog_flag(float_zero_div,infinity), - set_prolog_flag(float_undefined,nan), - - once(prolog_stack_property(global,min_free(Free))), % minimum free cells (8 bytes/cell) - (Free < 8196 -> set_prolog_stack(global,min_free(8196)) ; true), - - % create clpBNR specific flags - create_prolog_flag(clpBNR_iteration_limit,3000,[type(integer),keep(true)]), - create_prolog_flag(clpBNR_default_precision,6,[type(integer),keep(true)]), - create_prolog_flag(clpBNR_verbose,false,[type(boolean),keep(true)]), - - clpStatistics % initialize statistics -)). diff --git a/docs/Package/clpBNR-0.9.4/pack.pl b/pack.pl similarity index 100% rename from docs/Package/clpBNR-0.9.4/pack.pl rename to pack.pl diff --git a/src/clpBNR.pl b/prolog/clpBNR.pl old mode 100755 new mode 100644 similarity index 91% rename from src/clpBNR.pl rename to prolog/clpBNR.pl index f804a44..6a6e191 --- a/src/clpBNR.pl +++ b/prolog/clpBNR.pl @@ -91,6 +91,22 @@ :- style_check([-singleton, -discontiguous]). % :- discontiguous ... not reliable. +:- create_prolog_flag(clpBNR_iteration_limit,3000,[type(integer),keep(true)]). +:- create_prolog_flag(clpBNR_default_precision,6,[type(integer),keep(true)]). +:- create_prolog_flag(clpBNR_verbose,false,[type(boolean),keep(true)]). + +% Provide exception hooks for the global variables used by this library. +% This allows running the library in any thread. + +user:exception(undefined_global_variable, Var, retry) :- + init_gvar(Var, Value), + nb_setval(Var, Value). + +init_gvar('clpBNR:node_count', 0). +init_gvar('clpBNR:evalNode', 0). +init_gvar('clpBNR:evalNodeFail', 0). +init_gvar('clpBNR:iterations', 0). + % % SWIP optimise control - set flag to true for compiled arithmetic % @@ -121,8 +137,8 @@ g_incb(G) :- nb_getval(G,N), N1 is N+1, b_setval(G,N1). % undone on backtrack g_read(G,V) :- nb_getval(G,V). -:- discontiguous - clpBNR:clpStatistics/0, clpBNR:clpStatistic/1, +:- discontiguous + clpBNR:clpStatistics/0, clpBNR:clpStatistic/1, sandbox:safe_global_variable/1. sandbox:safe_global_variable('clpBNR:userTime'). @@ -149,7 +165,7 @@ clpStatistic(inferences(I)) :- statistics(inferences,I1), g_read('clpBNR:inferences',I0), I is I1-I0. -:- include(ia_primitives). % interval arithmetic relations via evalNode/4. +:- include(clpBNR/ia_primitives). % interval arithmetic relations via evalNode/4. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % @@ -172,7 +188,7 @@ % flags (list) abstraction get_interval_flags_(Int,Flags) :- get_attr(Int, clpBNR, interval(Type, Val, Nodelist, Flags)). - + set_interval_flags_(Int,Flags) :- % flags assumed to be ground so no copy required get_attr(Int, clpBNR, interval(Type, Val, Nodelist, _)), put_attr(Int, clpBNR, interval(Type, Val, Nodelist, Flags)). @@ -182,7 +198,7 @@ % universal_interval((-1.0Inf,1.0Inf)). -% Finite intervals - 64 bit IEEE reals, +% Finite intervals - 64 bit IEEE reals, finite_interval(real, (-1.0e+16,1.0e+16)). finite_interval(integer, (L,H)) :- %% SWIP: current_prolog_flag(bounded,false),!, % integers are unbounded, but use tagged limits for finite default @@ -196,7 +212,7 @@ % % non-backtracking set bounds for use with branch and bound % -nb_setbounds(Int, [L,U]) :- +nb_setbounds(Int, [L,U]) :- get_attr(Int, clpBNR, Def), Def = interval(_, Val, _, _), ^(Val,(L,U),NewVal), % new range is intersection (from ia_primitives) @@ -205,7 +221,7 @@ % % get current value % -getValue(Int, Val) :- +getValue(Int, Val) :- number(Int) -> Val=(Int,Int) % numeric point value ; get_attr(Int, clpBNR, interval(_, Val, _, _)). % interval, optimized for SWIP @@ -214,10 +230,10 @@ % set monitor action on an interval % watch(Int,Action) :- - atom(Action), + atom(Action), get_interval_flags_(Int,Flags), !, lists:subtract(Flags,[watch(_)],Flags1), % remove any previous setting - (Action = none -> true ; set_interval_flags_(Int,[watch(Action)|Flags1])). + (Action = none -> true ; set_interval_flags_(Int,[watch(Action)|Flags1])). watch(Ints,Action) :- is_list(Ints), watch_list_(Ints,Action). @@ -254,7 +270,7 @@ queue_nodes_(real,_,_,Nodes,Nodes). % type real - just use Nodes queue_nodes_(integer,(L,H),_,Nodes,Nodes) :- % type integer #1 with integral bounds - integral_(L), integral_(H), !. % run Nodes if bounds integers or infinities + integral_(L), integral_(H), !. % run Nodes if bounds integers or infinities queue_nodes_(integer,_,Int,_,[node(integral,_,0,$(Int))|_]). % type integer #2 % schedule an "integral" node to re-adjust bounds (Nodes run on adjustment) @@ -263,7 +279,7 @@ integral_(B) :- integer(B). % -% range(Int, Bounds) for compatability +% range(Int, Bounds) for compatability % % for interval(Int) and number(Int), check if value is (always) in specified range, unifying any vars with current value range(Int, [L,H]) :- getValue(Int, (IL,IH)), number(IL), !, % existing interval @@ -275,7 +291,7 @@ getValue(Int, (L,H)). % will bind any intput var's to values % -% domain(Int, Dom) for interval(Int) for compatability +% domain(Int, Dom) for interval(Int) for compatability % domain(Int, Dom) :- interval_object(Int, Type, Val, _), @@ -296,7 +312,7 @@ % % midpoint(Int, Wid) midpoint of an interval or numeric value % based on: -% Frédéric Goualard. How do you compute the midpoint of an interval?. +% Frédéric Goualard. How do you compute the midpoint of an interval?. % ACM Transactions on Mathematical Software, Association for Computing Machinery, 2014, % 40 (2), 10.1145/2493882. hal-00576641v1 % Exception, single infinite bound treated as largest finite FP value @@ -323,7 +339,7 @@ median_bound_(lo,L,FL), median_bound_(hi,H,FH), median_(FL,FH,Med), !. - + median_bound_(lo,B,FB) :- B=:=0, FB is nexttoward(B,1.0). median_bound_(lo,-1.0Inf,FB) :- FB is nexttoward(-1.0Inf,0.0). median_bound_(lo,B,FB) :- FB is roundtoward(float(B), to_negative). @@ -353,22 +369,22 @@ Ints::Dom :- is_list(Ints),!, intervals_(Ints,Dom). - -R::Dom :- var(R), var(Dom), !, % declare R = real(L,H), Note: R can be interval + +R::Dom :- var(R), var(Dom), !, % declare R = real(L,H), Note: R can be interval int_decl_(real,(_,_),R), domain(R,Dom). R::Dom :- var(Dom), !, % domain query (if interval(R) else fail) domain(R,Dom). % "domain" query, unify interval Type and Bounds -R::Dom :- % interval(R) or number(R) and nonvar(Dom) +R::Dom :- % interval(R) or number(R) and nonvar(Dom) Dom=..[Type|Bounds], (Bounds=[] -> Val=(_,_) ; Val=..[,|Bounds]), int_decl_(Type,Val,R). int_decl_(boolean,_,R) :- !, % boolean is integer; 0=false, 1=true, ignore any bounds. int_decl_(integer,(0,1),R). - + int_decl_(Type,(L,H),R) :- interval_object(R,IType,CVal,_NL), !, % already interval lower_bound_val_(Type,L,IL), % changing type,bounds? upper_bound_val_(Type,H,IH), @@ -381,7 +397,7 @@ lower_bound_val_(Type,L,IL), upper_bound_val_(Type,H,IH), IL= + (IL=IH -> R=IL ; % point range, can unify now put_attr(R, clpBNR, interval(Type, (IL,IH), _NL, [])) % attach clpBNR attribute ). @@ -422,7 +438,7 @@ upper_bound_val_(integer,H,IH), (IL=IH -> Int=IL % narrowed to point - ; (put_attr(Int,clpBNR,interval(integer,(IL,IH),NodeList,Flags)), % set Type (only case allowed) + ; (put_attr(Int,clpBNR,interval(integer,(IL,IH),NodeList,Flags)), % set Type (only case allowed) linkNodeList_(NodeList, Agenda, NewAgenda) ) ). @@ -470,7 +486,7 @@ mergeFlags_(Flags1,[],Flags1) :- !. mergeFlags_(Flags1,Flags2,Flags) :- lists:union(Flags1,Flags2,Flags). % ambiguous if both have same flag set to different values - + mergeNodes_([N],NodeList,NodeList) :- var(N),!. mergeNodes_([N|Ns],NodeList,[N|NewList]) :- N=node(Op,_,_,Ops), @@ -479,7 +495,7 @@ mergeNodes_([N|Ns],NodeList,NewList) :- mergeNodes_(Ns,NodeList,NewList). -notIn_([node(NOp,_,_,NOps)|Ns],Op,Ops) :- % equivalent node(Op, _, _,Ops) ==> failure +notIn_([node(NOp,_,_,NOps)|Ns],Op,Ops) :- % equivalent node(Op, _, _,Ops) ==> failure NOp==Op, NOps==Ops, % avoid unification of ops and operands (Op may be more than just an atom) !, fail. notIn_([N|Ns],Op,Ops) :- % keep searching @@ -509,16 +525,16 @@ addConstraints_(C,Agenda,NextAgenda), !, addConstraints_(Cs,NextAgenda,NewAgenda). -% low overhead version for internal use -constrain_(C) :- +% low overhead version for internal use +constrain_(C) :- buildConstraint_(C,T/T,Agenda), stable_(Agenda). - + buildConstraint_(C,Agenda,NewAgenda) :- debug_clpBNR_('Add ~p',{C}), build_(C, 1, boolean, Agenda, NewAgenda), !. -:- include(ia_simplify). % simplifies constraints to a hopefully more efficient equivalent +:- include(clpBNR/ia_simplify). % simplifies constraints to a hopefully more efficient equivalent % % build a node from an expression @@ -539,9 +555,9 @@ number(Num), !. build_(Exp, Num, VarType, Agenda, Agenda) :- % pre-compile ground Exp (including pi and e) ground(Exp), - safe_(Exp), % safe to evaluate + safe_(Exp), % safe to evaluate L is roundtoward(Exp,to_negative), % rounding only affects float evaluation - (float(L) % not precise -> interval + (float(L) % not precise -> interval -> (H is roundtoward(Exp,to_positive), int_decl_(VarType,(L,H),Num)) % if precise use L ; Num=L % else use precise value ), !. @@ -561,18 +577,18 @@ safe_([A|As]) :- !, safe_(A), safe_(As). -safe_(F) :- +safe_(F) :- current_arithmetic_function(F), % evaluable by is/2 F =.. [Op|Args], \+memberchk(Op,[**,sin,cos,tan,asin,acos,atan]), % unsafe operations due to rounding safe_(Args). -% a constraint must evaluate to a boolean +% a constraint must evaluate to a boolean constraint_(C) :- nonvar(C), C =..[Op|_], fmap_(Op,_,_,_,[boolean|_]), !. % map constraint operator to primitive/arity/types fmap_(+, add, ZXY, ZXY, [real,real,real]). -fmap_(-, add, [Z,X,Y], [X,Z,Y], [real,real,real]). % note subtract before minus +fmap_(-, add, [Z,X,Y], [X,Z,Y], [real,real,real]). % note subtract before minus fmap_(*, mul, ZXY, ZXY, [real,real,real]). fmap_(/, mul, [Z,X,Y], [X,Z,Y], [real,real,real]). fmap_(**, pow, ZXY, ZXY, [real,real,real]). @@ -612,7 +628,7 @@ % reverse map from Op and Args (used by "verbose" top level output to reverse compile constraints) remap_(Op,$(Z,X,Y),C) :- constraint_(Op), Z==1, !, % simplification for constraints - C=..[Op,X,Y]. + C=..[Op,X,Y]. remap_(Op,$(Z,X,Y),Z==C) :- !, C=..[Op,X,Y]. remap_(Op,$(Z,X),Z==C) :- @@ -642,7 +658,7 @@ g_read('clpBNR:node_count',C). % extend list with X -newmember([X|Xs],N) :- +newmember([X|Xs],N) :- nonvar(X), !, % not end of (indefinite) list newmember(Xs,N). newmember([N|_],N). % end of list @@ -669,7 +685,7 @@ sandbox:safe_global_variable('clpBNR:iterations'). clpStatistics :- - current_prolog_flag(clpBNR_iteration_limit,L), + current_prolog_flag(clpBNR_iteration_limit,L), g_assign('clpBNR:iterations',L), % set "low" water mark to upper limit fail. % backtrack to reset other statistics. @@ -680,7 +696,7 @@ % % Execute a node on the queue -% Note: "special" ops like instantiate and integral not counted as narrowing op in clpStatistics +% Note: "special" ops like instantiate and integral not counted as narrowing op in clpStatistics % % Comment out the following to enable Op tracing: goal_expansion(traceIntOp_(Op, Args, PrimArgs, New),true). @@ -693,9 +709,9 @@ getValue(YArg,YVal), evalNode(Op, P, $(ZVal,XVal,YVal), $(NZVal,NXVal,NYVal)), % can fail causing stable_ to fail => backtracking traceIntOp_(Op, [ZArg,XArg,YArg], [ZVal,XVal,YVal], [NZVal,NXVal,NYVal]), % in ia_utilities - updateValue_(ZVal, NZVal, ZArg, OpsLeft, Agenda, AgendaZ), - updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, AgendaZX), - updateValue_(YVal, NYVal, YArg, OpsLeft, AgendaZX, NewAgenda). + updateValue_(ZVal, NZVal, ZArg, OpsLeft, Agenda, AgendaZ), + updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, AgendaZX), + updateValue_(YVal, NYVal, YArg, OpsLeft, AgendaZX, NewAgenda). doNode_($(ZArg,XArg), Op, P, OpsLeft, Agenda, NewAgenda) :- % Arity 2 Op var(P), !, % check persistent bit @@ -713,7 +729,7 @@ getValue(Arg,Val), evalNode(Op, P, $(Val), $(NVal)), % can fail causing stable_ to fail => backtracking traceIntOp_(Op, [Arg], [Val], [NVal]), % in ia_utilities - updateValue_(Val, NVal, Arg, 1, Agenda,NewAgenda). % always update value regardless of OpsLeft limiter + updateValue_(Val, NVal, Arg, 1, Agenda,NewAgenda). % always update value regardless of OpsLeft limiter doNode_(Args, Op, p, _, Agenda, Agenda) :- % persistent bit "set", skip node and trim trim_ops_(Args). @@ -733,14 +749,14 @@ trim_op_(Op1). trim_op_(Arg) :- number(Arg), !. -trim_op_(Arg) :- +trim_op_(Arg) :- get_attr(Arg, clpBNR, Def), % an interval ? Def = interval(_, _, NList, _), trim_persistent_(NList,TrimList), % if trimmed list empty, set to a new unshared var to avoid cycles(?) on backtracking (var(TrimList) -> setarg(3,Def,_) ; setarg(3,Def,TrimList)). % write trimmed node list -trim_persistent_(T,T) :- var(T), !. % end of indefinite list +trim_persistent_(T,T) :- var(T), !. % end of indefinite list trim_persistent_([node(_,P,_,_)|Ns],TNs) :- nonvar(P), !, trim_persistent_(Ns,TNs). trim_persistent_([N|Ns],[N|TNs]) :- trim_persistent_(Ns,TNs). @@ -775,7 +791,7 @@ linkNode_(List/[X|NextTail], X, List/NextTail) :- % add to list setarg(3,X,1). % set linked bit -:- include(ia_utilities). % print,solve, etc. +:- include(clpBNR/ia_utilities). % print,solve, etc. % % Get all defined statistics @@ -787,39 +803,52 @@ clpBNR_version_("0.9.4"). -:- initialization(( - % restore "optimise" flag - (nb_current('clpBNR:temp',Opt) - -> (nb_delete('clpBNR:temp'), set_prolog_flag(optimise,Opt)) - ; true +check_features :- + ( current_prolog_flag(bounded,true) + -> print_message(error, clpBNR(bounded)) + ; true ), + ( current_prolog_flag(float_overflow,_) + -> true + ; print_message(error, clpBNR(no_IEEE)) + ). - clpBNR_version_(Version), format("*** clpBNR v~walpha ***\n\n",[Version]), - (current_prolog_flag(bounded,true) - -> write("Error: clpBNR requires unbounded integers and rationals.\n") - ; true - ), - (current_prolog_flag(float_overflow,_) - -> true - ; write("Error: clpBNR requires support for IEEE arithmetic.\n") - ), - - % Set required arithmetic flags +clp_set_prolog_flags :- set_prolog_flag(prefer_rationals, true), % enable rational arithmetic set_prolog_flag(max_rational_size, 16), % rational size in bytes before .. set_prolog_flag(max_rational_size_action, float), % conversion to float set_prolog_flag(float_overflow,infinity), % enable IEEE continuation values set_prolog_flag(float_zero_div,infinity), - set_prolog_flag(float_undefined,nan), - - once(prolog_stack_property(global,min_free(Free))), % minimum free cells (8 bytes/cell) - (Free < 8196 -> set_prolog_stack(global,min_free(8196)) ; true), - - % create clpBNR specific flags - create_prolog_flag(clpBNR_iteration_limit,3000,[type(integer),keep(true)]), - create_prolog_flag(clpBNR_default_precision,6,[type(integer),keep(true)]), - create_prolog_flag(clpBNR_verbose,false,[type(boolean),keep(true)]), - - clpStatistics % initialize statistics -)). + set_prolog_flag(float_undefined,nan). + +set_min_free(Amount) :- + ( prolog_stack_property(global,min_free(Free)) + -> ( Free < Amount + -> set_prolog_stack(global,min_free(Amount)) + ; true + ) + ; true + ). + +restore_optimise :- + ( nb_current('clpBNR:temp',Opt) + -> nb_delete('clpBNR:temp'), + set_prolog_flag(optimise,Opt) + ; true + ). + +finish_up :- + check_features, + clp_set_prolog_flags, + set_min_free(8196), + restore_optimise. + +:- initialization(finish_up, now). + +:- multifile prolog:message//1. + +prolog:message(clpBNR(bounded)) --> + [ 'clpBNR requires unbounded integers and rationals.'-[] ]. +prolog:message(clpBNR(no_IEEE)) --> + [ 'clpBNR requires support for IEEE arithmetic.'-[] ]. diff --git a/docs/Package/clpBNR-0.9.4/prolog/ia_primitives.pl b/prolog/clpBNR/ia_primitives.pl old mode 100755 new mode 100644 similarity index 100% rename from docs/Package/clpBNR-0.9.4/prolog/ia_primitives.pl rename to prolog/clpBNR/ia_primitives.pl diff --git a/docs/Package/clpBNR-0.9.4/prolog/ia_simplify.pl b/prolog/clpBNR/ia_simplify.pl similarity index 100% rename from docs/Package/clpBNR-0.9.4/prolog/ia_simplify.pl rename to prolog/clpBNR/ia_simplify.pl diff --git a/docs/Package/clpBNR-0.9.4/prolog/ia_utilities.pl b/prolog/clpBNR/ia_utilities.pl old mode 100755 new mode 100644 similarity index 100% rename from docs/Package/clpBNR-0.9.4/prolog/ia_utilities.pl rename to prolog/clpBNR/ia_utilities.pl diff --git a/src/ia_primitives.pl b/src/ia_primitives.pl deleted file mode 100755 index a265485..0000000 --- a/src/ia_primitives.pl +++ /dev/null @@ -1,838 +0,0 @@ -/* The MIT License (MIT) - * - * Copyright (c) 2019,2020 Rick Workman - * - * Permission is hereby granted, free of charge, to any person obtaining a copy - * of this software and associated documentation files (the "Software"), to deal - * in the Software without restriction, including without limitation the rights - * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the Software is - * furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included in all - * copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - * SOFTWARE. - */ - -% -% API function for executing primitives -% -% evalNode(+primitive_name, ?persistence_flag, +$(inputs..), -$(outputs..)) -% -evalNode(Op, P, Is, R) :- - g_inc('clpBNR:evalNode'), % count of primitive calls - narrowing_op(Op, P, Is, R), - !. -evalNode(Op, _, _, _):- - g_inc('clpBNR:evalNodeFail'), % count of primitive call failures - debugging(clpBNR, true), % fail immediately unless debug=true - current_node_(Node), - debug_clpBNR_('Fail ~w',Node), - fail. - -sandbox:safe_global_variable('clpBNR:evalNode'). -sandbox:safe_global_variable('clpBNR:evalNodeFail'). - -clpStatistics :- - g_assign('clpBNR:evalNode',0), - g_assign('clpBNR:evalNodeFail',0), - fail. % backtrack to reset other statistics. - -clpStatistic(narrowingOps(C)) :- g_read('clpBNR:evalNode',C). - -clpStatistic(narrowingFails(C)) :- g_read('clpBNR:evalNodeFail',C). - -% SWIP optimization for non_empty/2, replaces nexttoward function call with constant -goal_expansion(L < RHS, L < MaxFloat) :- - RHS == nexttoward( 1.0Inf,0), - current_prolog_flag(float_max,MaxFloat). % MaxFloat is nexttoward( 1.0Inf,0). -goal_expansion(H > RHS, H > NegMaxFloat) :- - RHS == nexttoward( -1.0Inf,0), - current_prolog_flag(float_max,MaxFloat), - NegMaxFloat is -MaxFloat. % NegMaxFloat is nexttoward( -1.0Inf,0). -% -% non-empty inteval test, L= L < nexttoward( 1.0Inf,0) ; true), - (float(H) -> H > nexttoward(-1.0Inf,0) ; true). - -% -% forces all intervals to boolean range, (optimize if already boolean) -% -booleanVal_((0,0),(0,0)). -booleanVal_((1,1),(1,1)). -booleanVal_((0,1),(0,1)). -booleanVal_(V,(0,1)):- ^(V,(0,1),(0,1)). % constrain non-booleans to (0,1) - -% -% interval category: non-negative (includes zero), non-positive, or split (contains 0) -% -intCase(p, (L,H)) :- L>=0, !. % non-negative -intCase(n, (L,H)) :- H=<0, !. % non-positive -intCase(s, (L,H)). % split - -% -% interval primitive functions -% X, Y and Z are intervals -% - -% Z := X ^ Y (intersection) -^((Xl,Xh), (Yl,Yh), (Zl,Zh)) :- - Zl is max(Xl, Yl), Zh is min(Xh, Yh), - non_empty(Zl,Zh). - -% -% wrap repeating interval onto a prime cylinder of width W and -% return projected interval and "mulipliers" to re-project -% Note: fails if interval wider than W so beware of outward rounding on ranges, e.g., -% range of [-pi/2,pi/2] on tan argument actually spans 3 cylinders (-1 to 1). -% -wrap_((Xl,Xh), W, (MXl,MXh), (Xpl,Xph)) :- % project onto cylinder from -W/2 to W/2 - MXl is round(Xl/W), - MXh is round(Xh/W), - MXh-MXl =< 1, Xh-Xl =< W, % MX check first to avoid overflow - Xpl is Xl - (MXl*W), Xph is Xh-(MXh*W). - -% -% unwrap projected interval back to original range -% -unwrap_((Xpl,Xph), W, (MXl,MXh), (Xl,Xh)) :- - Xl is Xpl+W*MXl, Xh is Xph+W*MXh. - -union_(X,[],X) :-!. % set union (including empty set) -union_([],Y,Y) :-!. -union_((Xl,Xh),(Yl,Yh),(Zl,Zh)) :- Zl is min(Xl,Yl), Zh is max(Xh,Yh). - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% Relational Operations (narrowing operations) -% -:- discontiguous clpBNR:narrowing_op/4. -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% integral (contract to integer bounds) - -narrowing_op(integral, P, $((L,H)), $((NL,NH))) :- - NL is ceiling(L), NH is floor(H), - NL= P=p ; true). % if a point, mark as persistent (makes a huge difference on some problems) - -narrowing_op(eq, p, $(Z, X, Y), $(NewZ, X, Y)) :- % persistent, X and Y don't intersect, Z is false - \+(^(X,Y,_)), !, - ^(Z, (0,0), NewZ). - -narrowing_op(eq, p, $(Z, (X,X), (Y,Y)), $(NewZ, (X,X), (Y,Y))) :- % if X and Y are necessarily equal, Z is true - X =:= Y, !, - ^(Z, (1,1), NewZ). - -narrowing_op(eq, _, $(Z,X,Y), $(NewZ,X,Y)) :- ^(Z,(0,1),NewZ). % else no change, but narrow Z to boolean - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==(X<>Y) % (Z boolean, X,Y integer) - -narrowing_op(ne, p, $(Z, X, Y), $(NewZ, X, Y)) :- % persistent: disjoint or necessarily equal - not_equal_(X,Y,Z1), !, - ^(Z, Z1, NewZ). - -narrowing_op(ne, _, $((1,1), X, Y), $((1,1), NewX, NewY)) :- % Z true, X/Y a point = bound of Y/X - ne_int_(X,Y,NewX), % narrow X if Y is a point and a bound of X (always succeeds but may not narrow) - ne_int_(Y,NewX,NewY), !. % narrow Y if NewX is a point and a bound of Y (always succeeds but may not narrow) - -narrowing_op(ne, _, $(Z,X,Y), $(NewZ,X,Y)) :- ^(Z,(0,1),NewZ). % none of the above, narrow Z to boolean - -not_equal_((Xl,Xh), (Yl,Yh), (1,1)) :- (Xh < Yl ; Yh < Xl). % X and Y disjoint, Z true -not_equal_((X,X), (Y,Y), (0,0)) :- X =:= Y. % pt(X)=:=pt(Y), Z false - -% Z <> X, where where Z and X are integer intervals (enforced elsewhere) -ne_int_((X,H), (X,X), (NewL,H)) :- !, % X is a point, and low bound of Z - NewL is X+1, NewL= le_true(X,Y,New,P) % Z true - ; (Z=(0,0) - -> le_false(X,Y,New,P) % Z false - ; le_bool(X,Y,New,P) % Z unknown - ) - ). -le_true(X, (Yl,Yh), $((1,1), (NXl,NXh), (NYl,NYh)), P) :- - ^(X, (-1.0Inf,Yh), (NXl,NXh)), % NewX := (Xl,Xh) ^(NI,Yh) - ^((Yl,Yh), (NXl,1.0Inf), (NYl,NYh)), % NewY := (Yl,Yh) ^(Xl,PI) - (NXh =< NYl -> P=p ; true). % will? always be true - -le_false(X, Y, $((0,0), NewX, NewY), P) :- - lt_true(Y,X,$(_,NewY,NewX),P). - -le_bool((Xl,Xh), (Yl,Yh), $(NewZ, (Xl,Xh), (Yl,Yh)), P) :- - (Yh < Xl - -> (NewZ=(0,0), P=p) % false persistant - ; (Xh =< Yl - -> (NewZ=(1,1), P=p) % true persistant - ; NewZ=(0,1) % indefinite boolean - ) - ). - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==(X lt_true(X,Y,New,P) % Z true - ; (Z=(0,0) - -> lt_false(X,Y,New,P) % Z false - ; lt_bool(X,Y,New,P) % Z unknown - ) - ). - -lt_true( X, (Yl,Yh), $((1,1), (NXl,NXh), (NYl,NYh)), P) :- - next_lt_(Yh,-1,YhD), % YhD is next downward value from Yh - ^(X, (-1.0Inf,YhD), (NXl,NXh)), % NewX := (Xl,Xh) ^(NInf,YhD) - next_lt_(NXl,1,NXlU), % NXlU is next upward value from NXl - ^((Yl,Yh), (NXlU,1.0Inf), (NYl,NYh)), % NewY := (Yl,Yh) ^(NXlU,PInf) - (NXh < NYl -> P=p ; true). % will? always be true - -lt_false(X, Y, $((0,0), NewX, NewY), P) :- - le_true(Y,X,$(_,NewY,NewX),P). - -lt_bool((Xl,Xh), (Yl,Yh), $(NewZ, (Xl,Xh), (Yl,Yh)), P) :- - (Yh =< Xl - -> (NewZ=(0,0), P=p) % false persistant - ; (Xh < Yl - -> (NewZ=(1,1), P=p) % true persistant - ; NewZ=(0,1) % indefinite boolean - ) - ). - -% Note: can't narrow an infinite bound, minimize change to bound -% Repeat, not sound on reals (uses nexttoward, missing values between floats) -next_lt_( 1.0Inf, _, 1.0Inf) :- !. -next_lt_(-1.0Inf, _, -1.0Inf) :- !. -next_lt_(V, -1, NV) :- NV is max(V-1,nexttoward(V,-1.0Inf)). % integers will get sorted out with `integral` -next_lt_(V, 1, NV) :- NV is min(V+1,nexttoward(V, 1.0Inf)). - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==(X<=Y) % inclusion, constrains X to be subinterval of Y (Z boolean) - -% Only two cases: either X and Y intersect or they don't. -narrowing_op(in, _, $(Z, X, Y), $(NewZ, NewX, Y)):- - ^(X,Y,NewX), !, % NewX is intersection of X and Y - ^(Z,(1,1),NewZ). - -narrowing_op(in, p, $((0,0), X, Y), $((0,0), NewX, Y)):- % persistent, X and Y don't intersect' - ^(Z,(0,0),NewZ). - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==X+Y - -narrowing_op(add, _, $((Zl,Zh), (Xl,Xh), (Yl,Yh)), $((NZl,NZh), (NXl,NXh), (NYl,NYh))) :- - NZl is max(Zl, roundtoward( Xl+Yl, to_negative)), % NewZ := Z ^ (X+Y), - NZh is min(Zh, roundtoward( Xh+Yh, to_positive)), - non_empty(NZl,NZh), - % Note: subtraction done by adding minus values so rounding mode consistent - % during any numeric type conversion. - NXl is max(Xl, roundtoward(NZl+(-Yh), to_negative)), % NewX := X ^ (NZ-Y), - NXh is min(Xh, roundtoward(NZh+(-Yl), to_positive)), - non_empty(NXl,NXh), - NYl is max(Yl, roundtoward(NZl+(-NXh), to_negative)), % NewY := Y ^ (NZ-NX). - NYh is min(Yh, roundtoward(NZh+(-NXl), to_positive)), - non_empty(NYl,NYh). - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==X*Y - -narrowing_op(mul, _, $(Z,X,Y), $(NewZ, NewX, NewY)) :- - mul_(X,Y,Z,NewZ), % NewZ := Z ^ (X*Y) - odiv_(NewZ,X,Y,Yp), % Yp := Y ^ (Z/X), - odiv_(NewZ,Yp,X,NewX), % NewX := X ^ (Z/Yp), - % if X narrowed it may be necessary to recalculate Y due to non-optimal ordering. - mul_redoY(Y,Yp,X,NewX,NewY,NewZ). - -mul_redoY(Y,Y,X,NewX,NewY,NewZ) :- % if initial Y narrowing didn't change (contain 0?) - X \= NewX, !, % and X did narrow - %/(NewZ,NewX,Y1), ^(Y,Y1,NewY). % recalculate Y with NewX - odiv_(NewZ,NewX,Y,NewY). % recalculate Y with NewX -mul_redoY(_,NewY,_,_,NewY,_). % else keep Y as NewY - -% NZ := Z ^ (X * Y) (multiply) -mul_(X, Y, Z, NZ) :- - intCase(Cx,X), - intCase(Cy,Y), - multCase(Cx,Cy,X,Y,Z,NZ), !. - -% NZ := Z ^ (X / Y) (odiv) -odiv_(X, Y, Z, NZ) :- - intCase(Cx,X), - intCase(Cy,Y), - odivCase(Cx,Cy,X,Y,Z,NZ). % !'s in odiv' - -% -% * cases ("Interval Arithmetic: from Principles to Implementation", Fig. 3) -% -%multCase(z,_, X, _, _, X) :- !. % X==0 -%multCase(_,z, _, Y, _, Y). % Y==0 - -multCase(p,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- - NZl is max(Zl,roundtoward(Xl*Yl,to_negative)), - NZh is min(Zh,roundtoward(Xh*Yh,to_positive)), - non_empty(NZl,NZh). -multCase(p,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- - NZl is max(Zl,roundtoward(Xh*Yl,to_negative)), - NZh is min(Zh,roundtoward(Xl*Yh,to_positive)), - non_empty(NZl,NZh). -multCase(n,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- - NZl is max(Zl,roundtoward(Xl*Yh,to_negative)), - NZh is min(Zh,roundtoward(Xh*Yl,to_positive)), - non_empty(NZl,NZh). -multCase(n,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- - NZl is max(Zl,roundtoward(Xh*Yh,to_negative)), - NZh is min(Zh,roundtoward(Xl*Yl,to_positive)), - non_empty(NZl,NZh). - -multCase(p,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- - NZl is max(Zl,roundtoward(Xh*Yl,to_negative)), - NZh is min(Zh,roundtoward(Xh*Yh,to_positive)), - non_empty(NZl,NZh). -multCase(n,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- - NZl is max(Zl,roundtoward(Xl*Yh,to_negative)), - NZh is min(Zh,roundtoward(Xl*Yl,to_positive)), - non_empty(NZl,NZh). -multCase(s,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- - NZl is max(Zl,roundtoward(Xl*Yh,to_negative)), - NZh is min(Zh,roundtoward(Xh*Yh,to_positive)), - non_empty(NZl,NZh). -multCase(s,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- - NZl is max(Zl,roundtoward(Xh*Yl,to_negative)), - NZh is min(Zh,roundtoward(Xl*Yl,to_positive)), - non_empty(NZl,NZh). - -multCase(s,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- - NZl is max(Zl,min(roundtoward(Xl*Yh,to_negative),roundtoward(Xh*Yl,to_negative))), - NZh is min(Zh,max(roundtoward(Xl*Yl,to_positive),roundtoward(Xh*Yh,to_positive))), - non_empty(NZl,NZh). - -% -% / cases ("Interval Arithmetic: from Principles to Implementation", Fig. 4) -% Tricky handling the "zero" cases - returns universal interval when no narowing possible: -% 1. denominator is zero or contains zero (z and s) -% 2. denominator is bounded by zero (p and n) and numerator contains zero (see numeric guards) -% Numeric guards check for 0/0, e.g., p,p -> Xh+Yh>0 - -odivCase(p,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :- - sign(Yl)+sign(Xl) > 0, !, % X > 0 or Y > 0 - NZl is max(Zl,roundtoward(Xl/Yh,to_negative)), - NZh is min(Zh,roundtoward(Xh/max(0.0,Yl),to_positive)), - non_empty(NZl,NZh). -odivCase(p,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :- - sign(Xl)-sign(Yh) > 0, !, % X > 0 or Y < 0 - NZl is max(Zl,roundtoward(Xh/min(-0.0,Yh),to_negative)), - NZh is min(Zh,roundtoward(Xl/Yl,to_positive)), - non_empty(NZl,NZh). -odivCase(p,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :- - Xl > 0, !, % X > 0 - ZIl is roundtoward(Xl/Yh,to_negative), - ZIh is roundtoward(Xl/Yl,to_positive), - (Zl > ZIh -> NZl is max(Zl,ZIl) ; NZl = Zl), % similar to ^/3 - (Zh < ZIl -> NZh is min(Zh,ZIh) ; NZh = Zh), - non_empty(NZl,NZh). - -odivCase(n,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):- - sign(Yl)-sign(Xh) > 0, !, % X < 0 or Y > 0 - NZl is max(Zl,roundtoward(Xl/max(0.0,Yl),to_negative)), - NZh is min(Zh,roundtoward(Xh/Yh,to_positive)), - non_empty(NZl,NZh). -odivCase(n,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :- - sign(Yh)+sign(Xh) < 0, !, % X < 0 or Y < 0 - NZl is max(Zl,roundtoward(Xh/Yl,to_negative)), - NZh is min(Zh,roundtoward(Xl/min(-0.0,Yh),to_positive)), - non_empty(NZl,NZh). -odivCase(n,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :- - Xh < 0, !, % X < 0 - ZIl is roundtoward(Xh/Yl,to_negative), - ZIh is roundtoward(Xh/Yh,to_positive), - (Zl > ZIh -> NZl is max(Zl,ZIl) ; NZl = Zl), % similar to ^/3 - (Zh < ZIl -> NZh is min(Zh,ZIh) ; NZh = Zh), - non_empty(NZl,NZh). - -odivCase(s,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :- - Yl > 0, !, % Y > 0 - NZl is max(Zl,roundtoward(Xl/Yl,to_negative)), - NZh is min(Zh,roundtoward(Xh/Yl,to_positive)), - non_empty(NZl,NZh). -odivCase(s,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :- - Yh < 0, !, % Y < 0 - NZl is max(Zl,roundtoward(Xh/Yh,to_negative)), - NZh is min(Zh,roundtoward(Xl/Yh,to_positive)), - non_empty(NZl,NZh). - -odivCase(_,_, _, _, Z, Z). % all other cases, no narrowing - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==min(X,Y) Z==max(X,Y) - -narrowing_op(min, _, $((Zl,Zh),(Xl,Xh),(Yl,Yh)), New) :- - Z1l is max(Zl,min(Xl,Yl)), % Z1 := Z ^ min(X,Y), - Z1h is min(Zh,min(Xh,Yh)), - minimax((Zl,1.0Inf), $((Z1l,Z1h),(Xl,Xh),(Yl,Yh)), New). - -narrowing_op(max, _, $((Zl,Zh),(Xl,Xh),(Yl,Yh)), New) :- - Z1l is max(Zl,max(Xl,Yl)), % Z1 := Z ^ max(X,Y), - Z1h is min(Zh,max(Xh,Yh)), - minimax((-1.0Inf,Zh), $((Z1l,Z1h),(Xl,Xh),(Yl,Yh)), New). - -minimax(_, $(Z,X,Y), $(New, X, New)) :- % Case 1, X not in Z1 - \+(^(Z,X,_)), !, % _ := Z1 \^ X, - ^(Y,Z,New). % New := Y ^ Z1. - -minimax(_, $(Z,X,Y), $(New, New, Y)) :- % Case 2, Y not in Z1 - \+(^(Z,Y,_)), !, % _ := Z1 \^ Y, - ^(X,Z,New). % New := X ^ Z1. - -minimax(Zpartial, $(Z,X,Y), $(Z, NewX, NewY)) :- % Case 3, overlap - ^(X,Zpartial,NewX), % NewX := X ^ Zpartial, - ^(Y,Zpartial,NewY). % NewY := Y ^ Zpartial. - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==abs(X) - -narrowing_op(abs, _, $(Z,X), $(NewZ, NewX)) :- - intCase(Cx,X), - absCase(Cx,X,Z,NewZ), !, - non_empty(NewZ), - intersect_regions_(X,NewZ,NewX), - non_empty(NewX). - -% -% intersect and join Z with positive and negative regions of X -% -intersect_regions_(Z,X,NewZ) :- - (^(Z,X,ZPos) -> true ; ZPos = []), % positive Z region - (-(X,Z,ZNeg) -> true ; ZNeg = []), % negative Z region - union_(ZPos,ZNeg,NewZ). - -% NZ := Z ^ -X (unary minus) --((Xl,Xh), (Zl,Zh), (NZl,NZh)) :- - NZl is max(Zl,-Xh), NZh is min(Zh,-Xl), - NZl =< NZh. % no infinity/overflow check required - -% -% abs(X) cases -% -absCase(p, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is max(Zl,Xl), NZh is min(Zh,Xh). -absCase(n, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is max(Zl,-Xh), NZh is min(Zh,-Xl). -absCase(s, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is max(Zl,0), NZh is min(Zh,max(-Xl,Xh)). - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z== -X - -narrowing_op(minus, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :- - NZl is max(Zl,-Xh), NZh is min(Zh,-Xl), % NewZ := Z ^ -X, no empty check required - NXl is max(Xl,-NZh), NXh is min(Xh,-NZl). % NewX := X ^ -Z, no empty check required - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z== sqrt(X) (Z is positive root of positive X) - -narrowing_op(sqrt, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :- - Xh>=0, Xl1 is max(0,Xl), % narrow X to positive range - NZl is max(max(0,Zl),roundtoward(sqrt(Xl1),to_negative)), - NZh is min(Zh,roundtoward(sqrt(Xh),to_positive)), - non_empty(NZl,NZh), - NXh is min(Xh,roundtoward(NZh*NZh,to_positive)), - NXl is max(Xl,-NXh), - non_empty(NXl,NXh). - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z== exp(X) - -narrowing_op(exp, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :- - NZl is max(Zl,roundtoward(exp(Xl),to_negative)), - NZh is min(Zh,roundtoward(exp(Xh),to_positive)), - non_empty(NZl,NZh), - NXl is max(Xl,roundtoward(log(NZl),to_negative)), - NXh is min(Xh,roundtoward(log(NZh),to_positive)), - non_empty(NXl,NXh). - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z== X**Y - -narrowing_op(pow, _, $(Z,X,(R,R)), $(NewZ, NewX, (R,R))) :- !, % Y = point exponent - pt_powr_(Z,X,R,NewZ), non_empty(NewZ), - Ri is 1/R, - pt_powr_(X,NewZ,Ri,NewX), non_empty(NewX). - -narrowing_op(pow, _, $(Z,(Xl,Xh),Y), $(NewZ, NewX, NewY)) :- % interval Y, X must be positive - XZl is max(0,Xl), Xp = (XZl,Xh), Xh >= 0, % X must be positive (>=0) - powr_(Z,Xp,Y,NewZ), non_empty(NewZ), - universal_interval(UI), odiv_((1,1),Y,UI,Yi), - powr_(Xp,NewZ,Yi,NewX), non_empty(NewX), - ln(NewZ,Ynum), ln(NewX,Yden), - odiv_(Ynum,Yden,Y,NewY). % NY := Y ^ log(NZ)/log(NX)*/ - -% requires conservative rounding using nexttoward -pt_powr_(Z,X,R,NewZ) :- - R < 0, !, % negative R - Rp is -R, - universal_interval(UI), - odiv_((1,1),Z,UI,Zi), - pt_powr_(Zi,X,Rp,NZi), - odiv_((1,1),NZi,Z,NewZ). -pt_powr_(Z,X,R,NewZ) :- - integer(R), R rem 2 =:= 0, !, % even integer R - Z = (Zl,Zh), X = (Xl,Xh), - (float(Xl) - -> Zl1 is max(0,nexttoward(roundtoward(Xl**R,to_negative),-1.0Inf)) - ; Zl1 is max(0,roundtoward(Xl**R,to_negative)) - ), - (float(Xh) - -> Zh1 is nexttoward(roundtoward(Xh**R,to_positive), 1.0Inf) - ; Zh1 is roundtoward(Xh**R,to_positive) - ), - ( sign(Xl)*sign(Xh) =< 0 - -> NZl is max(0,Zl) % X bounded by or includes 0 - ; NZl is max(min(Zl1,Zh1),Zl) % X doesn't include 0 - ), - NZh is min(max(Zl1,Zh1),Zh), - NewZ = (NZl,NZh). -pt_powr_(Z,X,R,NewZ) :- - rational(R), denominator(R) rem 2 =:= 0, !, % even rational R - Z = (Zl,Zh), X = (Xl,Xh), - (float(Xl) - -> Zl1 is max(0,nexttoward(roundtoward(Xl**R,to_negative),-1.0Inf)) - ; Zl1 is max(0,roundtoward(Xl**R,to_negative)) - ), - (float(Xh) - -> Zh1 is nexttoward(roundtoward(Xh**R,to_positive), 1.0Inf) - ; Zh1 is roundtoward(Xh**R,to_positive) - ), - Zpl is max(Zl,Zl1), Zph is min(Zh,Zh1), % positive case - (Zpl =< Zph -> Zps = (Zpl,Zph) ; Zps = []), - Znl is max(Zl,-Zh1), Znh is min(Zh,-Zl1), % negative case - (Znl =< Znh -> Zns = (Znl,Znh) ; Zns = []), - union_(Zps,Zns,NewZ). % union of positive and negative cases -pt_powr_((Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :- % all other cases - (float(Xl) - -> Zl1 is nexttoward(roundtoward(Xl**R,to_negative),-1.0Inf) - ; Zl1 is roundtoward(Xl**R,to_negative) - ), - (float(Xh) - -> Zh1 is nexttoward(roundtoward(Xh**R,to_positive), 1.0Inf) - ; Zh1 is roundtoward(Xh**R,to_positive) - ), - NZl is max(Zl,min(Zl1,Zh1)), NZh is min(Zh,max(Zl1,Zh1)). - -ln((Xl,Xh), (Zl,Zh)) :- - Zl is roundtoward(log(Xl),to_negative), % rounding dependable? - Zh is roundtoward(log(Xh),to_positive). - -powr_((Zl,Zh), (Xl,Xh), (Yl,Yh), (NZl,NZh)) :- % X assumed positive - Zll is max(0,nexttoward(roundtoward(Xl**Yl,to_negative),-1.0Inf)), - Zlh is nexttoward(roundtoward(Xl**Yh,to_positive), 1.0Inf), - Zhl is max(0,nexttoward(roundtoward(Xh**Yl,to_negative),-1.0Inf)), - Zhh is nexttoward(roundtoward(Xh**Yh,to_positive), 1.0Inf), - NZlmin is min(Zll, min(Zlh, min(Zhl, Zhh))), % min of all - NZl is max(Zl,NZlmin), - NZhmax is max(Zll, max(Zlh, max(Zhl, Zhh))), % max of all - NZh is min(Zh,NZhmax). - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z== sin(X) - -narrowing_op(sin, _, $(Z,X), $(NewZ, NewX)) :- - wrap_(X,2*pi,MX,Xp), !, % fails if X too wide - sin_(MX,Xp,Z,NMX,X1,NewZ), - unwrap_(X1,2*pi,NMX,X2), % project back to original cylinder - non_empty(NewZ), - ^(X,X2,NewX). -% ^(Z1,(-1,1),NewZ). % take care of any rounding beyond range - -narrowing_op(sin, _, $(Z,X), $(NewZ,X)) :- % no narrowing possible, just constrain Z - ^(Z,(-1,1),NewZ). - -sin_((MX,MX), X, Z, (MX,MX), NewX, NewZ) :- - !, % same cylinder, split into 3 interval convex sectors - Pi_2 is pi/2, - sin_sector_(lo, Pi_2, X, Z, NXlo, NZlo), - sin_sector_(mid,Pi_2, X, Z, NXmid, NZmid), - sin_sector_(hi, Pi_2, X, Z, NXhi, NZhi), - union_(NXlo,NXmid,NX1), union_(NX1,NXhi,NewX), - union_(NZlo,NZmid,NZ1), union_(NZ1,NZhi,NewZ). - -sin_((MXl,MXh), (Xl,Xh), Z, (NMXl,NMXh), NewX, NewZ) :- - % adjacent cylinders, - Pi is pi, MPi is -pi, - try_sin_((Xl,Pi), Z, NX1, NZ1,MXl,MXh,NMXl), - try_sin_((MPi,Xh),Z, NX2, NZ2,MXh,MXl,NMXh), - union_(NZ1,NZ2,NewZ), - union_(NX1,NX2,NewX). - -try_sin_(X,Z,NewX,NewZ,MXS,MXF,MXS) :- sin_((MXS,MXS), X, Z, _, NewX, NewZ),!. -try_sin_(X,Z,[],[],MXS,MXF,MXF). % if sin_ fails, return empty X interval for union - -sin_sector_(lo,Pi_2,(Xl,Xh),Z,(NXl,NXh),NewZ) :- - Xl =< -Pi_2, !, % X intersects lo sector, flip to mid - X1l is max(-pi-Xh,-Pi_2), X1h is -pi-Xl, - sin_prim_((X1l,X1h),Z,(NX1l,NX1h),NewZ), - NXl is -pi-NX1h, NXh is -pi-NX1l. - -sin_sector_(hi,Pi_2,(Xl,Xh),Z,(NXl,NXh),NewZ) :- - Xh >= Pi_2, !, % X intersects hi sector, flip to mid - X1l is pi-Xh, X1h is min(pi-Xl,Pi_2), - sin_prim_((X1l,X1h),Z,(NX1l,NX1h),NewZ), - NXl is pi-NX1h, NXh is pi-NX1l. - -sin_sector_(mid,Pi_2,(Xl,Xh),Z,NewX,NewZ) :- - Xl =< Pi_2, Xh >= -Pi_2, !, % mid sector, valid asin function range - X1l is max(Xl,-Pi_2), X1h is min(Xh, Pi_2), - sin_prim_((X1l,X1h),Z,NewX,NewZ). - -sin_sector_(_Any,_,_X,_Z,[],[]). - -% both sin and asin monotonic, increasing in range -pi/2 to pi/2 -sin_prim_((Xl,Xh),(Zl,Zh),(NXl,NXh),(NZl,NZh)) :- - % Note: can't depend on transcendentals to obey rounding - NZl is max(Zl,max(-1,nexttoward(sin(Xl),-1.0Inf))), - NZh is min(Zh,min( 1,nexttoward(sin(Xh), 1.0Inf))), - NXl is max(Xl,nexttoward(asin(NZl),-1.0Inf)), - NXh is min(Xh,nexttoward(asin(NZh), 1.0Inf)). - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z== cos(X) - -narrowing_op(cos, _, $(Z,X), $(NewZ, NewX)) :- - wrap_(X,2*pi,MX,Xp), !, % fails if X too wide - cos_(MX,Xp,Z,NMX,X1,NewZ), - unwrap_(X1,2*pi,NMX,X2), % project back to original cylinder - non_empty(NewZ), - ^(X,X2,NewX). -% ^(Z1,(-1,1),NewZ). % take care of any rounding beyond range - -narrowing_op(cos, _, $(Z,X), $(NewZ,X)) :- % no narrowing possible, just constrain Z - ^(Z,(-1,1),NewZ). - -cos_((MX,MX), X, Z, (MX,MX), NewX, NewZ) :- - !, % same cylinder, split into 2 interval convex sectors - cos_sector_(neg, X, Z, NXneg, NZneg), - cos_sector_(pos, X, Z, NXpos, NZpos), - union_(NZneg,NZpos,NewZ), - union_(NXneg,NXpos,NewX). - -cos_((MXl,MXh), (Xl,Xh), Z, (NMXl,NMXh), NewX, NewZ) :- - % adjacent cylinders, - Pi is pi, MPi is -pi, - try_cos_((Xl,Pi), Z, NX1, NZ1,MXl,MXh,NMXl), - try_cos_((MPi,Xh),Z, NX2, NZ2,MXh,MXl,NMXh), - union_(NZ1,NZ2,NewZ), - union_(NX1,NX2,NewX). - -try_cos_(X,Z,NewX,NewZ,MXS,MXF,MXS) :- cos_((MXS,MXS), X, Z, _, NewX, NewZ),!. -try_cos_(X,Z,[],[],MXS,MXF,MXF). % if cos_ fails, return empty X interval for union - -cos_sector_(neg,(Xl,Xh),Z,(NXl,NXh),NewZ) :- % Lo is -pi, Hi is 0 - Xl =< 0, !, - X1l is -Xh, X1h is max(0,-Xl), - cos_prim_((X1l,X1h),Z,(NX1l,NX1h),NewZ), - NXl is -NX1h, NXh is -NX1l. - -cos_sector_(pos,(Xl,Xh),Z,NewX,NewZ) :- % Lo is 0, Hi is pi - Xh >=0, !, - X1l is max(0,Xl), - cos_prim_((X1l,Xh),Z,NewX,NewZ). - -cos_sector_(_Any,_X,_Z,[],[]). - -% both cos and acos monotonic decreasing in range 0 to pi -cos_prim_((Xl,Xh),(Zl,Zh),(NXl,NXh),(NZl,NZh)) :- - % Note: can't depend on transcendentals to obey rounding - NZl is max(Zl,max(-1,nexttoward(cos(Xh),-1.0Inf))), - NZh is min(Zh,min( 1,nexttoward(cos(Xl), 1.0Inf))), - NXl is max(Xl,nexttoward(acos(NZh),-1.0Inf)), - NXh is min(Xh,nexttoward(acos(NZl), 1.0Inf)). - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z== tan(X) - -narrowing_op(tan, _, $(Z,X), $(NewZ, NewX)) :- - wrap_(X,pi,MX,Xp), !, % fails if X too wide - tan_(MX,Xp,Z,NMX,X1,NewZ), - unwrap_(X1,pi,NMX,X2), % project back to original cylinder - non_empty(NewZ), - ^(X,X2,NewX). - -narrowing_op(tan, _, ZX, ZX). % no narrowing possible, e.g., not same or adjacent cylinder. - -% both tan and atan monotonic, increasing in range -pi/2 to pi/2 -tan_((MX,MX), (Xl,Xh), (Zl,Zh), (MX,MX), (NXl,NXh), (NZl,NZh)) :- - !, % same cylinder - % Note: can't depend on transcendentals to obey rounding - NZl is max(Zl,nexttoward(tan(Xl),-1.0Inf)), - NZh is min(Zh,nexttoward(tan(Xh), 1.0Inf)), - non_empty(NZl,NZh), - NXl is max(Xl,nexttoward(atan(NZl),-1.0Inf)), - NXh is min(Xh,nexttoward(atan(NZh), 1.0Inf)), - non_empty(NXl,NXh). - -% not same cylinder, must span infinite discontinuity --> no change -%tan_(MX, X, Z, MX, X, Z). - -tan_((MXl,MXh), (Xl,Xh), Z, (NMXl,NMXh), NewX, NewZ) :- -% MXl is MXh-1, % adjacent cylinders - PiBy2 is pi/2, MPiBy2 is -PiBy2, - try_tan_((Xl,PiBy2), Z, NX1, NZ1,MXl,MXh,NMXl), - try_tan_((MPiBy2,Xh), Z, NX2, NZ2,MXh,MXl,NMXh), - union_(NZ1,NZ2,NewZ), % will fail if both empty - union_(NX1,NX2,NewX). - -try_tan_(X,Z,NewX,NewZ,MXS,MXF,MXS) :- tan_((MXS,MXS), X, Z, _, NewX, NewZ), !. -try_tan_(X,Z,[],[],MXS,MXF,MXF). % if tan_ fails, return empty X interval for union - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z== ~X boolean negation (Z and X boolean) - -narrowing_op(not, P, $(Z,X), $(NewZ, NewX)) :- - booleanVal_(Z,ZB), booleanVal_(X,XB), % constrain to boolean - notB_(XB,ZB,NewZ,P), - notB_(NewZ,XB,NewX,P). - %(NewZ=(B,B) -> P=p ; true). % if either was a point, both are now - -notB_((0,0), (_,1), (1,1), p). -notB_((1,1), (0,_), (0,0), p). -notB_((0,1), Z, Z, _). - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==X and Y boolean 'and' - -narrowing_op(and, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :- - booleanVal_(Z,ZB), booleanVal_(X,XB), booleanVal_(Y,YB), - andB_rel_(ZB,XB,YB, NewZ, NewX, NewY, P), !. - -andB_rel_(Z,(1,1),(1,1), NewZ,(1,1),(1,1), p) :- !, ^(Z,(1,1),NewZ). % all cases persistent (points) except last -andB_rel_(Z,(0,0),Y, NewZ,(0,0),Y, p) :- !, ^(Z,(0,0),NewZ). -andB_rel_(Z,X,(0,0), NewZ,X,(0,0), p) :- !, ^(Z,(0,0),NewZ). -andB_rel_((1,1),X,Y, (1,1),NewX,NewY, p) :- !, ^(X,(1,1),NewX), ^(Y,(1,1),NewY). -andB_rel_((0,0),X,(1,1), (0,0),NewX,(1,1), p) :- !, ^(X,(0,0),NewX). -andB_rel_((0,0),(1,1),Y, (0,0),(1,1),NewY, p) :- !, ^(Y,(0,0),NewY). -andB_rel_(Z,X,Y, Z,X,Y, P). % still indeterminate - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==X and Y boolean 'nand' - -narrowing_op(nand, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :- - booleanVal_(Z,ZB), booleanVal_(X,XB), booleanVal_(Y,YB), - nandB_rel_(ZB,XB,YB, NewZ, NewX, NewY, P), !. - -nandB_rel_(Z,(1,1),(1,1), NewZ,(1,1),(1,1), p) :- !, ^(Z,(0,0),NewZ). % all cases persistent except last -nandB_rel_(Z,(0,0),Y, NewZ,(0,0),Y, p) :- !, ^(Z,(1,1),NewZ). -nandB_rel_(Z,X,(0,0), NewZ,X,(0,0), p) :- !, ^(Z,(1,1),NewZ). -nandB_rel_((0,0),X,Y, (0,0),NewX,NewY, p) :- !, ^(X,(1,1),NewX), ^(Y,(1,1),NewY). -nandB_rel_((1,1),X,(1,1), (1,1),NewX,(1,1), p) :- !, ^(X,(0,0),NewX). -nandB_rel_((1,1),(1,1),Y, (1,1),(1,1),NewY, p) :- !, ^(Y,(0,0),NewY). -nandB_rel_(Z,X,Y, Z,X,Y, P). % still indeterminate - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==X or Y boolean 'or' - -narrowing_op(or, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :- - booleanVal_(Z,ZB), booleanVal_(X,XB), booleanVal_(Y,YB), - orB_rel_(ZB,XB,YB, NewZ, NewX, NewY, P), !. - -orB_rel_(Z,(0,0),(0,0), NewZ,(0,0),(0,0), p) :- !, ^(Z,(0,0),NewZ). % all cases persistent (points) except last -orB_rel_(Z,(1,1),Y, NewZ,(1,1),Y, p) :- !, ^(Z,(1,1),NewZ). -orB_rel_(Z,X,(1,1), NewZ,X,(1,1), p) :- !, ^(Z,(1,1),NewZ). -orB_rel_((0,0),X,Y, (0,0),NewX,NewY, p) :- !, ^(X,(0,0),NewX), ^(Y,(0,0),NewY). -orB_rel_((1,1),X,(0,0), (1,1),NewX,(0,0), p) :- !, ^(X,(1,1),NewX). -orB_rel_((1,1),(0,0),Y, (1,1),(0,0),NewY, p) :- !, ^(Y,(1,1),NewY). -orB_rel_(Z,X,Y, Z,X,Y, P). % still indeterminate - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==X nor Y boolean 'nor' - -narrowing_op(nor, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :- - booleanVal_(Z,ZB), booleanVal_(X,XB), booleanVal_(Y,YB), - norB_rel_(ZB,XB,YB, NewZ, NewX, NewY, P), !. - -norB_rel_(Z,(0,0),(0,0), NewZ,(0,0),(0,0), p) :- !, ^(Z,(1,1),NewZ). % all cases persistent (points) except last -norB_rel_(Z,(1,1),Y, NewZ,(1,1),Y, p) :- !, ^(Z,(0,0),NewZ). -norB_rel_(Z,X,(1,1), NewZ,X,(1,1), p) :- !, ^(Z,(0,0),NewZ). -norB_rel_((1,1),X,Y, (1,1),NewX,NewY, p) :- !, ^(X,(0,0),NewX), ^(Y,(0,0),NewY). -norB_rel_((0,0),X,(0,0), (0,0),NewX,(0,0), p) :- !, ^(X,(1,1),NewX). -norB_rel_((0,0),(0,0),Y, (0,0),(0,0),NewY, p) :- !, ^(Y,(1,1),NewY). -norB_rel_(Z,X,Y, Z,X,Y, P). % still indeterminate - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==X xor Y boolean 'xor' - -narrowing_op(xor, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :- - booleanVal_(Z,ZB), booleanVal_(X,XB), booleanVal_(Y,YB), - xorB_rel_(ZB,XB,YB, NewZ, NewX, NewY, P), !. - -xorB_rel_(Z,(B,B),(B,B), NewZ,(B,B),(B,B), p) :- !, ^(Z,(0,0),NewZ). % all cases persistent (points) except last -xorB_rel_(Z,(X,X),(Y,Y), NewZ,(X,X),(Y,Y), p) :- !, ^(Z,(1,1),NewZ). -xorB_rel_((B,B),X,(B,B), (B,B),NewX,(B,B), p) :- !, ^(X,(0,0),NewX). -xorB_rel_((Z,Z),X,(Y,Y), (Z,Z),NewX,(Y,Y), p) :- !, ^(X,(1,1),NewX). -xorB_rel_((B,B),(B,B),Y, (B,B),(B,B),NewY, p) :- !, ^(Y,(0,0),NewY). -xorB_rel_((Z,Z),(X,X),Y, (Z,Z),(X,X),NewY, p) :- !, ^(Y,(1,1),NewY). -xorB_rel_(Z,X,Y, Z,X,Y, P). % still indeterminate - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Z==X -> Y boolean 'implies' - -narrowing_op(imB, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :- - booleanVal_(Z,ZB), booleanVal_(X,XB), booleanVal_(Y,YB), - imB_rel_(ZB,XB,YB, NewZ, NewX, NewY, P), !. - -imB_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 -imB_rel_(Z,(0,0),Y, NewZ,(0,0),Y, p) :- !, ^(Z,(1,1),NewZ). % X=0, Z=1, persistant -imB_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 -imB_rel_(Z,X,(1,1), NewZ,X,(1,1), P) :- !, ^(Z,(1,1),NewZ). % Y=1, Z=1, persistant -imB_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 -imB_rel_(Z,X,Y, Z,X,Y, P). % still indeterminate diff --git a/src/ia_simplify.pl b/src/ia_simplify.pl deleted file mode 100644 index 2d8d9e7..0000000 --- a/src/ia_simplify.pl +++ /dev/null @@ -1,386 +0,0 @@ -/* The MIT License (MIT) - * - * Copyright (c) 2019,2020 Rick Workman - * - * Permission is hereby granted, free of charge, to any person obtaining a copy - * of this software and associated documentation files (the "Software"), to deal - * in the Software without restriction, including without limitation the rights - * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the Software is - * furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included in all - * copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - * SOFTWARE. - */ -% -% Expression variables are Prolog var's; numeric constants are precise. -% Note floats are treated like variables and no arithmetic is performed on them -% - -% negatable operators -negate_op_(+,-). -negate_op_(-,+). - -% power operators -invert_op_(*,/). -invert_op_(/,*). - -% signed multiplier and power values -sign_val_(+,C,C) :- C>=0. -sign_val_(-,C,N) :- N is -C. - -pwr_val_(*,C,C) :- C>=0. -pwr_val_(/,C,N) :- N is -C. - - -% utility to distribute A*(B1+B2) if they have variables in common or A is a constant -% former is safe for intervals (sub-distributive law), latter is valid -distribute_(C,B,Exp) :- - number(C), !, - B=..[AddOp,B1,B2], negate_op_(AddOp,_), % watch out for vars - DExp=..[AddOp,C*B1,C*B2], - simplify(DExp,Exp). -distribute_(A,B,Exp) :- % not always a good thing, e.g., X*(X-1) may be better than X**2-X - B=..[Op,B1,B2], negate_op_(Op,_), % watch out for vars - shared_vars_(A,B), !, % any shared vars - DExp=..[Op,A*B1,A*B2], - simplify(DExp,Exp). - -% utility for (in)equality reduction -normalize_(A,B,Op,ROp,Exp) :- - B==0, !, % RHS = 0 - n_build_(Op, A, 0, Exp). -normalize_(A,B,Op,ROp,Exp) :- - A==0, !, % LHS = 0 - n_build_(ROp, B, 0, Exp). -normalize_(A,B,Op,ROp,Exp) :- - occurs_in_(A,B), !, % RHS is shared var - n_build_(Op, A-B, 0, Exp). -normalize_(A,B,Op,ROp,Exp) :- - occurs_in_(B,A), !, % LHS is shared var - n_build_(ROp, B-A, 0, Exp). -normalize_(A,B,Op,ROp,Exp) :- - compound(A), compound(B), % LHS and RHS are expressions with shared vars - shared_vars_(A,B), !, - n_build_(Op, A-B, 0, Exp). -normalize_(A,B,Op,ROp,Exp) :- % everything else, leave LHS and RHS alone - simplify(B,RHS), - n_build_(Op, A, RHS, Exp). - -occurs_in_(Exp, V) :- var(V), term_variables(Exp,Vs), shared_var_(Vs,V). - -n_build_(Op, L, RHS, Exp) :- simplify(L,LHS), Exp=..[Op,LHS,RHS]. - -shared_vars_(A,B) :- - term_variables(A,AVs), - term_variables(B,BVs), - shared_vars_in_(AVs,BVs). - -shared_vars_in_([AV|_AVs],BVs) :- - shared_var_(BVs,AV), !. -shared_vars_in_([_AV|AVs],BVs) :- - shared_vars_in_(AVs,BVs). - -shared_var_([BV|_BVs],AV) :- - AV == BV, !. -shared_var_([_BV|BVs],AV) :- - shared_var_(BVs,AV). - - -% -% simplify/2 -% -simplify(A,A) :- var(A),!. - -simplify(C,C) :- rational(C),!. % includes integers - -% possible distribute -simplify(A*B,Exp) :- - compound(B), - distribute_(A,B,Exp), !. -simplify(A*B,Exp) :- - compound(A), - distribute_(B,A,Exp), !. - -% simplify equalities and inequalities -simplify(A==B,Exp) :- - normalize_(A,B,==,==,Exp), !. - -simplify(A==,Exp), !. - -simplify(A>=B,Exp) :- - normalize_(A,B,>=,=<,Exp), !. - -simplify(A,Exp), !. - -simplify(A>B,Exp) :- - normalize_(A,B,>,<,Exp), !. -/* -% simplify "cascaded" divisions A/B/C = (A/B)/C = A*C/B -simplify(A/B,Exp) :- -% compound(B), B=Bn/Bd, !, -% simplify(A*Bd/Bn,Exp). - compound(A), A=An/Ad, !, - simplify(An*B, E1), - simplify(E1/Ad,Exp). -*/ -% -% General purpose simplify -% -simplify(E,Exp) :- - collect_exp_(E, L/L, Es/[]), % collect terms in a difference list - collect_exp_items(Es,Is), % collect like terms - reduce_exp_items_(Is,Exp), % and reduce back to expression - !. - -% use difference list to collect terms of an expression -collect_exp_(A, List, Head/NewT) :- - var(A), !, - List=Head/[term(A,1)|NewT]. % separate to keep head unification cheap - -collect_exp_(C, List, Head/NewT) :- - rational(C), !, - List=Head/[C|NewT]. - -collect_exp_(A, List, Head/NewT) :- % floats as terms, can collect but not do arithmetic - float(A), !, - List=Head/[term(A,1)|NewT]. - -collect_exp_(A+B, List, NewList) :- - !, - left_exp_(A,AE), - collect_exp_(AE,List,ListA), - simplify(B,BT), - collect_exp_(BT,ListA,NewList). - -collect_exp_(A-B,List,NewList) :- - !, - left_exp_(A,AE), - collect_exp_(AE,List,ListA), - simplify(B,BT), - collect_neg_(BT,ListA,NewList). - -collect_exp_(-B,List,NewList) :- - simplify(B,BT), - collect_neg_(BT,List,NewList). - -collect_exp_(T, List/[Term|NewT], List/NewT) :- - simplify_term(T,Term). - -left_exp_(A,A) :- var(A), !. -left_exp_(A,A) :- A=..[AOp|_], negate_op_(AOp,_), !. -left_exp_(A,AE) :- simplify(A,AE). - -collect_neg_(BT,List/[N|NewT], List/NewT) :- - rational(BT), - N is -BT. % rational constant expressed as AT-BT -collect_neg_(BT,ListA/T,ListA/NewT) :- - collect_exp_(BT,P/P,ListB), - negate_term_(ListB,T/NewT). - -negate_term_(T/T,NewT/NewT) :- var(T). -negate_term_([term(V,N)|Es]/T,[term(V,NN)|NEs]/NewT) :- - NN is -N, - negate_term_(Es/T,NEs/NewT). -negate_term_([N|Es]/T,[NN|NEs]/NewT) :- - NN is -N, - negate_term_(Es/T,NEs/NewT). - -collect_exp_items([],[]). -collect_exp_items([E|Es],[NE|NEs]) :- - collect_exp_item_(Es,E,NE,ENxt), !, - collect_exp_items(ENxt,NEs). - -collect_exp_item_([],E,E,[]). -collect_exp_item_([V|Es],U,Eo,EsNxt) :- - rational(U), rational(V), % constant values must be precise to add - S is U+V, - collect_exp_item_(Es,S,Eo,EsNxt). -collect_exp_item_([term(V,N1)|Es],term(U,N2),Eo,EsNxt) :- - V==U, - N is N1+N2, - collect_exp_item_(Es,term(V,N),Eo,EsNxt). -collect_exp_item_([Ei|Es],E,Eo,[Ei|EsNxt]) :- % Note that imprecise floats are not added - collect_exp_item_(Es,E,Eo,EsNxt). - -reduce_exp_items_([V],V) :- % terminate: var - var(V), !. -reduce_exp_items_([C],C) :- % terminate: constant - rational(C), !. -reduce_exp_items_([term(V,N)],Exp) :- !, % terminate: single term(_..) - reduce_exp_items_([0,term(V,N)],Exp). -reduce_exp_items_([Exp],Exp). % terminate: final expression -reduce_exp_items_([T1,T2|Ts],Exp) :- % 2 or more terms, combine and continue - reduce_exp_item_(T1,Op1,Exp1), - reduce_exp_item_(T2,Op2,Exp2), - build_exp_(Exp1,Exp2,Op1,Op2,ExpN), !, % ExpN =.. [Op,Exp1,Exp2], !, - reduce_exp_items_([ExpN|Ts],Exp). - -build_exp_(Z,Exp2,Op1,+,Exp2) :- Z==0. -build_exp_(Z,Exp2,Op1,-,NExp2) :- Z==0, build_Nexp_(Exp2,NExp2). -build_exp_(Exp1,Z,+,Op2,Exp1) :- Z==0. -build_exp_(Exp1,Z,-,Op2,NExp1) :- Z==0, build_Nexp_(Exp1,NExp1). -build_exp_(Exp1,Exp2,+,+,Exp1+Exp2). -build_exp_(Exp1,Exp2,+,-,Exp1-Exp2). -build_exp_(Exp1,Exp2,-,+,Exp2-Exp1). -build_exp_(Exp1,Exp2,-,-,NExp1-Exp2) :- build_Nexp_(Exp1,NExp1). - -build_Nexp_(Exp, NExp) :- % constant - rational(Exp), - NExp is -Exp. -build_Nexp_(Exp, NExp) :- % * or / expression, find head - nonvar(Exp), Exp =.. [Op,A1,A2], invert_op_(Op,_), - build_Nexp_(A1,NA1), - NExp =.. [Op,NA1,A2]. -build_Nexp_(Exp, -Exp). % other expression or var - - -reduce_exp_item_(V, +, V) :- var(V). -reduce_exp_item_(term(V,0), +, 0). -reduce_exp_item_(term(V,1), +, V). -reduce_exp_item_(term(V,-1), -, V). -reduce_exp_item_(term(V,N), Op, T) :- mult_term_(V, N, Op, T). -reduce_exp_item_(R, Op, M) :- rational(R), val_sign_(R,Op,M). % constants -reduce_exp_item_(-Exp, -, Exp). -reduce_exp_item_(Exp, +, Exp). - -mult_term_(T, N, Op, M*T) :- var(T), val_sign_(N,Op,M). -mult_term_(T, N, Op, V) :- ground(V), T is V*N. -mult_term_(X*Y, N, Op, Exp*Y) :- mult_term_(X,N,Op,Exp). % maintain Op associativity -mult_term_(X/Y, N, Op, Exp/Y) :- mult_term_(X,N,Op,Exp). -mult_term_(T, N, Op, M*T) :- val_sign_(N,Op,M). - -val_sign_(V,Op,Vp) :- Vp is abs(V), (V<0 -> Op = - ; Op = +). - -% -% simplify a term to an "expression" of term structures and constants -% -simplify_term(T,Term) :- - collect_term_(T, L/L, Ts/[]), % collect in a difference list - collect_term_items_([1|Ts],[C|Is]), % ensure there's a constant multiplier and collect like terms - sort(0, @=<, Is, SIs), - reduce_term_items_([1|SIs],ITerm), % reduce back to expression - % this does constant folding if reduction resulted in a constant - (rational(ITerm) -> Term is ITerm*C ; Term = term(ITerm,C)), - !. - -collect_term_(A, List, Head/NewT) :- - var(A), !, - List=Head/[elem(A,1)|NewT]. - -collect_term_(A, List, Head/NewT) :- - rational(A), !, - List=Head/[A|NewT]. - -collect_term_(A, List, Head/NewT) :- - float(A), !, - List=Head/[elem(A,1)|NewT]. - -collect_term_(-A, List, ListA/NewT) :- % -Term as Term * -1 for reducing signs - simplify(A,AT), collect_term_(AT,List,ListA/[-1|NewT]). - -collect_term_(A**N, List, Head/NewT) :- % possible special case of user written element - simplify(N,NT), rational(NT), - simplify(A,AT), - simplify_pwr_(AT,NT,Term), - List=Head/[Term|NewT]. - -collect_term_(A*B,List,NewList) :- - !, - left_term_(A,AE), - collect_term_(AE,List,ListA), - simplify(B,BT), - collect_term_(BT,ListA,NewList). - -collect_term_(A/B,List,ListA/NewT) :- - !, - left_term_(A,AE), - simplify(B,BT), - collect_term_(AE,List,ListA/T), - collect_term_(BT,P/P,ListB), - invert_term_(ListB,T/NewT). - -collect_term_(E,List/[elem(Exp,1)|NewT], List/NewT) :- - E =.. [F|IArgs], - simplify_list_(IArgs,OArgs), - Exp =.. [F|OArgs]. - -% simplify_pwr_: NT rational, -simplify_pwr_(AT,NT,Term) :- rational(AT), !, % constant folding - Term is AT**NT. -simplify_pwr_(Exp**P,NT,elem(Exp,Pwr)) :- integer(NT), rational(P), !, % (X**P)**NT - Pwr is NT*P. -simplify_pwr_(AT,NT,elem(AT,NT)). - -left_term_(A,A) :- var(A), !. -left_term_(A,A) :- A=..[AOp|_], invert_op_(AOp,_), !. -left_term_(A,AE) :- simplify(A,AE). - -simplify_list_([],[]). -simplify_list_([I|IArgs],[O|OArgs]) :- - simplify(I,O), - simplify_list_(IArgs,OArgs). - -invert_term_(T/T,NewT/NewT) :- var(T),!. -invert_term_([elem(V,N)|Es]/T,[elem(V,NN)|NEs]/NewT) :- !, - NN is -N, - invert_term_(Es/T,NEs/NewT). -invert_term_([N|Es]/T,[NN|NEs]/NewT) :- - NN is 1/N, - invert_term_(Es/T,NEs/NewT). - -collect_term_items_([],[]). -collect_term_items_([E|Es],[NE|NEs]) :- - collect_term_item_(Es,E,NE,ENxt), - collect_term_items_(ENxt,NEs). - -collect_term_item_([],E,E,[]). -collect_term_item_([V|Es],U,Eo,EsNxt) :- - rational(U), rational(V), % precise for multiply - S is U*V, - collect_term_item_(Es,S,Eo,EsNxt). -collect_term_item_([elem(V,N1)|Es],elem(U,N2),Eo,EsNxt) :- - V==U, - N is N1+N2, - collect_term_item_(Es,elem(V,N),Eo,EsNxt). -collect_term_item_([Ei|Es],E,Eo,[Ei|EsNxt]) :- - collect_term_item_(Es,E,Eo,EsNxt). - -reduce_term_items_([Exp],Exp) :- % terminating variable - var(Exp), !. -reduce_term_items_([Exp],Exp) :- % terminating constant - rational(Exp), !. -reduce_term_items_([elem(V,N)],Exp) :- !, % terminating single elem(_..) - reduce_term_items_([1,elem(V,N)],Exp). -reduce_term_items_([Exp],Exp). % terminating final expression -reduce_term_items_([T1,T2|Ts],Exp) :- - reduce_term_item_(T1,Exp1,_), - reduce_term_item_(T2,Exp2,Op), - build_term_(Exp1,Exp2,Op,ExpN), % ExpN =.. [Op,Exp1,Exp2], - !, - reduce_term_items_([ExpN|Ts],Exp). - - -build_term_(E, C, Op, Exp) :- var(E), Exp =.. [Op,E,C]. -build_term_(1, A**B,/, A**NB) :- rational(B), NB is -B . -build_term_(1, Exp, *, Exp). -build_term_(1, Exp, /, 1/Exp). -build_term_(One/B,C, *, C/B) :- One==1. -build_term_(E1, E2, Op, Exp) :- Exp =.. [Op,E1,E2]. - -reduce_term_item_(V, V, *) :- var(V),!. % already reduced to var -reduce_term_item_(elem(_, 0), 1, *). -reduce_term_item_(elem(V, 1), V, *). -reduce_term_item_(elem(V,-1), V, /). -reduce_term_item_(elem(V, E), V**P, Op) :- P is abs(E), (E>0 -> Op= * ; Op= /). -reduce_term_item_(R, R, *). % catchall: constant or expression diff --git a/src/ia_utilities.pl b/src/ia_utilities.pl deleted file mode 100755 index ac77e09..0000000 --- a/src/ia_utilities.pl +++ /dev/null @@ -1,797 +0,0 @@ -/* The MIT License (MIT) - * - * Copyright (c) 2019,2020 Rick Workman - * - * Permission is hereby granted, free of charge, to any person obtaining a copy - * of this software and associated documentation files (the "Software"), to deal - * in the Software without restriction, including without limitation the rights - * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the Software is - * furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included in all - * copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - * SOFTWARE. - */ - -% -% compatability predicates -% - -% list filter -list(X) :- is_list(X). - -% -% print_interval(Term), print_interval(Stream,Term) -% prints Term to optional Stream with intervals expanded to domains -% uses format/3 so extended Stream options, e.g., atom(A), are supported -% -print_interval(Term) :- print_interval(current_output,Term). - -print_interval(Stream,Term) :- - copy_term_nat(Term,Out), % copy of Term without attributes for output - term_variables(Term,TVars), - term_variables(Out,OVars), - name_vars_(TVars,OVars,0), % name variables, intervals replaced by declarations - format(Stream,'~w',[Out]). - -name_vars_([],[],_). -name_vars_([TVar|TVars],[OVar|OVars],N) :- - (domain(TVar,Dom) - -> OVar = (V::Dom) - ; OVar = V - ), - number_chars(N,C),atom_chars(V,['V'|C]), - N1 is N+1, - name_vars_(TVars,OVars,N1). - -% -% SWIP attribute portray hook - used if write_attributes=portray -% -:- set_prolog_flag(write_attributes, portray). % generally useful - -attr_portray_hook(interval(Type,Val,N,F),_Int) :- - interval_domain_(Type,Val,Dom), - write(Dom). - -% -% SWIP hook for top level output (ignores flags) -% two modes : -% Verbose=true, output all interval domains and associated constraints -% Verbose=false, only output query vars with no constraints -% -project_attributes(QueryVars, _) :- % mark QueryVars for output when Verbose=false - flag_query_vars_(QueryVars). - -flag_query_vars_([]). -flag_query_vars_([QV|QueryVars]) :- - get_interval_flags_(QV,Flags), - set_interval_flags_(QV,[queryVar|Flags]), !, % add 'queryVar' to Flags - flag_query_vars_(QueryVars). -flag_query_vars_([QV|QueryVars]) :- - flag_query_vars_(QueryVars). - -attribute_goals(X) --> - {current_prolog_flag(clpBNR_verbose,Verbose), - domain_goals_(Verbose,X,Goals) - }, - Goals. - -domain_goals_(true,X,Cs) :- % Verbose=true, QueryVars and constraints output - interval_object(X, Type, Val, Nodes), !, - interval_domain_(Type, Val, Dom), - constraints_(Nodes,X,NCs), % reverse map to list of original constraint - to_comma_exp_(NCs,X::Dom,Cs). -domain_goals_(false,X,[X::Dom]) :- % Verbose=false, only QueryVars output - get_interval_flags_(X,Flags), - memberchk(queryVar,Flags), !, % set by project_attributes - interval_object(X, Type, Val, _), - intValue_(Val,Type,Dom). -domain_goals_(_,X,[]). % catchall but normally non-query attvar, Verbose=false - -constraints_([Node],_,[]) :- var(Node), !. % end of indefinite list -constraints_([node(COp,P,_,Args)|Nodes],X,[C|Cs]) :- - var(P), % skip persistent nodes (already in value) - term_variables(Args,[V|_]), % this ensures constraints only get output once - V==X, - fmap_(Op,COp,_,_,_), !, - remap_(Op,Args,C), - constraints_(Nodes,X,Cs). -constraints_([_|Nodes],X,Cs) :- - constraints_(Nodes,X,Cs). - -to_comma_exp_([],Decl,[Decl]). -to_comma_exp_([N|NCs],Decl,[(Decl,{Cs})]) :- - to_comma_cexp_(NCs,N,Cs). - -to_comma_cexp_([],In,In). -to_comma_cexp_([C|Cs],In,Out) :- - to_comma_cexp_(Cs,(In,C),Out). - -% -% Simplified output format for interval ranges -% -% 1. Any 'real' rational bounds output as floats (conservatively rounded). -% 2. Any subnormal bounds converted to zero. -% 3. Sufficiently narrow reals output in elipsis format. -% 4. Any real intervals with bounds zero or subnormal output as 0.0... -% 5. Query interval variables output as V = V::Domain -% 6. Subterm intervals as sub-terms printed as domain (Type(L,H)). -% -:- ( % to avoid quoting ellipsis postfix format, also increase answer depth a bit - current_prolog_flag(answer_write_options,Opts), - lists:subtract(Opts,[quoted(_)],NOpts0), - lists:subtract(NOpts0,[max_depth(D)],NOpts), - (var(D) -> ND=20 ; ND is max(D,20)), - set_prolog_flag(answer_write_options, [quoted(false),max_depth(ND)|NOpts]) - ). - -intValue_((0,1),integer,boolean). % boolean -intValue_((L,H),real,Out...) :- % virtual zero (zero or subnormal) - zero_float_(L), - zero_float_(H), !, - format(chars(Zero),"~16f",0.0), - string_chars(Out,[' '|Zero]). % space prefix -intValue_((L,H),real,Out...) :- % two floats with minimal leading match - float_chars_(L,LC), - float_chars_(H,HC), - sign_chars_(LC,HC,ULC,UHC,Codes/Match), - leading_zero(ULC,UHC,ZLC), - matching_(ZLC,UHC,Match,0,Dec,MLen), - MLen>Dec+1, !, % minimum of one matching digit after decimal point - string_codes(Out,Codes). -intValue_((L,H),Type,Dom) :- % default, just convert rationals - rational_fraction_(L,FL), - rational_fraction_(H,FH), - interval_domain_(Type, (FL,FH), Dom). - -zero_float_(B) :- float(B), float_class(B,C), (C=zero ; C=subnormal). - -float_chars_(B,Cs) :- - rational_fraction_(B,F), - float(F), float_class(F, normal), - format(chars(Cs),'~16f',F), % same length after '.' (pads with trailing 0's) - length(Cs,Len), Len=<32. % reverts to precise format if too long - -rational_fraction_(B,FB) :- - rational(B,_,D), D \= 1, !, % non-integer rational - FB is float(B). -rational_fraction_(F,0.0) :- - float(F), - float_class(F,subnormal), !. -rational_fraction_(B,B). - -sign_chars_(['-'|LC],['-'|HC],HC,LC,[' ','-'|T]/T) :- !. % negate, save sign and reverse bounds -sign_chars_(LC,HC,LC,HC,[' '|T]/T). % need a character that doesn't print: 2=STX (Start of Text) - -leading_zero(['9'|ULC],['1','0'|_],['0','9'|ULC]) :- !. -leading_zero(ULC,_,ULC). - -% Note: match should never be exact (L=H) so [] case not required -% matching_([],[],[],N,Dec,N). -matching_([C,'.'|LCs],[C,'.'|HCs],[C,'.'|Cs],N,Dec,Nout) :- !, % absorbing "." - digit_(C,_), - Dec is N+1, - N1 is N+2, - matching_(LCs,HCs,Cs,N1,Dec,Nout). -matching_([LC,'.',LC1|LCs],[HC,'.',HC1|HCs],[HC,'.'|Cs],N,Dec,Nout) :- !, % absorbing "." - digit_match_(LC,LC1,HC,HC1), - Dec is N+1, - N1 is N+2, - matching_([LC1|LCs],[HC1|HCs],Cs,N1,Dec,Nout). -matching_([C|LCs],[C|HCs],[C|Cs],N,Dec,Nout) :- !, % matching digit - digit_(C,_), - N1 is N+1, - matching_(LCs,HCs,Cs,N1,Dec,Nout). -matching_([LC,LC1|LCs],[HC,HC1|HCs],[HC|Cs],N,Dec,Nout) :- % match after rounding - digit_match_(LC,LC1,HC,HC1), !, - N1 is N+1, - matching_([LC1|LCs],[HC1|HCs],Cs,N1,Dec,Nout). -matching_(LCs,HCs,[],N,Dec,N) :- % non-matching after '.' - nonvar(Dec). - -digit_match_(LC,LC1,HC,HC1) :- % rounding test if first digits different - digit_(LC,L), digit_(LC1,L1), digit_(HC,H), digit_(HC1,H1), - H is (L+1) mod 10, - L1 >= 5, H1 < 5. - -% char test for properly formatted number -digit_(DC,D) :- atom_number(DC,D), integer(D), 0==pt(XfMid),V2), % Step 6., 7. & 8. for V2 - tree_insert(ZTree1,V2,ZTree2), % Step 10. for V2 - select_min(ZTree2,NxtY,ZTreeY), % Step 9. and 11., remove Y from tree - iterate_MS(Z,Xs,P,NxtY,ZTreeY). % Step 13. -iterate_MS(Z,Xs,P,Zl-(Zh,XVs),ZTree) :- % termination criteria (Step 12.) met - {Zl=K,!, - tree_insert(R, NK-NData, NewR). -tree_insert(n(L,R,K-(U,Data)), K-(NU,NData), n(NewL,R,K-(U,Data))) :- NU=U - tree_insert(R, NData, NewR). - -select_min(n(L,R,Min),Min,R) :- var(L), !, - nonvar(Min). % fail if tree was empty -select_min(n(L,R,KData),Min,n(NewL,R,KData)) :- - select_min(L,Min,NewL). - -% construct box from interval bounds -box_([],[]). -box_([X|Xs],[R|Rs]) :- - getValue(X,R), - box_(Xs,Rs). - -% find widest interval and its midpoint -widest_MS([X|Xs],[XV|XVs],Xf,XfMid) :- - widest1_MS(Xs,XVs,XV,X,Xf,XfMid). - -widest1_MS([],[],(L,H),Xf,Xf,XfMid) :- - midpoint_(L,H,XfMid), !. -widest1_MS([X|Xs],[(L,H)|XVs],(L0,H0),X0,Xf,XfMid) :- - (H-L) > (H0-L0), !, - widest1_MS(Xs,XVs,(L,H),X,Xf,XfMid). -widest1_MS([X|Xs],[XV|XVs],W0,X0,Xf,XfMid) :- - widest1_MS(Xs,XVs,W0,X0,Xf,XfMid). - -% Midpoint test, Step 11+: -% Discard all pairs (Z,z) from the list that satisfy F(C) geometric/arithmetic mean -% M is min(sqrt(abs(L)),abs(L)/2)*sign(L)+min(sqrt(abs(H)),abs(H)/2)*sign(H). - -% -% splitsolve(Int) - joint search on list of intervals -% simple split, e.g., no filtering on splutions, etc. -% -splitsolve(X) :- - current_prolog_flag(clpBNR_default_precision,P), - splitsolve(X,P). - -splitsolve(X,P) :- - number(X), !. % already a point value -splitsolve(X,P) :- - interval(X), !, % if single interval, put it into a list - Err is 10**(-P), - simplesolveall_([X],Err). -splitsolve(X,P) :- % assumed to be a list - flatten_list(X,[],XF), % flatten before iteration - Err is 10**(-P), - simplesolveall_(XF,Err). - -% flatten list(s) using difference lists -flatten_list(V,Tail,[V|Tail]) :- var(V), !. -flatten_list([],Tail,Tail) :- !. -flatten_list([H|T], Tail, List) :- !, - flatten_list(H, HTail, List), - flatten_list(T, Tail, HTail). -flatten_list(N,Tail,[N|Tail]). - -simplesolveall_(Xs,Err) :- - select_wide_(Xs,D1,X), - interval_object(X, Type, (Xl,Xh), _), - midpoint_(Xl,Xh,Xmid), - choice_generator_(Type,X,(Xl,Xh),Xmid,Err,Choices), - !, - Choices, % generate alternatives - simplesolveall_(Xs,Err). -simplesolveall_(Xs,Err). % terminated - -select_wide_([X],_,X) :- !. % select last remaining element -select_wide_([X1,X2|Xs],D1,X) :- % compare widths and discard one interval - (var(D1) -> delta(X1,D1) ; true), - delta(X2,D2), - (D1 >= D2 - -> select_wide_([X1|Xs],D1,X) - ; select_wide_([X2|Xs],D2,X) - ). - -% bisect_interval(real,X,[Xl,Xh],Xmid,Err,({X= true ; absolve_l(X,Type,Delta,1,Limit)), - (not(not(upper_bound(X))) -> true ; absolve_r(X,Type,Delta,1,Limit)). - -absolve([],_). % extends to lists -absolve([X|Xs],Lim):- absolve(X,Lim),!, absolve(Xs,Lim). - -delta_(integer,(L,U),D) :- D is U div 2 - L div 2. -delta_(real, (L,U),D) :- D is U/2 - L/2. - -absolve_l(X, Type, DL, NL, Limit):- NL=pt(Split)}, % so X must be > - absolve_l(X,Type, DL1, NL1, Limit). -absolve_l(_,_,_,_,_). % final result - -absolve_r(X, Type, DU, NU, Limit):- NU=pt(Split)}),!, - {X= NPt=SPt ; split_real_lo(X,(L,SPt),NPt,Err)). - -split_real_hi(X,(Pt,H),NPt,Err) :- % search upper range for a split point - splitinterval_real_((Pt,H),SPt,Err), !, - (X\=SPt -> NPt=SPt ; split_real_hi(X,(SPt,H),NPt,Err)). - - -% -% splitinterval_integer_ and splitinterval_real_ require ! at point of call. -% -splitinterval_integer_((L,H),0) :- - L < 0, H > 0. % contains 0 but not bounded by 0 -splitinterval_integer_((-1.0Inf,H),Pt) :- - Pt is H*10-5. % subtract 5 in case H is 0. (-5, -15, -105, -1005, ...) -splitinterval_integer_((L,-1.0Inf),Pt) :- - Pt is L*10+5. % add 5 in case L is 0. (5, 15, 105, 1005, ...) -splitinterval_integer_((L,H),Pt) :- - H-L >= 16, % don't split ranges smaller than 16 - Pt is (L div 2) + (H div 2). % avoid overflow - -splitinterval_real_((L,H),0,E) :- % if interval contains 0, split on (precisely) 0. - L < 0, H > 0, !, % cut in case error criteria fails - (H-L) > E. % fail if width is less than error criteria, overflow is OK - -splitinterval_real_((-1.0Inf,H),Pt,_) :- % neg. infinity to H=<0 - !, % if following overflows, split failed - Pt is H*10-1, % subtract 1 in case H is 0. (-1, -11, -101, -1001, ...) - Pt > -1.0Inf. % Pt must be finite -splitinterval_real_((L,1.0Inf),Pt,_) :- % L>=0 to pos. infinity - !, % if following overflows, split failed - Pt is L*10+1, - Pt < 1.0Inf. % Pt must be finite - -splitinterval_real_((L,H),Pt,E) :- % finite L,H, positive or negative but not split, Pt\=0. - \+ chk_small(L,H,E), % only split if not small - splitMean_(L,H,Pt), !. - -%splitinterval_real_([L,H],Pt,E) :- -% writeln('FAIL'([L,H],Pt,E)),fail. - -% (approx.) geometric mean of L and H (same sign) -splitMean_(L,H,M) :- L >=0, !, % positive range - (L=:=0 -> M is min(H/2, sqrt( H)) ; M is sqrt(L)*sqrt(H)). -splitMean_(L,H,M) :- % negative range - (H=:=0 -> M is max(L/2,-sqrt(-L)) ; M is -sqrt(-L)*sqrt(-H)). - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% partial_derivative(E, X, D) -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% Perform symbolic differentiation followed by simplify/2: D = dE/dX -% -% located in utilities as it's a common requirement for applying series expansion -% to improve convergance. Avoids necessity of exposing simplify/2. -% Note that the set of rules only covers functions supported by {}. -% -partial_derivative(F,X,DFX) :- - pd(F,X,DD), - simplify(DD,DFX). - % simplify(DD,DX), - % simplify(DX,DFX). % post re-simplify for insurance?? - -pd(Y,X,D) :- - var(Y), (X==Y -> D=1 ; D=0), !. - -pd(C,_,0) :- % DX = 0 - atomic(C), !. % atom, number, .. - -pd(-U,X,DX) :- % DX = -DU - pd(U,X,DU), - pd_minus(DU,DX). - -pd(U+V,X,DX) :- % DX = DU+DV - pd(U,X,DU), pd(V,X,DV), - pd_add(DU,DV,DX). - -pd(U-V,X,DX) :- % DX = DU-DV - pd(U,X,DU), pd(V,X,DV), - pd_sub(DU,DV,DX). - -pd(U*V,X,DX) :- % DX = V*DU+U*DV - pd(U,X,DU), pd(V,X,DV), - pd_mul(V,DU,VDU), - pd_mul(U,DV,UDV), - pd_add(VDU,UDV,DX). - -pd(U/V,X,DX) :- % DX = (V*DU-U*DV)/(V**2) - pd(U,X,DU), pd(V,X,DV), - pd_mul(V,DU,VDU), - pd_mul(U,DV,UDV), - pd_sub(VDU,UDV,DXN), - pd_pow(V,2,DXD), - pd_div(DXN,DXD,DX). - -pd(U**N,X,DX) :- % DX = (N*U**(N-1))*DU - pd(U,X,DU), - pd_sub(N,1,N1), %N1 is N-1, - pd_pow(U,N1,UN1), - pd_mul(N,UN1,DX1), - pd_mul(DX1,DU,DX). - -pd(log(U),X,DX) :- % DX = DU/U - pd(U,X,DU), - pd_div(DU,U,DX). - -pd(exp(U),X,DX) :- % DX = exp(U)*DU - pd(U,X,DU), - pd_exp(U,ExpU), - pd_mul(ExpU,DU,DX). - -pd(sqrt(U),X,DX) :- % DX = DU/(2*sqrt(U)) - pd(U,X,DU), - pd_sqrt(U,SqrtU), - pd_mul(2,SqrtU,DXD), - pd_div(DU,DXD,DX). - -pd(sin(U),X,DX) :- % DX = cos(U)*DU - pd(U,X,DU), - pd_cos(U,CosU), - pd_mul(CosU,DU,DX). - -pd(cos(U),X,DX) :- % DX = -sin(U)*DU - pd(U,X,DU), - pd_sin(U,SinU), - pd_mul(SinU,DU,DSU), - pd_minus(DSU,DX). - -pd(tan(U),X,DX) :- % DX = DU/cos(U)**2 - pd(U,X,DU), - pd_cos(U,CosU), - pd_pow(CosU,2,DXD), - pd_div(DU,DXD,DX). - -pd(asin(U),X,DX) :- % DX = DU/sqrt(1-U**2) - pd(U,X,DU), - pd_pow(U,2,USQ), - pd_sub(1,USQ,USQ1), - pd_sqrt(USQ1,DXD), - pd_div(DU,DXD,DX). - -pd(acos(U),X,DX) :- % DX = -DU/sqrt(1-U**2) = -D(asin(U)) - pd(asin(U),X,DasinX), - pd_minus(DasinX,DX). - -pd(atan(U),X,DX) :- % DX = DU/(1+U**2) - pd(U,X,DU), - pd_pow(U,2,USQ), - pd_add(1,USQ,USQ1), - pd_div(DU,USQ1,DX). - - -% optimizations -% also facilitates compiled arithmetic, avoids using catch/3 for instantiation errors -pd_minus(DU,DX) :- ground(DU) -> DX is -DU ; DX = -DU. - -pd_add(DU,DV,DV) :- DU==0, !. -pd_add(DU,DV,DU) :- DV==0, !. -pd_add(DU,DV,DX) :- ground((DU,DV)) -> DX is DU+DV ; DX = DU+DV. - -pd_sub(DU,DV,DX) :- DU==0, !, pd_minus(DV,DX). -pd_sub(DU,DV,DU) :- DV==0, !. -pd_sub(DU,DV,DX) :- ground((DU,DV)) -> DX is DU-DV ; DX = DU-DV. - -pd_mul(DU,DV,0) :- DU==0, !. -pd_mul(DU,DV,0) :- DV==0, !. -pd_mul(DU,DV,DV) :- DU==1, !. -pd_mul(DU,DV,DU) :- DV==1, !. -pd_mul(DU,DV,DX) :- DU== -1, !, pd_minus(DV,DX). -pd_mul(DU,DV,DX) :- DV== -1, !, pd_minus(DU,DX). -pd_mul(DU,DV,DX) :- ground((DU,DV)) -> DX is DU*DV ; DX = DU*DV. - -pd_div(DU,DV,0) :- DU==0, !. -pd_div(DU,DV,0) :- DV==0, !, fail. -pd_div(DU,DV,DU) :- DV==1, !. -pd_div(DU,DV,DU) :- DV== -1, !, pd_minus(DU,DX). -pd_div(DU,DV,DX) :- ground((DU,DV)) -> DX is DU/DV ; DX = DU/DV. - -pd_pow(DU,DV,0) :- DU==0, !. -pd_pow(DU,DV,1) :- DV==0, !. -pd_pow(DU,DV,1) :- DU==1, !. -pd_pow(DU,DV,DU) :- DV==1, !. -pd_pow(DU,DV,DX) :- ground((DU,DV)) -> DX is DU**DV ; DX = DU**DV. - -pd_exp(DU,DX) :- ground(DU) -> DX is exp(DU) ; DX = exp(DU). - -pd_sqrt(DU,DX) :- ground(DU) -> DX is sqrt(DU) ; DX = sqrt(DU). - -pd_sin(DU,DX) :- ground(DU) -> DX is sin(DU) ; DX = sin(DU). - -pd_cos(DU,DX) :- ground(DU) -> DX is cos(DU) ; DX = cos(DU). - -% -% safe declarations must be after the predicate definitions -% -% the following use meta-predicates -sandbox:safe_primitive(clpBNR:solve(_)). -sandbox:safe_primitive(clpBNR:solve(_,_)). -sandbox:safe_primitive(clpBNR:splitsolve(_)). -sandbox:safe_primitive(clpBNR:splitsolve(_,_)). -sandbox:safe_primitive(clpBNR:global_minimum(_,_)). -sandbox:safe_primitive(clpBNR:global_minimum(_,_,_)). -sandbox:safe_primitive(clpBNR:global_maximum(_,_)). -sandbox:safe_primitive(clpBNR:global_maximum(_,_,_)). -sandbox:safe_primitive(clpBNR:partial_derivative(_,_,_)).