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