1/*	The MIT License (MIT)
    2 *
    3 *	Copyright (c) 2019-2025 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 
   45/*simplify_eq_(Diff,B,Op,Exp) :-  compound(Diff), Diff = A1-A2, % LHS is difference and RHS =:=0
   46	number(B), B =:= 0, !,
   47	simplify_eq_(A1,A2,Op,Exp). */
   48simplify_eq_(A,B,Op,Exp) :-         % LHS and RHS are expressions with shared vars
   49	num_exp_vars_(A,AVs),         % only consider vars in arithmetic Ops
   50	num_exp_vars_(B,BVs),
   51	shared_vars_in_(AVs,BVs), !,  % so far just a test, no unification of non-local vars
   52	simplify(A-B,ES),             % normalize with non-zero value on RHS
   53	num_separate_(ES,LHS,RHS),
   54	(RHS =:= 0, LHS = EA-EB -> Exp =.. [Op,EA,EB] ; Exp=..[Op,LHS,RHS]).
   55simplify_eq_(A,B,Op,Exp) :-         % everything else, leave LHS and RHS alone
   56	simplify(A,AS),
   57	num_separate_(AS,LHS,AN),
   58	simplify(B,BS),
   59	num_separate_(BS,BS1,BN),
   60	num_rhs_(AN,BN,BS1,RHS),      % leave as symbolic, for later fuzzing as req'd
   61	Exp=..[Op,LHS,RHS].
   62
   63num_separate_(Exp,Exp,0) :- var(Exp), !.
   64num_separate_(A+B,Out,Num) :- !,
   65	(number(B) 
   66	 -> Num is -B,
   67	 	Out=A
   68	 ;	num_separate_(A,AOut,Num),
   69	 	Out=AOut+B
   70	).
   71num_separate_(A-B,Out,Num) :- !,
   72	(number(B) 
   73	 -> Num = B,
   74	 	Out=A
   75	 ;	num_separate_(A,AOut,Num),
   76	 	Out=AOut-B
   77	).
   78num_separate_(Exp,Exp,0).
   79
   80num_rhs_(N,N,E,E)   :- !.
   81num_rhs_(0,N,E,RHS) :- !, (E==0 -> RHS = -N ; RHS = E-N).
   82num_rhs_(N,0,E,RHS) :- !, (E==0 -> RHS = N ; RHS = E+N).
   83num_rhs_(L,R,E,RHS) :- N is L-R, (E==0 -> RHS = N ; RHS = E+N).
   84
   85num_exp_vars_(Exp,Vars) :-  % collect vars in arithmetic expressions, fails if none
   86	exp_vars_(Exp,V/V,Vs/[]),	
   87	(Vs=[_|_] -> term_variables(Vs,Vars)).  % removes any duplicates, fails if no vars 
   88
   89exp_vars_(X,ListIn,ListIn) :- atomic(X), !.  % includes []
   90exp_vars_(X,List/[X|NextTail],List/NextTail) :- var(X), !. 
   91exp_vars_([X|Xs],ListIn,ListOut) :- !,
   92	exp_vars_(X,ListIn,ListNxt),
   93	exp_vars_(Xs,ListNxt,ListOut).
   94exp_vars_(Exp,ListIn,ListOut) :-
   95	Exp =.. [Op|Xs],
   96	(memberchk(Op,[-,+,*,/,**]) -> exp_vars_(Xs,ListIn,ListOut) ; ListIn = ListOut).
   97
   98shared_vars_in_([AV|AVs],BVs) :-
   99	shared_var_(BVs,AV) -> true ; shared_vars_in_(AVs,BVs).
  100
  101shared_var_([BV|BVs],AV) :-
  102	AV == BV -> true ; shared_var_(BVs,AV).
  103
  104%
  105% simplify/2
  106%
  107simplify(A,A) :- var(A),!.
  108
  109simplify(N,N) :- number(N),!.
  110
  111% possible distribute
  112simplify(A*B,Exp) :-
  113	compound(B),
  114	distribute_(A,B,Exp), !.
  115simplify(A*B,Exp) :-
  116	compound(A),
  117	distribute_(B,A,Exp), !.
  118
  119% simplify equalities and inequalities
  120simplify(A==B,Exp) :-
  121	simplify_eq_(A,B,==,Exp), !.
  122
  123simplify(A=<B,Exp) :-
  124	simplify_eq_(A,B,=<,Exp), !.
  125
  126simplify(A>=B,Exp) :-
  127	simplify_eq_(A,B,>=,Exp), !.
  128
  129simplify(A<B,Exp) :-
  130	simplify_eq_(A,B,<,Exp), !.
  131
  132simplify(A>B,Exp) :-
  133	simplify_eq_(A,B,>,Exp), !.
  134
  135simplify(A is B,A is B) :- !.  % Don't bother simplifying `is`
  136/*
  137% simplify "cascaded" divisions A/B/C = (A/B)/C = A*C/B
  138simplify(A/B,Exp) :- 
  139%	compound(B), B=Bn/Bd, !,
  140%	simplify(A*Bd/Bn,Exp).
  141	compound(A), A=An/Ad, !,
  142	simplify(An*B, E1),
  143	simplify(E1/Ad,Exp).
  144*/
  145%
  146% General purpose simplify
  147%
  148simplify(E,Exp) :-
  149	collect_exp_(E, L/L, Es/[]),  % collect terms in a difference list
  150	collect_exp_items(Es,Is),     % collect like terms
  151	reduce_exp_items_(Is,Exp),    % and reduce back to expression
  152	!.
  153	
  154% use difference list to collect terms of an expression
  155collect_exp_(A, List, Head/NewT) :-
  156	var(A), !,
  157	List=Head/[term(A,1)|NewT].  % separate to keep head unification cheap
  158	
  159collect_exp_(C, List, Head/NewT) :-
  160	preciseNumber(C), !,               % precise numeric value
  161	List=Head/[C|NewT].
  162
  163collect_exp_(A, List, Head/NewT) :-    % floats as terms, can collect but not do arithmetic
  164	float(A), !,
  165	List=Head/[term(A,1)|NewT].
  166
  167collect_exp_(A+B, List, NewList) :-
  168	!,
  169	left_exp_(A,AE),
  170	collect_exp_(AE,List,ListA),
  171	simplify(B,BT), 
  172	collect_exp_(BT,ListA,NewList).
  173
  174collect_exp_(A-B,List,NewList) :-
  175	!,
  176	left_exp_(A,AE),
  177	collect_exp_(AE,List,ListA),
  178	simplify(B,BT),
  179	collect_neg_(BT,ListA,NewList).
  180	
  181collect_exp_(-B,List,NewList) :-
  182	simplify(B,BT),
  183	collect_neg_(BT,List,NewList).
  184
  185collect_exp_(T, List/[Term|NewT], List/NewT) :-
  186	simplify_term(T,Term).
  187
  188left_exp_(A,A)  :- var(A), !.
  189left_exp_(A,A)  :- functor(A,+,2), !.
  190left_exp_(A,A)  :- functor(A,-,2), !.
  191left_exp_(A,AE) :- simplify(A,AE).
  192
  193collect_neg_(BT,List/[N|NewT], List/NewT) :-
  194	rational(BT),
  195	N is -BT.  % rational constant expressed as AT-BT
  196collect_neg_(BT,ListA/T,ListA/NewT) :-
  197	collect_exp_(BT,P/P,ListB),
  198	negate_term_(ListB,T/NewT).
  199
  200negate_term_(T/T,NewT/NewT) :- var(T).
  201negate_term_([term(V,N)|Es]/T,[term(V,NN)|NEs]/NewT) :-
  202	NN is -N,
  203	negate_term_(Es/T,NEs/NewT).
  204negate_term_([N|Es]/T,[NN|NEs]/NewT) :-
  205	NN is -N,
  206	negate_term_(Es/T,NEs/NewT).
  207
  208collect_exp_items(Es,NEs) :-
  209	sort(0,@>=,Es,Sorted),  % num_separate_/3 depends on sort order
  210	collect_exp_items_(Sorted,NEs).
  211
  212% assume items sorted
  213collect_exp_items_([],[]).
  214collect_exp_items_([U],[U]) :- !.
  215collect_exp_items_([U,V|Es],EsNxt) :-
  216	number(U), number(V), !,                % add precise constant values
  217	S is U+V, 
  218	collect_exp_items_([S|Es],EsNxt).
  219collect_exp_items_([term(V,N1),term(U,N2)|Es],EsNxt) :-
  220	V==U,
  221	(ground(V) -> V =\= 1.0Inf ; true), !,  % infinities don't obey algebraic rules
  222	N is N1+N2,
  223	collect_exp_items_([term(V,N)|Es],EsNxt).
  224collect_exp_items_([U,V|Es],[U|EsNxt]) :-   % Note that imprecise floats are not added
  225	collect_exp_items_([V|Es],EsNxt).
  226
  227reduce_exp_items_([T1,T2|Ts],Exp) :-        % 2 or more terms, combine and continue
  228	reduce_exp_item_(T1,Op1,Exp1),
  229	reduce_exp_item_(T2,Op2,Exp2),
  230	build_exp_(Exp1,Exp2,Op1,Op2,ExpN), !,  % ExpN =.. [Op,Exp1,Exp2], !,
  231	reduce_exp_items_([ExpN|Ts],Exp).
  232reduce_exp_items_([E],Exp) :-               % terminating case, 1 item left
  233	(compound(E), functor(E,term,2)
  234	 -> reduce_exp_items_([0,E],Exp)        % terminate single term(_..)
  235	 ;  Exp = E                             % already reduced to var, constant or expression
  236	).
  237
  238build_exp_(Z,Exp2,_Op1,Op2,NExp2) :- number(Z), Z =:= 0,
  239	(Op2 == '-' -> build_Nexp_(Exp2,NExp2) ;  NExp2 = Exp2).
  240build_exp_(Exp1,Z,Op1,_Op2,NExp1) :- number(Z), Z =:= 0,  
  241	(Op1 == '-' -> build_Nexp_(Exp1,NExp1) ; NExp1 = Exp1).
  242build_exp_(Exp1,Exp2,+,+,Exp1+Exp2).
  243build_exp_(Exp1,Exp2,+,-,Exp1-Exp2).
  244build_exp_(Exp1,Exp2,-,+,Exp2-Exp1).
  245build_exp_(Exp1,Exp2,-,-,NExp1-Exp2) :- build_Nexp_(Exp1,NExp1).
  246
  247build_Nexp_(Exp, NExp) :-                                   % constant
  248	rational(Exp),
  249	NExp is -Exp.
  250build_Nexp_(Exp, NExp) :-                                   % * or / expression, find head
  251	nonvar(Exp), Exp =.. [Op,A1,A2], (Op == * ; Op == /),
  252	build_Nexp_(A1,NA1),
  253	NExp =.. [Op,NA1,A2].
  254build_Nexp_(Exp, -Exp).                                     % other expression or var
  255
  256reduce_exp_item_(V,           +, V)    :- var(V).           % variables
  257reduce_exp_item_(R,           +, R)    :- rational(R).      % constants
  258reduce_exp_item_(term(V,1),   +, V).
  259reduce_exp_item_(term(V,-1),  -, V).
  260reduce_exp_item_(term(V,0),   +, 0)    :- finite_(V).       % reduces to +0 if V finite
  261reduce_exp_item_(term(V,N),   Op, T)   :- mult_term_(V, N, Op, T).
  262reduce_exp_item_(-Exp,        -, Exp).
  263reduce_exp_item_(Exp,         +, Exp).
  264
  265finite_(V) :- 
  266	(ground(V)
  267%	 -> abs(V) < 1.0Inf  % ground expression, evaluates to finite value
  268	 -> catch(abs(V) =\= 1.0Inf,_Err,true)  % ground expression, does not evaluate to infinite value
  269	 ;  interval_object(V, _Type, (LB,UB), _),  % interval with finite bounds
  270	    -1.0Inf < LB, UB < 1.0Inf
  271	).
  272
  273mult_term_(T,   N, Op, M*T)   :- var(T), val_sign_(N,Op,M).
  274mult_term_(X*Y, N, Op, Exp*Y) :- mult_term_(X,N,Op,Exp).     % maintain Op associativity
  275mult_term_(X/Y, N, Op, Exp/Y) :- mult_term_(X,N,Op,Exp).
  276mult_term_(T,   N, Op, M*T)   :- val_sign_(N,Op,M).
  277
  278val_sign_(V,Op,Vp) :- Vp is abs(V), (V<0 -> Op = (-) ; Op = (+)).
  279
  280%
  281% simplify a term to an "expression" of term structures and constants
  282%
  283simplify_term(T,Term) :-
  284	collect_term_(T, L/L, Ts/[]),    % collect in a difference list
  285	collect_term_items_(Ts,Is),      % collect like terms
  286	% this does constant folding if reduction resulted in a constant
  287	reduce_term_items_(Is,ITerm),    % reduce back to expression
  288	term_final_(ITerm,Term),
  289	!.
  290
  291term_final_(Term,Term) :- rational(Term), !.
  292term_final_(Term,term(T,C)) :- compound(Term), Term = C*T, rational(C), !.
  293term_final_(Term,term(Term,1)).
  294
  295collect_term_(A, List, Head/NewT) :-
  296	var(A), !,
  297	List=Head/[elem(A,1)|NewT].
  298	
  299collect_term_(A, List, Head/NewT) :-
  300	rational(A), !,
  301	List=Head/[A|NewT].
  302
  303collect_term_(A, List, Head/NewT) :-
  304	float(A), !,
  305	List=Head/[elem(A,1)|NewT].
  306
  307collect_term_(-A, List, ListA/NewT) :- % -Term as Term * -1 for reducing signs
  308	simplify(A,AT), collect_term_(AT,List,ListA/[-1|NewT]).
  309
  310collect_term_(A**N, List, Head/NewT) :-  % possible special case of user written element
  311	simplify(N,NT), rational(NT),
  312	simplify(A,AT),
  313	simplify_pwr_(AT,NT,Term),
  314	List=Head/[Term|NewT].
  315
  316collect_term_(A*B,List,NewList) :-
  317	!,
  318	left_term_(A,AE),
  319	collect_term_(AE,List,ListA),
  320	simplify(B,BT),
  321	collect_term_(BT,ListA,NewList).
  322	
  323collect_term_(A/B,List,ListA/NewT) :-
  324	!,
  325	left_term_(A,AE),
  326	simplify(B,BT),
  327	collect_term_(AE,List,ListA/T),
  328	collect_term_(BT,P/P,ListB),
  329	invert_term_(ListB,T/NewT).
  330
  331collect_term_(E,List/[elem(Exp,1)|NewT], List/NewT) :-
  332	E =.. [F|IArgs],
  333	simplify_list_(IArgs,OArgs),
  334	Exp =.. [F|OArgs].
  335
  336% simplify_pwr_: NT rational, 
  337simplify_pwr_(AT,NT,elem(Exp,Pwr)) :-
  338	compound(AT), AT=Exp**P, rational(P), !,  % (X**P)**NT, NT and P both rational
  339	Pwr is NT*P.
  340simplify_pwr_(AT,NT,elem(AT,NT)).
  341
  342left_term_(A,A)  :- var(A), !.
  343left_term_(A,A)  :- functor(A,*,2), !.
  344left_term_(A,A)  :- functor(A,/,2), !.
  345left_term_(A,AE) :- simplify(A,AE).
  346
  347simplify_list_([],[]).
  348simplify_list_([I|IArgs],[O|OArgs]) :-
  349	simplify(I,O),
  350	simplify_list_(IArgs,OArgs).
  351
  352invert_term_(T/T,NewT/NewT) :- var(T),!.
  353invert_term_([elem(V,N)|Es]/T,[elem(V,NN)|NEs]/NewT) :-  !,
  354	NN is -N,
  355	invert_term_(Es/T,NEs/NewT).
  356invert_term_([N|Es]/T,[NN|NEs]/NewT) :-
  357	NN is 1/N,
  358	invert_term_(Es/T,NEs/NewT).
  359	
  360collect_term_items_([],[]).
  361collect_term_items_(Es,NEs) :-
  362	msort([1|Es],Sorted),                  % ensure initial constant multiplier and sort
  363	collect_term_item_(Sorted,NEs).
  364
  365collect_term_item_([],[]).
  366collect_term_item_([U],[U]) :- !.
  367collect_term_item_([U,V|Es],EsNxt) :-
  368	rational(U), rational(V), !,           % precise for multiply
  369	S is U*V,
  370	collect_term_item_([S|Es],EsNxt).
  371collect_term_item_([elem(U,N1),elem(V,N2)|Es],EsNxt) :-
  372	V==U,
  373	(ground(V) -> V =\= 1.0Inf ; true), !, % infinities don't obey algebraic rules
  374	N is N1+N2,
  375	collect_term_item_([elem(U,N)|Es],EsNxt).
  376collect_term_item_([U,V|Es],[U|EsNxt]) :-
  377	collect_term_item_([V|Es],EsNxt).
  378
  379reduce_term_items_([Exp],Exp) :-           % terminating variable
  380	var(Exp), !.
  381reduce_term_items_([Exp],Exp) :-           % terminating constant
  382	rational(Exp), !.
  383reduce_term_items_([elem(V,N)],Exp) :- !,  % terminating single elem(_..)
  384	reduce_term_items_([1,elem(V,N)],Exp).
  385reduce_term_items_([Exp],Exp).             % terminating final expression
  386reduce_term_items_([T1,T2|Ts],Exp) :-
  387	reduce_term_item_(T1,Exp1,_),
  388	reduce_term_item_(T2,Exp2,Op),
  389	build_term_(Exp1,Exp2,Op,ExpN),        % ExpN =.. [Op,Exp1,Exp2],
  390	!,
  391	reduce_term_items_([ExpN|Ts],Exp).
  392
  393% optimize some cases involving ±1
  394build_term_(C, E, /, Exp) :- nonvar(E), E=A**B, one_term_(C,A**NB,Exp), rational(B), NB is -B.
  395build_term_(E, C, _, Exp) :- nonvar(E), E=C1/B, one_term_(C1,C/B,Exp).
  396build_term_(C, E, *, Exp) :- one_term_(C, E, Exp).
  397build_term_(E, C, _, Exp) :- one_term_(C, E, Exp). 
  398build_term_(E1,E2,Op,Exp) :- Exp =.. [Op,E1,E2].
  399
  400one_term_(One, Exp0, Exp) :-
  401	number(One), C is float(One),
  402	( C =  1.0 -> Exp =  Exp0
  403	;(C = -1.0 -> Exp = -Exp0)
  404	).
  405
  406reduce_term_item_(V,          V, *)     :- var(V),!.  % already reduced to var
  407reduce_term_item_(elem(_, 0), 1, *).
  408reduce_term_item_(elem(V, 1), V, *).
  409reduce_term_item_(elem(V,-1), V, /).
  410reduce_term_item_(elem(V, E), V**P, Op) :- P is abs(E), (E>0 -> Op= * ; Op= /).
  411reduce_term_item_(R,          R, *).  % catchall: constant or expression