share/algebra/solver/EnglishSource/SOLVER.TEX

   1 \begin{verbatim}\r
   2 /******************************************************************************/
   3 /*                                                                            */
   4 /* SOLVER - THE NEXT GENERATION                                               */
   5 /*                                                                            */
   6 /* Copyright (C) 2000 : Eckhard Hennig, Ralf Sommer                           */
   7 /* This library is free software; you can redistribute it and/or modify it    */
   8 /* under the terms of the GNU Library General Public License as published     */
   9 /* by the Free Software Foundation; either version 2 of the License, or (at   */
  10 /* your option) any later version.                                            */
  11 /*                                                                            */
  12 /* This library is distributed in the hope that it will be useful, but        */
  13 /* WITHOUT any WARRANTY; without even the implied warranty of                 */
  14 /* MERCHANTABILITY or FITNESS for A PARTICULAR PURPOSE.  See the GNU          */
  15 /* Library General Public License for more details.                           */
  16 /*                                                                            */
  17 /* You should have received a copy of the GNU Library General Public          */
  18 /* License along with this library; if not, write to the Free Software        */
  19 /* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA                  */
  20 /******************************************************************************/
  21 /* Credits: This program is based on many ideas from Henning Trispel's work   */
  22 /*          on the EASY-Solver in 1991.                                       */
  23 /*                                                                            */
  24 /******************************************************************************/
  25 /* Author(s)    : Eckhard Hennig, Ralf Sommer                                 */
  26 /* Project start: 19.01.1994                                                  */
  27 /* Completed    : 16.07.1994                                                  */
  28 /* last change  : 29.06.1995                                                  */
  29 /* Time         : 09:26                                                       */
  30 /******************************************************************************/
  31 /* Changes      : ||||| ||||| ||||| |||                                       */
  32 /******************************************************************************/
  33 /* Modified by  : Dan Stanger dan.stanger@ieee.org to work under maxima       */
  34 /******************************************************************************/
  35 load( "solver/linsolve.mac" )$
  36 load( "solver/slvrtbox.mac" )$
  37 load( "solver/slvrmsgs.mac" )$
  38 load( "solver/misc.mac" )$
  39
  40 put( 'SOLVER, 1, 'Version )$
  41
  42 SetVersion(
  43   /* KEY    = */ 'SOLVER,
  44   'MODULE      = "SOLVER",
  45   'DESCRIPTION = "Symbolic solver for parametric systems of equations.",
  46   'AUTHORS     = "Eckhard Hennig, Ralf Sommer",
  47   'DATE        = "19.01.1994",
  48   'LASTCHANGE  = "29.06.1995",
  49   'TIME        = "09:26",
  50   'PLAN        = "Add a-priori transforms"
  51 )$
  52
  53
  54 /******************************************************************************/
  55 /* last change: 29.06.1995                                                    */
  56 /* Time       : 09:26                                                         */
  57 /* By         : Eckhard Hennig                                                */
  58 /* Description: Option variable Solver_Valuate_All_Nonlin_Vars added.         */
  59 /*              Bug in ValuationSolver removed: rhs variables of solutions    */
  60 /*              were not added to the list of variables.                      */
  61 /******************************************************************************/
  62 /* last change: 28.05.1995                                                    */
  63 /* Time       : 11:59                                                         */
  64 /* By         : Eckhard Hennig                                                */
  65 /* Description: Solver break test added.                                      */
  66 /******************************************************************************/
  67 /* last change: 18.05.1995                                                    */
  68 /* Time       : 16:10                                                         */
  69 /* By         : Eckhard Hennig                                                */
  70 /* Description: Name conflict with AI keyword PARAMS removed.                 */
  71 /******************************************************************************/
  72 /* last change: 30.01.1995                                                    */
  73 /* Time       : 21.42                                                         */
  74 /* By         : Eckhard Hennig                                                */
  75 /* Description: Removal of map( 'num, ... ) has revealed some side effects on */
  76 /*              the list of equations in ValuationSolver. Bug repaired by     */
  77 /*              inserting two additional copylist calls.                      */
  78 /*              Arguments to both load commands above now written in lower-   */
  79 /*              case letters to avoid problems with UNIX versions.            */
  80 /******************************************************************************/
  81 /* last change: 24.01.1995                                                    */
  82 /* Time       : 16.33                                                         */
  83 /* By         : Eckhard Hennig                                                */
  84 /* Description: mapping 'num to expressions caused some invalid solutions,    */
  85 /*              therefore removed. Version property added.                    */
  86 /******************************************************************************/
  87 /* last change: 19.01.1995                                                    */
  88 /* Time       : 13.04                                                         */
  89 /* By         : Eckhard Hennig, Ralf Sommer                                   */
  90 /* Description: Bug in ImmediateAssignments corrected, Function AppendImmed   */
  91 /*              removed.                                                      */
  92 /******************************************************************************/
  93 /* last change: 19.01.1995                                                    */
  94 /* Time       : 11.28                                                         */
  95 /* By         : Eckhard Hennig, Ralf Sommer                                   */
  96 /* Description: SLVRTBOX and SLVRMSGS now autoloading.                        */
  97 /******************************************************************************/
  98 /* last change: 17.01.1995                                                    */
  99 /* Time       : 20.15                                                         */
 100 /* By         : Eckhard Hennig, Ralf Sommer                                   */
 101 /* Description: Option variables modified (underscores inserted according to  */
 102 /*              R. Petti's suggestions)                                       */
 103 /*              Function SetValuation added.                                  */
 104 /******************************************************************************/
 105 /* last change: 25.10.1994                                                    */
 106 /* Time       : 12.46                                                         */
 107 /* By         : Eckhard Hennig                                                */
 108 /* Description: errorMsg renamed to ErrMsg.                                   */
 109 /******************************************************************************/
 110 /* last change: 12.09.1994                                                    */
 111 /* Time       : 11.50                                                         */
 112 /* By         : Eckhard Hennig                                                */
 113 /* Description: errcatch wrapped around final fullratsimp.                    */
 114 /******************************************************************************/
 115 /* last change: 05.09.1994                                                    */
 116 /* Time       : 15.59                                                         */
 117 /* By         : Eckhard Hennig                                                */
 118 /* Description: mode_identity's inserted.                                     */
 119 /******************************************************************************/
 120 /* last change: 05.09.1994                                                    */
 121 /* Time       : 15.21                                                         */
 122 /* By         : Eckhard Hennig                                                */
 123 /* Description: Number of linear equations is now correctly displayed.        */
 124 /*              Bug in DumpToFile corrected.                                  */
 125 /*              Numbers of solution sets are now printed by the postprocessor.*/
 126 /******************************************************************************/
 127 /* last change: 30.08.1994                                                    */
 128 /* Time       : 13.45                                                         */
 129 /* By         : Eckhard Hennig                                                */
 130 /* Description: Single inconsistent solution paths are no longer returned.    */
 131 /******************************************************************************/
 132 /* last change: 29.08.1994                                                    */
 133 /* Time       : 17.16                                                         */
 134 /* By         : Eckhard Hennig                                                */
 135 /* Description: Slight modification of Immediate Assignment Solver.           */
 136 /*              Valuation property for sqrt added.                            */
 137 /*              New option variable SolverDefaultValuation.                   */
 138 /*              Change in Valuation. Now making use of the above option var.  */
 139 /******************************************************************************/
 140 /* last change: 26.08.1994                                                    */
 141 /* Time       : 15.29                                                         */
 142 /* By         : Eckhard Hennig                                                */
 143 /* Description: Bug in Linear Solver removed.                                 */
 144 /******************************************************************************/
 145 /* last change: 25.08.1994                                                    */
 146 /* Time       : 12.41                                                         */
 147 /* By         : Eckhard Hennig                                                */
 148 /* Description: New option variable SolverRatSimpSols.                        */
 149 /*              Linear Solver now calls the consistency check.                */
 150 /*              Capabilities of the transforms in ValuationSolver strongly    */
 151 /*              enhanced.                                                     */
 152 /******************************************************************************/
 153 /* last change: 24.08.1994                                                    */
 154 /* Time       : 23.14                                                         */
 155 /* By         : Eckhard Hennig                                                */
 156 /* Description: Bug fixed in postprocessor: compound expressions are now cor- */
 157 /*              rectly evaluated.                                             */
 158 /*              Linear Solver now checks for remaining equations & variables, */
 159 /*              and removes those "linear" eqs and vars which do not really   */
 160 /*              belong to the "true" linear subsystem.                        */
 161 /******************************************************************************/
 162 /* last change: 11.08.1994                                                    */
 163 /* Time       : 19.09                                                         */
 164 /* By         : Eckhard Hennig                                                */
 165 /* Description: Variable SolverDelEq2VarPref replaced by function variable    */
 166 /*              SolverDelEq. Heuristic algorithm for linear equation          */
 167 /*              extraction improved.                                          */
 168 /******************************************************************************/
 169 /* last change: 19.07.1994                                                    */
 170 /* Time       : 18.01                                                         */
 171 /* By         : Eckhard Hennig                                                */
 172 /* Description: Bug in post processor corrected.                              */
 173 /******************************************************************************/
 174
 175
 176 /******************************************************************************/
 177 /* Global variables for Solver                                                */
 178 /******************************************************************************/
 179
 180 /******************************************************************************/
 181 /* If Solver_Immed_Assign is true then the Solver searches the equations      */
 182 /* for immediate assignments of the form variable = constant and immediately  */
 183 /* inserts these constraints into the remaining equations.                    */
 184 /******************************************************************************/
 185
 186 define_variable( Solver_Immed_Assign, true, boolean )$
 187
 188
 189 /******************************************************************************/
 190 /* Solver_Repeat_Immed controls whether the search for immediate assignments  */
 191 /* is performed repeatedly until no more of them are found.                   */
 192 /******************************************************************************/
 193
 194 define_variable( Solver_Repeat_Immed, true, boolean )$
 195
 196
 197 /******************************************************************************/
 198 /* Solver_Subst_Powers controls whether the Solver substitutes powers of a    */
 199 /* variable by new symbols in one of the following cases:                     */
 200 /* 1. var^n appears raised to exactly one power n                             */
 201 /* 2. for all var^m in the equations : m=n*k , k integer                      */
 202 /******************************************************************************/
 203
 204 define_variable( Solver_Subst_Powers, false, boolean )$
 205
 206
 207 /******************************************************************************/
 208 /* Solver_Incons_Params controls whether the Solver terminates when a non-    */
 209 /* trivial equation containing only parameters is encountered. for example,   */
 210 /* if A and B are defined as parameters and the solution process yields       */
 211 /* A = B^2  then the Solver stops if Solver_Incons_Params is = 'BREAK. If set */
 212 /* to 'ASK the Solver asks whether  A - B^2  is zero and continues if it is.  */
 213 /* If set to 'IGNORE the Solver quietly assumes that the expression is zero   */
 214 /* if it does not directly contradict with any of the assumptions made before.*/
 215 /******************************************************************************/
 216
 217 define_variable( Solver_Incons_Params, 'ASK, any_check )$
 218
 219 put(
 220   'Solver_Incons_Params,
 221   lambda( [ x ],
 222     if not member( x, [ 'ASK, 'BREAK, 'IGNORE ] ) then
 223       ErrorHandler("InvIncPar", x, 'Fatal )
 224   ),
 225   'value_check
 226 )$
 227
 228
 229 /******************************************************************************/
 230 /* Solver_Linear controls whether the Solver tries to find and solve blocks   */
 231 /* of linear equations before submitting the remaining equations to the       */
 232 /* heuristic valuation solver. This is useful if the equations to be solved   */
 233 /* are known to contain a large linear part.                                  */
 234 /******************************************************************************/
 235
 236 define_variable( Solver_Linear, true, boolean )$
 237
 238
 239 /******************************************************************************/
 240 /* If Solver_Repeat_Linear is true then the LinearSolver will be called       */
 241 /* repeatedly until no more linear equations are found. If false then         */
 242 /* LinearSolver will be called only once.                                     */
 243 /******************************************************************************/
 244
 245 define_variable( Solver_Repeat_Linear, true, boolean )$
 246
 247
 248 /******************************************************************************/
 249 /* If Solver_Find_All_Linear_Vars is true then LinearSolver will try to find  */
 250 /* linear equations with respect to all available variables and solve these   */
 251 /* equations simultaneously. If false then LinearSolver will only search for  */
 252 /* linear equations with respect to the variables passed over in the function */
 253 /* call.                                                                      */
 254 /******************************************************************************/
 255
 256 define_variable( Solver_Find_All_Linear_Vars, true, boolean )$
 257
 258
 259 /******************************************************************************/
 260 /* Solver_Assumptions contains constraints on the parameters which should be  */
 261 /* checked by the user after the termination of Solver. These constraints     */
 262 /* result from the parameter consistency check the behavior of which is       */
 263 /* controlled by the setting of the option variable Solver_Incons_Params.     */
 264 /* Any numerical solution of the equations obtained by assigning numerical    */
 265 /* values to symbolic parameters should be checked for consistency with all   */
 266 /* expressions in Solver_Assumptions.                                         */
 267 /******************************************************************************/
 268
 269 define_variable( Solver_Assumptions, [], list )$
 270
 271
 272 /******************************************************************************/
 273 /* Solver_Del_Eq holds the name of a function which controls the behavior of  */
 274 /* the heuristic search algorithm which extract linear equations from the     */
 275 /* entire system of equations.                                                */
 276 /******************************************************************************/
 277
 278 define_variable( Solver_Del_Eq, 'MakeSquareLinearBlocks, any )$
 279
 280
 281 /******************************************************************************/
 282 /* If Solver_Valuate_All_Nonlin_Vars is true then the ValuationSolver will    */
 283 /* valuate the equations w.r.t. all remaining variables and not only w.r.t    */
 284 /* variables which are being searched for at the current step.                */
 285 /******************************************************************************/
 286
 287 define_variable( Solver_Valuate_All_Nonlin_Vars, false, boolean )$
 288
 289
 290 /******************************************************************************/
 291 /* Solver_Valuation_Strategy holds the name of the equation valuation         */
 292 /* strategy called by the valuation solver to determine the order by which    */
 293 /* the equations are to be solved.                                            */
 294 /******************************************************************************/
 295
 296 define_variable( Solver_Valuation_Strategy, 'MinVarPathsFirst, any_check )$
 297
 298 put(
 299   'Solver_Valuation_Strategy,
 300   lambda(
 301     [ x ],
 302     if not FunctionP( x ) then
 303       ErrorHandler( "UndefStrat", x, 'Fatal )
 304   ),
 305   'value_check
 306 )$
 307
 308
 309 /******************************************************************************/
 310 /* Solver_Default_Valuation contains the default valuation for arithmetic     */
 311 /* operators. Whenever an operator is encountered for which no valuation has  */
 312 /* been defined by SetProp( <operator>, 'Valuation, <valuation> ) this value  */
 313 /* is taken for the formula complexity calculations.                          */
 314 /******************************************************************************/
 315
 316 define_variable( Solver_Default_Valuation, 10, fixnum )$
 317
 318
 319 /******************************************************************************/
 320 /* Solver_Max_Len_Val_Order limits the length of the list of candidates for   */
 321 /* the solve calls in ValuationSolver. A low value will usually increase the  */
 322 /* efficiency of the valuation solver since, in general, the first or second  */
 323 /* attempt to solve an equation (hopefully) succeeds.                         */
 324 /******************************************************************************/
 325
 326 define_variable( Solver_Max_Len_Val_Order, 5, fixnum )$
 327
 328
 329 /******************************************************************************/
 330 /* Solver_Transforms is a list containing the names of functions which can be */
 331 /* applied to an equation after a failed solve call. These functions must     */
 332 /* take three arguments: the equation to be transformed, the variable to be   */
 333 /* solved for, and a list of (probably implicit) solutions the Solver has     */
 334 /* already found for the equation.                                            */
 335 /******************************************************************************/
 336
 337 define_variable( Solver_Transforms, [], list )$
 338
 339
 340 /******************************************************************************/
 341 /* If Solver_Postprocess is set to false no postprocessing of the results will*/
 342 /* be done. Instead, the solutions are displayed in the internal hierarchical */
 343 /* list format. Useful for debugging purposes.                                */
 344 /******************************************************************************/
 345
 346 define_variable( Solver_Postprocess, true, boolean )$
 347
 348
 349 /******************************************************************************/
 350 /* Solver_Backsubst controls the output format of Solver. If Solver_Backsubst */
 351 /* is true then the result will be displayed with fully evaluated right-hand  */
 352 /* sides for each variable. If the option variable is set to false then the   */
 353 /* right-hand sides of the solutions may still contain references to some     */
 354 /* of the other variables which have been solved for.                         */
 355 /******************************************************************************/
 356
 357 define_variable( Solver_Backsubst, true, boolean )$
 358
 359
 360 /******************************************************************************/
 361 /* With Solver_Disp_All_Sols set to true all solutions will be displayed      */
 362 /* including those for the variables which have been solved for in the        */
 363 /* solution process but have not been explicitly asked for.                   */
 364 /******************************************************************************/
 365
 366 define_variable( Solver_Disp_All_Sols, false, boolean )$
 367
 368
 369 /******************************************************************************/
 370 /* If Solver_RatSimp_Sols is true then the Solver Postprocessor will          */
 371 /* fullratsimp the solutions before returning them.                           */
 372 /******************************************************************************/
 373
 374 define_variable( Solver_RatSimp_Sols, true, boolean )$
 375
 376
 377 /******************************************************************************/
 378 /* If Solver_Dump_To_File is true then the ValuationSolver writes the         */
 379 /* solutions and yet unsolved equations to a file after each iteration. This  */
 380 /* might help if the Solver crashes.                                          */
 381 /******************************************************************************/
 382
 383 define_variable( Solver_Dump_To_File, false, boolean )$
 384
 385
 386 /******************************************************************************/
 387 /* Solver_Dump_File contains the name of the file to which the dump is        */
 388 /* written.                                                                   */
 389 /******************************************************************************/
 390
 391 define_variable( Solver_Dump_File, "SOLVER.DMP", any )$
 392
 393
 394 /******************************************************************************/
 395 /* Solver_Break_Test holds the name of a function which is called immediately */
 396 /* before attempting to solve an equation. This function then decides whether */
 397 /* Solver should try to solve the equation or whether it should stop, e.g.    */
 398 /* because the problem has become too complex. The arguments passed to the    */
 399 /* Solver_Break_Test are 1. the equation, 2. the variable, 3. the valuation.  */
 400 /* The Solver stops if the Solver_Break_Test returns true, it continues if    */
 401 /* the return value is false.                                                 */
 402 /******************************************************************************/
 403
 404 define_variable( Solver_Break_Test, 'SolverJustDoIt, any )$
 405
 406
 407 /******************************************************************************/
 408 /* Solver, main program                                                       */
 409 /******************************************************************************/
 410
 411 Solver( Equations, [ SolverParams ] ) := (
 412
 413   mode_declare(
 414     [ Equations, SolverParams ], list
 415   ),
 416
 417   block(
 418
 419     [
 420       Variables,       /* List of variables to be solved for            */
 421       UserVars,        /* List of variables specified by user           */
 422       Parameters,      /* List of symbols to be used as parameters      */
 423       Expressions,     /* List of compound expressions to be solved for */
 424       PowerSubst,      /* List of substituted symbols for powers        */
 425
 426       Solutions,       /* List of solutions found by Solver             */
 427       RemainingEqs,    /* List of remaining equations                   */
 428
 429       /* Assumptions made in linear solver should be local  */
 430       Solver_Assumptions : [],
 431
 432       Active
 433     ],
 434
 435     mode_declare(
 436       [
 437         Variables, UserVars, Parameters, Expressions, PowerSubst,
 438         Solutions
 439       ], list,
 440       Active, boolean
 441     ),
 442
 443     /* No assumptions to start with */ErrorHandlerSolver_Assumptions : [],
 444
 445     /* Initialize list of solutions */
 446     Solutions : [],
 447
 448     /* Do all necessary preprocessing */
 449     map(
 450       lambda([x,y],x::y),
 451       [
 452         'Equations, 'SolverParams, 'Variables, 'Parameters,
 453         'Expressions, 'PowerSubst, 'UserVars
 454       ],
 455       SetupSolver( Equations, SolverParams )
 456     ),
 457
 458     if MsgLevel = 'DEBUG then
 459       display(
 460         Equations, SolverParams, Variables, Parameters, Expressions,
 461         PowerSubst, UserVars, Solutions
 462       ),
 463
 464     block(
 465       [],
 466
 467       /* Search for and apply immediate assignments */
 468
 469       Active : Solver_Immed_Assign,
 470       while Active do (
 471         map(
 472           lambda([x,y],x::y),
 473           [ 'Active, 'Solutions, 'Equations, 'Variables ],
 474           ImmediateAssignments( Solutions, Equations, Variables, Parameters )
 475         ),
 476
 477         Active : Active and Solver_Repeat_Immed,
 478
 479         if MsgLevel = 'DEBUG then
 480           display(
 481             Equations, SolverParams, Variables, Parameters, Expressions,
 482             PowerSubst, UserVars, Solutions
 483           )
 484       ),
 485
 486       if Empty( Equations ) or Empty( Variables ) then
 487         return( false ),
 488
 489       Equations : map(
 490         lambda(
 491           [ Eq ],
 492           fullratsimp( lhs( Eq ) - rhs( Eq ) )
 493         ),
 494         Equations
 495       ),
 496
 497       /* Find and solve linear equations */
 498
 499       Active : Solver_Linear,
 500       while Active do (
 501         map(
 502           lambda([x,y],x::y),
 503           [ 'Active, 'Solutions, 'Equations, 'Variables ],
 504           LinearSolver( Solutions, Equations, Variables, Parameters )
 505         ),
 506
 507         Active : Active and Solver_Repeat_Linear,
 508
 509         if MsgLevel = 'DEBUG then
 510           display(
 511             Equations, SolverParams, Variables, Parameters, Expressions,
 512             PowerSubst, UserVars, Solutions
 513           )
 514       ),
 515
 516       if Empty( Equations ) or Empty( Variables ) then
 517         return( false ),
 518
 519
 520       if Solver_Valuate_All_Nonlin_Vars then
 521         Variables : Union(
 522           SetDifference( listofvars( Equations ), Parameters ),
 523           Variables
 524         ),
 525
 526       /* apply valuation strategies to solve the nonlinear equations. */
 527       map(
 528         lambda([x,y],x::y),
 529         [ 'Active, 'Solutions, 'Equations, 'Variables ],
 530         ValuationSolver( Solutions, Equations, Variables, Parameters )
 531       )
 532
 533     ), /* END block */
 534
 535
 536     /* Return the solutions and the unsolved equations */
 537
 538     if Solver_Postprocess then
 539       return( PostProcess( Solutions, UserVars, Expressions, PowerSubst ) )
 540     else
 541       return( Solutions )
 542   )
 543 )$
 544
 545
 546 /******************************************************************************/
 547 /* TerminateSolver terminates the Solver.                                     */
 548 /******************************************************************************/
 549
 550 TerminateSolver() := error( ErrMsg["SolvrTerm"] )$
 551
 552
 553 /******************************************************************************/
 554 /* SetupSolver preprocesses the equations and optional parameters before      */
 555 /* submitting them to the Solver. The equations are checked if they are       */
 556 /* really equations, and all equations of the form  NUMBER = NUMBER  or       */
 557 /* f( PARAMETERS ) = g( PARAMETERS )  are checked for consistency and then    */
 558 /* dropped.                                                                   */
 559 /******************************************************************************/
 560
 561 SetupSolver( Equations, SolverParams ) := (
 562
 563   mode_declare(
 564     [ Equations, SolverParams ], list
 565   ),
 566
 567   block(
 568
 569     [
 570        i, AllVars, Var, SubstSym, Power,
 571        Variables, Parameters, Expressions,
 572        PowerSubst, UserVars
 573     ],
 574
 575     mode_declare(
 576       i, fixnum,
 577       [
 578         AllVars, Power, Variables, Parameters, Expressions,
 579         PowerSubst, UserVars
 580       ], list,
 581       [ Var, SubstSym ], any
 582     ),
 583
 584
 585     Expressions : [],
 586     PowerSubst  : [],
 587     Parameters  : [],
 588
 589
 590     /* Make sure that Equations is a list. Abort if it is not. */
 591
 592     if not listp( Equations ) then
 593       ErrorHandler( "EqsNotLst", Equations, 'Fatal ),
 594
 595     /* delete all entries from Equations which are not equations */
 596
 597     Equations : sublist( Equations, 'EquationP ),
 598
 599     /* Convert all floating-point numbers to rational numbers. If this was */
 600     /* not done then rounding errors could fool the consistency check.     */
 601
 602     Equations : map( 'rat, Equations ),
 603
 604     /* Process the optional arguments to Solver */
 605
 606     if not Empty( SolverParams ) then (
 607
 608       Variables : SolverParams[1],
 609
 610       /* Make sure that Variables is a list. Abort if it is not. */
 611
 612       if not listp( Variables ) then
 613         ErrorHandler( "VarNotLst", Variables, 'Fatal ),
 614
 615       /* delete multiple occurrences of identical symbols */
 616
 617       Variables : Setify( Variables ),
 618
 619       PrintMsg( 'DETAIL, SolverMsg["VarsAre"], Variables ),
 620
 621       if length( SolverParams ) > 1 then (
 622
 623         Parameters : SolverParams[2],
 624
 625         /* Make sure that Parameters is a list. Abort if it is not. */
 626
 627         if not listp( Parameters ) then
 628           ErrorHandler( "ParNotLst", Parameters, 'Fatal ),
 629
 630         /* delete multiple occurrences of identical symbols */
 631
 632         Parameters : Setify( Parameters ),
 633
 634         PrintMsg( 'DETAIL, SolverMsg["ParsAre"], Parameters )
 635       )
 636
 637     ),
 638
 639
 640     /* Check if Variables and Parameters are disjoint sets of symbols */
 641
 642     if not DisjointP( Variables, Parameters ) then
 643       ErrorHandler(
 644         "VarParConfl", Intersection( Variables, Parameters ), 'Fatal
 645       ),
 646
 647
 648     /* Make a list of all variables to be solved for. This list may contain */
 649     /* more symbols than Variables when either no SolverParams have been    */
 650     /* given or when the user has specified only a subset of all existing   */
 651     /* symbols.                                                             */
 652
 653     AllVars : SetDifference( listofvars( Equations ), Parameters ),
 654
 655     /* solve for all available variables if Variables is empty. */
 656
 657     if Empty( Variables ) then
 658       Variables : AllVars,
 659
 660
 661     /* Check all those equations which contain only equations of the form  */
 662     /* NUMBER = NUMBER or which contain only parameters and none of the    */
 663     /* variables of interest for consistency. Remove consistent equations  */
 664     /* and store the assumptions made in the list Solver_Assumptions.      */
 665
 666     Equations : ParamConsistency( Equations, Parameters ),
 667
 668     /* Abort if no equations are left after the above step. */
 669
 670     if Empty( Equations ) or Empty( Variables ) then (
 671       PrintMsg( 'SHORT, SolverMsg["NoEqOrVar"] ),
 672       return( [ [], SolverParams, Variables, Parameters, [], [], [] ] )
 673     ),
 674
 675     /* If there are compound expressions to be solved for then the solver     */
 676     /* tries to solve for the variables contained in them first.              */
 677     /* Subsequently the expressions are rebuilt from the solutions of         */
 678     /* these variables and the parameters.                                    */
 679
 680     for i thru length( Variables ) do
 681
 682       /* Search the variable list for non-atomic expressions. */
 683
 684       if not atom( Var : Variables[i] ) then (
 685
 686         /* append expression to the expression list */
 687         Expressions : endcons( Var, Expressions ),
 688
 689         /* Insert the variables in the expression into the list of variables  */
 690         /* but keep the parameters out.                                       */
 691         Variables[i] : SetDifference( listofvars( Var ), Parameters ),
 692
 693         PrintMsg( 'SHORT, SolverMsg["TrySolve4"], Variables[i] ),
 694         PrintMsg( 'SHORT, SolverMsg["Solve4Exp"], Var )
 695       ),
 696
 697     /* Make a list of all the variables the user actually wants to know */
 698     UserVars : sublist( Variables, 'atom ),
 699
 700     /* Flatten the list of variables and make it a set. */
 701     Variables : Setify( Flatten( Variables ) ),
 702
 703
 704     /* If Solver_Subst_Powers is true then substitute powers by new symbols */
 705
 706     if Solver_Subst_Powers then (
 707
 708       PrintMsg( 'SHORT, SolverMsg["SubstPwrs"] ),
 709
 710       for Var in AllVars do (
 711
 712         /* get all powers of Var in Equations. Convert negative powers to */
 713         /* positive ones.                                                 */
 714         Power : Setify( abs( ListOfPowers( Equations, Var ) ) ),
 715
 716         /* If there is more than one power of Var then substitute each  */
 717         /* var^m only if all m are integer multiples of the lowest      */
 718         /* power > 0. The check is done by examining the modulus of all */
 719         /* powers with respect to the lowest power.                     */
 720
 721         if
 722           not member( 'false, map( 'integerp, Power ) )
 723           and
 724           Power[1] # 1
 725         then
 726
 727           if ( length( Power ) = 1 ) or
 728
 729             block(
 730               [ Modulus : Power[1] ],
 731               not member(
 732                 'false,
 733                 map( 'ZeroP, totaldisrep( rat( rest( Power ) ) ) )
 734               )
 735             )
 736
 737           then (
 738             /* Make a new symbol for var^power */
 739             SubstSym : concat( Var, "^", Power[1] ),
 740
 741             /* Store a reference to the original term in an assoc list */
 742             PowerSubst : endcons( SubstSym = Var^Power[1], PowerSubst ),
 743
 744             /* Substitute the new symbol for the original term */
 745             Equations  : ratsubst( SubstSym, Var^Power[1], Equations ),
 746             Variables  : subst( SubstSym, Var, Variables ),
 747
 748             /* Notify user */
 749             PrintMsg( 'DETAIL, SolverMsg["subst"], Var^Power )
 750           )
 751
 752       ) /* END for Var in AllVars */
 753
 754     ), /* END if Solver_Subst_Powers */
 755
 756     return(
 757       [
 758         Equations, SolverParams, Variables, Parameters, Expressions,
 759         PowerSubst, UserVars
 760       ]
 761     )
 762   )
 763 )$
 764
 765
 766 /******************************************************************************/
 767 /* ParamConsistency checks equations of the form NUMBER = NUMBER and equa-    */
 768 /* tions which contain only parameters for consistency. If the system cannot  */
 769 /* determine whether a parametric expression is zero then the user is optio-  */
 770 /* nally asked to supply the required information. The assumptions made are   */
 771 /* stored in the list Solver_Assumptions.                                     */
 772 /******************************************************************************/
 773
 774 ParamConsistency( Eqs, Pars, [ Action ] ) := (
 775
 776   mode_declare(
 777     [ Eqs, Pars ], list,
 778     Action, any
 779   ),
 780
 781   block(
 782
 783     [ Eq, lhsminusrhs, i, consistent ],
 784
 785     mode_declare(
 786       [ Eq, lhsminusrhs ], any,
 787       i, fixnum,
 788       consistent, boolean
 789     ),
 790
 791     PrintMsg( 'SHORT, SolverMsg["ConsChk"] ),
 792
 793     if Empty( Action ) then
 794       Action : 'BREAK
 795     else
 796       Action : Action[1],
 797
 798     consistent : true,
 799
 800     for i thru length( Eqs ) do (
 801
 802       Eq : Eqs[i],
 803
 804       /* Does the equation contain only numbers or parameters? */
 805       if Empty(
 806         SetDifference( listofvars( Eq ), Pars )
 807       )
 808       then (
 809
 810         /* If so, check for consistency */
 811         lhsminusrhs :  expand( lhs( Eq ) - rhs( Eq ) ),
 812
 813         /* Test if difference of both rhs's is zero */
 814         if lhsminusrhs # 0 then
 815
 816           if SolverAssumeZero( lhsminusrhs ) then
 817             PrintMsg( 'SHORT, SolverMsg["Assum"], lhsminusrhs = 0 )
 818           else
 819             if Action = 'BREAK then (
 820               /* Abort if difference is non-zero */
 821               PrintMsg( 'SHORT, SolverMsg["Incons"], lhs( Eq ) = rhs( Eq ) ),
 822               TerminateSolver()
 823             )
 824             else
 825               return( consistent : false ),
 826
 827         /* Kill the now redundant equation */
 828         if Action = 'BREAK then
 829           Eqs[i] : []
 830
 831       ) /* END if Empty */
 832
 833     ), /* END for i */
 834
 835     if consistent then
 836       PrintMsg( 'SHORT, SolverMsg["NoneFnd"] ),
 837
 838     if Action = 'BREAK then
 839       return( delete( [], Eqs ) )
 840     else
 841       return( consistent )
 842   )
 843 )$
 844
 845
 846 /******************************************************************************/
 847 /* SolverAssumeZero checks whether Expression is (assumed to be) equal to     */
 848 /* zero. If it isn't, the function returns false or asks the user for his     */
 849 /* decision. New assumptions are appended to the list Solver_Assumptions.     */
 850 /******************************************************************************/
 851
 852 SolverAssumeZero( Expression ) := (
 853
 854   mode_declare(
 855     Expression, any
 856   ),
 857
 858   block(
 859
 860     [ AssumptionExists, i ],
 861
 862     mode_declare(
 863       i, fixnum,
 864       AssumptionExists, boolean
 865     ),
 866
 867     /* Return false immediately if Expression is a number # 0 or if  */
 868     /* Solver_Incons_Params is set to 'BREAK.                          */
 869
 870     if
 871       Empty( listofvars( Expression ) )
 872       or ( Solver_Incons_Params = 'BREAK )
 873     then
 874       return( false ),
 875
 876     /* Do a simple check to find out whether assumption already exists:    */
 877     /* Zero is substituted for Expression in the stored assumption. If the */
 878     /* result is zero then the assumption already exists.                  */
 879
 880     AssumptionExists : false,
 881     for i thru length( Solver_Assumptions ) while not AssumptionExists do
 882       if
 883         fullratsimp(
 884           ratsubst( 0, Expression, lhs( Solver_Assumptions[i] ) )
 885         ) = 0
 886       then
 887         AssumptionExists : true,
 888
 889     /* If the expression is not yet assumed to be equal to zero, ask the user */
 890     /* to decide whether it is.                                               */
 891
 892     if not AssumptionExists then
 893
 894       if
 895         ( Solver_Incons_Params = 'ASK )
 896         and
 897         ( AskZeroNonzero( Expression ) # 'zero )
 898       then
 899         return( false )
 900       else
 901         /* If difference is equal to zero or assumed to be so then store    */
 902         /* the constraint in the global list Solver_Assumptions. So the user*/
 903         /* has access to all assumptions made during the solution process   */
 904         /* and can check any numerical solutions for consistency with the   */
 905         /* assumptions.                                                     */
 906         Solver_Assumptions : endcons(  Expression = 0, Solver_Assumptions )
 907
 908     else
 909       PrintMsg( 'DETAIL, SolverMsg["AssmFnd"], Expression = 0 ),
 910
 911     return( true )
 912   )
 913 )$
 914
 915 AskZeroNonzero( Expression ) :=
 916   if equal( x, 0 ) = true then 'zero
 917   else block(
 918     [s : 'pnz],
 919     while s#'zero and s#'nonzero do (
 920       s : read( "Is", Expression, "zero or nonzero?" )
 921     ),
 922     s
 923   )$
 924
 925 /******************************************************************************/
 926 /* ListOfPowers returns the list of powers # 0 of a variable in a set of      */
 927 /* equations.                                                                 */
 928 /******************************************************************************/
 929
 930 ListOfPowers( Eqs, Var ) := (
 931
 932   mode_declare(
 933     Eqs, list,
 934     Var, any
 935   ),
 936
 937   delete(
 938     0,
 939     apply(
 940       'Union,
 941       map(
 942         lambda(
 943           [ Eq ],
 944           Powers( expand( lhs( Eq ) - rhs( Eq ) ), Var )
 945         ),
 946         Eqs
 947       )
 948     )
 949   )
 950 )$
 951
 952
 953 /******************************************************************************/
 954 /* ImmediateAssignments directly applies all immediate assignments of the     */
 955 /* form  var = rhs...  before the actual Solver is called.                    */
 956 /******************************************************************************/
 957
 958 ImmediateAssignments( Solutions, RemainingEqs, Variables, Parameters ) := (
 959
 960   mode_declare(
 961     [ Solutions, RemainingEqs, Variables, Parameters ], list
 962   ),
 963
 964   block(
 965
 966     [ i, Vars, AssignmentMade, Left, Right ],
 967
 968     mode_declare(
 969       i, fixnum,
 970       Vars, list,
 971       AssignmentMade, boolean,
 972       [ Left, Right ], any
 973     ),
 974
 975     PrintMsg( 'SHORT, SolverMsg["SrchImmed"] ),
 976
 977     AssignmentMade : false,
 978
 979     for i thru length( RemainingEqs ) do block(
 980       [],
 981
 982       Left  : lhs( RemainingEqs[i] ),
 983       Right : rhs( RemainingEqs[i] ),
 984
 985       /* Scan the equations for simple assignments of the form  X = Expr or   */
 986       /* Expr = X  and apply this assignment only if X is not a parameter and */
 987       /* if Expr contains only numbers or parameters.                         */
 988       /* Remark: Equations containing only parameters have already been remo- */
 989       /* ved from the equation list in SetupSolver.                           */
 990
 991       if symbolp( Left ) then (
 992
 993         Vars : listofvars( Right ),
 994
 995         /* Check if lhs is an isolated variable and make sure that there are */
 996         /* only parameters on the right-hand side.                           */
 997
 998         if freeof( Left, Right )
 999           and Empty( SetDifference( Vars, Parameters ) )
1000         then (
1001
1002           if assoc( Left, Solutions ) = false then (
1003             PrintMsg(
1004               'DETAIL,
1005               SolverMsg["Assign"], totaldisrep( Left = Right )
1006             ),
1007
1008             if member( Left, map( lhs, Solutions ) ) then (
1009               if assoc( Left, Solutions ) # Right then (
1010                 PrintMsg( 'SHORT, SolverMsg["Incons"], Left = Right ),
1011                 TerminateSolver()
1012               )
1013             )
1014             else
1015               Solutions : endcons( Left = Right, Solutions ),
1016             RemainingEqs[i] : [],
1017
1018             AssignmentMade : true
1019           ),
1020
1021           /* This return prevents Macsyma from entering the next if statement */
1022           /* which must only be executed when the current outer if statement  */
1023           /* has not been entered.                                            */
1024           return( 'DONE )
1025         )
1026       ),
1027
1028       if symbolp( Right ) then (
1029
1030         Vars : listofvars( Left ),
1031
1032         if freeof( Right, Left )
1033           and Empty( SetDifference( Vars, Parameters ) )
1034         then (
1035
1036           if assoc( Right, Solutions ) = false then (
1037             PrintMsg(
1038               'DETAIL,
1039               SolverMsg["Assign"], totaldisrep( Right = Left )
1040             ),
1041
1042             if member( Right, map( lhs, Solutions ) ) then (
1043               if assoc( Right, Solutions ) # Left then (
1044                 PrintMsg( 'SHORT, SolverMsg["Incons"], Left = Right ),
1045                 TerminateSolver()
1046               )
1047             )
1048             else
1049               Solutions : endcons( Right = Left, Solutions ),
1050             RemainingEqs[i] : [],
1051
1052             AssignmentMade : true
1053           )
1054
1055         )
1056       )
1057     ), /* END for i thru length */
1058
1059     if AssignmentMade then (
1060
1061       /* delete the used equations from the list */
1062       RemainingEqs : delete( [], RemainingEqs ),
1063
1064       if not Empty( RemainingEqs ) then (
1065
1066         /* Remove all the variables from the working list which have been */
1067         /* determined by an immediate assignment.                         */
1068         Variables : SetDifference( Variables, map( 'lhs, Solutions ) ),
1069
1070         /* Evaluate the remaining equations with the constraints */
1071         RemainingEqs : ev( RemainingEqs, Solutions ),
1072
1073         /* Do a parameter consistency check */
1074         RemainingEqs : ParamConsistency( RemainingEqs, Parameters )
1075       )
1076
1077       else
1078         PrintMsg( 'SHORT, SolverMsg["NoEqTerm"] )
1079
1080     )
1081     else
1082       PrintMsg( 'SHORT, SolverMsg["NoImmed"] ),
1083
1084     /* END if AssignmentMade */
1085
1086     return( [ AssignmentMade, Solutions, RemainingEqs, Variables ] )
1087   )
1088 )$
1089
1090
1091 /******************************************************************************/
1092 /* LinearSolver extracts blocks of linear equations from an arbitrary system  */
1093 /* of equations by a heuristic searching strategy and solves the linear block */
1094 /* if there is one.                                                           */
1095 /******************************************************************************/
1096
1097 LinearSolver( Solutions, Equations, Variables, Parameters ) := (
1098
1099   mode_declare(
1100     [ Solutions, Equations, Variables, Parameters ], list
1101   ),
1102
1103   block(
1104
1105     [
1106       CoeffMatrix, ValuationMatrix, ActiveVars,
1107       LinEqNos, LinVarNos, NewVars, LinSolVars,
1108       LinearEqs, LinearVars, LinearSolutions,
1109       i, j, me, mv, NumVars, NumEqs,
1110       MaxValVar, MaxValEq,
1111       EqValuation, VarValuation,
1112       rhsExpressions,
1113
1114       linsolvewarn             : false,
1115       linsolve_params          : false,
1116       Solve_Inconsistent_Error : false,
1117
1118       EqWasLast : false
1119     ],
1120
1121     mode_declare(
1122       [ CoeffMatrix, ValuationMatrix ], any,
1123       [
1124         LinEqNos, LinVarNos, ActiveVars, LinearEqs, LinearVars, LinSolVars,
1125         EqValuation, VarValuation, LinearSolutions, NewVars, rhsExpressions
1126       ], list,
1127       [ i, j, me, mv, NumVars, NumEqs, MaxValVar, MaxValEq ], fixnum
1128     ),
1129
1130     if Empty( Equations ) then (
1131
1132       if Empty( Variables ) then
1133         PrintMsg( 'SHORT, SolverMsg["AllSolved"] )
1134       else
1135         PrintMsg( 'SHORT, SolverMsg["NoEqLeft"], Variables ),
1136
1137       return( [ false, Solutions, Equations, Variables ] )
1138     )
1139     else if Empty( Variables ) then (
1140       PrintMsg( 'SHORT, SolverMsg["EqLeft"] ),
1141       return( [ false, Solutions, Equations, Variables ] )
1142     ),
1143
1144     PrintMsg( 'SHORT, SolverMsg["SrchLinEq"] ),
1145
1146     /* Make a list of all remaining variables */
1147     NewVars : listofvars( Equations ),
1148
1149     if Solver_Find_All_Linear_Vars then (
1150
1151       /* The following rather weird commands put the variables into the order */
1152       /* [ variables to be currently solved for, other variables ]. If the    */
1153       /* current variables are linear variables then it is more likely that   */
1154       /* they will have explicit solutions if they are located in the left    */
1155       /* half of the system of equations. Otherwise the will more likely be   */
1156       /* used as parameters of the null space (if there is one).              */
1157
1158       ActiveVars : append(
1159         Intersection( Variables, NewVars ),
1160         SetDifference( NewVars, append( Variables, Parameters ) )
1161       ),
1162       LinearEqs : Equations,
1163       Equations : []
1164     )
1165     else (
1166       ActiveVars : Intersection( Variables, NewVars ),
1167       LinearEqs : [],
1168       for i thru length( Equations ) do
1169         if
1170           not Empty( Intersection( listofvars( Equations[i] ), ActiveVars ) )
1171         then (
1172           LinearEqs : endcons( Equations[i], LinearEqs ),
1173           Equations[i] : []
1174         ),
1175       Equations : delete( [], Equations )
1176     ),
1177
1178     PrintMsg( 'SHORT, SolverMsg["wrt"], ActiveVars ),
1179
1180     /* Set up the complete coefficient matrix w.r.t. all variables */
1181     CoeffMatrix : ComplCoeffMatrix( LinearEqs, ActiveVars ),
1182
1183     /* The valuation matrix contains a 1 at each Position where a variable */
1184     /* appears in a nonlinear form, and 0's otherwise.                     */
1185
1186     ValuationMatrix : matrixmap(
1187       lambda(
1188         [ x ],
1189         if x = 'false then
1190           1
1191         else
1192           0
1193       ),
1194       CoeffMatrix
1195     ),
1196
1197     /* EqValuation contains the number of nonlinear variables for each eq. */
1198     EqValuation : map(
1199       lambda( [ Row ], apply( "+", Row ) ),
1200       ListMatrix( ValuationMatrix )
1201     ),
1202
1203     /* VarValuation contains for all vars the number of equations in which */
1204     /* var[i] appears in a nonlinear form.                                 */
1205     VarValuation : map(
1206       lambda( [ Row ], apply( "+", Row ) ),
1207       ListMatrix( transpose( ValuationMatrix ) )
1208     ),
1209
1210     LinEqNos  : makelist( i, i, 1, NumEqs : length( EqValuation ) ),
1211     LinVarNos : makelist( i, i, 1, NumVars : length( VarValuation ) ),
1212
1213
1214     /* Remove nonlinear equations and/or variables until only a linear    */
1215     /* block remains, i.e. for all i, j VarValuation[i] = 0 and           */
1216     /* EqValuation[j] = 0.                                                */
1217
1218     while
1219       ( apply( "+", VarValuation ) # 0 )
1220       and
1221       ( apply( "+", EqValuation ) # 0 )
1222     do (
1223
1224       /* Determine maximum equation valuation and number of corresponding */
1225       /* equation.                                                        */
1226
1227       me : 0,
1228       MaxValEq : -1,
1229       for i thru NumEqs do
1230         if EqValuation[i] > MaxValEq then (
1231           me : i,
1232           MaxValEq : mode_identity( fixnum, EqValuation[i] )
1233         ),
1234
1235       /* Determine maximum variable valuation and number of corresponding */
1236       /* variable.                                                        */
1237
1238       mv : 0,
1239       MaxValVar : -1,
1240       for j thru NumVars do
1241         if VarValuation[j] > MaxValVar then (
1242           mv : j,
1243           MaxValVar : mode_identity( fixnum, VarValuation[j] )
1244         ),
1245
1246
1247       if apply( Solver_Del_Eq, [ MaxValEq, MaxValVar ] ) then (
1248         i : me,
1249
1250         for j thru NumVars do (
1251           VarValuation[j] : VarValuation[j] - ValuationMatrix[i, j],
1252           ValuationMatrix[i, j] : 0
1253         ),
1254
1255         /* Mark equation as deleted */
1256         EqValuation[i] : 0,
1257         LinEqNos[i]    : 0
1258       )
1259       else (
1260         j : mv,
1261
1262         for i thru NumEqs do (
1263           EqValuation[i] : EqValuation[i] - ValuationMatrix[i, j],
1264           ValuationMatrix[i, j] : 0
1265         ),
1266
1267         /* Mark variable as deleted */
1268         VarValuation[j] : 0,
1269         LinVarNos[j]    : 0
1270       )
1271
1272     ), /* END while */
1273
1274     /* Make list of linear equations */
1275     LinEqNos  : delete( 0, LinEqNos ),
1276
1277     /* Make list of linear variables */
1278     LinVarNos : delete( 0, LinVarNos ),
1279
1280     if Empty( LinEqNos ) or Empty( LinVarNos ) then (
1281       PrintMsg( 'SHORT, SolverMsg["NoLinEqs"] ),
1282       return( [ false, Solutions, append( LinearEqs, Equations ), Variables ] )
1283     ),
1284
1285     /* Extract linear equations. append all nonlinear equations to Equations */
1286     /* again.                                                                */
1287     for i thru length( LinearEqs ) do
1288       if not member( i, LinEqNos ) then (
1289         Equations : endcons( LinearEqs[i], Equations ),
1290         LinearEqs[i] : []
1291       ),
1292
1293     LinearEqs : delete( [], LinearEqs ),
1294
1295
1296     /* Extract linear variables. Since the extraction of linear equations may */
1297     /* also have removed linear variables (the coefficients are now 0) it is  */
1298     /* necessary to intersect the set of the linear variables with the set of */
1299     /* those variables which actually appear in the linear equations.         */
1300
1301     LinearVars : Intersection(
1302       map( lambda( [ i ], ActiveVars[i] ), LinVarNos ),
1303       listofvars( LinearEqs )
1304     ),
1305
1306     /* By analogy, the same applies to the linear equations. Thus, keep only  */
1307     /* those equations which still contain any of the linear variables.       */
1308
1309     for i thru length( LinearEqs ) do
1310       if DisjointP( listofvars( LinearEqs[i] ), LinearVars ) then (
1311         Equations : endcons( LinearEqs[i], Equations ),
1312         LinearEqs[i] : []
1313       ),
1314
1315     LinearEqs : delete( [], LinearEqs ),
1316
1317     /* Return if no linear equations are left */
1318
1319     if Empty( LinearVars ) or Empty( LinearEqs ) then (
1320       PrintMsg( 'SHORT, SolverMsg["NoLinEqs"] ),
1321       return( [ false, Solutions, append( LinearEqs, Equations ), Variables ] )
1322     ),
1323
1324     PrintMsg(
1325       'SHORT,
1326       SolverMsg["Found"], length( LinearEqs ), SolverMsg["LinEqs"],
1327       length( LinearVars ), SolverMsg["LinVars"]
1328     ),
1329     PrintMsg( 'SHORT,  SolverMsg["VarsAre"],  LinearVars ),
1330     PrintMsg( 'DETAIL, SolverMsg["EqsAre"],   LinearEqs ),
1331     PrintMsg( 'SHORT,  SolverMsg["SolvLinEq"] ),
1332
1333     rhsExpressions : [],
1334     Solve_Inconsistent_Eqn_Nos : [ 0 ],
1335
1336     while not Empty( Solve_Inconsistent_Eqn_Nos ) do (
1337
1338       /* solve the linear equations */
1339       LinearSolutions : LinsolveM( LinearEqs , LinearVars ),
1340
1341       /* Check for "inconsistent" equations. */
1342       if not Empty( Solve_Inconsistent_Eqn_Nos ) then (
1343
1344         PrintMsg( 'DEBUG, SolverMsg["Incons"], Solve_Inconsistent_Eqn_Nos ),
1345
1346         /* Remove "inconsistent" equations */
1347         for i in Solve_Inconsistent_Eqn_Nos  do
1348             LinearEqs[i] : [],
1349
1350         LinearEqs : delete( [], LinearEqs ),
1351
1352         /* append rhs = 0 to Solver_Assumptions if rhs contains only */
1353         /* parameters.                                               */
1354         rhsExpressions : ParamConsistency(
1355           Solve_Inconsistent_Terms, Parameters
1356         )
1357
1358       )
1359
1360     ),
1361
1362     /* append rhs's which have led to "inconsistencies" but still   */
1363     /* contain variables to the list of equations.                  */
1364     Equations : append( Equations, rhsExpressions ),
1365
1366     PrintMsg( 'DETAIL, SolverMsg["Solutions"], LinearSolutions ),
1367
1368     /* Insert the solutions from linsolve into the remaining equations */
1369     Equations : ParamConsistency(
1370         fullratsimp( ev( Equations, LinearSolutions ) ),
1371         Parameters
1372     ),
1373
1374     /* append the linear solutions to the list of solutions */
1375     Solutions : append( Solutions, LinearSolutions ),
1376
1377     /* append all variables to the working list which appear on the rhs's of  */
1378     /* the linear solutions. delete all variables which have been solved for. */
1379
1380     /* Linear variables for which a solution has been obtained. */
1381     LinSolVars : map( 'lhs, LinearSolutions ),
1382
1383     /* Linear variables which are free parameters of the null space. */
1384     LinearVars : SetDifference( LinearVars, LinSolVars ),
1385
1386     /* append all those variables which are parameters of the null space of */
1387     /* the linear equations and which do not appear in the remaining        */
1388     /* equations to the list of parameters.                                 */
1389
1390     for Var in LinearVars do
1391       if freeof( Var, Equations ) then (
1392         PrintMsg( 'SHORT, SolverMsg["FreeVar2Par"], Var ),
1393         Parameters : endcons( Var, Parameters ),
1394         Variables  : delete( Var, Variables )
1395       ),
1396
1397     NewVars : SetDifference(
1398       listofvars( map( 'rhs, LinearSolutions ) ),
1399       Parameters
1400     ),
1401
1402     Variables : Union(
1403       SetDifference( Variables, LinSolVars ),
1404       NewVars
1405     ),
1406
1407     return( [ true, Solutions, Equations, Variables ] )
1408   )
1409 )$
1410
1411
1412 /******************************************************************************/
1413 /* The following strategies decide whether the linear solver should delete a  */
1414 /* nonlinear equation or a nonlinear variable from the system while searching */
1415 /* for linear subblocks of equations.                                         */
1416 /******************************************************************************/
1417
1418 define_variable( EqWasLast, false, boolean )$
1419
1420 MakeSquareLinearBlocks( ValEq, ValVar ) := (
1421
1422   mode_declare(
1423     [ ValEq, ValVar ], fixnum
1424   ),
1425
1426   if ValEq = ValVar then
1427     EqWasLast : not EqWasLast
1428   else
1429     if ValEq > ValVar then
1430       EqWasLast : true
1431     else
1432       EqWasLast : false
1433 )$
1434
1435
1436 DelEqBeforeVar( ValEq, ValVar ) := (
1437
1438   mode_declare(
1439     [ ValEq, ValVar ], fixnum
1440   ),
1441
1442   if ValEq >= ValVar then
1443     true
1444   else
1445     false
1446 )$
1447
1448
1449 /******************************************************************************/
1450 /* ComplCoeffMatrix returns a matrix whose row size is equal to the number of */
1451 /* equations and whose column size is equal to the number of variables. The   */
1452 /* entry at Position [i,j] is RatCoeff( equation[i], variable[j] ) if         */
1453 /* equation[i] is linear w.r.t. variable[j] and 'false if equation[i] is      */
1454 /* nonlinear w.r.t. variable[j].                                              */
1455 /******************************************************************************/
1456
1457 ComplCoeffMatrix( Eqs, ActiveVars ) := (
1458
1459   mode_declare(
1460     [ Eqs, ActiveVars ], list
1461   ),
1462
1463   block(
1464
1465
1466     apply(
1467
1468       'matrix,
1469
1470       /* for each equation do */
1471       map(
1472
1473         lambda(
1474           [ Eq ],
1475
1476           /* for each variable do */
1477           map(
1478             lambda(
1479               [ Var ],
1480               block(
1481                 [ rc ],
1482                 rc  : LinCoeff( Eq, Var ),
1483
1484                 /* Return the RatCoeff only if it contains none of the active */
1485                 /* variables or if the equation doesn't contain var at all.   */
1486
1487                 if
1488                   ( ( rc # 0 ) and DisjointP( listofvars( rc ), ActiveVars ) )
1489                   or
1490                   freeof( Var, Eq )
1491                 then
1492                   rc
1493                 else
1494                   false
1495               )
1496             ),
1497
1498             ActiveVars
1499           ) /* END lambda( [ Var ] ) */
1500
1501         ), /* END lambda( [ Eq ] ) */
1502
1503
1504         /* map target: Transform all equations into homogeneous form. */
1505
1506         map(
1507           lambda(
1508             [ Eq ],
1509             fullratsimp( expand( lhs( Eq ) - rhs( Eq ) ) )
1510           ),
1511           Eqs
1512         )
1513
1514
1515       ) /* END map( lambda( [ Eq ] ) ) */
1516
1517     ) /* END apply */
1518   )
1519 )$
1520
1521
1522 /******************************************************************************/
1523 /* LinCoeff returns the linear coefficient of Var within Eq if Var appears    */
1524 /* raised to the first power only.                                            */
1525 /******************************************************************************/
1526
1527 LinCoeff( Eq, Var ) := (
1528
1529   mode_declare(
1530     [ Eq, Var ], any
1531   ),
1532
1533   block(
1534
1535     [ BCoeff ],
1536
1537     mode_declare(
1538       BCoeff, list
1539     ),
1540
1541     if ListOfPowers( [ Eq ], Var ) = [ 1 ] then (
1542       BCoeff : bothcoeff( Eq, Var ),
1543       if freeof( Var, second( BCoeff ) ) then
1544         return( first( BCoeff ) )
1545     ),
1546
1547     return( 0 )
1548   )
1549 )$
1550
1551
1552 /******************************************************************************/
1553 /* ValuationSolver                                                            */
1554 /******************************************************************************/
1555
1556 ValuationSolver( Solutions, Equations, Variables, Parameters ) := (
1557
1558   mode_declare(
1559     [ Solutions, Equations, Variables, Parameters ], list
1560   ),
1561
1562   block(
1563
1564     [
1565       VarPaths, ValMatrix, Eq, Var, Trans, TempEq, TransEq,
1566       SolveOrder, SolveInfo, Transform, Solution, SolCheck,
1567       Status, Solved, Failed, UniqueSol, TryToSolve, CheckSol,
1568       i, k
1569     ],
1570
1571     mode_declare(
1572       [ VarPaths, ValMatrix, Eq, Var, Trans, TempEq, TransEq ], any,
1573       [ SolveOrder, SolveInfo, Transform, Solution, SolCheck ], list,
1574       [ Status, Solved, Failed, UniqueSol, TryToSolve, CheckSol ], boolean,
1575       [ i, k ], fixnum
1576     ),
1577
1578
1579     UniqueSol : true,
1580
1581     LOOP,
1582
1583     PrintMsg( 'SHORT, SolverMsg["Chk4RemEq"] ),
1584
1585     if Empty( Equations ) then (
1586
1587       if Empty( Variables ) then
1588         PrintMsg( 'SHORT, SolverMsg["AllSolved"] )
1589       else
1590         PrintMsg( 'SHORT, SolverMsg["NoEqLeft"], Variables ),
1591
1592       Status : false
1593
1594     )
1595     else if Empty( Variables ) then (
1596       PrintMsg( 'SHORT, SolverMsg["EqLeft"] ),
1597       Status : false
1598     )
1599     else (
1600
1601       PrintMsg(
1602         'SHORT,
1603         length( Equations ), SolverMsg["Eqs"],
1604         length( Variables ), SolverMsg["Vars"]
1605       ),
1606       PrintMsg( 'DETAIL, SolverMsg["VarsAre"], Variables ),
1607       PrintMsg( 'DEBUG, SolverMsg["EqsAre"], Equations ),
1608
1609       /* Dump solutions and remaining equations to file if requested. */
1610       if Solver_Dump_To_File then
1611         DumpToFile( Solutions, Equations, Variables ),
1612
1613       if ( length( Variables ) = 1 ) and ( length( Equations ) = 1 ) then
1614         SolveOrder : [ [ 1, 1, "(irrelevant)" ] ]
1615
1616       else (
1617         PrintMsg( 'SHORT, SolverMsg["ValStrat"] ),
1618
1619         /* Set up the valuation matrices. */
1620         VarPaths  : OccurrenceMatrix( Equations, Variables ),
1621         ValMatrix : ValuationMatrix( Equations, Variables ),
1622
1623         /* Determine an order by which the equations should be solved. */
1624         SolveOrder : apply( Solver_Valuation_Strategy, [ VarPaths, ValMatrix ] )
1625       ),
1626
1627       Solved : false,
1628
1629       unless Solved or Empty( SolveOrder ) do (
1630
1631         SolveInfo : Pop( SolveOrder ),
1632         if listp(SolveInfo[1]) then SolveInfo : first(SolveInfo),
1633         k   : mode_identity( fixnum, first( SolveInfo ) ),
1634         Eq  : part( Equations, k ),
1635         Var : Variables[ second( SolveInfo ) ],
1636
1637         PrintMsg(
1638           'SHORT,
1639           SolverMsg["TrySolveEq"], k, SolverMsg["ForVar"], Var
1640         ),
1641         PrintMsg( 'SHORT, SolverMsg["Valuation"], third( SolveInfo ) ),
1642         PrintMsg( 'DETAIL, SolverMsg["EqIs"], Eq = 0 ),
1643
1644         TryToSolve : true,
1645         CheckSol   : true,
1646         Transform  : copylist( Solver_Transforms ),
1647
1648         /* Do the solver break test to check whether it is worth */
1649         /* attempting to solve the equation at all.              */
1650         Failed : apply( Solver_Break_Test, [Eq, Var, third( SolveInfo )] ),
1651
1652         unless Solved or Failed do (
1653
1654           /* Try to solve the selected equation */
1655           if TryToSolve then
1656             Solution : solve( Eq, Var ),
1657
1658           /* Check if the equation was solved correctly */
1659           if CheckSol then (
1660             PrintMsg( 'SHORT, SolverMsg["CheckSol"] ),
1661             SolCheck : SolutionOK( Solution, Var ),
1662
1663             PrintMsg( 'DETAIL, SolverMsg["Solutions"], Solution )
1664           )
1665           else
1666             SolCheck : [ false ],
1667
1668           /* All solutions OK? */
1669           if member( true, SolCheck ) then (
1670             PrintMsg( 'SHORT, SolverMsg["SolOK"] ),
1671             Solved : true
1672           )
1673
1674           /* If not, apply transformations */
1675           else (
1676             if CheckSol then
1677               PrintMsg( 'SHORT, SolverMsg["SolNotOK"] ),
1678
1679             /* Give up if no transformations are left */
1680             if Empty( Transform ) then (
1681               PrintMsg( 'SHORT, SolverMsg["GiveUp"] ),
1682               Failed : true
1683             )
1684
1685             else (
1686               /* Retrieve one transformation function */
1687               Trans : Pop( Transform ),
1688               PrintMsg( 'SHORT,  SolverMsg["AppTrans"], Trans ),
1689
1690               /* and apply it to the equation, the variable, and the solution */
1691               TransEq : apply( Trans, [ Eq, Var, Solution ] ),
1692
1693               /* The transformation should return an equation as its function */
1694               /* value. However, if no reasonable transformation of the       */
1695               /* equation was possible then the SOLVE function should not be  */
1696               /* tried again. Hence, to signal a failure, the transformation  */
1697               /* must return an empty list, which will instruct the Solver to */
1698               /* try the next transformation instead. In addition, the        */
1699               /* transformation may itself take care of solving the equation. */
1700               /* It must then return a list of solutions:                     */
1701               /*   [ var = solution_1, var = solution_2, ... ]                */
1702
1703               /* Did the transformation fail? */
1704               if TransEq = [] then (
1705                 PrintMsg( 'SHORT, SolverMsg["TransFail"] ),
1706
1707                 /* Instruct the Solver to try the next transformation */
1708                 TryToSolve : false,
1709                 CheckSol   : false
1710               )
1711
1712               /* Did it solve the equation by itself? */
1713               else if listp( TransEq ) then (
1714                 PrintMsg( 'SHORT, SolverMsg["TransSolv"] ),
1715                 Solution : TransEq,
1716
1717                 /* Instruct the Solver not to call SOLVE again */
1718                 TryToSolve : false,
1719                 CheckSol   : true
1720               )
1721
1722               /* Transformation thinks it has succeeded, so try again */
1723               else (
1724                 Eq : TransEq,
1725                 PrintMsg( 'DETAIL, SolverMsg["ResTrans"], Eq = 0 ),
1726                 PrintMsg( 'SHORT,  SolverMsg["RetryTrans"] ),
1727
1728                 TryToSolve : true,
1729                 CheckSol   : true
1730               )
1731
1732             ) /* END if Empty( Transform ) else */
1733
1734           ) /* if member( true, SolCheck ) else */
1735
1736         ) /* END unless Solved or Failed */
1737
1738       ), /* END unless Solved or Empty( SolveOrder ) */
1739
1740       if Solved then (
1741         if length( Solution ) > 1 then
1742           PrintMsg( 'SHORT, SolverMsg["NotUnique"] ),
1743
1744         if member( false, SolCheck ) then
1745           PrintMsg( 'SHORT, SolverMsg["SolsLost"] ),
1746
1747         /* Remove solved equation from list of equations. Store it in TempEq  */
1748         /* so it can be appended to Equations again if the consistency check  */
1749         /* fails.                                                             */
1750         TempEq : part( Equations, k ),
1751         Equations : delete( [], Set_Element( Equations, k, [] ) ),
1752
1753         /* Check solutions for consistency with remaining equations. */
1754         for i thru length( Solution ) do (
1755
1756           if part( SolCheck, i ) then (
1757
1758             PrintMsg(
1759               'DETAIL, SolverMsg["Solution"], i, SolverMsg["ForVar"], Var
1760             ),
1761
1762             if not ParamConsistency(
1763               fullratsimp( ev( Equations, part( Solution, i ) ) ),
1764               Parameters,
1765               'CONTINUE
1766             ) then (
1767               PrintMsg( 'SHORT, SolverMsg["Contradict"], part( Solution, i ) ),
1768               Set_Element( Solution, i, 'INCONSISTENT_PATH )
1769             )
1770
1771           )
1772           else (
1773             PrintMsg( 'DETAIL, SolverMsg["Dropped"], part( Solution, i ) ),
1774             Set_Element( Solution, i, [] )
1775           )
1776
1777         ),
1778
1779         /* delete all implicit or empty solutions */
1780         Solution : delete( [], Solution ),
1781
1782         if Empty( Solution ) then (
1783           PrintMsg( 'SHORT, SolverMsg["NoValidSol"], Var ),
1784
1785           Equations : endcons( TempEq, Equations ),
1786           Solved : false
1787         )
1788
1789         /* Check if there are any consistent solutions */
1790         else if not member(
1791           true,
1792           map(
1793             lambda(
1794               [x],
1795               if x = 'INCONSISTENT_PATH then
1796                 false
1797               else
1798                 true
1799             ),
1800             Solution
1801           )
1802         )
1803         then (
1804           PrintMsg( 'SHORT, SolverMsg["NoConsSol"], Var ),
1805
1806           Solutions : endcons( 'INCONSISTENT_PATH, Solutions ),
1807           Solved : false
1808         )
1809
1810         /* append consistent solutions to the solution list */
1811         else (
1812           PrintMsg( 'DETAIL, SolverMsg["ConsSol"], Var, ":", Solution ),
1813
1814           Variables : delete( Var, Variables ),
1815
1816           /* If the solution is unique or if there's only one consistent */
1817           /* solution then ...                                           */
1818           if length( Solution ) = 1 then (
1819
1820             /* ... store it, insert it into the remaining equations, and */
1821             /* add its rhs variables to the list of unknowns.            */
1822             Solutions : append( Solutions, Solution ),
1823             Equations : copylist(
1824               fullratsimp( ev( Equations, Solution ) )
1825             ),
1826             Variables : Union(
1827               Variables,
1828               SetDifference(
1829                 listofvars( rhs( first( Solution ) ) ),
1830                 Parameters
1831               )
1832             )
1833           )
1834
1835           else /* length( Solution ) > 1, call ValuationSolver recursively */
1836
1837             block(
1838               [ MultipleSolutions, RSolutions, RVars, REqs, Sol, Stat ],
1839
1840               mode_declare(
1841                 [ MultipleSolutions, RSolutions, RVars, REqs ], list,
1842                 Sol, any,
1843                 Stat, boolean
1844               ),
1845
1846               MultipleSolutions : [],
1847
1848               for Sol in Solution do (
1849
1850                 map(
1851                   lambda([x,y], x::y),
1852                   [ 'Stat, 'RSolutions, 'REqs, 'RVars ],
1853                   ValuationSolver(
1854                     [ Sol ],
1855                     copylist( fullratsimp( ev( Equations, Sol ) ) ),
1856                     Union(
1857                       Variables,
1858                       SetDifference( listofvars( rhs( Sol ) ), Parameters )
1859                     ),
1860                     Parameters
1861                   )
1862                 ),
1863
1864                 MultipleSolutions : endcons( RSolutions, MultipleSolutions )
1865
1866               ), /* END for Sol */
1867
1868               Solutions : endcons( MultipleSolutions, Solutions ),
1869
1870               UniqueSol : false
1871             ) /* END block */
1872
1873         ) /* END if Empty( Solution ) */
1874
1875       )
1876
1877       else (
1878
1879         /* append remaining equations to the solutions if no further */
1880         /* solutions could be determined.                            */
1881
1882         for e in Equations do (
1883           if is( equal( e, 0 ) ) # 'unknown then (
1884             PrintMsg( SHORT, "Inconsistent equation", e = 0),
1885             TerminateSolver()
1886           )
1887         ),
1888
1889         Solutions : endcons( [ Equations ], Solutions )
1890       ), /* END if Solved */
1891
1892       Status : Solved
1893     ),
1894
1895     if Status and UniqueSol then
1896       go( LOOP )
1897     else
1898       return( [ Status, Solutions, Equations, Variables ] )
1899   )
1900 )$
1901
1902
1903 /******************************************************************************/
1904 /* SolutionOK checks whether the result of a call to the solve function is    */
1905 /* indeed a solution of the form  var = expression_free_of_var.               */
1906 /******************************************************************************/
1907
1908 SolutionOK( Solution, Var ) := (
1909
1910   mode_declare(
1911     [ Solution, Var ], any
1912   ),
1913
1914   if listp( Solution ) then
1915
1916     /* List of solutions must not be empty. */
1917     if Empty( Solution ) then
1918       [ false ]
1919
1920     else
1921
1922       /* Check if the lhs of each solution is equal to var and make sure */
1923       /* that var does not appear on the rhs's.                          */
1924       map(
1925         lambda(
1926           [ Sol ],
1927           ( lhs( Sol ) = Var ) and freeof( Var, rhs( Sol ) )
1928         ),
1929         Solution
1930       )
1931
1932   else
1933     [ false ]
1934
1935 )$
1936
1937
1938 /******************************************************************************/
1939 /* ValuationMatrix generates a matrix of valuations with respect to each      */
1940 /* equation and each variable.                                                */
1941 /******************************************************************************/
1942
1943 ValuationMatrix( Equations, Variables ) := (
1944
1945   mode_declare(
1946     [ Equations, Variables ], list
1947   ),
1948
1949   genmatrix(
1950     lambda( [ i, j ], Valuation( Equations[i], Variables[j] ) ),
1951     length( Equations ), length( Variables )
1952   )
1953 )$
1954
1955
1956 /******************************************************************************/
1957 /* Operator valuation factors for expression valuation.                       */
1958 /******************************************************************************/
1959
1960 SetProp( 'sin,   'Valuation, 10 )$
1961 SetProp( 'cos,   'Valuation, 10 )$
1962 SetProp( 'tan,   'Valuation, 10 )$
1963 SetProp( 'asin,  'Valuation, 12 )$
1964 SetProp( 'acos,  'Valuation, 12 )$
1965 SetProp( 'atan,  'Valuation, 12 )$
1966 SetProp( 'sinh,  'Valuation, 12 )$
1967 SetProp( 'cosh,  'Valuation, 12 )$
1968 SetProp( 'tanh,  'Valuation, 12 )$
1969 SetProp( 'asinh, 'Valuation, 12 )$
1970 SetProp( 'acosh, 'Valuation, 12 )$
1971 SetProp( 'atanh, 'Valuation, 12 )$
1972 SetProp( "+",    'Valuation,  1 )$
1973 SetProp( "-",    'Valuation,  1 )$
1974 SetProp( "*",    'Valuation,  4 )$
1975 SetProp( "/",    'Valuation,  4 )$
1976 SetProp( "^",    'Valuation, 10 )$
1977 SetProp( 'sqrt,  'Valuation, 10 )$
1978 SetProp( 'exp,   'Valuation, 10 )$
1979 SetProp( 'log,   'Valuation, 10 )$
1980
1981
1982 /******************************************************************************/
1983 /* With SetValuation, the operator valuation factors can be redefined.        */
1984 /******************************************************************************/
1985
1986 SetValuation( Operator, Valuation ) :=
1987   SetProp( Operator, 'Valuation, Valuation )$
1988
1989
1990 /******************************************************************************/
1991 /* Valuation measures the complexity of an expression with respect to Var by  */
1992 /* weighting the operator tree representation of Expr.                        */
1993 /******************************************************************************/
1994
1995 Valuation( Expr, Var ) := (
1996
1997   mode_declare(
1998     [ Expr, Var ], any
1999   ),
2000
2001
2002   block(
2003
2004     [ OpFactor ],
2005
2006     mode_declare(
2007       OpFactor, fixnum
2008     ),
2009
2010     /* Return zero if Expr does not contain Var. */
2011     if freeof( Var, Expr ) then
2012       return( 0 )
2013
2014     else
2015       /* Return 1 if Expr is an atom, i.e. Expr = Var. */
2016       if atom( Expr ) then
2017         return( 1 )
2018
2019       /* If Expr is an algebraic expression then retrieve the valuation    */
2020       /* factor associated with the operator of Expr and recursively apply */
2021       /* the valuation function to each subexpression of Expr.             */
2022       else (
2023
2024         if (
2025             OpFactor : mode_identity( fixnum, get( op( Expr ), 'Valuation ) )
2026            ) = false
2027         then
2028           OpFactor : Solver_Default_Valuation,
2029
2030         return(
2031           OpFactor * apply(
2032             "+",
2033             map(
2034               lambda( [ SubExpr ], Valuation( SubExpr, Var ) ),
2035               substpart( "[", Expr, 0 )
2036             )
2037           )
2038         )
2039
2040       ) /* END if atom */
2041
2042   )
2043 )$
2044
2045
2046 /******************************************************************************/
2047 /* OccurrenceMatrix sets up a matrix in which the number of occurrences of     */
2048 /* each variable in each equation is counted.                                 */
2049 /******************************************************************************/
2050
2051 OccurrenceMatrix( Equations, Variables ) := (
2052
2053   mode_declare(
2054     [ Equations, Variables ], list
2055   ),
2056
2057   genmatrix(
2058     lambda( [ i, j ], Occurences( Equations[i], Variables[j] ) ),
2059     length( Equations ), length( Variables )
2060   )
2061 )$
2062
2063
2064 /******************************************************************************/
2065 /* Occurences counts the number of occurences of Var in Expr, i.e. the number */
2066 /* of paths to distinct occurrences of the atom Var in the internal tree      */
2067 /* representation of Expr.                                                    */
2068 /******************************************************************************/
2069
2070 Occurences( Expr, Var ) := (
2071
2072   mode_declare(
2073     [ Expr, Var ], any
2074   ),
2075
2076   if atom( Expr ) then
2077
2078     if Expr = Var then
2079       1
2080     else
2081       0
2082
2083   else
2084     apply(
2085       "+",
2086       map(
2087         lambda( [ SubExpr], Occurences( SubExpr, Var ) ),
2088         substpart( "[", Expr, 0)
2089       )
2090     )
2091 )$
2092
2093
2094 /******************************************************************************/
2095 /* MinVarPathsFirst tries to find variables which can be easily isolated.     */
2096 /* These are variables which appear only once in an entire expression tree    */
2097 /* (= 1 in OccurrenceMatrix).                                                 */
2098 /******************************************************************************/
2099
2100 MinVarPathsFirst( OccMat, ValMat ) := (
2101
2102   mode_declare(
2103     [ OccMat, ValMat ], any
2104   ),
2105
2106   block(
2107
2108     [
2109       SolveOrder, SolveOrder1, SumVarPaths,
2110       i, j, v, ne, nv
2111     ],
2112
2113     mode_declare(
2114       [ SolverOrder, SolveOrder1, SumVarPaths ], list,
2115       [ i, j, v, ne, nv, Function( RowSize, ColSize, Position ) ], fixnum
2116     ),
2117
2118     SolveOrder : [],
2119
2120     ne : RowSize( OccMat ),
2121     nv : ColSize( OccMat ),
2122
2123     SumVarPaths : map(
2124       lambda( [ Row ], apply( "+", Row ) ),
2125       ListMatrix( OccMat )
2126     ),
2127
2128     /* Search for equations which contain only one variable in one path. */
2129     for i thru length( SumVarPaths ) do
2130       if SumVarPaths[i] = 1 then (
2131         SolveOrder : endcons(
2132           [ i, j : Position( 1, OccMat[i] ), ValMat[i, j] ],
2133           SolveOrder
2134         ),
2135         /* Mark eq/var Position as used. */
2136         OccMat[i, j] : 0,
2137         ValMat[i, j] : 0
2138       ),
2139
2140     /* Sort SolveOrder by least valuation. */
2141     if not Empty( SolveOrder ) then
2142       SolveOrder : SortSolveOrder( SolveOrder ),
2143
2144     /* Find all variables with only one path in the expression tree. */
2145     SolveOrder1 : [],
2146     for i thru ne do
2147       for j thru nv do
2148         if OccMat[i, j] = 1 then (
2149           SolveOrder1 : endcons( [ i, j, ValMat[i, j] ], SolveOrder1 ),
2150           OccMat[i, j] : 0,
2151           ValMat[i, j] : 0
2152         ),
2153
2154     /* Sort variables by least valuation. */
2155     if not Empty( SolveOrder1 ) then
2156       SolveOrder : append(
2157         SolveOrder,
2158         SortSolveOrder( SolveOrder1 )
2159       ),
2160
2161     /* append additional candidates if necessary. */
2162     if length( SolveOrder ) < Solver_Max_Len_Val_Order then (
2163
2164       SolveOrder1 : [],
2165       for i thru ne do
2166         for j thru nv do
2167           if ( v : mode_identity( fixnum, ValMat[i, j] ) ) # 0 then
2168             SolveOrder1 : endcons( [ i, j, v ], SolveOrder1 ),
2169
2170       SolveOrder : append(
2171         SolveOrder,
2172         SortSolveOrder( SolveOrder1 )
2173       ),
2174
2175       /* Return only as many candidates as given by Solver_Max_Len_Val_Order */
2176       if ( i : length( SolveOrder ) ) > Solver_Max_Len_Val_Order then
2177         SolveOrder : rest(
2178           SolveOrder, Solver_Max_Len_Val_Order - i
2179         )
2180     ),
2181
2182     return( SolveOrder )
2183   )
2184 )$
2185
2186
2187 /******************************************************************************/
2188 /* PostProcess does all the postprocessing needed to display the results.     */
2189 /* This includes expansion of the solution list hierarchies, backsubsitution, */
2190 /* and extraction of the variables which the user explicitly asked for.       */
2191 /******************************************************************************/
2192
2193 PostProcess( Solutions, UserVars, Expressions, PowerSubst ) := (
2194
2195   mode_declare(
2196     [ Solutions, UserVars, Expressions, PowerSubst ], list
2197   ),
2198
2199   block(
2200     [
2201       SolSet, UsrSolSet, UserSolutions, UnsolvedEqs, InternalSols,
2202       Var, TempVar, EvalVar,
2203       i
2204     ],
2205
2206     mode_declare(
2207       [ SolSet, UsrSolSet, UserSolutions, UnsolvedEqs, InternalSols ], list,
2208       [ Var, TempVar, EvalVar ], any,
2209       i, fixnum
2210     ),
2211
2212
2213     PrintMsg( 'SHORT, SolverMsg["PostPr"] ),
2214
2215     /* first of all, Flatten the solution list hierarchy and drop all */
2216     /* inconsistent solution paths.                                   */
2217     Solutions : sublist(
2218       ExpandSolutionHierarchy( Solutions ),
2219       lambda( [Set], last( Set ) # 'INCONSISTENT_PATH )
2220     ),
2221
2222     /* Return an empty list if no consistent solution paths are left. */
2223     if Empty( Solutions ) then
2224       return( [] ),
2225
2226     /* Do the backsubstitutions. */
2227     if not Empty( Solutions ) then (
2228
2229       UserSolutions : [],
2230       i : 0,
2231
2232       for SolSet in Solutions do (
2233
2234         i : i + 1,
2235         PrintMsg( 'SHORT, SolverMsg["SolSet"], i ),
2236
2237         UsrSolSet   : [],
2238
2239         /* Extract the unsolved equations */
2240         UnsolvedEqs : sublist( SolSet, lambda( [x], not EquationP( x ) ) ),
2241
2242         /* and the solutions. */
2243         SolSet : sublist( SolSet, 'EquationP ),
2244
2245         /* If no complete backsubstitution is requested then variables on the */
2246         /* right-hand sides of the solutions will only be substituted if they */
2247         /* do not belong to the variables specified in the command line.      */
2248         if not Solver_Backsubst then
2249           InternalSols : sublist(
2250             SolSet,
2251             lambda( [x], not member( lhs( x ), UserVars ) )
2252           ),
2253
2254         /* Evaluate all variables and expressions with the solutions. */
2255         for Var in append( UserVars, Expressions ) do (
2256
2257           EvalVar : if member( Var, map(lhs, SolSet) ) or not atom( Var ) then
2258
2259             /* There may be errors when indeterminate expressions are */
2260             /* encountered.                                           */
2261
2262             errcatch(
2263               if Solver_Backsubst then
2264                 ev( Var, SolSet, infeval )
2265               else (
2266                 TempVar : ev( Var, SolSet ),
2267                 ev( TempVar, InternalSols, infeval )
2268               )
2269             )
2270
2271           else
2272             [],
2273
2274           if EvalVar # [] then
2275             UsrSolSet : endcons( Var = EvalVar[1], UsrSolSet )
2276           else
2277             PrintMsg( 'SHORT, SolverMsg["NoSol"], Var )
2278         ),
2279
2280
2281         /* append the unsolved equations. */
2282         if not Empty( UnsolvedEqs ) then
2283           UsrSolSet : endcons( UnsolvedEqs, UsrSolSet ),
2284
2285         if Solver_RatSimp_Sols then
2286           UsrSolSet : errcatch( fullratsimp( UsrSolSet ) )
2287         else
2288           UsrSolSet : [UsrSolSet],
2289
2290         if UsrSolSet = [] then
2291           PrintMsg( 'SHORT, SolverMsg["SolSetDrp"] )
2292         else
2293           UserSolutions : endcons( UsrSolSet[1], UserSolutions )
2294
2295       ), /* END for SolSet */
2296
2297       if Solver_Dump_To_File then
2298         DumpToFile( UserSolutions, [], [] )
2299
2300     ), /* END if not Empty( Solutions ) */
2301
2302
2303     return( UserSolutions )
2304
2305   ) /* END block */
2306 )$
2307
2308
2309 /******************************************************************************/
2310 /* ExpandSolutionHierarchy transforms the hierarchically structured list of   */
2311 /* solutions into a list of flat lists of solutions.                          */
2312 /******************************************************************************/
2313
2314 ExpandSolutionHierarchy( Solutions ) := (
2315
2316   mode_declare(
2317     Solutions, list
2318   ),
2319
2320   block(
2321     [ FlatSolutions ],
2322
2323     mode_declare(
2324       FlatSolutions, list
2325     ),
2326
2327     if length( Solutions ) = 0 then return ( [] )
2328     /* listp = true indicates an additional recursion level */
2329     else if listp( last( Solutions ) ) then (
2330
2331       FlatSolutions : rest( Solutions, -1 ),
2332
2333       return(
2334
2335         map(
2336           lambda( [ x ], append( FlatSolutions, x ) ),
2337
2338           apply(
2339             'append,
2340             map( 'ExpandSolutionHierarchy, last( Solutions ) )
2341           )
2342         )
2343
2344       )
2345
2346     )
2347
2348     else
2349       return( [ Solutions ] )
2350
2351   )
2352 )$
2353
2354
2355 /******************************************************************************/
2356 /* DumpToFile dumps the current set of solutions, equations and variables to  */
2357 /* the file <Solver_Dump_File>.                                               */
2358 /******************************************************************************/
2359
2360 DumpToFile( Sols, Eqs, Vars ) := (
2361
2362   mode_declare(
2363     [ Sols, Eqs, Vars ], list
2364   ),
2365
2366   block(
2367     [ Solutions, Equations, Variables ],
2368
2369     PrintMsg( 'SHORT, SolverMsg["Dump"], Solver_Dump_File ),
2370
2371     apply(
2372       'StringOut,
2373       [
2374         Solver_Dump_File,
2375         'Solutions = Sols,
2376         'Equations = Eqs,
2377         'Variables = Vars
2378       ]
2379     )
2380   )
2381 )$
2382
2383 tma():=trace(
2384 TerminateSolver,
2385 SetupSolver,
2386 ParamConsistency,
2387 SolverAssumeZero,
2388 ListOfPowers,
2389 ImmediateAssignments,
2390 LinearSolver,
2391 MakeSquareLinearBlocks,
2392 DelEqBeforeVar,
2393 ComplCoeffMatrix,
2394 LinCoeff,
2395 ValuationSolver,
2396 SolutionOK,
2397 ValuationMatrix,
2398 SetValuation,
2399 Valuation,
2400 OccurrenceMatrix,
2401 Occurences,
2402 MinVarPathsFirst,
2403 PostProcess,
2404 ExpandSolutionHierarchy,
2405 DumpToFile)$
2406
2407 \end{verbatim}\r