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