1% BEGIN LICENSE BLOCK 2% Version: CMPL 1.1 3% 4% The contents of this file are subject to the Cisco-style Mozilla Public 5% License Version 1.1 (the "License"); you may not use this file except 6% in compliance with the License. You may obtain a copy of the License 7% at www.eclipse-clp.org/license. 8% 9% Software distributed under the License is distributed on an "AS IS" 10% basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See 11% the License for the specific language governing rights and limitations 12% under the License. 13% 14% The Original Code is The ECLiPSe Constraint Logic Programming System. 15% The Initial Developer of the Original Code is Cisco Systems, Inc. 16% Portions created by the Initial Developer are 17% Copyright (C) 1997-2006 Cisco Systems, Inc. All Rights Reserved. 18% 19% Contributor(s): Joachim Schimpf, IC-Parc 20% 21% END LICENSE BLOCK 22% ---------------------------------------------------------------------- 23% System: ECLiPSe Constraint Logic Programming System 24% Version: $Id: linearize.pl,v 1.2 2013/06/15 14:44:38 jschimpf Exp $ 25% 26% Description: Expression simplifier 27% 28% Authors: J.Schimpf, IC-Parc 29% 30% 31% A linear expression is normalised into a list (sum) of the form 32% 33% [C0*1, C1*X1, C2*X2, ...] 34% 35% where Ci are numbers and Xi are distinct variables. 36% The first (constant) term is always present, Ci (i>=1) are nonzero. 37% 38% linearize/3 converts a general arithmetic expression into a 39% normalised linear form. The nonlinear parts are extracted using 40% auxiliary variables, and returned as Residue, which is a list of 41% Aux=NonLinExpr constraints. 42% 43% linearize/3 understands all arithmetic expressions plus 44% 45% lin(Lin) embedded normalised linear expression 46% 47% 48% linrenorm/2 re-normalises a linear normal form expression after variable 49% bindings, unifications. 50% 51% delinearize/2 converts back into normal arithmetic expression, but 52% that may often not be needed since linearize/3 understands lin(Lin). 53% 54% TODO: 55% -- E1*E2-clause of linearize/7 always makes a choicepoint... 56% 57% ---------------------------------------------------------------------- 58 59 60% ---------------------------------------------------------------------- 61:- module(linearize). 62% ---------------------------------------------------------------------- 63 64:- comment(categories, ["Data Structures"]). 65:- comment(summary, "Normalizers for arithmetic expressions"). 66:- comment(author, "Joachim Schimpf, IC-Parc"). 67:- comment(date, "$Date: 2013/06/15 14:44:38 $"). 68:- comment(copyright, "Cisco Systems, Inc."). 69 70:- export 71 linearize/3, 72 linrenorm/2, 73 delinearize/2. 74 75:- comment(linearize/3, [ 76 summary:"Split and arithmetic expression into linear and nonlinear parts", 77 amode:linearize(?,-,-), 78 args:["Expression":"Arithmetic expression with constants and variables", 79 "Linear":"Normalized lists of monomials", 80 "Residue":"Residual nonlinear components in the form AuxVar=Expr"], 81 see_also:[polynorm/3,linrenorm/2,delinearize/2], 82 desc:html("\ 83 The linear component of expression is normalised into a list 84 (sum) of monomials of the form 85<PRE> 86 [C0*1, C1*X1, C2*X2, ...] 87</PRE> 88 where Ci are numbers and Xi are distinct variables. 89 The first (constant) term is always present, Ci (i>=1) are nonzero. 90<P> 91 linearize/3 converts a general arithmetic expression into a 92 normalised linear form. The nonlinear parts are extracted using 93 auxiliary variables, and returned as Residue, which is a list of 94 Aux=NonLinExpr constraints. 95<P> 96 linearize/3 understands all arithmetic expressions plus 97<PRE> 98 lin(Lin) 99</PRE> 100 where Lin is an already normalised linear expression. All variables 101 within Expression (which are free at linearization time) are taken 102 to be numerical variables. If you intend to have variables which can 103 be bound to symbolic expressions rather than numbers, they must be 104 wrapped into an eval/1 functor. 105 "), 106 eg:" 107 ?- linearize(3*X-7*Y+2*(X+Y), L, R). 108 X = X 109 Y = Y 110 L = [0 * 1, 5 * X, -5 * Y] 111 R = [] 112 yes. 113 114 ?- linearize(12*3, L, R). 115 L = [36 * 1] 116 R = [] 117 yes. 118 119 ?- linearize(X*(3+7*Y)-X, L, R). 120 Y = Y 121 X = X 122 L = [0 * 1, 1 * _308, 2 * X] 123 R = [_308 = X * (7 * Y)] 124 yes. 125 " 126]). 127 128:- comment(linrenorm/2, [ 129 summary:"Renormalize a linear form", 130 amode:linrenorm(+,-), 131 args:["LinOld":"Possibly denormal lists of monomials", 132 "LinNew":"Normalized lists of monomials"], 133 see_also:[linearize/3,delinearize/2], 134 desc:html("\ 135 The normal form of linear expressions is 136<PRE> 137 [C0*1, C1*X1, C2*X2, ...] 138</PRE> 139 where Ci are numbers and Xi are distinct variables. 140 The first (constant) term is always present, Ci (i>=1) are nonzero. 141<P> 142 Such a form can become denormalized due to unifications 143 (instantiation or variable-variable aliasing). This predicate 144 renormalizes it. Note that variables may only become instantiated 145 to numbers! 146 "), 147 eg:" 148 ?- linearize(3*X-7*Y+2*(X+Y), L1, R), writeln(L1), 149 Y = 3, 150 linrenorm(L1,L2), writeln(L2). 151 152 [0 * 1, 5 * X, -5 * Y] 153 [-15 * 1, 5 * X] 154 " 155]). 156 157:- comment(delinearize/2, [ 158 summary:"Convert a linear form back to a standard arithmetic expression", 159 amode:delinearize(?,-), 160 args:[ 161 "Linear":"Normalized lists of monomials", 162 "Expression":"Arithmetic expression with constants and variables"], 163 see_also:[linrenorm/2,linearize/3], 164 eg:" 165 ?- linearize(3*X-7*Y+2*(X+Y), L, R), delinearize(L, E). 166 X = X 167 Y = Y 168 R = [] 169 L = [0 * 1, 5 * X, -5 * Y] 170 E = 5 * X - 5 * Y 171 yes. 172 " 173]). 174 175% ---------------------------------------------------------------------- 176 177:-tool(linearize/3,linearize_body/4). 178 179linearize_body(Expr, Lin, Res, Mod) :- 180 linearize(Expr, 0, C, L, [], Res, [], Mod), 181 Lin0 = [C*1|L], 182 linrenorm(Lin0, Lin). 183 184 185% linearize(Expr, CstIn, CstOut, LinOut, LinIn, Res, ResT, Module). 186% Split arithmetic expression Expr into constant, linear and residue part. 187 188:- mode linearize(?, +, -, -, +, -, +, +). 189linearize(X, C, C, [1*X|L], L, R, R, _Mod) :- 190 var(X), !. 191linearize(K, C0, C, L, L, R, R, _Mod) :- 192 number(K), !, 193 ( C0 = 0 -> 194 % Avoid adding things to zero as we may lose the sign of 195 % a -0.0__-0.0 breal in doing so 196 C = K 197 ; 198 ( K = 0 -> 199 C = C0 200 ; 201 C is C0+K 202 ) 203 ). 204linearize(E, C0, C, L, L0, R, R0, Mod) :- E = subscript(Array,Index), !, 205 ( nonground(Index) -> 206 C = C0, L = [1*Aux|L0], R = [Aux=E|R0] 207 ; 208 subscript(Array,Index,Elem), 209 linearize(Elem, C0, C, L, L0, R, R0, Mod) 210 ). 211linearize(lin(Lin), C0, C, L, L0, R, R, _Mod) :- !, 212 Lin = [K*1|L1], 213 C is C0+K, 214 difflistify(L1, L, L0). 215linearize(+E, C0, C, L, L0, R, R0, Mod) :- !, 216 linearize(E, C0, C, L, L0, R, R0, Mod). 217linearize(-E, C0, C, L, L0, R, R0, Mod) :- !, 218 linearize((-1)*E, C0, C, L, L0, R, R0, Mod). 219linearize(E1+E2, C0, C, L, L0, R, R0, Mod) :- !, 220 linearize(E1, C0, C1, L1, L0, R1, R0, Mod), 221 linearize(E2, C1, C, L, L1, R, R1, Mod). 222linearize(sum(List), C0, C, L, L0, R, R0, Mod) :- !, 223 linearize_sum(List, C0, C, L, L0, R, R0, Mod). 224linearize(E1-E2, C0, C, L, L0, R, R0, Mod) :- !, 225 linearize(E1, C0, C1, L1, L0, R1, R0, Mod), 226 linearize((-1)*E2, C1, C, L, L1, R, R1, Mod). 227linearize(E1*E2, C0, C, L, L0, R, R0, Mod) :- !, 228 ( make_prod_sum(E1, E2, ProdSum) -> 229 linearize_sum(ProdSum, C0, C, L, L0, R, R0,Mod) 230 ; 231 linearize(E1, 0, E1C, E1Lin, [], ProdR1, ProdR0, Mod), 232 linearize(E2, 0, E2C, E2Lin, [], ProdR, ProdR1, Mod), 233 EC is E1C*E2C, 234 ( C0 = 0 -> 235 % if the accumulated constant is an integer 0 then 236 % do not add it to the newly calculated sum, doing 237 % this preserves the potential -0.0__-0.0 (negative 238 % zero) result of the multiplication which would 239 % otherwise be lost 240 C is EC 241 ; 242 ( EC = 0 -> 243 C = C0 244 ; 245 C is C0 + EC 246 ) 247 ), 248 lin_mult(E1C, E2Lin, L1, L0), 249 lin_mult(E2C, E1Lin, L2, L1), 250 ( E1Lin==[] -> 251 L=L2, R=ProdR, R0=ProdR0 252 ; E2Lin==[] -> 253 L=L2, R=ProdR, R0=ProdR0 254 ; 255 delinearize(E1Lin,E1LinExpr), 256 delinearize(E2Lin,E2LinExpr), 257 ProdR0 = [], 258 resubstitute_free_nonlinears(ProdR,L2,R0,R1), 259% E1LinExpr = lin([0*1|E1Lin]), 260% E2LinExpr = lin([0*1|E2Lin]), 261 L = [1*Aux|L2], 262 % E1LinExpr*E2LinExpr could be simplified 263 % and possibly a constant factored out. 264 R = [Aux=E1LinExpr*E2LinExpr|R1] 265 ) 266 ). 267linearize(E, C0, C, L, L0, R, R0, Mod) :- E = eval(Expr), !, 268 ( var(Expr) -> 269 C = C0, L = [1*Aux|L0], R = [Aux=E|R0] 270 ; 271 linearize(Expr, C0, C, L, L0, R, R0, Mod) 272 ). 273linearize(E, C, C, [1*Aux|L0], L0, [Aux=E|R0], R0, _Mod) :- 274 nonground(E), !. 275%linearize(E, C0, C, L, L, R, R, Mod) :- 276% (K is E)@Mod, C is C0+K. % ground(E) 277linearize(E, C0, C, L, L0, R, R0, Mod) :- 278 %% ground(E) 279 (subcall((K is E)@Mod,[]) -> 280 %% function returned without delaying goals 281 C is C0+K, 282 L=L0, 283 R=R0 284 ; 285 %% function left delayed goals, so return it to the caller 286 C0=C, 287 L=[1*Aux|L0], 288 R=[Aux=E|R0] 289 ). 290 291 %% Suceeds when the variable V is present in the linear expression 292 var_exists_in([_C*Y | Tail],V):- 293 ( Y==V -> 294 true 295 ; 296 var_exists_in(Tail,V) 297 ). 298 299 %% 300 %% resubstitute([Temp=SubExp|..],+Lin, R0, R) 301 %% 302 %% Reconstructs any non-linear expressions which had sub- 303 %% expressions replaced by temporaries (as given in the first 304 %% argument). The caveat is that these bindings CANNOT be 305 %% replaced when the temporary variable exists in the normailised 306 %% linear expression Lin (as this would break the simple [1* 307 %% Const, Const*Var, Const*Var...] format). 308 %% 309 resubstitute_free_nonlinears([],_Lin,R,R). 310 resubstitute_free_nonlinears([LHS=RHS|More], Lin, R0, R) :- 311 (var_exists_in(Lin,LHS) -> 312 %% var=expr bindings which refer to variables in the 313 %% linear expression must be returned 314 R1=[LHS=RHS|R0] 315 ; 316 %% Only substitue bindings that are NOT present in the 317 %% Linear expression 318 LHS=RHS,R1=R0 319 ), 320 resubstitute_free_nonlinears(More,Lin,R1,R). 321 322 :- mode linearize_sum(?, +, -, -, +, -, +, +). 323 linearize_sum(X, C, C, [1*Aux|L], L, [Aux=sum(X)|R], R, _Mod) :- 324 var(X), !. 325 linearize_sum([], C, C, L, L, R, R, _Mod) :- !. 326 linearize_sum([E|Es], C0, C, L, L0, R, R0, Mod) :- !, 327 linearize(E, C0, C1, L1, L0, R1, R0, Mod), 328 linearize_sum(Es, C1, C, L, L1, R, R1, Mod). 329 linearize_sum(List, C0, C, L, L0, R, R0, Mod) :- List = subscript(Array,Index), !, 330 ( nonvar(Array), ground(Index) -> 331 subscript(Array,Index,Elems), 332 linearize_sum(Elems, C0, C, L, L0, R, R0, Mod) 333 ; 334 C = C0, L = [1*Aux|L0], R = [Aux=sum(List)|R0] 335 ). 336 linearize_sum(X, C, C, L, L, R, R, _Mod) :- 337 error(5, _ is sum(X)). 338 339 340 % lin_mult(+Cst, +LinIn, -LinOut, +LinOutT). 341 lin_mult(K, CXs, CKXs, CKXs0) :- 342 % `not K =\= 0' is like `K =:= 0' except it fails for bounded reals 343 % which span 0 rather than setting up a delayed goal. 344 ( not K =\= 0 -> CKXs=CKXs0 345 ; lin_mult1(K, CXs, CKXs, CKXs0) ). 346 347 :- mode lin_mult1(+,+,-,+). 348 lin_mult1(_, [], CKXs, CKXs). 349 lin_mult1(K, [C*X|CXs], [CK*X|CKXs1], CKXs) :- 350 CK is C*K, 351 lin_mult1(K, CXs, CKXs1, CKXs). 352 353 make_prod_sum([], [], Sum) :- -?-> Sum=[]. 354 make_prod_sum([X|Xs], [Y|Ys], Sum) :- -?-> 355 Sum = [X*Y|XYs], 356 make_prod_sum(Xs, Ys, XYs). 357 358 359% ---------------------------------------------------------------------- 360% (re)normalise a linear term (see above) by collecting constants, 361% collapsing multiple occurrences of variables and eliminating zeros. 362% ---------------------------------------------------------------------- 363 364linrenorm(Lin, LinNorm) :- 365 sort(2, >=, Lin, LinSorted), 366 sepia_kernel:collapse_linear(LinSorted, LinNorm). 367 368 369% ---------------------------------------------------------------------- 370% Normal form -> standard expressions 371% ---------------------------------------------------------------------- 372 373delinearize(Norm, Expr) :- 374 delinearize(Norm, 0, Expr). 375 376delinearize([], Expr, Expr). 377delinearize([F*X|Ts], Expr0, Expr) :- 378 ( Expr0 == 0 -> 379 Fact = F, Expr1 = Term 380 ; sgn(F, -1) -> 381 Fact is -F, Expr1 = Expr0 - Term 382 ; 383 Fact = F, Expr1 = Expr0 + Term 384 ), 385 ( nonvar(X) -> 386 Term is Fact * X 387 ; not Fact =\= 1 -> 388 Term = X 389 ; 390 Term = Fact*X 391 ), 392 delinearize(Ts, Expr1, Expr). 393 394 395% ---------------------------------------------------------------------- 396 397difflistify([], DXs, DXs). 398difflistify([X|Xs], [X|DXs], DXs0) :- 399 difflistify(Xs, DXs, DXs0). 400 401 402 403% ---------------------------------------------------------------------- 404% Normalizer for polynomials 405% 406% A monomial is a list of constants and variables and represents 407% the product of all the list elements. 408% In a normalised monomial, the list is sorted, the first element is 409% the (only) constant and the others are variables. 410% 411% A normalised polynomial is represented as a list of lists of 412% normalised monomials. The sublists represent groups of monomials 413% of the same degree in ascending order. If there are no monomials 414% for a certain degree, the list element is missing: 415% 416% [ConstantMonos, LinearMonos, QuadraticMonos, CubicMonos, ...] 417% 418% In a normalised polynomial, all monomials are normalised and 419% all monomials with identical variables are merged. 420% 421% Non-polynomial components are factored out by introducing an 422% auxiliary variable and adding a term Aux=NonPolyExpr to the 423% NonPoly result list. 424% 425% We use an intermediate form called sum_of_prods which is just a 426% list (sum) of lists (products of numbers and vars) and not further 427% normalised. 428% ---------------------------------------------------------------------- 429 430:- export polynorm/3. 431:- export polyrenorm/2. 432:- export polydenorm/2. 433 434:- comment(polynorm/3, [ 435 summary:"Extracts and normalises the polynomial part of an arithmetic expression", 436 amode:polynorm(?,-,-), 437 args:["Expression":"Arithmetic expression with constants and variables", 438 "NormPoly":"Normalized polynomial", 439 "Residue":"Residual nonpolynomial components in the form AuxVar=Expr"], 440 see_also:[polydenorm/2,polyrenorm/2,linearize/3], 441 desc:html("\ 442 This predicate converts a general arithmetic expression into a 443 normalized polynomial representation and possibly a nonpolynomial 444 residue. The normalized polynomial representation is a follows: 445<P> 446 A <EM>monomial</EM> is a list of constants and variables and represents 447 the product of all the list elements. 448 In a <EM>normalised monomial</EM>, the list is sorted, the first element is 449 the (only) constant and the others are variables. 450<P> 451 A <EM>normalised polynomial</EM> is represented as a list of lists of 452 normalised monomials. The sublists represent groups of monomials 453 of the same degree in ascending order. If there are no monomials 454 for a certain degree, the list element is missing: 455<PRE> 456 [ConstantMonos, LinearMonos, QuadraticMonos, CubicMonos, ...] 457</PRE> 458 In a normalised polynomial, all monomials are normalised and 459 all monomials with identical variables are merged. 460<P> 461 Non-polynomial components are factored out by introducing an 462 auxiliary variable and adding a term Aux=NonPolyExpr to the 463 Residue result list. All variables within Expression (which 464 are free at normalization time) are taken to be numerical 465 variables. If you intend to have variables which can be bound 466 to symbolic expressions rather than number, they must be 467 wrapped into an eval/1 functor. 468 "), 469 eg:" 470 ?- polynorm(2*5 + 3*(X+5*Y+7)*Z, Poly, Res). 471 X = X 472 Y = Y 473 Z = Z 474 Poly = [[[10]], [[21, Z]], [[3, X, Z], [15, Y, Z]]] 475 Res = [] 476 yes. 477 478 ?- polynorm(3*(X+Y), Poly, Res). 479 X = X 480 Y = Y 481 Poly = [[[3, X], [3, Y]]] 482 Res = [] 483 yes. 484 485 ?- polynorm(3, Poly, Res). 486 Poly = [[[3]]] 487 Res = [] 488 yes. 489 " 490]). 491 492:- comment(polyrenorm/2, [ 493 summary:"Renormalize a polynomial form", 494 amode:polyrenorm(+,-), 495 args:["PolyOld":"Possibly denormal polynomial form", 496 "PolyNew":"Normalized polynomial form"], 497 see_also:[polynorm/3,polydenorm/2], 498 desc:html("\ 499 See polynorm/3 for the definition of the polynomial form. 500 Such a form can become denormalized due to unifications 501 (instantiation or variable-variable aliasing). This predicate 502 renormalizes it. 503 "), 504 eg:" 505 ?- polynorm(3*(X+Y), Poly1, []), writeln(Poly1), 506 Y = 3, 507 polyrenorm(Poly1, Poly2), writeln(Poly2). 508 509 [[[3, X], [3, Y]]] 510 [[[9]], [[3, X]]] 511 " 512]). 513 514:- comment(polydenorm/2, [ 515 summary:"Convert a polynomial form back to a standard arithmetic expression", 516 amode:polydenorm(?,-), 517 args:[ 518 "NormPoly":"Normalized polynomial form", 519 "Expression":"Arithmetic expression with constants and variables"], 520 see_also:[polyrenorm/2,polynorm/3], 521 eg:" 522 ?- polynorm(2*5 + 3*(X+5*Y+7)*Z, Poly, []), polydenorm(Poly, Expr). 523 X = X 524 Y = Y 525 Z = Z 526 Poly = [[[10]], [[21, Z]], [[3, X, Z], [15, Y, Z]]] 527 Expr = 10 + 21 * Z + 3 * X * Z + 15 * Y * Z 528 yes. 529 " 530]). 531 532polynorm(Expr, NormPoly, NonPoly) :- 533 sum_of_prods(Expr, SumOfProds, [], NonPoly, []), 534 sum_of_prods_to_poly(SumOfProds, NormPoly). 535 536 537 % renormalise a polynomial that may have become 538 % denormal by instantiation and bindings 539polyrenorm(FormerNormPoly, NormPoly) :- 540 ( 541 foreach(Group, FormerNormPoly), 542 fromto(MessyPoly, MessyPoly1, MessyPoly2, []) 543 do 544 append(Group, MessyPoly2, MessyPoly1) 545 ), 546 sum_of_prods_to_poly(MessyPoly, NormPoly). 547 548% ---------------------------------------------------------------------- 549 550:- export quadnorm/6. 551 552:- comment(quadnorm/6, [ 553 summary:"Extracts constant, linear and quadratic part of an arithmetic expression", 554 amode:quadnorm(?,-,-,-,-,-), 555 args:["Expression":"Arithmetic expression with constants and variables", 556 "Const":"Variable or number", 557 "Linear":"Variable or normalized linear polynomial", 558 "Quadratic":"Variable or normalized quadratic polynomial", 559 "PolyRes":"Variable or normalized superquadratic polynomial", 560 "Residue":"Residual nonpolynomial components in the form AuxVar=Expr"], 561 see_also:[polynorm/3,linearize/3], 562 desc:html("\ 563 This predicate is a simplified interface to polynorm/3 for the case 564 where one is only interested in linear and quadratic components. 565 See polynorm/3 for details. 566 "), 567 eg:" 568 ?- quadnorm(2*5 + 3*(X+5*Y+7)*Z, Const, Lin, Quad, Poly, Res). 569 X = X 570 Y = Y 571 Z = Z 572 Const = 10 573 Lin = [[21, Z]] 574 Quad = [[3, X, Z], [15, Y, Z]] 575 Poly = [] 576 Res = [] 577 yes. 578 " 579]). 580 581quadnorm(Expr, Const, Lin, Quad, RestPoly, NonPoly) :- 582 polynorm(Expr, Poly, NonPoly), 583 extract_const(Poly, Const, Poly1), 584 extract_lin(Poly1, Lin, Poly2), 585 extract_quad(Poly2, Quad, RestPoly). 586 587 extract_const([[[Const]]|RestPoly], Const, RestPoly) :- !. 588 extract_const(RestPoly, 0, RestPoly). 589 590 extract_lin([CXs|RestPoly], CXs, RestPoly) :- CXs = [[_,_]|_], !. 591 extract_lin(RestPoly, [], RestPoly). 592 593 extract_quad([CXYs|RestPoly], CXYs, RestPoly) :- CXYs = [[_,_,_]|_], !. 594 extract_quad(RestPoly, [], RestPoly). 595 596% ---------------------------------------------------------------------- 597 598 % sum_of_prods(?Expr, -SumOfProds, +..., -NonPoly, +...) 599 % Flatten an expression into a (additive) list of 600 % (multiplicative) lists of numbers and variables. 601 % variables and constants are mixed and may be repeated. 602 % No arithmetic operations are performed in here! 603:- mode sum_of_prods(?, -, +, -, +). 604sum_of_prods(X, [[X]|P0], P0, R, R) :- (var(X);number(X)), !. 605sum_of_prods(E1+E2, P, P0, R, R0) :- !, 606 sum_of_prods(E1, P, P1, R, R1), 607 sum_of_prods(E2, P1, P0, R1, R0). 608sum_of_prods(E1-E2, P, P0, R, R0) :- !, 609 sum_of_prods(E1, P, P1, R, R1), 610 sum_of_prods((-1)*E2, P1, P0, R1, R0). 611sum_of_prods(+E, P, P0, R, R0) :- !, 612 sum_of_prods(E, P, P0, R, R0). 613sum_of_prods(-E, P, P0, R, R0) :- !, 614 sum_of_prods(-1*E, P, P0, R, R0). 615sum_of_prods(sum(Es), P, P0, R, R0) :- !, 616 sum_list(Es, P, P0, R, R0). 617sum_of_prods(EX*EY, P, P0, R, R0) :- !, 618 ( list_times_list(EX, EY, P, P0, R, R0) -> 619 true 620 ; 621 sum_of_prods(EX, PX, [], R, R1), 622 sum_of_prods(EY, PY, [], R1, R0), 623 poly_times_poly(PX, PY, P, P0) 624 ). 625sum_of_prods(EX^N, P, P0, R, R0) :- integer(N), !, 626 sum_of_prods(EX, PX, [], R, R0), 627 poly_power(PX, N, P, P0). 628sum_of_prods(lin(Lin), P, P0, R, R0) :- !, 629 R = R0, 630 ( foreach(C*V, Lin), fromto(P, [[C,V]|P1], P1, P0) do true ). 631sum_of_prods(E, P, P0, R, R0) :- E = eval(Expr), !, 632 ( var(Expr) -> 633 P = [[Aux]|P0], R = [Aux=E|R0] 634 ; 635 sum_of_prods(Expr, P, P0, R, R0) 636 ). 637sum_of_prods(E, P, P0, R, R0) :- E = subscript(Array,Index), !, 638 ( nonground(Index) -> 639 P = [[Aux]|P0], R = [Aux=E|R0] 640 ; 641 subscript(Array,Index,Elem), 642 sum_of_prods(Elem, P, P0, R, R0) 643 ). 644sum_of_prods(E, [[Aux]|P0], P0, [Aux=E|R0], R0) :- 645 nonground(E), !. 646sum_of_prods(E, [[C]|P0], P0, R, R) :- 647 % ground(E), 648 C is E. 649 650 poly_times_poly([], _MonoYs, P, P). 651 poly_times_poly([MonoX|MonoXs], MonoYs, P, P0) :- 652 ( 653 foreach(MonoY, MonoYs), 654 fromto(P, [MonoXY|MonoXYs], MonoXYs, P1), 655 param(MonoX) 656 do 657 append(MonoX, MonoY, MonoXY) % MonoX*MonoY 658 ), 659 poly_times_poly(MonoXs, MonoYs, P1, P0). 660 661 sum_list(List, P, P0, R, R0) :- var(List), !, 662 P = [[Aux]|P0], R = [Aux=sum(List)|R0]. 663 sum_list([], P, P0, R, R0) :- !, P=P0, R=R0. 664 sum_list([X|Xs], P, P0, R, R0) :- !, 665 sum_of_prods(X, P, P1, R, R1), 666 sum_list(Xs, P1, P0, R1, R0). 667 sum_list(List, P, P0, R, R0) :- List = subscript(Array,Index), !, 668 ( nonground(Index) -> 669 P = [[Aux]|P0], R = [Aux=sum(List)|R0] 670 ; 671 subscript(Array, Index, Elems), 672 sum_list(Elems, P, P0, R, R0) 673 ). 674 sum_list(X, _P, _P0, _R, _R0) :- 675 error(5, _ is sum(X)). 676 677 list_times_list([], [], P, P0, R, R0) ?- P=P0, R=R0. 678 list_times_list([X|Xs], [Y|Ys], P, P0, R, R0) ?- 679 sum_of_prods(X, PX, [], R, R2), 680 sum_of_prods(Y, PY, [], R2, R1), 681 poly_times_poly(PX, PY, P, P1), 682 list_times_list(Xs, Ys, P1, P0, R1, R0). 683 684 poly_power(Poly, 1, P, P0) :- !, 685 append(Poly, P0, P). 686 poly_power(Poly, N, P, P0) :- 687 N1 is N//2, N2 is (N+1)//2, % N =:= N1+N2 688 poly_power(Poly, N1, P1, []), 689 poly_power(Poly, N2, P2, []), 690 poly_times_poly(P1, P2, P, P0). 691 692%---------------------------------------------------------------------- 693 694sum_of_prods_to_poly(MessyPoly, NormPoly) :- 695 ( 696 foreach(Xs,MessyPoly), 697 foreach([Cnew|VarXs],NormMonos) 698 do 699 sort(0, >=, Xs, XsSorted), % constants first, then variables 700 separate_constants(XsSorted, 1, Cnew, VarXs) 701 ), 702 split_by_degree(NormMonos, NonUniqueMonosByDegree), 703 ( foreach(NonUniqueMonosSameDegree,NonUniqueMonosByDegree), 704 foreach(NormMonosSameDegree, NormPoly) 705 do 706 sort(2, =<, NonUniqueMonosSameDegree, SortedMonos), 707 collapse_monos(SortedMonos, NormMonosSameDegree) 708 ). 709 710 711% separate_constants(+XsSorted, +Cin, -Cout, -VarXs) 712% XsSorted is a list with numbers at the beginning, 713% which are all multiplied and returned as Cout, 714% the rest list is returned as VarXs 715separate_constants([], C, C, []). 716separate_constants(XsSorted, C, Cnew, VarXs) :- 717 XsSorted = [X|Xs], 718 ( number(X) -> 719 C1 is C*X, 720 separate_constants(Xs, C1, Cnew, VarXs) 721 ; 722 Cnew = C, VarXs = XsSorted 723 ). 724 725degree_sort(Ms, SMs) :- 726 ( foreach(M,Ms), foreach(N-M,NMs) do length(M,N) ), 727 keysort(NMs, SNMs), 728 ( foreach(M,SMs), foreach(_-M,SNMs) do true ). 729 730split_by_degree(Ms, SMs) :- 731 ( foreach(M,Ms), foreach(N-M,NMs) do length(M,N) ), 732 keysort(NMs, SNMs), 733 group_same_key_values(SNMs, GSNMs), 734 ( foreach(_-SM,GSNMs), foreach(SM,SMs) do true ). 735% strip_key_and_fill_gaps(GSNMs, SMs, 1). 736 737group_same_key_values([], []). 738group_same_key_values([K-V|List], [K-[V|KVs]|GroupedList]) :- 739 group_same_key_values(List, K, KVs, GroupedList). 740 741 group_same_key_values([], _, [], []). 742 group_same_key_values([K-V|List], K, [V|KVs], GroupedList) :- !, 743 group_same_key_values(List, K, KVs, GroupedList). 744 group_same_key_values([K-V|List], _K, [], [K-[V|KVs]|GroupedList]) :- 745 group_same_key_values(List, K, KVs, GroupedList). 746 747 strip_key_and_fill_gaps([], [], _N). 748 strip_key_and_fill_gaps(KVsKVs, SVs, N) :- 749 KVsKVs = [K-Vs|KVs], 750 N1 is N+1, 751 ( K > N -> 752 SVs = [[]|SVs0], % missing degree - insert a filler 753 strip_key_and_fill_gaps(KVsKVs, SVs0, N1) 754 ; 755 SVs = [Vs|SVs0], 756 strip_key_and_fill_gaps(KVs, SVs0, N1) 757 ). 758 759 760collapse_monos([], []). 761collapse_monos([Mono|Monos], NormPoly) :- 762 collapse_monos(Mono, Monos, NormPoly). 763 764 collapse_monos(M, [], [M]). 765 collapse_monos(Mono1, [Mono2|Monos], NormPoly) :- 766 Mono2 = [C2|X2s], 767 Mono1 = [C1|X1s], 768 ( X1s == X2s -> 769 C is C1+C2, 770 collapse_monos([C|X1s], Monos, NormPoly) 771 ; 772 NormPoly = [Mono1|NormPoly1], 773 collapse_monos(Mono2, Monos, NormPoly1) 774 ). 775 776 777%---------------------------------------------------------------------- 778 779polydenorm([], 0). 780polydenorm([[]|Degrees], Expr) :- !, 781 polydenorm(Degrees, Expr). 782polydenorm([[[C0|Xs0]|RestDegree0]|Degrees], Expr) :- 783 list_to_times(Xs0, C0, Prod0), 784 ( 785 foreach(Degree, [RestDegree0|Degrees]), 786 fromto(Prod0, Expr1, Expr4, Expr) 787 do 788 ( 789 foreach([C|Xs], Degree), 790 fromto(Expr1, Expr2, Expr3, Expr4) 791 do 792 % The following test avoids a potential delayed goal if C 793 % is a bounded real. 794 ( not C > 0 -> 795 CNeg is -C, 796 list_to_times(Xs, CNeg, Prod), 797 Expr3 = (Expr2 - Prod) 798 ; 799 list_to_times(Xs, C, Prod), 800 Expr3 = (Expr2 + Prod) 801 ) 802 ) 803 ). 804 805 806 list_to_times([], X0, X0). 807 list_to_times([X|Xs], X0, Expr) :- 808 ( X0 == 1 -> 809 list_to_times(Xs, X, Expr) 810 ; 811 list_to_times(Xs, X0*X, Expr) 812 ). 813 814 815 816% ---------------------------------------------------------------------- 817% Normalizer for arbitrary arithmetic expressions 818% 819% Given an expression, returns the generalised normalisation which is defined 820% to be a record with the following fields 821% 822% Note: There are many possible 'normal' forms of an expression, this function 823% returns just one of them. 824% ---------------------------------------------------------------------- 825:-export normalize/3. 826normalize( _Expr, _Lin, _Constraints ) :- 827 true. 828%---------------------------------------------------------------------- 829% Simplifies an expression, by reducing the number of variable references 830%---------------------------------------------------------------------- 831:-export simplify/2. 832simplify(E,E):- 833 var(E),!. 834simplify(C0-C1,C2):- 835 number(C0),number(C1),C2 is C0-C1,!. 836simplify(C0/C1,C2):- 837 number(C0),number(C1),C2 is C0/C1,!. 838 839% Addition 840simplify(C0+C1,C2):- 841 number(C0),number(C1),C2 is C0+C1,!. 842simplify(C+V,lin([C*1,1*V])):- 843 number(C),var(V),!. 844simplify(V+C,lin([C*1,1*V])):- 845 number(C),var(V),!. 846simplify(V0+V1,lin([0*1,1*V0,1*V1])):- 847 var(V0),var(V1),!. 848simplify(C+lin([K0*1|Tail]),lin([K1*1|Tail])):- 849 number(C),K1 is K0 + C,!. 850simplify(V+lin([K0*1|Tail1]),lin([K0*1|Tail2])):- 851 var(V), 852 simplify_merge([1*V],Tail1,Tail2), 853 !. 854simplify(lin([K0*1|Tail])+C,lin([K1*1|Tail])):- 855 number(C),K1 is K0 + C,!. 856simplify(lin([K0*1|Tail1])+V,lin([K0*1|Tail2])):- 857 var(V), 858 simplify_merge([1*V],Tail1,Tail2), 859 !. 860simplify(lin([K0*1|Tail0])+lin([K1*1|Tail1]),lin([K2*1|Tail2])):- 861 K2 is K0 + K1, 862 simplify_merge(Tail0,Tail1,Tail2), 863 !. 864simplify(E1+E2,ExprOut):- 865 simplify(E1,E3), 866 simplify(E2,E4), 867 simplify(E3+E4,ExprOut),!. 868 869% multiplication 870simplify(C0*C1,C2):- 871 number(C0),number(C1),C2 is C0*C1,!. 872simplify(C*V,lin([0*1,C*V])):- 873 number(C),var(V),!. 874simplify(V*C,lin([0*1,C*V])):- 875 number(C),var(V),!. 876simplify(V0*V1,poly([V0,V1])):- 877 var(V0),var(V1),!. 878simplify(C*lin(Lin0), lin(Lin1)):- 879 number(C), 880 (foreach(C0*V0,Lin0), param(C), fromto(Lin1,Out,In,[]) do 881 C1 is C0*C, 882 Out = [C1*V0|In] 883 ),!. 884simplify(E1*E2,ExprOut):- 885 simplify(E1,E3), 886 writeq(E1),write('simp-to->'),writeq(E3),nl, 887 simplify(E2,E4), 888 writeq(E2),write('simp-to->'),writeq(E4),nl, 889 simplify(E3*E4,ExprOut), 890 writeq(E3*E4),write('simp-to->'),writeq(ExprOut),nl, 891 !. 892 893% catch all 894simplify(E, E):-!. 895 896 897% merge lists of ordered terms together maintaining the order 898simplify_merge([],Acc,Acc):-!. 899simplify_merge(Acc,[],Acc):-!. 900simplify_merge([C0*V0|T0],[C1*V1|T1],[C2*V2|Acc0]):- 901 (V0 == V1 -> 902 C2 is C0+C1, 903 V2 = V0, 904 simplify_merge(T0,T1,Acc0) 905 ; 906 (V0 @< V1 -> 907 C2 is C0, 908 V2 = V0, 909 simplify_merge(T0,[C1*V1|T1],Acc0) 910 ; 911 C2 is C1, 912 V2 = V1, 913 simplify_merge([C0*V0|T0],T1,Acc0) 914 ) 915 ). 916 917%---------------------------------------------------------------------- 918% Test code 919%---------------------------------------------------------------------- 920 921 922%ex(3). 923%ex(X). 924%ex(3*X). 925%ex(X*3). 926%ex(4*3). 927%ex(3*X+4*Y). 928%ex(3*(X+Y)). 929 930%ex(3*(X+5*Y)). 931/* 932ex(3*(X+5*Y+7)*Z). 933ex(2*5 + 3*(X+5*Y+7)*Z). 934ex(2*5 + 3*(X+5*(Z+7)/(3*X + 4*Y))*Z). 935 936test :- 937 ex(Expr0), 938 nl, 939 nl, 940 writeq(Expr0),nl, 941 simplify(Expr0,Expr), 942 write('- simplifies to ->'), write(Expr), nl, 943 %polynorm(Expr, Norm, _NonPoly), 944 %writeq(Norm), nl, 945 %polydenorm(Norm, NormExpr), 946 %writeq(NormExpr), nl, 947 fail. 948*/ 949%:-test. 950 951