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