ginac/matrix.cpp

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