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