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