ginac/matrix.cpp

   1 /** @file matrix.cpp
   2  *
   3  *  Implementation of symbolic matrices */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2000 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 <algorithm>
  24 #include <stdexcept>
  25
  26 #include "matrix.h"
  27 #include "archive.h"
  28 #include "utils.h"
  29 #include "debugmsg.h"
  30
  31 #ifndef NO_GINAC_NAMESPACE
  32 namespace GiNaC {
  33 #endif // ndef NO_GINAC_NAMESPACE
  34
  35 GINAC_IMPLEMENT_REGISTERED_CLASS(matrix, basic)
  36
  37 //////////
  38 // default constructor, destructor, copy constructor, assignment operator
  39 // and helpers:
  40 //////////
  41
  42 // public
  43
  44 /** Default ctor.  Initializes to 1 x 1-dimensional zero-matrix. */
  45 matrix::matrix()
  46     : inherited(TINFO_matrix), row(1), col(1)
  47 {
  48     debugmsg("matrix default constructor",LOGLEVEL_CONSTRUCT);
  49     m.push_back(_ex0());
  50 }
  51
  52 matrix::~matrix()
  53 {
  54     debugmsg("matrix destructor",LOGLEVEL_DESTRUCT);
  55 }
  56
  57 matrix::matrix(matrix const & other)
  58 {
  59     debugmsg("matrix copy constructor",LOGLEVEL_CONSTRUCT);
  60     copy(other);
  61 }
  62
  63 matrix const & matrix::operator=(matrix const & other)
  64 {
  65     debugmsg("matrix operator=",LOGLEVEL_ASSIGNMENT);
  66     if (this != &other) {
  67         destroy(1);
  68         copy(other);
  69     }
  70     return *this;
  71 }
  72
  73 // protected
  74
  75 void matrix::copy(matrix const & other)
  76 {
  77     inherited::copy(other);
  78     row=other.row;
  79     col=other.col;
  80     m=other.m;  // use STL's vector copying
  81 }
  82
  83 void matrix::destroy(bool call_parent)
  84 {
  85     if (call_parent) inherited::destroy(call_parent);
  86 }
  87
  88 //////////
  89 // other constructors
  90 //////////
  91
  92 // public
  93
  94 /** Very common ctor.  Initializes to r x c-dimensional zero-matrix.
  95  *
  96  *  @param r number of rows
  97  *  @param c number of cols */
  98 matrix::matrix(unsigned r, unsigned c)
  99     : inherited(TINFO_matrix), row(r), col(c)
 100 {
 101     debugmsg("matrix constructor from unsigned,unsigned",LOGLEVEL_CONSTRUCT);
 102     m.resize(r*c, _ex0());
 103 }
 104
 105 // protected
 106
 107 /** Ctor from representation, for internal use only. */
 108 matrix::matrix(unsigned r, unsigned c, exvector const & m2)
 109     : inherited(TINFO_matrix), row(r), col(c), m(m2)
 110 {
 111     debugmsg("matrix constructor from unsigned,unsigned,exvector",LOGLEVEL_CONSTRUCT);
 112 }
 113
 114 //////////
 115 // archiving
 116 //////////
 117
 118 /** Construct object from archive_node. */
 119 matrix::matrix(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 120 {
 121     debugmsg("matrix constructor from archive_node", LOGLEVEL_CONSTRUCT);
 122     if (!(n.find_unsigned("row", row)) || !(n.find_unsigned("col", col)))
 123         throw (std::runtime_error("unknown matrix dimensions in archive"));
 124     m.reserve(row * col);
 125     for (unsigned int i=0; true; i++) {
 126         ex e;
 127         if (n.find_ex("m", e, sym_lst, i))
 128             m.push_back(e);
 129         else
 130             break;
 131     }
 132 }
 133
 134 /** Unarchive the object. */
 135 ex matrix::unarchive(const archive_node &n, const lst &sym_lst)
 136 {
 137     return (new matrix(n, sym_lst))->setflag(status_flags::dynallocated);
 138 }
 139
 140 /** Archive the object. */
 141 void matrix::archive(archive_node &n) const
 142 {
 143     inherited::archive(n);
 144     n.add_unsigned("row", row);
 145     n.add_unsigned("col", col);
 146     exvector::const_iterator i = m.begin(), iend = m.end();
 147     while (i != iend) {
 148         n.add_ex("m", *i);
 149         i++;
 150     }
 151 }
 152
 153 //////////
 154 // functions overriding virtual functions from bases classes
 155 //////////
 156
 157 // public
 158
 159 basic * matrix::duplicate() const
 160 {
 161     debugmsg("matrix duplicate",LOGLEVEL_DUPLICATE);
 162     return new matrix(*this);
 163 }
 164
 165 void matrix::print(ostream & os, unsigned upper_precedence) const
 166 {
 167     debugmsg("matrix print",LOGLEVEL_PRINT);
 168     os << "[[ ";
 169     for (unsigned r=0; r<row-1; ++r) {
 170         os << "[[";
 171         for (unsigned c=0; c<col-1; ++c) {
 172             os << m[r*col+c] << ",";
 173         }
 174         os << m[col*(r+1)-1] << "]], ";
 175     }
 176     os << "[[";
 177     for (unsigned c=0; c<col-1; ++c) {
 178         os << m[(row-1)*col+c] << ",";
 179     }
 180     os << m[row*col-1] << "]] ]]";
 181 }
 182
 183 void matrix::printraw(ostream & os) const
 184 {
 185     debugmsg("matrix printraw",LOGLEVEL_PRINT);
 186     os << "matrix(" << row << "," << col <<",";
 187     for (unsigned r=0; r<row-1; ++r) {
 188         os << "(";
 189         for (unsigned c=0; c<col-1; ++c) {
 190             os << m[r*col+c] << ",";
 191         }
 192         os << m[col*(r-1)-1] << "),";
 193     }
 194     os << "(";
 195     for (unsigned c=0; c<col-1; ++c) {
 196         os << m[(row-1)*col+c] << ",";
 197     }
 198     os << m[row*col-1] << "))";
 199 }
 200
 201 /** nops is defined to be rows x columns. */
 202 unsigned matrix::nops() const
 203 {
 204     return row*col;
 205 }
 206
 207 /** returns matrix entry at position (i/col, i%col). */
 208 ex & matrix::let_op(int const i)
 209 {
 210     return m[i];
 211 }
 212
 213 /** expands the elements of a matrix entry by entry. */
 214 ex matrix::expand(unsigned options) const
 215 {
 216     exvector tmp(row*col);
 217     for (unsigned i=0; i<row*col; ++i) {
 218         tmp[i]=m[i].expand(options);
 219     }
 220     return matrix(row, col, tmp);
 221 }
 222
 223 /** Search ocurrences.  A matrix 'has' an expression if it is the expression
 224  *  itself or one of the elements 'has' it. */
 225 bool matrix::has(ex const & other) const
 226 {
 227     GINAC_ASSERT(other.bp!=0);
 228
 229     // tautology: it is the expression itself
 230     if (is_equal(*other.bp)) return true;
 231
 232     // search all the elements
 233     for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) {
 234         if ((*r).has(other)) return true;
 235     }
 236     return false;
 237 }
 238
 239 /** evaluate matrix entry by entry. */
 240 ex matrix::eval(int level) const
 241 {
 242     debugmsg("matrix eval",LOGLEVEL_MEMBER_FUNCTION);
 243
 244     // check if we have to do anything at all
 245     if ((level==1)&&(flags & status_flags::evaluated)) {
 246         return *this;
 247     }
 248
 249     // emergency break
 250     if (level == -max_recursion_level) {
 251         throw (std::runtime_error("matrix::eval(): recursion limit exceeded"));
 252     }
 253
 254     // eval() entry by entry
 255     exvector m2(row*col);
 256     --level;
 257     for (unsigned r=0; r<row; ++r) {
 258         for (unsigned c=0; c<col; ++c) {
 259             m2[r*col+c] = m[r*col+c].eval(level);
 260         }
 261     }
 262
 263     return (new matrix(row, col, m2))->setflag(status_flags::dynallocated |
 264                                                status_flags::evaluated );
 265 }
 266
 267 /** evaluate matrix numerically entry by entry. */
 268 ex matrix::evalf(int level) const
 269 {
 270     debugmsg("matrix evalf",LOGLEVEL_MEMBER_FUNCTION);
 271
 272     // check if we have to do anything at all
 273     if (level==1) {
 274         return *this;
 275     }
 276
 277     // emergency break
 278     if (level == -max_recursion_level) {
 279         throw (std::runtime_error("matrix::evalf(): recursion limit exceeded"));
 280     }
 281
 282     // evalf() entry by entry
 283     exvector m2(row*col);
 284     --level;
 285     for (unsigned r=0; r<row; ++r) {
 286         for (unsigned c=0; c<col; ++c) {
 287             m2[r*col+c] = m[r*col+c].evalf(level);
 288         }
 289     }
 290     return matrix(row, col, m2);
 291 }
 292
 293 // protected
 294
 295 int matrix::compare_same_type(basic const & other) const
 296 {
 297     GINAC_ASSERT(is_exactly_of_type(other, matrix));
 298     matrix const & o=static_cast<matrix &>(const_cast<basic &>(other));
 299
 300     // compare number of rows
 301     if (row != o.rows()) {
 302         return row < o.rows() ? -1 : 1;
 303     }
 304
 305     // compare number of columns
 306     if (col != o.cols()) {
 307         return col < o.cols() ? -1 : 1;
 308     }
 309
 310     // equal number of rows and columns, compare individual elements
 311     int cmpval;
 312     for (unsigned r=0; r<row; ++r) {
 313         for (unsigned c=0; c<col; ++c) {
 314             cmpval=((*this)(r,c)).compare(o(r,c));
 315             if (cmpval!=0) return cmpval;
 316         }
 317     }
 318     // all elements are equal => matrices are equal;
 319     return 0;
 320 }
 321
 322 //////////
 323 // non-virtual functions in this class
 324 //////////
 325
 326 // public
 327
 328 /** Sum of matrices.
 329  *
 330  *  @exception logic_error (incompatible matrices) */
 331 matrix matrix::add(matrix const & other) const
 332 {
 333     if (col != other.col || row != other.row) {
 334         throw (std::logic_error("matrix::add(): incompatible matrices"));
 335     }
 336
 337     exvector sum(this->m);
 338     exvector::iterator i;
 339     exvector::const_iterator ci;
 340     for (i=sum.begin(), ci=other.m.begin();
 341          i!=sum.end();
 342          ++i, ++ci) {
 343         (*i) += (*ci);
 344     }
 345     return matrix(row,col,sum);
 346 }
 347
 348 /** Difference of matrices.
 349  *
 350  *  @exception logic_error (incompatible matrices) */
 351 matrix matrix::sub(matrix const & other) const
 352 {
 353     if (col != other.col || row != other.row) {
 354         throw (std::logic_error("matrix::sub(): incompatible matrices"));
 355     }
 356
 357     exvector dif(this->m);
 358     exvector::iterator i;
 359     exvector::const_iterator ci;
 360     for (i=dif.begin(), ci=other.m.begin();
 361          i!=dif.end();
 362          ++i, ++ci) {
 363         (*i) -= (*ci);
 364     }
 365     return matrix(row,col,dif);
 366 }
 367
 368 /** Product of matrices.
 369  *
 370  *  @exception logic_error (incompatible matrices) */
 371 matrix matrix::mul(matrix const & other) const
 372 {
 373     if (col != other.row) {
 374         throw (std::logic_error("matrix::mul(): incompatible matrices"));
 375     }
 376
 377     exvector prod(row*other.col);
 378     for (unsigned i=0; i<row; ++i) {
 379         for (unsigned j=0; j<other.col; ++j) {
 380             for (unsigned l=0; l<col; ++l) {
 381                 prod[i*other.col+j] += m[i*col+l] * other.m[l*other.col+j];
 382             }
 383         }
 384     }
 385     return matrix(row, other.col, prod);
 386 }
 387
 388 /** operator() to access elements.
 389  *
 390  *  @param ro row of element
 391  *  @param co column of element
 392  *  @exception range_error (index out of range) */
 393 ex const & matrix::operator() (unsigned ro, unsigned co) const
 394 {
 395     if (ro<0 || ro>=row || co<0 || co>=col) {
 396         throw (std::range_error("matrix::operator(): index out of range"));
 397     }
 398
 399     return m[ro*col+co];
 400 }
 401
 402 /** Set individual elements manually.
 403  *
 404  *  @exception range_error (index out of range) */
 405 matrix & matrix::set(unsigned ro, unsigned co, ex value)
 406 {
 407     if (ro<0 || ro>=row || co<0 || co>=col) {
 408         throw (std::range_error("matrix::set(): index out of range"));
 409     }
 410
 411     ensure_if_modifiable();
 412     m[ro*col+co]=value;
 413     return *this;
 414 }
 415
 416 /** Transposed of an m x n matrix, producing a new n x m matrix object that
 417  *  represents the transposed. */
 418 matrix matrix::transpose(void) const
 419 {
 420     exvector trans(col*row);
 421
 422     for (unsigned r=0; r<col; ++r) {
 423         for (unsigned c=0; c<row; ++c) {
 424             trans[r*row+c] = m[c*col+r];
 425         }
 426     }
 427     return matrix(col,row,trans);
 428 }
 429
 430 /* Determiant of purely numeric matrix, using pivoting. This routine is only
 431  * called internally by matrix::determinant(). */
 432 ex determinant_numeric(const matrix & M)
 433 {
 434     GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
 435     matrix tmp(M);
 436     ex det=_ex1();
 437     ex piv;
 438
 439     for (unsigned r1=0; r1<M.rows(); ++r1) {
 440         int indx = tmp.pivot(r1);
 441         if (indx == -1) {
 442             return _ex0();
 443         }
 444         if (indx != 0) {
 445             det *= _ex_1();
 446         }
 447         det = det * tmp.m[r1*M.cols()+r1];
 448         for (unsigned r2=r1+1; r2<M.rows(); ++r2) {
 449             piv = tmp.m[r2*M.cols()+r1] / tmp.m[r1*M.cols()+r1];
 450             for (unsigned c=r1+1; c<M.cols(); c++) {
 451                 tmp.m[r2*M.cols()+c] -= piv * tmp.m[r1*M.cols()+c];
 452             }
 453         }
 454     }
 455     return det;
 456 }
 457
 458 // Compute the sign of a permutation of a vector of things, used internally
 459 // by determinant_symbolic_perm() where it is instantiated for int.
 460 template <class T>
 461 int permutation_sign(vector<T> s)
 462 {
 463     if (s.size() < 2)
 464         return 0;
 465     int sigma=1;
 466     for (typename vector<T>::iterator i=s.begin(); i!=s.end()-1; ++i) {
 467         for (typename vector<T>::iterator j=i+1; j!=s.end(); ++j) {
 468             if (*i == *j)
 469                 return 0;
 470             if (*i > *j) {
 471                 iter_swap(i,j);
 472                 sigma = -sigma;
 473             }
 474         }
 475     }
 476     return sigma;
 477 }
 478
 479 /** Determinant built by application of the full permutation group. This
 480  *  routine is only called internally by matrix::determinant(). */
 481 ex determinant_symbolic_perm(const matrix & M)
 482 {
 483     GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
 484
 485     if (M.rows()==1) {  // speed things up
 486         return M(0,0);
 487     }
 488
 489     ex det;
 490     ex term;
 491     vector<unsigned> sigma(M.cols());
 492     for (unsigned i=0; i<M.cols(); ++i) sigma[i]=i;
 493
 494     do {
 495         term = M(sigma[0],0);
 496         for (unsigned i=1; i<M.cols(); ++i) term *= M(sigma[i],i);
 497         det += permutation_sign(sigma)*term;
 498     } while (next_permutation(sigma.begin(), sigma.end()));
 499
 500     return det;
 501 }
 502
 503 /** Recursive determiant for small matrices having at least one symbolic entry.
 504  *  This algorithm is also known as Laplace-expansion. This routine is only
 505  *  called internally by matrix::determinant(). */
 506 ex determinant_symbolic_minor(const matrix & M)
 507 {
 508     GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
 509
 510     if (M.rows()==1) {  // end of recursion
 511         return M(0,0);
 512     }
 513     if (M.rows()==2) {  // speed things up
 514         return (M(0,0)*M(1,1)-
 515                 M(1,0)*M(0,1));
 516     }
 517     if (M.rows()==3) {  // speed things up even a little more
 518         return ((M(2,1)*M(0,2)-M(2,2)*M(0,1))*M(1,0)+
 519                 (M(1,2)*M(0,1)-M(1,1)*M(0,2))*M(2,0)+
 520                 (M(2,2)*M(1,1)-M(2,1)*M(1,2))*M(0,0));
 521     }
 522
 523     ex det;
 524     matrix minorM(M.rows()-1,M.cols()-1);
 525     for (unsigned r1=0; r1<M.rows(); ++r1) {
 526         // assemble the minor matrix
 527         for (unsigned r=0; r<minorM.rows(); ++r) {
 528             for (unsigned c=0; c<minorM.cols(); ++c) {
 529                 if (r<r1) {
 530                     minorM.set(r,c,M(r,c+1));
 531                 } else {
 532                     minorM.set(r,c,M(r+1,c+1));
 533                 }
 534             }
 535         }
 536         // recurse down
 537         if (r1%2) {
 538             det -= M(r1,0) * determinant_symbolic_minor(minorM);
 539         } else {
 540             det += M(r1,0) * determinant_symbolic_minor(minorM);
 541         }
 542     }
 543     return det;
 544 }
 545
 546 /*  Leverrier algorithm for large matrices having at least one symbolic entry.
 547  *  This routine is only called internally by matrix::determinant(). The
 548  *  algorithm is deemed bad for symbolic matrices since it returns expressions
 549  *  that are very hard to canonicalize. */
 550 /*ex determinant_symbolic_leverrier(const matrix & M)
 551  *{
 552  *    GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
 553  *
 554  *    matrix B(M);
 555  *    matrix I(M.row, M.col);
 556  *    ex c=B.trace();
 557  *    for (unsigned i=1; i<M.row; ++i) {
 558  *        for (unsigned j=0; j<M.row; ++j)
 559  *            I.m[j*M.col+j] = c;
 560  *        B = M.mul(B.sub(I));
 561  *        c = B.trace()/ex(i+1);
 562  *    }
 563  *    if (M.row%2) {
 564  *        return c;
 565  *    } else {
 566  *        return -c;
 567  *    }
 568  *}*/
 569
 570 /** Determinant of square matrix.  This routine doesn't actually calculate the
 571  *  determinant, it only implements some heuristics about which algorithm to
 572  *  call.  When the parameter for normalization is explicitly turned off this
 573  *  method does not normalize its result at the end, which might imply that
 574  *  the symbolic 2x2 matrix [[a/(a-b),1],[b/(a-b),1]] is not immediatly
 575  *  recognized to be unity.  (This is Mathematica's default behaviour, it
 576  *  should be used with care.)
 577  *
 578  *  @param     normalized may be set to false if no normalization of the
 579  *             result is desired (i.e. to force Mathematica behavior, Maple
 580  *             does normalize the result).
 581  *  @return    the determinant as a new expression
 582  *  @exception logic_error (matrix not square) */
 583 ex matrix::determinant(bool normalized) const
 584 {
 585     if (row != col) {
 586         throw (std::logic_error("matrix::determinant(): matrix not square"));
 587     }
 588
 589     // check, if there are non-numeric entries in the matrix:
 590     for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) {
 591         if (!(*r).info(info_flags::numeric)) {
 592             if (normalized) {
 593                 return determinant_symbolic_minor(*this).normal();
 594             } else {
 595                 return determinant_symbolic_perm(*this);
 596             }
 597         }
 598     }
 599     // if it turns out that all elements are numeric
 600     return determinant_numeric(*this);
 601 }
 602
 603 /** Trace of a matrix.
 604  *
 605  *  @return    the sum of diagonal elements
 606  *  @exception logic_error (matrix not square) */
 607 ex matrix::trace(void) const
 608 {
 609     if (row != col) {
 610         throw (std::logic_error("matrix::trace(): matrix not square"));
 611     }
 612
 613     ex tr;
 614     for (unsigned r=0; r<col; ++r) {
 615         tr += m[r*col+r];
 616     }
 617     return tr;
 618 }
 619
 620 /** Characteristic Polynomial.  The characteristic polynomial of a matrix M is
 621  *  defined as the determiant of (M - lambda * 1) where 1 stands for the unit
 622  *  matrix of the same dimension as M.  This method returns the characteristic
 623  *  polynomial as a new expression.
 624  *
 625  *  @return    characteristic polynomial as new expression
 626  *  @exception logic_error (matrix not square)
 627  *  @see       matrix::determinant() */
 628 ex matrix::charpoly(ex const & lambda) const
 629 {
 630     if (row != col) {
 631         throw (std::logic_error("matrix::charpoly(): matrix not square"));
 632     }
 633
 634     matrix M(*this);
 635     for (unsigned r=0; r<col; ++r) {
 636         M.m[r*col+r] -= lambda;
 637     }
 638     return (M.determinant());
 639 }
 640
 641 /** Inverse of this matrix.
 642  *
 643  *  @return    the inverted matrix
 644  *  @exception logic_error (matrix not square)
 645  *  @exception runtime_error (singular matrix) */
 646 matrix matrix::inverse(void) const
 647 {
 648     if (row != col) {
 649         throw (std::logic_error("matrix::inverse(): matrix not square"));
 650     }
 651
 652     matrix tmp(row,col);
 653     // set tmp to the unit matrix
 654     for (unsigned i=0; i<col; ++i) {
 655         tmp.m[i*col+i] = _ex1();
 656     }
 657     // create a copy of this matrix
 658     matrix cpy(*this);
 659     for (unsigned r1=0; r1<row; ++r1) {
 660         int indx = cpy.pivot(r1);
 661         if (indx == -1) {
 662             throw (std::runtime_error("matrix::inverse(): singular matrix"));
 663         }
 664         if (indx != 0) {  // swap rows r and indx of matrix tmp
 665             for (unsigned i=0; i<col; ++i) {
 666                 tmp.m[r1*col+i].swap(tmp.m[indx*col+i]);
 667             }
 668         }
 669         ex a1 = cpy.m[r1*col+r1];
 670         for (unsigned c=0; c<col; ++c) {
 671             cpy.m[r1*col+c] /= a1;
 672             tmp.m[r1*col+c] /= a1;
 673         }
 674         for (unsigned r2=0; r2<row; ++r2) {
 675             if (r2 != r1) {
 676                 ex a2 = cpy.m[r2*col+r1];
 677                 for (unsigned c=0; c<col; ++c) {
 678                     cpy.m[r2*col+c] -= a2 * cpy.m[r1*col+c];
 679                     tmp.m[r2*col+c] -= a2 * tmp.m[r1*col+c];
 680                 }
 681             }
 682         }
 683     }
 684     return tmp;
 685 }
 686
 687 void matrix::ffe_swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2)
 688 {
 689     ensure_if_modifiable();
 690
 691     ex tmp=ffe_get(r1,c1);
 692     ffe_set(r1,c1,ffe_get(r2,c2));
 693     ffe_set(r2,c2,tmp);
 694 }
 695
 696 void matrix::ffe_set(unsigned r, unsigned c, ex e)
 697 {
 698     set(r-1,c-1,e);
 699 }
 700
 701 ex matrix::ffe_get(unsigned r, unsigned c) const
 702 {
 703     return operator()(r-1,c-1);
 704 }
 705
 706 /** Solve a set of equations for an m x n matrix by fraction-free Gaussian
 707  *  elimination. Based on algorithm 9.1 from 'Algorithms for Computer Algebra'
 708  *  by Keith O. Geddes et al.
 709  *
 710  *  @param vars n x p matrix
 711  *  @param rhs m x p matrix
 712  *  @exception logic_error (incompatible matrices)
 713  *  @exception runtime_error (singular matrix) */
 714 matrix matrix::fraction_free_elim(matrix const & vars,
 715                                   matrix const & rhs) const
 716 {
 717     if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col)) {
 718         throw (std::logic_error("matrix::solve(): incompatible matrices"));
 719     }
 720
 721     matrix a(*this); // make a copy of the matrix
 722     matrix b(rhs);     // make a copy of the rhs vector
 723
 724     // given an m x n matrix a, reduce it to upper echelon form
 725     unsigned m=a.row;
 726     unsigned n=a.col;
 727     int sign=1;
 728     ex divisor=1;
 729     unsigned r=1;
 730
 731     // eliminate below row r, with pivot in column k
 732     for (unsigned k=1; (k<=n)&&(r<=m); ++k) {
 733         // find a nonzero pivot
 734         unsigned p;
 735         for (p=r; (p<=m)&&(a.ffe_get(p,k).is_equal(_ex0())); ++p) {}
 736         // pivot is in row p
 737         if (p<=m) {
 738             if (p!=r) {
 739                 // switch rows p and r
 740                 for (unsigned j=k; j<=n; ++j) {
 741                     a.ffe_swap(p,j,r,j);
 742                 }
 743                 b.ffe_swap(p,1,r,1);
 744                 // keep track of sign changes due to row exchange
 745                 sign=-sign;
 746             }
 747             for (unsigned i=r+1; i<=m; ++i) {
 748                 for (unsigned j=k+1; j<=n; ++j) {
 749                     a.ffe_set(i,j,(a.ffe_get(r,k)*a.ffe_get(i,j)
 750                                   -a.ffe_get(r,j)*a.ffe_get(i,k))/divisor);
 751                     a.ffe_set(i,j,a.ffe_get(i,j).normal() /*.normal() */ );
 752                 }
 753                 b.ffe_set(i,1,(a.ffe_get(r,k)*b.ffe_get(i,1)
 754                               -b.ffe_get(r,1)*a.ffe_get(i,k))/divisor);
 755                 b.ffe_set(i,1,b.ffe_get(i,1).normal() /*.normal() */ );
 756                 a.ffe_set(i,k,0);
 757             }
 758             divisor=a.ffe_get(r,k);
 759             r++;
 760         }
 761     }
 762     // optionally compute the determinant for square or augmented matrices
 763     // if (r==m+1) { det=sign*divisor; } else { det=0; }
 764
 765     /*
 766     for (unsigned r=1; r<=m; ++r) {
 767         for (unsigned c=1; c<=n; ++c) {
 768             cout << a.ffe_get(r,c) << "\t";
 769         }
 770         cout << " | " <<  b.ffe_get(r,1) << endl;
 771     }
 772     */
 773
 774 #ifdef DO_GINAC_ASSERT
 775     // test if we really have an upper echelon matrix
 776     int zero_in_last_row=-1;
 777     for (unsigned r=1; r<=m; ++r) {
 778         int zero_in_this_row=0;
 779         for (unsigned c=1; c<=n; ++c) {
 780             if (a.ffe_get(r,c).is_equal(_ex0())) {
 781                zero_in_this_row++;
 782             } else {
 783                 break;
 784             }
 785         }
 786         GINAC_ASSERT((zero_in_this_row>zero_in_last_row)||(zero_in_this_row=n));
 787         zero_in_last_row=zero_in_this_row;
 788     }
 789 #endif // def DO_GINAC_ASSERT
 790
 791     // assemble solution
 792     matrix sol(n,1);
 793     unsigned last_assigned_sol=n+1;
 794     for (unsigned r=m; r>0; --r) {
 795         unsigned first_non_zero=1;
 796         while ((first_non_zero<=n)&&(a.ffe_get(r,first_non_zero).is_zero())) {
 797             first_non_zero++;
 798         }
 799         if (first_non_zero>n) {
 800             // row consists only of zeroes, corresponding rhs must be 0 as well
 801             if (!b.ffe_get(r,1).is_zero()) {
 802                 throw (std::runtime_error("matrix::fraction_free_elim(): singular matrix"));
 803             }
 804         } else {
 805             // assign solutions for vars between first_non_zero+1 and
 806             // last_assigned_sol-1: free parameters
 807             for (unsigned c=first_non_zero+1; c<=last_assigned_sol-1; ++c) {
 808                 sol.ffe_set(c,1,vars.ffe_get(c,1));
 809             }
 810             ex e=b.ffe_get(r,1);
 811             for (unsigned c=first_non_zero+1; c<=n; ++c) {
 812                 e=e-a.ffe_get(r,c)*sol.ffe_get(c,1);
 813             }
 814             sol.ffe_set(first_non_zero,1,
 815                         (e/a.ffe_get(r,first_non_zero)).normal());
 816             last_assigned_sol=first_non_zero;
 817         }
 818     }
 819     // assign solutions for vars between 1 and
 820     // last_assigned_sol-1: free parameters
 821     for (unsigned c=1; c<=last_assigned_sol-1; ++c) {
 822         sol.ffe_set(c,1,vars.ffe_get(c,1));
 823     }
 824
 825     /*
 826     for (unsigned c=1; c<=n; ++c) {
 827         cout << vars.ffe_get(c,1) << "->" << sol.ffe_get(c,1) << endl;
 828     }
 829     */
 830
 831 #ifdef DO_GINAC_ASSERT
 832     // test solution with echelon matrix
 833     for (unsigned r=1; r<=m; ++r) {
 834         ex e=0;
 835         for (unsigned c=1; c<=n; ++c) {
 836             e=e+a.ffe_get(r,c)*sol.ffe_get(c,1);
 837         }
 838         if (!(e-b.ffe_get(r,1)).normal().is_zero()) {
 839             cout << "e=" << e;
 840             cout << "b.ffe_get(" << r<<",1)=" << b.ffe_get(r,1) << endl;
 841             cout << "diff=" << (e-b.ffe_get(r,1)).normal() << endl;
 842         }
 843         GINAC_ASSERT((e-b.ffe_get(r,1)).normal().is_zero());
 844     }
 845
 846     // test solution with original matrix
 847     for (unsigned r=1; r<=m; ++r) {
 848         ex e=0;
 849         for (unsigned c=1; c<=n; ++c) {
 850             e=e+ffe_get(r,c)*sol.ffe_get(c,1);
 851         }
 852         try {
 853         if (!(e-rhs.ffe_get(r,1)).normal().is_zero()) {
 854             cout << "e=" << e << endl;
 855             e.printtree(cout);
 856             ex en=e.normal();
 857             cout << "e.normal()=" << en << endl;
 858             en.printtree(cout);
 859             cout << "rhs.ffe_get(" << r<<",1)=" << rhs.ffe_get(r,1) << endl;
 860             cout << "diff=" << (e-rhs.ffe_get(r,1)).normal() << endl;
 861         }
 862         } catch (...) {
 863             ex xxx=e-rhs.ffe_get(r,1);
 864             cerr << "xxx=" << xxx << endl << endl;
 865         }
 866         GINAC_ASSERT((e-rhs.ffe_get(r,1)).normal().is_zero());
 867     }
 868 #endif // def DO_GINAC_ASSERT
 869
 870     return sol;
 871 }
 872
 873 /** Solve simultaneous set of equations. */
 874 matrix matrix::solve(matrix const & v) const
 875 {
 876     if (!(row == col && col == v.row)) {
 877         throw (std::logic_error("matrix::solve(): incompatible matrices"));
 878     }
 879
 880     // build the extended matrix of *this with v attached to the right
 881     matrix tmp(row,col+v.col);
 882     for (unsigned r=0; r<row; ++r) {
 883         for (unsigned c=0; c<col; ++c) {
 884             tmp.m[r*tmp.col+c] = m[r*col+c];
 885         }
 886         for (unsigned c=0; c<v.col; ++c) {
 887             tmp.m[r*tmp.col+c+col] = v.m[r*v.col+c];
 888         }
 889     }
 890     for (unsigned r1=0; r1<row; ++r1) {
 891         int indx = tmp.pivot(r1);
 892         if (indx == -1) {
 893             throw (std::runtime_error("matrix::solve(): singular matrix"));
 894         }
 895         for (unsigned c=r1; c<tmp.col; ++c) {
 896             tmp.m[r1*tmp.col+c] /= tmp.m[r1*tmp.col+r1];
 897         }
 898         for (unsigned r2=r1+1; r2<row; ++r2) {
 899             for (unsigned c=r1; c<tmp.col; ++c) {
 900                 tmp.m[r2*tmp.col+c]
 901                     -= tmp.m[r2*tmp.col+r1] * tmp.m[r1*tmp.col+c];
 902             }
 903         }
 904     }
 905
 906     // assemble the solution matrix
 907     exvector sol(v.row*v.col);
 908     for (unsigned c=0; c<v.col; ++c) {
 909         for (unsigned r=col-1; r>=0; --r) {
 910             sol[r*v.col+c] = tmp[r*tmp.col+c];
 911             for (unsigned i=r+1; i<col; ++i) {
 912                 sol[r*v.col+c]
 913                     -= tmp[r*tmp.col+i] * sol[i*v.col+c];
 914             }
 915         }
 916     }
 917     return matrix(v.row, v.col, sol);
 918 }
 919
 920 // protected
 921
 922 /** Partial pivoting method.
 923  *  Usual pivoting returns the index to the element with the largest absolute
 924  *  value and swaps the current row with the one where the element was found.
 925  *  Here it does the same with the first non-zero element. (This works fine,
 926  *  but may be far from optimal for numerics.) */
 927 int matrix::pivot(unsigned ro)
 928 {
 929     unsigned k=ro;
 930
 931     for (unsigned r=ro; r<row; ++r) {
 932         if (!m[r*col+ro].is_zero()) {
 933             k = r;
 934             break;
 935         }
 936     }
 937     if (m[k*col+ro].is_zero()) {
 938         return -1;
 939     }
 940     if (k!=ro) {  // swap rows
 941         for (unsigned c=0; c<col; ++c) {
 942             m[k*col+c].swap(m[ro*col+c]);
 943         }
 944         return k;
 945     }
 946     return 0;
 947 }
 948
 949 //////////
 950 // global constants
 951 //////////
 952
 953 const matrix some_matrix;
 954 type_info const & typeid_matrix=typeid(some_matrix);
 955
 956 #ifndef NO_GINAC_NAMESPACE
 957 } // namespace GiNaC
 958 #endif // ndef NO_GINAC_NAMESPACE