1/*** Functionality replaced by module(simplify_constraint) - remove in future release ***/
    2
    3/*	The MIT License (MIT)
    4 *
    5 *	Copyright (c) 2019-2025 Rick Workman
    6 *
    7 *	Permission is hereby granted, free of charge, to any person obtaining a copy
    8 *	of this software and associated documentation files (the "Software"), to deal
    9 *	in the Software without restriction, including without limitation the rights
   10 *	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
   11 *	copies of the Software, and to permit persons to whom the Software is
   12 *	furnished to do so, subject to the following conditions:
   13 *
   14 *	The above copyright notice and this permission notice shall be included in all
   15 *	copies or substantial portions of the Software.
   16 *
   17 *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   18 *	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   19 *	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   20 *	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   21 *	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
   22 *	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   23 *	SOFTWARE.
   24 */
   25%
   26% Expression variables are Prolog var's; numeric constants are precise.
   27% Note floats are treated like variables and no arithmetic is performed on them
   28%
   29
   30% signed multiplier and power values
   31sign_val_(+,C,C) :- C>=0.
   32sign_val_(-,C,N) :- N is -C.
   33
   34pwr_val_(*,C,C) :- C>=0.
   35pwr_val_(/,C,N) :- N is -C.
   36
   37
   38% utility to distribute A*(B1+B2) if they have variables in common or A is a constant
   39% former is safe for intervals (sub-distributive law), latter is valid
   40distribute_(C,B,Exp) :-
   41	number(C), !,
   42	B=..[AddOp,B1,B2], (AddOp == '+' ; AddOp == '-' ),  % watch out for vars
   43	DExp=..[AddOp,C*B1,C*B2],
   44	simplify(DExp,Exp).
   45
   46% utility for (in)equality reduction 
   47/*simplify_eq_(Diff,B,Op,Exp) :-  compound(Diff), Diff = A1-A2, % LHS is difference and RHS =:=0
   48	number(B), B =:= 0, !,
   49	simplify_eq_(A1,A2,Op,Exp). */
   50simplify_eq_(A,B,Op,Exp) :-         % LHS and RHS are expressions with shared vars
   51	num_exp_vars_(A,AVs),         % only consider vars in arithmetic Ops
   52	num_exp_vars_(B,BVs),
   53	shared_vars_in_(AVs,BVs), !,  % so far just a test, no unification of non-local vars
   54	simplify(A-B,ES),             % normalize with non-zero value on RHS
   55	num_separate_(ES,LHS,RHS),
   56	(RHS =:= 0, LHS = EA-EB -> Exp =.. [Op,EA,EB] ; Exp=..[Op,LHS,RHS]).
   57simplify_eq_(A,B,Op,Exp) :-         % everything else, leave LHS and RHS alone
   58	simplify(A,AS),
   59	num_separate_(AS,LHS,AN),
   60	simplify(B,BS),
   61	num_separate_(BS,BS1,BN),
   62	num_rhs_(AN,BN,BS1,RHS),      % leave as symbolic, for later fuzzing as req'd
   63	Exp=..[Op,LHS,RHS].
   64
   65num_separate_(Exp,Exp,0) :- var(Exp), !.
   66num_separate_(A+B,Out,Num) :- !,
   67	(number(B) 
   68	 -> Num is -B,
   69	 	Out=A
   70	 ;	num_separate_(A,AOut,Num),
   71	 	Out=AOut+B
   72	).
   73num_separate_(A-B,Out,Num) :- !,
   74	(number(B) 
   75	 -> Num = B,
   76	 	Out=A
   77	 ;	num_separate_(A,AOut,Num),
   78	 	Out=AOut-B
   79	).
   80num_separate_(Exp,Exp,0).
   81
   82num_rhs_(N,N,E,E)   :- !.
   83num_rhs_(0,N,E,RHS) :- !, (E==0 -> RHS = -N ; RHS = E-N).
   84num_rhs_(N,0,E,RHS) :- !, (E==0 -> RHS = N ; RHS = E+N).
   85num_rhs_(L,R,E,RHS) :- N is L-R, (E==0 -> RHS = N ; RHS = E+N).
   86
   87num_exp_vars_(Exp,Vars) :-  % collect vars in arithmetic expressions, fails if none
   88	exp_vars_(Exp,V/V,Vs/[]),	
   89	(Vs=[_|_] -> term_variables(Vs,Vars)).  % removes any duplicates, fails if no vars 
   90
   91exp_vars_(X,ListIn,ListIn) :- atomic(X), !.  % includes []
   92exp_vars_(X,List/[X|NextTail],List/NextTail) :- var(X), !. 
   93exp_vars_([X|Xs],ListIn,ListOut) :- !,
   94	exp_vars_(X,ListIn,ListNxt),
   95	exp_vars_(Xs,ListNxt,ListOut).
   96exp_vars_(Exp,ListIn,ListOut) :-
   97	Exp =.. [Op|Xs],
   98	(memberchk(Op,[-,+,*,/,**]) -> exp_vars_(Xs,ListIn,ListOut) ; ListIn = ListOut).
   99
  100shared_vars_in_([AV|AVs],BVs) :-
  101	shared_var_(BVs,AV) -> true ; shared_vars_in_(AVs,BVs).
  102
  103shared_var_([BV|BVs],AV) :-
  104	AV == BV -> true ; shared_var_(BVs,AV).
  105
  106%
  107% simplify/2
  108%
  109simplify(A,A) :- var(A),!.
  110
  111simplify(N,N) :- number(N),!.
  112
  113% possible distribute
  114simplify(A*B,Exp) :-
  115	compound(B),
  116	distribute_(A,B,Exp), !.
  117simplify(A*B,Exp) :-
  118	compound(A),
  119	distribute_(B,A,Exp), !.
  120
  121% simplify equalities and inequalities
  122simplify(A==B,Exp) :-
  123	simplify_eq_(A,B,==,Exp), !.
  124
  125simplify(A=<B,Exp) :-
  126	simplify_eq_(A,B,=<,Exp), !.
  127
  128simplify(A>=B,Exp) :-
  129	simplify_eq_(A,B,>=,Exp), !.
  130
  131simplify(A<B,Exp) :-
  132	simplify_eq_(A,B,<,Exp), !.
  133
  134simplify(A>B,Exp) :-
  135	simplify_eq_(A,B,>,Exp), !.
  136
  137simplify(A is B,A is B) :- !.  % Don't bother simplifying `is`
  138/*
  139% simplify "cascaded" divisions A/B/C = (A/B)/C = A*C/B
  140simplify(A/B,Exp) :- 
  141%	compound(B), B=Bn/Bd, !,
  142%	simplify(A*Bd/Bn,Exp).
  143	compound(A), A=An/Ad, !,
  144	simplify(An*B, E1),
  145	simplify(E1/Ad,Exp).
  146*/
  147%
  148% General purpose simplify
  149%
  150simplify(E,Exp) :-
  151	collect_exp_(E, L/L, Es/[]),  % collect terms in a difference list
  152	collect_exp_items(Es,Is),     % collect like terms
  153	reduce_exp_items_(Is,Exp),    % and reduce back to expression
  154	!.
  155	
  156% use difference list to collect terms of an expression
  157collect_exp_(A, List, Head/NewT) :-
  158	var(A), !,
  159	List=Head/[term(A,1)|NewT].  % separate to keep head unification cheap
  160	
  161collect_exp_(C, List, Head/NewT) :-
  162	preciseNumber(C), !,               % precise numeric value
  163	List=Head/[C|NewT].
  164
  165collect_exp_(A, List, Head/NewT) :-    % floats as terms, can collect but not do arithmetic
  166	float(A), !,
  167	List=Head/[term(A,1)|NewT].
  168
  169collect_exp_(A+B, List, NewList) :-
  170	!,
  171	left_exp_(A,AE),
  172	collect_exp_(AE,List,ListA),
  173	simplify(B,BT), 
  174	collect_exp_(BT,ListA,NewList).
  175
  176collect_exp_(A-B,List,NewList) :-
  177	!,
  178	left_exp_(A,AE),
  179	collect_exp_(AE,List,ListA),
  180	simplify(B,BT),
  181	collect_neg_(BT,ListA,NewList).
  182	
  183collect_exp_(-B,List,NewList) :-
  184	simplify(B,BT),
  185	collect_neg_(BT,List,NewList).
  186
  187collect_exp_(T, List/[Term|NewT], List/NewT) :-
  188	simplify_term(T,Term).
  189
  190left_exp_(A,A)  :- var(A), !.
  191left_exp_(A,A)  :- functor(A,+,2), !.
  192left_exp_(A,A)  :- functor(A,-,2), !.
  193left_exp_(A,AE) :- simplify(A,AE).
  194
  195collect_neg_(BT,List/[N|NewT], List/NewT) :-
  196	rational(BT),
  197	N is -BT.  % rational constant expressed as AT-BT
  198collect_neg_(BT,ListA/T,ListA/NewT) :-
  199	collect_exp_(BT,P/P,ListB),
  200	negate_term_(ListB,T/NewT).
  201
  202negate_term_(T/T,NewT/NewT) :- var(T).
  203negate_term_([term(V,N)|Es]/T,[term(V,NN)|NEs]/NewT) :-
  204	NN is -N,
  205	negate_term_(Es/T,NEs/NewT).
  206negate_term_([N|Es]/T,[NN|NEs]/NewT) :-
  207	NN is -N,
  208	negate_term_(Es/T,NEs/NewT).
  209
  210collect_exp_items(Es,NEs) :-
  211	sort(0,@>=,Es,Sorted),  % num_separate_/3 depends on sort order
  212	collect_exp_items_(Sorted,NEs).
  213
  214% assume items sorted
  215collect_exp_items_([],[]).
  216collect_exp_items_([U],[U]) :- !.
  217collect_exp_items_([U,V|Es],EsNxt) :-
  218	number(U), number(V), !,                % add precise constant values
  219	S is U+V, 
  220	collect_exp_items_([S|Es],EsNxt).
  221collect_exp_items_([term(V,N1),term(U,N2)|Es],EsNxt) :-
  222	V==U,
  223	(ground(V), \+ string(V) -> V =\= 1.0Inf ; true), !,  % infinities don't obey algebraic rules
  224	N is N1+N2,
  225	collect_exp_items_([term(V,N)|Es],EsNxt).
  226collect_exp_items_([U,V|Es],[U|EsNxt]) :-   % Note that imprecise floats are not added
  227	collect_exp_items_([V|Es],EsNxt).
  228
  229reduce_exp_items_([T1,T2|Ts],Exp) :-        % 2 or more terms, combine and continue
  230	reduce_exp_item_(T1,Op1,Exp1),
  231	reduce_exp_item_(T2,Op2,Exp2),
  232	build_exp_(Exp1,Exp2,Op1,Op2,ExpN), !,  % ExpN =.. [Op,Exp1,Exp2], !,
  233	reduce_exp_items_([ExpN|Ts],Exp).
  234reduce_exp_items_([E],Exp) :-               % terminating case, 1 item left
  235	(compound(E), functor(E,term,2)
  236	 -> reduce_exp_items_([0,E],Exp)        % terminate single term(_..)
  237	 ;  Exp = E                             % already reduced to var, constant or expression
  238	).
  239
  240build_exp_(Z,Exp2,_Op1,Op2,NExp2) :- number(Z), Z =:= 0,
  241	(Op2 == '-' -> build_Nexp_(Exp2,NExp2) ;  NExp2 = Exp2).
  242build_exp_(Exp1,Z,Op1,_Op2,NExp1) :- number(Z), Z =:= 0,  
  243	(Op1 == '-' -> build_Nexp_(Exp1,NExp1) ; NExp1 = Exp1).
  244build_exp_(Exp1,Exp2,+,+,Exp1+Exp2).
  245build_exp_(Exp1,Exp2,+,-,Exp1-Exp2).
  246build_exp_(Exp1,Exp2,-,+,Exp2-Exp1).
  247build_exp_(Exp1,Exp2,-,-,NExp1-Exp2) :- build_Nexp_(Exp1,NExp1).
  248
  249build_Nexp_(Exp, NExp) :-                                   % constant
  250	rational(Exp),
  251	NExp is -Exp.
  252build_Nexp_(Exp, NExp) :-                                   % * or / expression, find head
  253	nonvar(Exp), Exp =.. [Op,A1,A2], (Op == * ; Op == /),
  254	build_Nexp_(A1,NA1),
  255	NExp =.. [Op,NA1,A2].
  256build_Nexp_(Exp, -Exp).                                     % other expression or var
  257
  258reduce_exp_item_(V,           +, V)    :- var(V).           % variables
  259reduce_exp_item_(R,           +, R)    :- rational(R).      % constants
  260reduce_exp_item_(term(V,1),   +, V).
  261reduce_exp_item_(term(V,-1),  -, V).
  262reduce_exp_item_(term(V,0),   +, 0)    :- finite_(V).       % reduces to +0 if V finite
  263reduce_exp_item_(term(V,N),   Op, T)   :- mult_term_(V, N, Op, T).
  264reduce_exp_item_(-Exp,        -, Exp).
  265reduce_exp_item_(Exp,         +, Exp).
  266
  267finite_(V) :- 
  268	(ground(V)
  269%	 -> abs(V) < 1.0Inf  % ground expression, evaluates to finite value
  270	 -> catch(abs(V) =\= 1.0Inf,_Err,true)  % ground expression, does not evaluate to infinite value
  271	 ;  interval_object(V, _Type, (LB,UB), _),  % interval with finite bounds
  272	    -1.0Inf < LB, UB < 1.0Inf
  273	).
  274
  275mult_term_(T,   N, Op, M*T)   :- var(T), val_sign_(N,Op,M).
  276mult_term_(X*Y, N, Op, Exp*Y) :- mult_term_(X,N,Op,Exp).     % maintain Op associativity
  277mult_term_(X/Y, N, Op, Exp/Y) :- mult_term_(X,N,Op,Exp).
  278mult_term_(T,   N, Op, M*T)   :- val_sign_(N,Op,M).
  279
  280val_sign_(V,Op,Vp) :- Vp is abs(V), (V<0 -> Op = (-) ; Op = (+)).
  281
  282%
  283% simplify a term to an "expression" of term structures and constants
  284%
  285simplify_term(T,Term) :-
  286	collect_term_(T, L/L, Ts/[]),    % collect in a difference list
  287	collect_term_items_(Ts,Is),      % collect like terms
  288	% this does constant folding if reduction resulted in a constant
  289	reduce_term_items_(Is,ITerm),    % reduce back to expression
  290	term_final_(ITerm,Term),
  291	!.
  292
  293term_final_(Term,Term) :- rational(Term), !.
  294term_final_(Term,term(T,C)) :- compound(Term), Term = C*T, rational(C), !.
  295term_final_(Term,term(Term,1)).
  296
  297collect_term_(A, List, Head/NewT) :-
  298	var(A), !,
  299	List=Head/[elem(A,1)|NewT].
  300	
  301collect_term_(A, List, Head/NewT) :-
  302	rational(A), !,
  303	List=Head/[A|NewT].
  304
  305collect_term_(A, List, Head/NewT) :-
  306	float(A), !,
  307	List=Head/[elem(A,1)|NewT].
  308
  309collect_term_(-A, List, ListA/NewT) :- % -Term as Term * -1 for reducing signs
  310	simplify(A,AT), collect_term_(AT,List,ListA/[-1|NewT]).
  311
  312collect_term_(A**N, List, Head/NewT) :-  % possible special case of user written element
  313	simplify(N,NT), rational(NT),
  314	simplify(A,AT),
  315	simplify_pwr_(AT,NT,Term),
  316	List=Head/[Term|NewT].
  317
  318collect_term_(A*B,List,NewList) :-
  319	!,
  320	left_term_(A,AE),
  321	collect_term_(AE,List,ListA),
  322	simplify(B,BT),
  323	collect_term_(BT,ListA,NewList).
  324	
  325collect_term_(A/B,List,ListA/NewT) :-
  326	!,
  327	left_term_(A,AE),
  328	simplify(B,BT),
  329	collect_term_(AE,List,ListA/T),
  330	collect_term_(BT,P/P,ListB),
  331	invert_term_(ListB,T/NewT).
  332
  333collect_term_(E,List/[elem(Exp,1)|NewT], List/NewT) :-
  334	E =.. [F|IArgs],
  335	simplify_list_(IArgs,OArgs),
  336	Exp =.. [F|OArgs].
  337
  338% simplify_pwr_: NT rational, 
  339simplify_pwr_(AT,NT,elem(Exp,Pwr)) :-
  340	compound(AT), AT=Exp**P, rational(P), !,  % (X**P)**NT, NT and P both rational
  341	Pwr is NT*P.
  342simplify_pwr_(AT,NT,elem(AT,NT)).
  343
  344left_term_(A,A)  :- var(A), !.
  345left_term_(A,A)  :- functor(A,*,2), !.
  346left_term_(A,A)  :- functor(A,/,2), !.
  347left_term_(A,AE) :- simplify(A,AE).
  348
  349simplify_list_([],[]).
  350simplify_list_([I|IArgs],[O|OArgs]) :-
  351	simplify(I,O),
  352	simplify_list_(IArgs,OArgs).
  353
  354invert_term_(T/T,NewT/NewT) :- var(T),!.
  355invert_term_([elem(V,N)|Es]/T,[elem(V,NN)|NEs]/NewT) :-  !,
  356	NN is -N,
  357	invert_term_(Es/T,NEs/NewT).
  358invert_term_([N|Es]/T,[NN|NEs]/NewT) :-
  359	NN is 1/N,
  360	invert_term_(Es/T,NEs/NewT).
  361	
  362collect_term_items_([],[]).
  363collect_term_items_(Es,NEs) :-
  364	msort([1|Es],Sorted),                  % ensure initial constant multiplier and sort
  365	collect_term_item_(Sorted,NEs).
  366
  367collect_term_item_([],[]).
  368collect_term_item_([U],[U]) :- !.
  369collect_term_item_([U,V|Es],EsNxt) :-
  370	rational(U), rational(V), !,           % precise for multiply
  371	S is U*V,
  372	collect_term_item_([S|Es],EsNxt).
  373collect_term_item_([elem(U,N1),elem(V,N2)|Es],EsNxt) :-
  374	V==U,
  375	(ground(V), \+string(V) -> V =\= 1.0Inf ; true), !,  % infinities don't obey algebraic rules
  376	N is N1+N2,
  377	collect_term_item_([elem(U,N)|Es],EsNxt).
  378collect_term_item_([U,V|Es],[U|EsNxt]) :-
  379	collect_term_item_([V|Es],EsNxt).
  380
  381reduce_term_items_([Exp],Exp) :-           % terminating variable
  382	var(Exp), !.
  383reduce_term_items_([Exp],Exp) :-           % terminating constant
  384	rational(Exp), !.
  385reduce_term_items_([elem(V,N)],Exp) :- !,  % terminating single elem(_..)
  386	reduce_term_items_([1,elem(V,N)],Exp).
  387reduce_term_items_([Exp],Exp).             % terminating final expression
  388reduce_term_items_([T1,T2|Ts],Exp) :-
  389	reduce_term_item_(T1,Exp1,_),
  390	reduce_term_item_(T2,Exp2,Op),
  391	build_term_(Exp1,Exp2,Op,ExpN),        % ExpN =.. [Op,Exp1,Exp2],
  392	!,
  393	reduce_term_items_([ExpN|Ts],Exp).
  394
  395% optimize some cases involving ±1
  396build_term_(C, E, /, Exp) :- nonvar(E), E=A**B, one_term_(C,A**NB,Exp), rational(B), NB is -B.
  397build_term_(E, C, _, Exp) :- nonvar(E), E=C1/B, one_term_(C1,C/B,Exp).
  398build_term_(C, E, *, Exp) :- one_term_(C, E, Exp).
  399build_term_(E, C, _, Exp) :- one_term_(C, E, Exp). 
  400build_term_(E1,E2,Op,Exp) :- Exp =.. [Op,E1,E2].
  401
  402one_term_(One, Exp0, Exp) :-
  403	number(One), C is float(One),
  404	( C =  1.0 -> Exp =  Exp0
  405	;(C = -1.0 -> Exp = -Exp0)
  406	).
  407
  408reduce_term_item_(V,          V, *)     :- var(V),!.  % already reduced to var
  409reduce_term_item_(elem(_, 0), 1, *).
  410reduce_term_item_(elem(V, 1), V, *).
  411reduce_term_item_(elem(V,-1), V, /).
  412reduce_term_item_(elem(V, E), V**P, Op) :- P is abs(E), (E>0 -> Op= * ; Op= /).
  413reduce_term_item_(R,          R, *).  % catchall: constant or expression