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