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