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