check/exam_lsolve.cpp

   1 /** @file exam_lsolve.cpp
   2  *
   3  *  These exams test solving small linear systems of symbolic equations. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  21  */
  22
  23 #include "exams.h"
  24
  25 static unsigned exam_lsolve1()
  26 {
  27         // A trivial example.
  28         unsigned result = 0;
  29         symbol x("x");
  30         ex eq, aux;
  31
  32         eq = (3*x+5 == numeric(8));
  33         aux = lsolve(eq, x);
  34         if (aux != 1) {
  35                 ++result;
  36                 clog << "solution of 3*x+5==8 erroneously returned "
  37                      << aux << endl;
  38         }
  39
  40         return result;
  41 }
  42
  43 static unsigned exam_lsolve2a()
  44 {
  45         // An example from the Maple online help.
  46         unsigned result = 0;
  47         symbol a("a"), b("b"), x("x"), y("y");
  48         lst eqns, vars;
  49         ex sol;
  50
  51         // Create the linear system [a*x+b*y==3,x-y==b]...
  52         eqns.append(a*x+b*y==3).append(x-y==b);
  53         // ...to be solved for [x,y]...
  54         vars.append(x).append(y);
  55         // ...and solve it:
  56         sol = lsolve(eqns, vars);
  57         ex sol_x = sol.op(0).rhs();  // rhs of solution for first variable (x)
  58         ex sol_y = sol.op(1).rhs();  // rhs of solution for second variable (y)
  59
  60         // It should have returned [x==(3+b^2)/(a+b),y==(3-a*b)/(a+b)]
  61         if (!normal(sol_x - (3+pow(b,2))/(a+b)).is_zero() ||
  62                 !normal(sol_y - (3-a*b)/(a+b)).is_zero()) {
  63                 ++result;
  64                 clog << "solution of the system " << eqns << " for " << vars
  65                      << " erroneously returned " << sol << endl;
  66         }
  67
  68         return result;
  69 }
  70
  71 static unsigned exam_lsolve2b()
  72 {
  73         // A boring example from Mathematica's online help.
  74         unsigned result = 0;
  75         symbol x("x"), y("y");
  76         lst eqns, vars;
  77         ex sol;
  78
  79         // Create the linear system [3*x+y==7,2*x-5*y==8]...
  80         eqns.append(3*x+y==7).append(2*x-5*y==8);
  81         // ...to be solved for [x,y]...
  82         vars.append(x).append(y);
  83         // ...and solve it:
  84         sol = lsolve(eqns, vars);
  85         ex sol_x = sol.op(0).rhs();  // rhs of solution for first variable (x)
  86         ex sol_y = sol.op(1).rhs();  // rhs of solution for second variable (y)
  87
  88         // It should have returned [x==43/17,y==-10/17]
  89         if ((sol_x != numeric(43,17)) ||
  90                 (sol_y != numeric(-10,17))) {
  91                 ++result;
  92                 clog << "solution of the system " << eqns << " for " << vars
  93                      << " erroneously returned " << sol << endl;
  94         }
  95
  96         return result;
  97 }
  98
  99 static unsigned exam_lsolve2c()
 100 {
 101         // A more interesting example from the Maple online help.
 102         unsigned result = 0;
 103         symbol x("x"), y("y");
 104         lst eqns, vars;
 105         ex sol;
 106
 107         // Create the linear system [I*x+y==1,I*x-y==2]...
 108         eqns.append(I*x+y==1).append(I*x-y==2);
 109         // ...to be solved for [x,y]...
 110         vars.append(x).append(y);
 111         // ...and solve it:
 112         sol = lsolve(eqns, vars);
 113         ex sol_x = sol.op(0).rhs();  // rhs of solution for first variable (x)
 114         ex sol_y = sol.op(1).rhs();  // rhs of solution for second variable (y)
 115
 116         // It should have returned [x==-3/2*I,y==-1/2]
 117         if ((sol_x != numeric(-3,2)*I) ||
 118                 (sol_y != numeric(-1,2))) {
 119                 ++result;
 120                 clog << "solution of the system " << eqns << " for " << vars
 121                      << " erroneously returned " << sol << endl;
 122         }
 123
 124         return result;
 125 }
 126
 127 static unsigned exam_lsolve2S()
 128 {
 129         // A degenerate example that went wrong in GiNaC 0.6.2.
 130         unsigned result = 0;
 131         symbol x("x"), y("y"), t("t");
 132         lst eqns, vars;
 133         ex sol;
 134
 135         // Create the linear system [0*x+0*y==0,0*x+1*y==t]...
 136         eqns.append(0*x+0*y==0).append(0*x+1*y==t);
 137         // ...to be solved for [x,y]...
 138         vars.append(x).append(y);
 139         // ...and solve it:
 140         sol = lsolve(eqns, vars);
 141         ex sol_x = sol.op(0).rhs();  // rhs of solution for first variable (x)
 142         ex sol_y = sol.op(1).rhs();  // rhs of solution for second variable (y)
 143
 144         // It should have returned [x==x,y==t]
 145         if ((sol_x != x) ||
 146                 (sol_y != t)) {
 147                 ++result;
 148                 clog << "solution of the system " << eqns << " for " << vars
 149                      << " erroneously returned " << sol << endl;
 150         }
 151
 152         return result;
 153 }
 154
 155 static unsigned exam_lsolve3S()
 156 {
 157         // A degenerate example that went wrong while trying to improve elimination
 158         unsigned result = 0;
 159         symbol b("b"), c("c");
 160         symbol x("x"), y("y"), z("z");
 161         lst eqns, vars;
 162         ex sol;
 163
 164         // Create the linear system [y+z==b,-y+z==c] with one additional row...
 165         eqns.append(ex(0)==ex(0)).append(b==z+y).append(c==z-y);
 166         // ...to be solved for [x,y,z]...
 167         vars.append(x).append(y).append(z);
 168         // ...and solve it:
 169         sol = lsolve(eqns, vars);
 170         ex sol_x = sol.op(0).rhs();  // rhs of solution for first variable (x)
 171         ex sol_y = sol.op(1).rhs();  // rhs of solution for second variable (y)
 172         ex sol_z = sol.op(2).rhs();  // rhs of solution for third variable (z)
 173
 174         // It should have returned [x==x,y==t,]
 175         if ((sol_x != x) ||
 176                 (sol_y != (b-c)/2) ||
 177                 (sol_z != (b+c)/2)) {
 178                 ++result;
 179                 clog << "solution of the system " << eqns << " for " << vars
 180                      << " erroneously returned " << sol << endl;
 181         }
 182
 183         return result;
 184 }
 185
 186 unsigned exam_lsolve()
 187 {
 188         unsigned result = 0;
 189
 190         cout << "examining linear solve" << flush;
 191         clog << "----------linear solve:" << endl;
 192
 193         result += exam_lsolve1();  cout << '.' << flush;
 194         result += exam_lsolve2a();  cout << '.' << flush;
 195         result += exam_lsolve2b();  cout << '.' << flush;
 196         result += exam_lsolve2c();  cout << '.' << flush;
 197         result += exam_lsolve2S();  cout << '.' << flush;
 198         result += exam_lsolve3S();  cout << '.' << flush;
 199
 200         if (!result) {
 201                 cout << " passed " << endl;
 202                 clog << "(no output)" << endl;
 203         } else {
 204                 cout << " failed " << endl;
 205         }
 206
 207         return result;
 208 }