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