1/*
    2
    3    Part of CLP(R) (Constraint Logic Programming over Reals)
    4
    5    Author:        Leslie De Koninck
    6    E-mail:        Leslie.DeKoninck@cs.kuleuven.be
    7    WWW:           http://www.swi-prolog.org
    8		   http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09
    9    Copyright (C): 2004, K.U. Leuven and
   10		   1992-1995, Austrian Research Institute for
   11		              Artificial Intelligence (OFAI),
   12			      Vienna, Austria
   13
   14    This software is part of Leslie De Koninck's master thesis, supervised
   15    by Bart Demoen and daily advisor Tom Schrijvers.  It is based on CLP(Q,R)
   16    by Christian Holzbaur for SICStus Prolog and distributed under the
   17    license details below with permission from all mentioned authors.
   18
   19    This program is free software; you can redistribute it and/or
   20    modify it under the terms of the GNU General Public License
   21    as published by the Free Software Foundation; either version 2
   22    of the License, or (at your option) any later version.
   23
   24    This program is distributed in the hope that it will be useful,
   25    but WITHOUT ANY WARRANTY; without even the implied warranty of
   26    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   27    GNU General Public License for more details.
   28
   29    You should have received a copy of the GNU Lesser General Public
   30    License along with this library; if not, write to the Free Software
   31    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
   32
   33    As a special exception, if you link this library with other files,
   34    compiled with a Free Software compiler, to produce an executable, this
   35    library does not by itself cause the resulting executable to be covered
   36    by the GNU General Public License. This exception does not however
   37    invalidate any other reasons why the executable file might be covered by
   38    the GNU General Public License.
   39*/
   40
   41
   42:- module(nf_r,
   43	[
   44	    {}/1,
   45	    nf/2,
   46	    entailed/1,
   47	    split/3,
   48	    repair/2,
   49	    nf_constant/2,
   50	    wait_linear/3,
   51	    nf2term/2
   52	]).   53
   54:- use_module('../clpqr/geler',
   55	[
   56	    geler/3
   57	]).   58:- use_module(bv_r,
   59	[
   60	    export_binding/2,
   61	    log_deref/4,
   62	    solve/1,
   63	    'solve_<'/1,
   64	    'solve_=<'/1,
   65	    'solve_=\\='/1
   66	]).   67:- use_module(ineq_r,
   68	[
   69	    ineq_one/4,
   70	    ineq_one_s_p_0/1,
   71	    ineq_one_s_n_0/1,
   72	    ineq_one_n_p_0/1,
   73	    ineq_one_n_n_0/1
   74	]).   75:- use_module(store_r,
   76	[
   77	    add_linear_11/3,
   78	    normalize_scalar/2
   79	]).   80:- use_module('../clpqr/highlight', []).   81:- use_module(library(error), [instantiation_error/1, type_error/2]).   82
   83
   84goal_expansion(geler(X,Y),geler(clpr,X,Y)).
   85
   86% -------------------------------------------------------------------------
   87
   88% {Constraint}
   89%
   90% Adds the constraint Constraint to the constraint store.
   91%
   92% First rule is to prevent binding with other rules when a variable is input
   93% Constraints are converted to normal form and if necessary, submitted to the linear
   94% equality/inequality solver (bv + ineq) or to the non-linear store (geler)
   95
   96{Rel} :-
   97	var(Rel),
   98	!,
   99	instantiation_error(Rel).
  100{R,Rs} :-
  101	!,
  102	{R},{Rs}.
  103{R;Rs} :-
  104	!,
  105	({R};{Rs}). % for entailment checking
  106{L < R} :-
  107	!,
  108	nf(L-R,Nf),
  109	submit_lt(Nf).
  110{L > R} :-
  111	!,
  112	nf(R-L,Nf),
  113	submit_lt(Nf).
  114{L =< R} :-
  115	!,
  116	nf(L-R,Nf),
  117	submit_le( Nf).
  118{<=(L,R)} :-
  119	!,
  120	nf(L-R,Nf),
  121	submit_le(Nf).
  122{L >= R} :-
  123	!,
  124	nf(R-L,Nf),
  125	submit_le(Nf).
  126{L =\= R} :-
  127	!,
  128	nf(L-R,Nf),
  129	submit_ne(Nf).
  130{L =:= R} :-
  131	!,
  132	nf(L-R,Nf),
  133	submit_eq(Nf).
  134{L = R} :-
  135	!,
  136	nf(L-R,Nf),
  137	submit_eq(Nf).
  138{Rel} :-
  139	type_error(clpr_constraint, Rel).
  140
  141% entailed(C)
  142%
  143% s -> c = ~s v c = ~(s /\ ~c)
  144% where s is the store and c is the constraint for which
  145% we want to know whether it is entailed.
  146% C is negated and added to the store. If this fails, then c is entailed by s
  147
  148entailed(C) :-
  149	negate(C,Cn),
  150	\+ {Cn}.
  151
  152% negate(C,Res).
  153%
  154% Res is the negation of constraint C
  155% first rule is to prevent binding with other rules when a variable is input
  156
  157negate(Rel,_) :-
  158	var(Rel),
  159	!,
  160	throw(instantiation_error(entailed(Rel),1)).
  161negate((A,B),(Na;Nb)) :-
  162	!,
  163	negate(A,Na),
  164	negate(B,Nb).
  165negate((A;B),(Na,Nb)) :-
  166	!,
  167	negate(A,Na),
  168	negate(B,Nb).
  169negate(A<B,A>=B) :- !.
  170negate(A>B,A=<B) :- !.
  171negate(A=<B,A>B) :- !.
  172negate(A>=B,A<B) :- !.
  173negate(A=:=B,A=\=B) :- !.
  174negate(A=B,A=\=B) :- !.
  175negate(A=\=B,A=:=B) :- !.
  176negate(Rel,_) :-
  177	type_error(clpr_constraint, Rel).
  178
  179% submit_eq(Nf)
  180%
  181% Submits the equality Nf = 0 to the constraint store, where Nf is in normal form.
  182% The following cases may apply:
  183% a) Nf = []
  184% b) Nf = [A]
  185%	b1) A = k
  186%	b2) invertible(A)
  187%	b3) linear -> A = 0
  188%	b4) nonlinear -> geler
  189% c) Nf=[A,B|Rest]
  190%	c1) A=k
  191%		c11) (B=c*X^+1 or B=c*X^-1), Rest=[] -> B=-k/c or B=-c/k
  192%		c12) invertible(A,B)
  193%		c13) linear(B|Rest)
  194%		c14) geler
  195%	c2) linear(Nf)
  196%	c3) nonlinear -> geler
  197
  198submit_eq([]).	% trivial success: case a
  199submit_eq([T|Ts]) :-
  200	submit_eq(Ts,T).
  201submit_eq([],A) :- submit_eq_b(A).	% case b
  202submit_eq([B|Bs],A) :- submit_eq_c(A,B,Bs).	% case c
  203
  204% submit_eq_b(A)
  205%
  206% Handles case b of submit_eq/1
  207
  208% case b1: A is a constant (non-zero)
  209submit_eq_b(v(_,[])) :-
  210	!,
  211	fail.
  212% case b2/b3: A is n*X^P => X = 0
  213submit_eq_b(v(_,[X^P])) :-
  214	var(X),
  215	P > 0,
  216	!,
  217	export_binding(X,0.0).
  218% case b2: non-linear is invertible: NL(X) = 0 => X - inv(NL)(0) = 0
  219submit_eq_b(v(_,[NL^1])) :-
  220	nonvar(NL),
  221	nl_invertible(NL,X,0.0,Inv),
  222	!,
  223	nf(-Inv,S),
  224	nf_add(X,S,New),
  225	submit_eq(New).
  226% case b4: A is non-linear and not invertible => submit equality to geler
  227submit_eq_b(Term) :-
  228	term_variables(Term,Vs),
  229	geler(Vs,nf_r:resubmit_eq([Term])).
  230
  231% submit_eq_c(A,B,Rest)
  232%
  233% Handles case c of submit_eq/1
  234
  235% case c1: A is a constant
  236submit_eq_c(v(I,[]),B,Rest) :-
  237	!,
  238	submit_eq_c1(Rest,B,I).
  239% case c2: A,B and Rest are linear
  240submit_eq_c(A,B,Rest) :-	% c2
  241	A = v(_,[X^1]),
  242	var(X),
  243	B = v(_,[Y^1]),
  244	var(Y),
  245	linear(Rest),
  246	!,
  247	Hom = [A,B|Rest],
  248	% 'solve_='(Hom).
  249	nf_length(Hom,0,Len),
  250	log_deref(Len,Hom,[],HomD),
  251	solve(HomD).
  252% case c3: A, B or Rest is non-linear => geler
  253submit_eq_c(A,B,Rest) :-
  254	Norm = [A,B|Rest],
  255	term_variables(Norm,Vs),
  256	geler(Vs,nf_r:resubmit_eq(Norm)).
  257
  258% submit_eq_c1(Rest,B,K)
  259%
  260% Handles case c1 of submit_eq/1
  261
  262% case c11a:
  263% i+kX^p=0 if p is an odd integer
  264% special case: one solution if P is negativ but also for a negative X
  265submit_eq_c1([], v(K,[X^P]), I) :-
  266    var(X),
  267    P =\= 0,
  268    0 > (-I/K),
  269    integer(P) =:= P,
  270	1 =:= integer(P) mod 2,
  271	!,
  272	Val is -((I/K) ** (1/P)),
  273	export_binding(X,Val).
  274% case c11b:
  275% i+kX^p=0 for p =\= 0, integer(P) =:= P
  276% special case: generate 2 solutions if p is a positive even integer
  277submit_eq_c1([], v(K,[X^P]), I) :-
  278    var(X),
  279    P =\= 0,
  280    0 =< (-I/K),
  281    integer(P) =:= P,
  282	0 =:= integer(P) mod 2,
  283	!,
  284	Val is (-I/K) ** (1/P),
  285	(	export_binding(X,Val)
  286	;
  287		ValNeg is -Val,
  288		export_binding(X, ValNeg)
  289    ).
  290% case c11c:
  291% i+kX^p=0 for p =\= 0, 0 =< (-I/K)
  292submit_eq_c1([], v(K,[X^P]), I) :-
  293    var(X),
  294    P =\= 0,
  295    0 =< (-I/K),
  296	!,
  297	Val is (-I/K) ** (1/P),
  298	export_binding(X,Val).
  299% case c11d: fail if var(X) and none of the above.
  300submit_eq_c1([], v(_K,[X^_P]), _I) :-
  301    var(X),
  302    !,
  303    fail.
  304% case c11e: fail for { -25 = _X^2.5 } and { -25 = _X^(-2.5) } and may be others!
  305%			 if you uncomment this case { -25 = _X^2.5 } throw an error(evaluation_error(undefined))
  306%			 and { -25 = _X^(-2.5) } succeed with an unbound X
  307submit_eq_c1([], v(K,[X^P]), I) :-
  308    nonvar(X),
  309    X = exp(_,_),   % TLS added 03/12
  310    1 =:= abs(P),
  311    0 >= I,
  312    0 >= K,
  313    !,
  314    fail.
  315% case c12: non-linear, invertible: cNL(X)^1+k=0 => inv(NL)(-k/c) = 0 ;
  316%				    cNL(X)^-1+k=0 => inv(NL)(-c/k) = 0
  317submit_eq_c1([],v(K,[NL^P]),I) :-
  318	nonvar(NL),
  319	(   P =:=  1,
  320	    Y is -I/K
  321	;   P =:= -1,
  322	    Y is -K/I
  323	),
  324	nl_invertible(NL,X,Y,Inv),
  325	!,
  326	nf(-Inv,S),
  327	nf_add(X,S,New),
  328	submit_eq(New).
  329% case c13: linear: X + Y + Z + c = 0 =>
  330submit_eq_c1(Rest,B,I) :-
  331	B = v(_,[Y^1]),
  332	var(Y),
  333	linear(Rest),
  334	!,
  335	% 'solve_='( [v(I,[]),B|Rest]).
  336	Hom = [B|Rest],
  337	nf_length(Hom,0,Len),
  338	normalize_scalar(I,Nonvar),
  339	log_deref(Len,Hom,[],HomD),
  340	add_linear_11(Nonvar,HomD,LinD),
  341	solve(LinD).
  342% case c14: other cases => geler
  343submit_eq_c1(Rest,B,I) :-
  344	Norm = [v(I,[]),B|Rest],
  345	term_variables(Norm,Vs),
  346	geler(Vs,nf_r:resubmit_eq(Norm)).
  347
  348% -----------------------------------------------------------------------
  349
  350% submit_lt(Nf)
  351%
  352% Submits the inequality Nf<0 to the constraint store, where Nf is in normal form.
  353
  354% 0 < 0 => fail
  355submit_lt([]) :- fail.
  356% A + B < 0
  357submit_lt([A|As]) :- submit_lt(As,A).
  358
  359% submit_lt(As,A)
  360%
  361% Does what submit_lt/1 does where Nf = [A|As]
  362
  363% v(K,P) < 0
  364submit_lt([],v(K,P)) :- submit_lt_b(P,K).
  365% A + B + Bs < 0
  366submit_lt([B|Bs],A) :- submit_lt_c(Bs,A,B).
  367
  368% submit_lt_b(P,K)
  369%
  370% Does what submit_lt/2 does where A = [v(K,P)] and As = []
  371
  372% c < 0
  373submit_lt_b([],I) :-
  374	!,
  375	I < -1.0e-10.
  376% cX^1 < 0 : if c < 0 then X > 0, else X < 0
  377submit_lt_b([X^1],K) :-
  378	var(X),
  379	!,
  380	(   K > 1.0e-10
  381	->  ineq_one_s_p_0(X)	% X is strictly negative
  382	;   ineq_one_s_n_0(X)	% X is strictly positive
  383	).
  384% non-linear => geler
  385submit_lt_b(P,K) :-
  386	term_variables(P,Vs),
  387	geler(Vs,nf_r:resubmit_lt([v(K,P)])).
  388
  389% submit_lt_c(Bs,A,B)
  390%
  391% Does what submit_lt/2 does where As = [B|Bs].
  392
  393%  c + kX < 0 => kX < c
  394submit_lt_c([],A,B) :-
  395	A = v(I,[]),
  396	B = v(K,[Y^1]),
  397	var(Y),
  398	!,
  399	ineq_one(strict,Y,K,I).
  400% linear < 0 => solve, non-linear < 0 => geler
  401submit_lt_c(Rest,A,B) :-
  402	Norm = [A,B|Rest],
  403	(   linear(Norm)
  404	->  'solve_<'(Norm)
  405	;   term_variables(Norm,Vs),
  406	    geler(Vs,nf_r:resubmit_lt(Norm))
  407	).
  408
  409% submit_le(Nf)
  410%
  411% Submits the inequality Nf =< 0 to the constraint store, where Nf is in normal form.
  412% See also submit_lt/1
  413
  414% 0 =< 0 => success
  415submit_le([]).
  416% A + B =< 0
  417submit_le([A|As]) :- submit_le(As,A).
  418
  419% submit_le(As,A)
  420%
  421% See submit_lt/2. This handles less or equal.
  422
  423% v(K,P) =< 0
  424submit_le([],v(K,P)) :- submit_le_b(P,K).
  425% A + B + Bs =< 0
  426submit_le([B|Bs],A) :- submit_le_c(Bs,A,B).
  427
  428% submit_le_b(P,K)
  429%
  430% See submit_lt_b/2. This handles less or equal.
  431
  432% c =< 0
  433submit_le_b([],I) :-
  434	!,
  435	I < 1.0e-10.
  436% cX^1 =< 0: if c < 0 then X >= 0, else X =< 0
  437submit_le_b([X^1],K) :-
  438	var(X),
  439	!,
  440	(   K > 1.0e-10
  441	->  ineq_one_n_p_0(X)	% X is non-strictly negative
  442	;   ineq_one_n_n_0(X)	% X is non-strictly positive
  443	).
  444%  cX^P =< 0 => geler
  445submit_le_b(P,K) :-
  446	term_variables(P,Vs),
  447	geler(Vs,nf_r:resubmit_le([v(K,P)])).
  448
  449% submit_le_c(Bs,A,B)
  450%
  451% See submit_lt_c/3. This handles less or equal.
  452
  453% c + kX^1 =< 0 => kX =< 0
  454submit_le_c([],A,B) :-
  455	A = v(I,[]),
  456	B = v(K,[Y^1]),
  457	var(Y),
  458	!,
  459	ineq_one(nonstrict,Y,K,I).
  460% A, B & Rest are linear => solve, otherwise => geler
  461submit_le_c(Rest,A,B) :-
  462	Norm = [A,B|Rest],
  463	(   linear(Norm)
  464	->  'solve_=<'(Norm)
  465	;   term_variables(Norm,Vs),
  466	    geler(Vs,nf_r:resubmit_le(Norm))
  467	).
  468
  469% submit_ne(Nf)
  470%
  471% Submits the inequality Nf =\= 0 to the constraint store, where Nf is in normal form.
  472% if Nf is a constant => check constant = 0, else if Nf is linear => solve else => geler
  473
  474submit_ne(Norm1) :-
  475	(   nf_constant(Norm1,K)
  476	->  \+ (K >= -1.0e-10, K =< 1.0e-10) % K =\= 0
  477	;   linear(Norm1)
  478	->  'solve_=\\='(Norm1)
  479	;   term_variables(Norm1,Vs),
  480	    geler(Vs,nf_r:resubmit_ne(Norm1))
  481	).
  482
  483% linear(A)
  484%
  485% succeeds when A is linear: all elements are of the form v(_,[]) or v(_,[X^1])
  486
  487linear([]).
  488linear(v(_,Ps)) :- linear_ps(Ps).
  489linear([A|As]) :-
  490	linear(A),
  491	linear(As).
  492
  493% linear_ps(A)
  494%
  495% Succeeds when A = V^1 with V a variable.
  496% This reflects the linearity of v(_,A).
  497
  498linear_ps([]).
  499linear_ps([V^1]) :- var(V).	% excludes sin(_), ...
  500
  501%
  502% Goal delays until Term gets linear.
  503% At this time, Var will be bound to the normalform of Term.
  504%
  505:- meta_predicate wait_linear( ?, ?, :).  506%
  507wait_linear(Term,Var,Goal) :-
  508	nf(Term,Nf),
  509	(   linear(Nf)
  510	->  Var = Nf,
  511	    call(Goal)
  512	;   term_variables(Nf,Vars),
  513	    geler(Vars,nf_r:wait_linear_retry(Nf,Var,Goal))
  514	).
  515%
  516% geler clients
  517%
  518resubmit_eq(N) :-
  519	repair(N,Norm),
  520	submit_eq(Norm).
  521resubmit_lt(N) :-
  522	repair(N,Norm),
  523	submit_lt(Norm).
  524resubmit_le(N) :-
  525	repair(N,Norm),
  526	submit_le(Norm).
  527resubmit_ne(N) :-
  528	repair(N,Norm),
  529	submit_ne(Norm).
  530wait_linear_retry(Nf0,Var,Goal) :-
  531	repair(Nf0,Nf),
  532	(   linear(Nf)
  533	->  Var = Nf,
  534	    call(Goal)
  535	;   term_variables(Nf,Vars),
  536	    geler(Vars,nf_r:wait_linear_retry(Nf,Var,Goal))
  537	).
  538% -----------------------------------------------------------------------
  539
  540% nl_invertible(F,X,Y,Res)
  541%
  542% Res is the evaluation of the inverse of nonlinear function F in variable X
  543% where X is Y
  544
  545nl_invertible(sin(X),X,Y,Res) :- Res is asin(Y).
  546nl_invertible(cos(X),X,Y,Res) :- Res is acos(Y).
  547nl_invertible(tan(X),X,Y,Res) :- Res is atan(Y).
  548nl_invertible(exp(B,C),X,A,Res) :-
  549	(   nf_constant(B,Kb)
  550	->  A > 1.0e-10,
  551	    Kb > 1.0e-10,
  552	    TestKb is Kb - 1.0, % Kb =\= 1.0
  553	    \+ (TestKb >= -1.0e-10, TestKb =< 1.0e-10),
  554	    X = C, % note delayed unification
  555	    Res is log(A)/log(Kb)
  556	;   nf_constant(C,Kc),
  557	    \+ (A >= -1.0e-10, A =< 1.0e-10), % A =\= 0
  558	    Kc > 1.0e-10, % Kc > 0
  559	    X = B, % note delayed unification
  560	    Res is A**(1.0/Kc)
  561	).
  562
  563% -----------------------------------------------------------------------
  564
  565% nf(Exp,Nf)
  566%
  567% Returns in Nf, the normal form of expression Exp
  568%
  569% v(A,[B^C,D^E|...]) means A*B^C*D^E*... where A is a scalar (number)
  570% v(A,[]) means scalar A
  571
  572% variable X => 1*X^1
  573nf(X,Norm) :-
  574	var(X),
  575	!,
  576	Norm = [v(1.0,[X^1])].
  577nf(X,Norm) :-
  578	number(X),
  579	!,
  580	nf_number(X,Norm).
  581%
  582nf(#(Const),Norm) :-
  583	monash_constant(Const,Value),
  584	!,
  585	Norm = [v(Value,[])].
  586%
  587nf(-A,Norm) :-
  588	!,
  589	nf(A,An),
  590	nf_mul_factor(v(-1.0,[]),An,Norm).
  591nf(+A,Norm) :-
  592	!,
  593	nf(A,Norm).
  594%
  595nf(A+B,Norm) :-
  596	!,
  597	nf(A,An),
  598	nf(B,Bn),
  599	nf_add(An,Bn,Norm).
  600nf(A-B,Norm) :-
  601	!,
  602	nf(A,An),
  603	nf(-B,Bn),
  604	nf_add(An,Bn,Norm).
  605%
  606nf(A*B,Norm) :-
  607	!,
  608	nf(A,An),
  609	nf(B,Bn),
  610	nf_mul(An,Bn,Norm).
  611nf(A/B,Norm) :-
  612	!,
  613	nf(A,An),
  614	nf(B,Bn),
  615	nf_div(Bn,An,Norm).
  616% non-linear function, one argument: Term = f(Arg) equals f'(Sa1) = Skel
  617nf(Term,Norm) :-
  618	nonlin_1(Term,Arg,Skel,Sa1),
  619	!,
  620	nf(Arg,An),
  621	nf_nonlin_1(Skel,An,Sa1,Norm).
  622% non-linear function, two arguments: Term = f(A1,A2) equals f'(Sa1,Sa2) = Skel
  623nf(Term,Norm) :-
  624	nonlin_2(Term,A1,A2,Skel,Sa1,Sa2),
  625	!,
  626	nf(A1,A1n),
  627	nf(A2,A2n),
  628	nf_nonlin_2(Skel,A1n,A2n,Sa1,Sa2,Norm).
  629%
  630nf(Term,_) :-
  631	type_error(clpr_expression, Term).
  632
  633% nf_number(N,Res)
  634%
  635% If N is a number, N is normalized
  636
  637nf_number(N,Res) :-
  638	number(N),
  639	(   (N >= -1.0e-10, N =< 1.0e-10) % N =:= 0
  640	->  Res = []
  641	;   Res = [v(N,[])]
  642	).
  643
  644nonlin_1(abs(X),X,abs(Y),Y).
  645nonlin_1(sin(X),X,sin(Y),Y).
  646nonlin_1(cos(X),X,cos(Y),Y).
  647nonlin_1(tan(X),X,tan(Y),Y).
  648nonlin_2(min(A,B),A,B,min(X,Y),X,Y).
  649nonlin_2(max(A,B),A,B,max(X,Y),X,Y).
  650nonlin_2(exp(A,B),A,B,exp(X,Y),X,Y).
  651nonlin_2(pow(A,B),A,B,exp(X,Y),X,Y).	% pow->exp
  652nonlin_2(A^B,A,B,exp(X,Y),X,Y).
  653
  654nf_nonlin_1(Skel,An,S1,Norm) :-
  655	(   nf_constant(An,S1)
  656	->  nl_eval(Skel,Res),
  657	    nf_number(Res,Norm)
  658	;   S1 = An,
  659	    Norm = [v(1.0,[Skel^1])]).
  660nf_nonlin_2(Skel,A1n,A2n,S1,S2,Norm) :-
  661	(   nf_constant(A1n,S1),
  662	    nf_constant(A2n,S2)
  663	->  nl_eval(Skel,Res),
  664	    nf_number(Res,Norm)
  665	;   Skel=exp(_,_),
  666	    nf_constant(A2n,Exp),
  667	    integerp(Exp,I)
  668	->  nf_power(I,A1n,Norm)
  669	;   S1 = A1n,
  670	    S2 = A2n,
  671	    Norm = [v(1.0,[Skel^1])]
  672	).
  673
  674% evaluates non-linear functions in one variable where the variable is bound
  675nl_eval(abs(X),R) :- R is abs(X).
  676nl_eval(sin(X),R) :- R is sin(X).
  677nl_eval(cos(X),R) :- R is cos(X).
  678nl_eval(tan(X),R) :- R is tan(X).
  679% evaluates non-linear functions in two variables where both variables are
  680% bound
  681nl_eval(min(X,Y),R) :- R is min(X,Y).
  682nl_eval(max(X,Y),R) :- R is max(X,Y).
  683nl_eval(exp(X,Y),R) :- R is X**Y.
  684
  685monash_constant(X,_) :-
  686	var(X),
  687	!,
  688	fail.
  689monash_constant(p,3.14259265).
  690monash_constant(pi,3.14259265).
  691monash_constant(e,2.71828182).
  692monash_constant(zero,1.0e-10).
  693
  694%
  695% check if a Nf consists of just a constant
  696%
  697
  698nf_constant([],0.0).
  699nf_constant([v(K,[])],K).
  700
  701% split(NF,SNF,C)
  702%
  703% splits a normalform expression NF into two parts:
  704% - a constant term C (which might be 0)
  705% - the homogene part of the expression
  706%
  707% this method depends on the polynf ordering, i.e. [] < [X^1] ...
  708
  709split([],[],0.0).
  710split([First|T],H,I) :-
  711	(   First = v(I,[])
  712	->  H = T
  713	;   I = 0.0,
  714	    H = [First|T]
  715	).
  716
  717% nf_add(A,B,C): merges two normalized additions into a new normalized addition
  718%
  719% a normalized addition is one where the terms are ordered, e.g. X^1 < Y^1, X^1 < X^2 etc.
  720% terms in the same variable with the same exponent are added,
  721% e.g. when A contains v(5,[X^1]) and B contains v(4,[X^1]) then C contains v(9,[X^1]).
  722
  723nf_add([],Bs,Bs).
  724nf_add([A|As],Bs,Cs) :- nf_add(Bs,A,As,Cs).
  725
  726nf_add([],A,As,Cs) :- Cs = [A|As].
  727nf_add([B|Bs],A,As,Cs) :-
  728	A = v(Ka,Pa),
  729	B = v(Kb,Pb),
  730	compare(Rel,Pa,Pb),
  731	nf_add_case(Rel,A,As,Cs,B,Bs,Ka,Kb,Pa).
  732
  733% nf_add_case(Rel,A,As,Cs,B,Bs,Ka,Kb,Pa)
  734%
  735% merges sorted lists [A|As] and [B|Bs] into new sorted list Cs
  736% A = v(Ka,Pa) and B = v(Kb,_)
  737% Rel is the ordering relation (<, > or =) between A and B.
  738% when Rel is =, Ka and Kb are added to form a new scalar for Pa
  739nf_add_case(<,A,As,Cs,B,Bs,_,_,_) :-
  740	Cs = [A|Rest],
  741	nf_add(As,B,Bs,Rest).
  742nf_add_case(>,A,As,Cs,B,Bs,_,_,_) :-
  743	Cs = [B|Rest],
  744	nf_add(Bs,A,As,Rest).
  745nf_add_case(=,_,As,Cs,_,Bs,Ka,Kb,Pa) :-
  746	Kc is Ka + Kb,
  747	(   (Kc >= -1.0e-10, Kc =< 1.0e-10) % Kc =:= 0.0
  748	->  nf_add(As,Bs,Cs)
  749	;   Cs = [v(Kc,Pa)|Rest],
  750	    nf_add(As,Bs,Rest)
  751	).
  752
  753nf_mul(A,B,Res) :-
  754	nf_length(A,0,LenA),
  755	nf_length(B,0,LenB),
  756	nf_mul_log(LenA,A,[],LenB,B,Res).
  757
  758nf_mul_log(0,As,As,_,_,[]) :- !.
  759nf_mul_log(1,[A|As],As,Lb,B,R) :-
  760	!,
  761	nf_mul_factor_log(Lb,B,[],A,R).
  762nf_mul_log(2,[A1,A2|As],As,Lb,B,R) :-
  763	!,
  764	nf_mul_factor_log(Lb,B,[],A1,A1b),
  765	nf_mul_factor_log(Lb,B,[],A2,A2b),
  766	nf_add(A1b,A2b,R).
  767nf_mul_log(N,A0,A2,Lb,B,R) :-
  768	P is N>>1,
  769	Q is N-P,
  770	nf_mul_log(P,A0,A1,Lb,B,Rp),
  771	nf_mul_log(Q,A1,A2,Lb,B,Rq),
  772	nf_add(Rp,Rq,R).
  773
  774
  775% nf_add_2: does the same thing as nf_add, but only has 2 elements to combine.
  776nf_add_2(Af,Bf,Res) :-		%  unfold: nf_add([Af],[Bf],Res).
  777	Af = v(Ka,Pa),
  778	Bf = v(Kb,Pb),
  779	compare(Rel,Pa,Pb),
  780	nf_add_2_case(Rel,Af,Bf,Res,Ka,Kb,Pa).
  781
  782nf_add_2_case(<,Af,Bf,[Af,Bf],_,_,_).
  783nf_add_2_case(>,Af,Bf,[Bf,Af],_,_,_).
  784nf_add_2_case(=,_, _,Res,Ka,Kb,Pa) :-
  785	Kc is Ka + Kb,
  786	(   (Kc >= -1.0e-10, Kc =< 1.0e-10) % Kc =:= 0
  787	->  Res = []
  788	;   Res = [v(Kc,Pa)]
  789	).
  790
  791% nf_mul_k(A,B,C)
  792%
  793% C is the result of the multiplication of each element of A (of the form v(_,_)) with scalar B (which shouldn't be 0)
  794nf_mul_k([],_,[]).
  795nf_mul_k([v(I,P)|Vs],K,[v(Ki,P)|Vks]) :-
  796	Ki is K*I,
  797	nf_mul_k(Vs,K,Vks).
  798
  799% nf_mul_factor(A,Sum,Res)
  800%
  801% multiplies each element of the list Sum with factor A which is of the form v(_,_)
  802% and puts the result in the sorted list Res.
  803nf_mul_factor(v(K,[]),Sum,Res) :-
  804	!,
  805	nf_mul_k(Sum,K,Res).
  806nf_mul_factor(F,Sum,Res) :-
  807	nf_length(Sum,0,Len),
  808	nf_mul_factor_log(Len,Sum,[],F,Res).
  809
  810% nf_mul_factor_log(Len,[Sum|SumTail],SumTail,F,Res)
  811%
  812% multiplies each element of Sum with F and puts the result in the sorted list Res
  813% Len is the length of Sum
  814% Sum is split logarithmically each step
  815
  816nf_mul_factor_log(0,As,As,_,[]) :- !.
  817nf_mul_factor_log(1,[A|As],As,F,[R]) :-
  818	!,
  819	mult(A,F,R).
  820nf_mul_factor_log(2,[A,B|As],As,F,Res) :-
  821	!,
  822	mult(A,F,Af),
  823	mult(B,F,Bf),
  824	nf_add_2(Af,Bf,Res).
  825nf_mul_factor_log(N,A0,A2,F,R) :-
  826	P is N>>1, % P is rounded(N/2)
  827	Q is N-P,
  828	nf_mul_factor_log(P,A0,A1,F,Rp),
  829	nf_mul_factor_log(Q,A1,A2,F,Rq),
  830	nf_add(Rp,Rq,R).
  831
  832% mult(A,B,C)
  833%
  834% multiplies A and B into C each of the form v(_,_)
  835
  836mult(v(Ka,La),v(Kb,Lb),v(Kc,Lc)) :-
  837	Kc is Ka*Kb,
  838	pmerge(La,Lb,Lc).
  839
  840% pmerge(A,B,C)
  841%
  842% multiplies A and B into sorted C, where each is of the form of the second argument of v(_,_)
  843
  844pmerge([],Bs,Bs).
  845pmerge([A|As],Bs,Cs) :- pmerge(Bs,A,As,Cs).
  846
  847pmerge([],A,As,Res) :- Res = [A|As].
  848pmerge([B|Bs],A,As,Res) :-
  849	A = Xa^Ka,
  850	B = Xb^Kb,
  851	compare(R,Xa,Xb),
  852	pmerge_case(R,A,As,Res,B,Bs,Ka,Kb,Xa).
  853
  854% pmerge_case(Rel,A,As,Res,B,Bs,Ka,Kb,Xa)
  855%
  856% multiplies and sorts [A|As] with [B|Bs] into Res where each is of the form of
  857% the second argument of v(_,_)
  858%
  859% A is Xa^Ka and B is Xb^Kb, Rel is ordening relation between Xa and Xb
  860
  861pmerge_case(<,A,As,Res,B,Bs,_,_,_) :-
  862	Res = [A|Tail],
  863	pmerge(As,B,Bs,Tail).
  864pmerge_case(>,A,As,Res,B,Bs,_,_,_) :-
  865	Res = [B|Tail],
  866	pmerge(Bs,A,As,Tail).
  867pmerge_case(=,_,As,Res,_,Bs,Ka,Kb,Xa) :-
  868	Kc is Ka + Kb,
  869	(   Kc =:= 0
  870	->  pmerge(As,Bs,Res)
  871	;   Res = [Xa^Kc|Tail],
  872	    pmerge(As,Bs,Tail)
  873	).
  874
  875% nf_div(Factor,In,Out)
  876%
  877% Out is the result of the division of each element in In (which is of the form v(_,_)) by Factor.
  878
  879% division by zero
  880nf_div([],_,_) :-
  881	!,
  882	zero_division.
  883% division by v(K,P) => multiplication by v(1/K,P^-1)
  884nf_div([v(K,P)],Sum,Res) :-
  885	!,
  886	Ki is 1.0/K,
  887	mult_exp(P,-1,Pi),
  888	nf_mul_factor(v(Ki,Pi),Sum,Res).
  889nf_div(D,A,[v(1.0,[(A/D)^1])]).
  890
  891% zero_division
  892%
  893% called when a division by zero is performed
  894zero_division :- fail.	% raise_exception(_) ?
  895
  896% mult_exp(In,Factor,Out)
  897%
  898% Out is the result of the multiplication of the exponents of the elements in In
  899% (which are of the form X^Exp by Factor.
  900mult_exp([],_,[]).
  901mult_exp([X^P|Xs],K,[X^I|Tail]) :-
  902	I is K*P,
  903	mult_exp(Xs,K,Tail).
  904%
  905% raise to integer powers
  906%
  907% | ?- time({(1+X+Y+Z)^15=0}). (sicstus, try with SWI)
  908% Timing 00:00:02.610	  2.610   iterative
  909% Timing 00:00:00.660	  0.660   binomial
  910nf_power(N,Sum,Norm) :-
  911	integer(N),
  912	compare(Rel,N,0),
  913	(   Rel = (<)
  914	->  Pn is -N,
  915	    % nf_power_pos(Pn,Sum,Inorm),
  916	    binom(Sum,Pn,Inorm),
  917	    nf_div(Inorm,[v(1.0,[])],Norm)
  918	;   Rel = (>)
  919	->  % nf_power_pos(N,Sum,Norm)
  920	    binom(Sum,N,Norm)
  921	;   Rel = (=)
  922	->  % 0^0 is indeterminate but we say 1
  923	    Norm = [v(1.0,[])]
  924	).
  925%
  926% N>0
  927%
  928% iterative method: X^N = X*(X^N-1)
  929nf_power_pos(1,Sum,Norm) :-
  930	!,
  931	Sum = Norm.
  932nf_power_pos(N,Sum,Norm) :-
  933	N1 is N-1,
  934	nf_power_pos(N1,Sum,Pn1),
  935	nf_mul(Sum,Pn1,Norm).
  936%
  937% N>0
  938%
  939% binomial method
  940binom(Sum,1,Power) :-
  941	!,
  942	Power = Sum.
  943binom([],_,[]).
  944binom([A|Bs],N,Power) :-
  945	(   Bs = []
  946	->  nf_power_factor(A,N,Ap),
  947	    Power = [Ap]
  948	;   Bs = [_|_]
  949	->  factor_powers(N,A,v(1.0,[]),Pas),
  950	    sum_powers(N,Bs,[v(1.0,[])],Pbs,[]),
  951	    combine_powers(Pas,Pbs,0,N,1,[],Power)
  952	).
  953
  954combine_powers([],[],_,_,_,Pi,Pi).
  955combine_powers([A|As],[B|Bs],L,R,C,Pi,Po) :-
  956	nf_mul(A,B,Ab),
  957	nf_mul_k(Ab,C,Abc),
  958	nf_add(Abc,Pi,Pii),
  959	L1 is L+1,
  960	R1 is R-1,
  961	C1 is C*R//L1,
  962	combine_powers(As,Bs,L1,R1,C1,Pii,Po).
  963
  964nf_power_factor(v(K,P),N,v(Kn,Pn)) :-
  965	Kn is K**N,
  966	mult_exp(P,N,Pn).
  967
  968factor_powers(0,_,Prev,[[Prev]]) :- !.
  969factor_powers(N,F,Prev,[[Prev]|Ps]) :-
  970	N1 is N-1,
  971	mult(Prev,F,Next),
  972	factor_powers(N1,F,Next,Ps).
  973sum_powers(0,_,Prev,[Prev|Lt],Lt) :- !.
  974sum_powers(N,S,Prev,L0,Lt) :-
  975	N1 is N-1,
  976	nf_mul(S,Prev,Next),
  977	sum_powers(N1,S,Next,L0,[Prev|Lt]).
  978
  979% ------------------------------------------------------------------------------
  980repair(Sum,Norm) :-
  981	nf_length(Sum,0,Len),
  982	repair_log(Len,Sum,[],Norm).
  983repair_log(0,As,As,[]) :- !.
  984repair_log(1,[v(Ka,Pa)|As],As,R) :-
  985	!,
  986	repair_term(Ka,Pa,R).
  987repair_log(2,[v(Ka,Pa),v(Kb,Pb)|As],As,R) :-
  988	!,
  989	repair_term(Ka,Pa,Ar),
  990	repair_term(Kb,Pb,Br),
  991	nf_add(Ar,Br,R).
  992repair_log(N,A0,A2,R) :-
  993	P is N>>1,
  994	Q is N-P,
  995	repair_log(P,A0,A1,Rp),
  996	repair_log(Q,A1,A2,Rq),
  997	nf_add(Rp,Rq,R).
  998
  999repair_term(K,P,Norm) :-
 1000	length(P,Len),
 1001	repair_p_log(Len,P,[],Pr,[v(1.0,[])],Sum),
 1002	nf_mul_factor(v(K,Pr),Sum,Norm).
 1003
 1004repair_p_log(0,Ps,Ps,[],L0,L0) :- !.
 1005repair_p_log(1,[X^P|Ps],Ps,R,L0,L1) :-
 1006	!,
 1007	repair_p(X,P,R,L0,L1).
 1008repair_p_log(2,[X^Px,Y^Py|Ps],Ps,R,L0,L2) :-
 1009	!,
 1010	repair_p(X,Px,Rx,L0,L1),
 1011	repair_p(Y,Py,Ry,L1,L2),
 1012	pmerge(Rx,Ry,R).
 1013repair_p_log(N,P0,P2,R,L0,L2) :-
 1014	P is N>>1,
 1015	Q is N-P,
 1016	repair_p_log(P,P0,P1,Rp,L0,L1),
 1017	repair_p_log(Q,P1,P2,Rq,L1,L2),
 1018	pmerge(Rp,Rq,R).
 1019
 1020repair_p(Term,P,[Term^P],L0,L0) :- var(Term).
 1021repair_p(Term,P,[],L0,L1) :-
 1022	nonvar(Term),
 1023	repair_p_one(Term,TermN),
 1024	nf_power(P,TermN,TermNP),
 1025	nf_mul(TermNP,L0,L1).
 1026%
 1027% An undigested term a/b is distinguished from an
 1028% digested one by the fact that its arguments are
 1029% digested -> cuts after repair of args!
 1030%
 1031repair_p_one(Term,TermN) :-
 1032	nf_number(Term,TermN),	% freq. shortcut for nf/2 case below
 1033	!.
 1034repair_p_one(A1/A2,TermN) :-
 1035	repair(A1,A1n),
 1036	repair(A2,A2n),
 1037	!,
 1038	nf_div(A2n,A1n,TermN).
 1039repair_p_one(Term,TermN) :-
 1040	nonlin_1(Term,Arg,Skel,Sa),
 1041	repair(Arg,An),
 1042	!,
 1043	nf_nonlin_1(Skel,An,Sa,TermN).
 1044repair_p_one(Term,TermN) :-
 1045	nonlin_2(Term,A1,A2,Skel,Sa1,Sa2),
 1046	repair(A1,A1n),
 1047	repair(A2,A2n),
 1048	!,
 1049	nf_nonlin_2(Skel,A1n,A2n,Sa1,Sa2,TermN).
 1050repair_p_one(Term,TermN) :-
 1051	nf(Term,TermN).
 1052
 1053nf_length([],Li,Li).
 1054nf_length([_|R],Li,Lo) :-
 1055	Lii is Li+1,
 1056	nf_length(R,Lii,Lo).
 1057% ------------------------------------------------------------------------------
 1058% nf2term(NF,Term)
 1059%
 1060% transforms a normal form into a readable term
 1061
 1062% empty normal form = 0
 1063nf2term([],0.0).
 1064% term is first element (+ next elements)
 1065nf2term([F|Fs],T) :-
 1066	f02t(F,T0),	% first element
 1067	yfx(Fs,T0,T).	% next elements
 1068
 1069yfx([],T0,T0).
 1070yfx([F|Fs],T0,TN) :-
 1071	fn2t(F,Ft,Op),
 1072	T1 =.. [Op,T0,Ft],
 1073	yfx(Fs,T1,TN).
 1074
 1075% f02t(v(K,P),T)
 1076%
 1077% transforms the first element of the normal form (something of the form v(K,P))
 1078% into a readable term
 1079f02t(v(K,P),T) :-
 1080	(   % just a constant
 1081	    P = []
 1082	->  T = K
 1083	;   TestK is K - 1.0, % K =:= 1
 1084	    (TestK >= -1.0e-10, TestK =< 1.0e-10)
 1085	->  p2term(P,T)
 1086	;   TestK is K + 1.0, % K =:= -1
 1087	    (TestK >= -1.0e-10, TestK =< 1.0e-10)
 1088	->  T = -Pt,
 1089	    p2term(P,Pt)
 1090	;   T = K*Pt,
 1091	    p2term(P,Pt)
 1092	).
 1093
 1094% f02t(v(K,P),T,Op)
 1095%
 1096% transforms a next element of the normal form (something of the form v(K,P))
 1097% into a readable term
 1098fn2t(v(K,P),Term,Op) :-
 1099	(   TestK is K - 1.0, % K =:= 1
 1100	    (TestK >= -1.0e-10, TestK =< 1.0e-10)
 1101	->  Term = Pt,
 1102	    Op = +
 1103	;   TestK is K + 1.0, % K =:= -1
 1104	    (TestK >= -1.0e-10, TestK =< 1.0e-10)
 1105	->  Term = Pt,
 1106	    Op = -
 1107	;   K < -1.0e-10 % K < 0
 1108	->  Kf is -K,
 1109	    Term = Kf*Pt,
 1110	    Op = -
 1111	;   % K > 0
 1112	    Term = K*Pt,
 1113	    Op = +
 1114	),
 1115	p2term(P,Pt).
 1116
 1117% transforms the P part in v(_,P) into a readable term
 1118p2term([X^P|Xs],Term) :-
 1119	(   Xs = []
 1120	->  pe2term(X,Xt),
 1121	    exp2term(P,Xt,Term)
 1122	;   Xs = [_|_]
 1123	->  Term = Xst*Xtp,
 1124	    pe2term(X,Xt),
 1125	    exp2term(P,Xt,Xtp),
 1126	    p2term(Xs,Xst)
 1127	).
 1128
 1129%
 1130exp2term(1,X,X) :- !.
 1131exp2term(-1,X,1.0/X) :- !.
 1132exp2term(P,X,Term) :-
 1133	% Term = exp(X,Pn)
 1134	Term = X^P.
 1135
 1136pe2term(X,Term) :-
 1137	var(X),
 1138	Term = X.
 1139pe2term(X,Term) :-
 1140	nonvar(X),
 1141	X =.. [F|Args],
 1142	pe2term_args(Args,Argst),
 1143	Term =.. [F|Argst].
 1144
 1145pe2term_args([],[]).
 1146pe2term_args([A|As],[T|Ts]) :-
 1147	nf2term(A,T),
 1148	pe2term_args(As,Ts).
 1149
 1150% transg(Goal,[OutList|OutListTail],OutListTail)
 1151%
 1152% puts the equalities and inequalities that are implied by the elements in Goal
 1153% in the difference list OutList
 1154%
 1155% called by geler.pl for project.pl
 1156
 1157transg(resubmit_eq(Nf)) -->
 1158	{
 1159	    nf2term([],Z),
 1160	    nf2term(Nf,Term)
 1161	},
 1162	[clpr:{Term=Z}].
 1163transg(resubmit_lt(Nf)) -->
 1164	{
 1165	    nf2term([],Z),
 1166	    nf2term(Nf,Term)
 1167	},
 1168	[clpr:{Term<Z}].
 1169transg(resubmit_le(Nf)) -->
 1170	{
 1171	    nf2term([],Z),
 1172	    nf2term(Nf,Term)
 1173	},
 1174	[clpr:{Term=<Z}].
 1175transg(resubmit_ne(Nf)) -->
 1176	{
 1177	    nf2term([],Z),
 1178	    nf2term(Nf,Term)
 1179	},
 1180	[clpr:{Term=\=Z}].
 1181transg(wait_linear_retry(Nf,Res,Goal)) -->
 1182	{
 1183	    nf2term(Nf,Term)
 1184	},
 1185	[clpr:{Term=Res},Goal].
 1186
 1187integerp(X) :-
 1188	floor(X)=:=X.
 1189
 1190integerp(X,I) :-
 1191	floor(X)=:=X,
 1192	I is integer(X).
 1193
 1194		 /*******************************
 1195		 *	       SANDBOX		*
 1196		 *******************************/
 1197:- multifile
 1198	sandbox:safe_primitive/1. 1199
 1200sandbox:safe_primitive(nf_r:{_}).
 1201sandbox:safe_primitive(nf_r:entailed(_))