GiNaC  1.8.3
matrix.cpp
Go to the documentation of this file.
1 
5 /*
6  * GiNaC Copyright (C) 1999-2022 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21  */
22 
23 #include "matrix.h"
24 #include "numeric.h"
25 #include "lst.h"
26 #include "idx.h"
27 #include "indexed.h"
28 #include "add.h"
29 #include "power.h"
30 #include "symbol.h"
31 #include "operators.h"
32 #include "normal.h"
33 #include "archive.h"
34 #include "utils.h"
35 
36 #include <algorithm>
37 #include <iostream>
38 #include <map>
39 #include <sstream>
40 #include <stdexcept>
41 #include <string>
42 
43 namespace GiNaC {
44 
47  print_func<print_latex>(&matrix::do_print_latex).
48  print_func<print_tree>(&matrix::do_print_tree).
49  print_func<print_python_repr>(&matrix::do_print_python_repr))
50 
51 // default constructor
54 
56 matrix::matrix() : row(1), col(1), m(1, _ex0)
57 {
59 }
60 
62 // other constructors
64 
65 // public
66 
71 matrix::matrix(unsigned r, unsigned c) : row(r), col(c), m(r*c, _ex0)
72 {
74 }
75 
80 matrix::matrix(unsigned r, unsigned c, const lst & l)
81  : row(r), col(c), m(r*c, _ex0)
82 {
84 
85  size_t i = 0;
86  for (auto & it : l) {
87  size_t x = i % c;
88  size_t y = i / c;
89  if (y >= r)
90  break; // matrix smaller than list: throw away excessive elements
91  m[y*c+x] = it;
92  ++i;
93  }
94 }
95 
99 matrix::matrix(std::initializer_list<std::initializer_list<ex>> l)
100  : row(l.size()), col(l.begin()->size())
101 {
103 
104  m.reserve(row*col);
105  for (const auto & r : l) {
106  unsigned c = 0;
107  for (const auto & e : r) {
108  m.push_back(e);
109  ++c;
110  }
111  if (c != col)
112  throw std::invalid_argument("matrix::matrix{{}}: wrong dimension");
113  }
114 }
115 
116 // protected
117 
119 matrix::matrix(unsigned r, unsigned c, const exvector & m2)
120  : row(r), col(c), m(m2)
121 {
123 }
124 matrix::matrix(unsigned r, unsigned c, exvector && m2)
125  : row(r), col(c), m(std::move(m2))
126 {
128 }
129 
131 // archiving
133 
134 void matrix::read_archive(const archive_node &n, lst &sym_lst)
135 {
136  inherited::read_archive(n, sym_lst);
137 
138  if (!(n.find_unsigned("row", row)) || !(n.find_unsigned("col", col)))
139  throw (std::runtime_error("unknown matrix dimensions in archive"));
140  m.reserve(row * col);
141  // XXX: default ctor inserts a zero element, we need to erase it here.
142  m.pop_back();
143  auto range = n.find_property_range("m", "m");
144  for (auto i=range.begin; i != range.end; ++i) {
145  ex e;
146  n.find_ex_by_loc(i, e, sym_lst);
147  m.emplace_back(e);
148  }
149 }
151 
153 {
154  inherited::archive(n);
155  n.add_unsigned("row", row);
156  n.add_unsigned("col", col);
157  for (auto & i : m) {
158  n.add_ex("m", i);
159  }
160 }
161 
163 // functions overriding virtual functions from base classes
165 
166 // public
167 
168 void matrix::print_elements(const print_context & c, const char *row_start, const char *row_end, const char *row_sep, const char *col_sep) const
169 {
170  for (unsigned ro=0; ro<row; ++ro) {
171  c.s << row_start;
172  for (unsigned co=0; co<col; ++co) {
173  m[ro*col+co].print(c);
174  if (co < col-1)
175  c.s << col_sep;
176  else
177  c.s << row_end;
178  }
179  if (ro < row-1)
180  c.s << row_sep;
181  }
182 }
183 
184 void matrix::do_print(const print_context & c, unsigned level) const
185 {
186  c.s << "[";
187  print_elements(c, "[", "]", ",", ",");
188  c.s << "]";
189 }
190 
191 void matrix::do_print_latex(const print_latex & c, unsigned level) const
192 {
193  c.s << "\\left(\\begin{array}{" << std::string(col,'c') << "}";
194  print_elements(c, "", "", "\\\\", "&");
195  c.s << "\\end{array}\\right)";
196 }
197 
198 void matrix::do_print_python_repr(const print_python_repr & c, unsigned level) const
199 {
200  c.s << class_name() << '(';
201  print_elements(c, "[", "]", ",", ",");
202  c.s << ')';
203 }
204 
206 size_t matrix::nops() const
207 {
208  return static_cast<size_t>(row) * static_cast<size_t>(col);
209 }
210 
212 ex matrix::op(size_t i) const
213 {
214  GINAC_ASSERT(i<nops());
215 
216  return m[i];
217 }
218 
220 ex & matrix::let_op(size_t i)
221 {
222  GINAC_ASSERT(i<nops());
223 
225  return m[i];
226 }
227 
228 ex matrix::subs(const exmap & mp, unsigned options) const
229 {
230  exvector m2(row * col);
231  for (unsigned r=0; r<row; ++r)
232  for (unsigned c=0; c<col; ++c)
233  m2[r*col+c] = m[r*col+c].subs(mp, options);
234 
235  return matrix(row, col, std::move(m2)).subs_one_level(mp, options);
236 }
237 
240 {
241  std::unique_ptr<exvector> ev(nullptr);
242  for (auto i=m.begin(); i!=m.end(); ++i) {
243  ex x = i->conjugate();
244  if (ev) {
245  ev->push_back(x);
246  continue;
247  }
248  if (are_ex_trivially_equal(x, *i)) {
249  continue;
250  }
251  ev.reset(new exvector);
252  ev->reserve(m.size());
253  for (auto j=m.begin(); j!=i; ++j) {
254  ev->push_back(*j);
255  }
256  ev->push_back(x);
257  }
258  if (ev) {
259  return matrix(row, col, std::move(*ev));
260  }
261  return *this;
262 }
263 
265 {
266  exvector v;
267  v.reserve(m.size());
268  for (auto & i : m)
269  v.push_back(i.real_part());
270  return matrix(row, col, std::move(v));
271 }
272 
274 {
275  exvector v;
276  v.reserve(m.size());
277  for (auto & i : m)
278  v.push_back(i.imag_part());
279  return matrix(row, col, std::move(v));
280 }
281 
282 // protected
283 
284 int matrix::compare_same_type(const basic & other) const
285 {
286  GINAC_ASSERT(is_exactly_a<matrix>(other));
287  const matrix &o = static_cast<const matrix &>(other);
288 
289  // compare number of rows
290  if (row != o.rows())
291  return row < o.rows() ? -1 : 1;
292 
293  // compare number of columns
294  if (col != o.cols())
295  return col < o.cols() ? -1 : 1;
296 
297  // equal number of rows and columns, compare individual elements
298  int cmpval;
299  for (unsigned r=0; r<row; ++r) {
300  for (unsigned c=0; c<col; ++c) {
301  cmpval = ((*this)(r,c)).compare(o(r,c));
302  if (cmpval!=0) return cmpval;
303  }
304  }
305  // all elements are equal => matrices are equal;
306  return 0;
307 }
308 
309 bool matrix::match_same_type(const basic & other) const
310 {
311  GINAC_ASSERT(is_exactly_a<matrix>(other));
312  const matrix & o = static_cast<const matrix &>(other);
313 
314  // The number of rows and columns must be the same. This is necessary to
315  // prevent a 2x3 matrix from matching a 3x2 one.
316  return row == o.rows() && col == o.cols();
317 }
318 
320 ex matrix::eval_indexed(const basic & i) const
321 {
322  GINAC_ASSERT(is_a<indexed>(i));
323  GINAC_ASSERT(is_a<matrix>(i.op(0)));
324 
325  bool all_indices_unsigned = static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint);
326 
327  // Check indices
328  if (i.nops() == 2) {
329 
330  // One index, must be one-dimensional vector
331  if (row != 1 && col != 1)
332  throw (std::runtime_error("matrix::eval_indexed(): vector must have exactly 1 index"));
333 
334  const idx & i1 = ex_to<idx>(i.op(1));
335 
336  if (col == 1) {
337 
338  // Column vector
339  if (!i1.get_dim().is_equal(row))
340  throw (std::runtime_error("matrix::eval_indexed(): dimension of index must match number of vector elements"));
341 
342  // Index numeric -> return vector element
343  if (all_indices_unsigned) {
344  unsigned n1 = ex_to<numeric>(i1.get_value()).to_int();
345  if (n1 >= row)
346  throw (std::runtime_error("matrix::eval_indexed(): value of index exceeds number of vector elements"));
347  return (*this)(n1, 0);
348  }
349 
350  } else {
351 
352  // Row vector
353  if (!i1.get_dim().is_equal(col))
354  throw (std::runtime_error("matrix::eval_indexed(): dimension of index must match number of vector elements"));
355 
356  // Index numeric -> return vector element
357  if (all_indices_unsigned) {
358  unsigned n1 = ex_to<numeric>(i1.get_value()).to_int();
359  if (n1 >= col)
360  throw (std::runtime_error("matrix::eval_indexed(): value of index exceeds number of vector elements"));
361  return (*this)(0, n1);
362  }
363  }
364 
365  } else if (i.nops() == 3) {
366 
367  // Two indices
368  const idx & i1 = ex_to<idx>(i.op(1));
369  const idx & i2 = ex_to<idx>(i.op(2));
370 
371  if (!i1.get_dim().is_equal(row))
372  throw (std::runtime_error("matrix::eval_indexed(): dimension of first index must match number of rows"));
373  if (!i2.get_dim().is_equal(col))
374  throw (std::runtime_error("matrix::eval_indexed(): dimension of second index must match number of columns"));
375 
376  // Pair of dummy indices -> compute trace
377  if (is_dummy_pair(i1, i2))
378  return trace();
379 
380  // Both indices numeric -> return matrix element
381  if (all_indices_unsigned) {
382  unsigned n1 = ex_to<numeric>(i1.get_value()).to_int(), n2 = ex_to<numeric>(i2.get_value()).to_int();
383  if (n1 >= row)
384  throw (std::runtime_error("matrix::eval_indexed(): value of first index exceeds number of rows"));
385  if (n2 >= col)
386  throw (std::runtime_error("matrix::eval_indexed(): value of second index exceeds number of columns"));
387  return (*this)(n1, n2);
388  }
389 
390  } else
391  throw (std::runtime_error("matrix::eval_indexed(): matrix must have exactly 2 indices"));
392 
393  return i.hold();
394 }
395 
397 ex matrix::add_indexed(const ex & self, const ex & other) const
398 {
399  GINAC_ASSERT(is_a<indexed>(self));
400  GINAC_ASSERT(is_a<matrix>(self.op(0)));
401  GINAC_ASSERT(is_a<indexed>(other));
402  GINAC_ASSERT(self.nops() == 2 || self.nops() == 3);
403 
404  // Only add two matrices
405  if (is_a<matrix>(other.op(0))) {
406  GINAC_ASSERT(other.nops() == 2 || other.nops() == 3);
407 
408  const matrix &self_matrix = ex_to<matrix>(self.op(0));
409  const matrix &other_matrix = ex_to<matrix>(other.op(0));
410 
411  if (self.nops() == 2 && other.nops() == 2) { // vector + vector
412 
413  if (self_matrix.row == other_matrix.row)
414  return indexed(self_matrix.add(other_matrix), self.op(1));
415  else if (self_matrix.row == other_matrix.col)
416  return indexed(self_matrix.add(other_matrix.transpose()), self.op(1));
417 
418  } else if (self.nops() == 3 && other.nops() == 3) { // matrix + matrix
419 
420  if (self.op(1).is_equal(other.op(1)) && self.op(2).is_equal(other.op(2)))
421  return indexed(self_matrix.add(other_matrix), self.op(1), self.op(2));
422  else if (self.op(1).is_equal(other.op(2)) && self.op(2).is_equal(other.op(1)))
423  return indexed(self_matrix.add(other_matrix.transpose()), self.op(1), self.op(2));
424 
425  }
426  }
427 
428  // Don't know what to do, return unevaluated sum
429  return self + other;
430 }
431 
433 ex matrix::scalar_mul_indexed(const ex & self, const numeric & other) const
434 {
435  GINAC_ASSERT(is_a<indexed>(self));
436  GINAC_ASSERT(is_a<matrix>(self.op(0)));
437  GINAC_ASSERT(self.nops() == 2 || self.nops() == 3);
438 
439  const matrix &self_matrix = ex_to<matrix>(self.op(0));
440 
441  if (self.nops() == 2)
442  return indexed(self_matrix.mul(other), self.op(1));
443  else // self.nops() == 3
444  return indexed(self_matrix.mul(other), self.op(1), self.op(2));
445 }
446 
448 bool matrix::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
449 {
450  GINAC_ASSERT(is_a<indexed>(*self));
451  GINAC_ASSERT(is_a<indexed>(*other));
452  GINAC_ASSERT(self->nops() == 2 || self->nops() == 3);
453  GINAC_ASSERT(is_a<matrix>(self->op(0)));
454 
455  // Only contract with other matrices
456  if (!is_a<matrix>(other->op(0)))
457  return false;
458 
459  GINAC_ASSERT(other->nops() == 2 || other->nops() == 3);
460 
461  const matrix &self_matrix = ex_to<matrix>(self->op(0));
462  const matrix &other_matrix = ex_to<matrix>(other->op(0));
463 
464  if (self->nops() == 2) {
465 
466  if (other->nops() == 2) { // vector * vector (scalar product)
467 
468  if (self_matrix.col == 1) {
469  if (other_matrix.col == 1) {
470  // Column vector * column vector, transpose first vector
471  *self = self_matrix.transpose().mul(other_matrix)(0, 0);
472  } else {
473  // Column vector * row vector, swap factors
474  *self = other_matrix.mul(self_matrix)(0, 0);
475  }
476  } else {
477  if (other_matrix.col == 1) {
478  // Row vector * column vector, perfect
479  *self = self_matrix.mul(other_matrix)(0, 0);
480  } else {
481  // Row vector * row vector, transpose second vector
482  *self = self_matrix.mul(other_matrix.transpose())(0, 0);
483  }
484  }
485  *other = _ex1;
486  return true;
487 
488  } else { // vector * matrix
489 
490  // B_i * A_ij = (B*A)_j (B is row vector)
491  if (is_dummy_pair(self->op(1), other->op(1))) {
492  if (self_matrix.row == 1)
493  *self = indexed(self_matrix.mul(other_matrix), other->op(2));
494  else
495  *self = indexed(self_matrix.transpose().mul(other_matrix), other->op(2));
496  *other = _ex1;
497  return true;
498  }
499 
500  // B_j * A_ij = (A*B)_i (B is column vector)
501  if (is_dummy_pair(self->op(1), other->op(2))) {
502  if (self_matrix.col == 1)
503  *self = indexed(other_matrix.mul(self_matrix), other->op(1));
504  else
505  *self = indexed(other_matrix.mul(self_matrix.transpose()), other->op(1));
506  *other = _ex1;
507  return true;
508  }
509  }
510 
511  } else if (other->nops() == 3) { // matrix * matrix
512 
513  // A_ij * B_jk = (A*B)_ik
514  if (is_dummy_pair(self->op(2), other->op(1))) {
515  *self = indexed(self_matrix.mul(other_matrix), self->op(1), other->op(2));
516  *other = _ex1;
517  return true;
518  }
519 
520  // A_ij * B_kj = (A*Btrans)_ik
521  if (is_dummy_pair(self->op(2), other->op(2))) {
522  *self = indexed(self_matrix.mul(other_matrix.transpose()), self->op(1), other->op(1));
523  *other = _ex1;
524  return true;
525  }
526 
527  // A_ji * B_jk = (Atrans*B)_ik
528  if (is_dummy_pair(self->op(1), other->op(1))) {
529  *self = indexed(self_matrix.transpose().mul(other_matrix), self->op(2), other->op(2));
530  *other = _ex1;
531  return true;
532  }
533 
534  // A_ji * B_kj = (B*A)_ki
535  if (is_dummy_pair(self->op(1), other->op(2))) {
536  *self = indexed(other_matrix.mul(self_matrix), other->op(1), self->op(2));
537  *other = _ex1;
538  return true;
539  }
540  }
541 
542  return false;
543 }
544 
545 
547 // non-virtual functions in this class
549 
550 // public
551 
555 matrix matrix::add(const matrix & other) const
556 {
557  if (col != other.col || row != other.row)
558  throw std::logic_error("matrix::add(): incompatible matrices");
559 
560  exvector sum(this->m);
561  auto ci = other.m.begin();
562  for (auto & i : sum)
563  i += *ci++;
564 
565  return matrix(row, col, std::move(sum));
566 }
567 
568 
572 matrix matrix::sub(const matrix & other) const
573 {
574  if (col != other.col || row != other.row)
575  throw std::logic_error("matrix::sub(): incompatible matrices");
576 
577  exvector dif(this->m);
578  auto ci = other.m.begin();
579  for (auto & i : dif)
580  i -= *ci++;
581 
582  return matrix(row, col, std::move(dif));
583 }
584 
585 
589 matrix matrix::mul(const matrix & other) const
590 {
591  if (this->cols() != other.rows())
592  throw std::logic_error("matrix::mul(): incompatible matrices");
593 
594  exvector prod(this->rows()*other.cols());
595 
596  for (unsigned r1=0; r1<this->rows(); ++r1) {
597  for (unsigned c=0; c<this->cols(); ++c) {
598  // Quick test: can we shortcut?
599  if (m[r1*col+c].is_zero())
600  continue;
601  for (unsigned r2=0; r2<other.cols(); ++r2)
602  prod[r1*other.col+r2] += (m[r1*col+c] * other.m[c*other.col+r2]);
603  }
604  }
605  return matrix(row, other.col, std::move(prod));
606 }
607 
608 
610 matrix matrix::mul(const numeric & other) const
611 {
612  exvector prod(row * col);
613 
614  for (unsigned r=0; r<row; ++r)
615  for (unsigned c=0; c<col; ++c)
616  prod[r*col+c] = m[r*col+c] * other;
617 
618  return matrix(row, col, std::move(prod));
619 }
620 
621 
623 matrix matrix::mul_scalar(const ex & other) const
624 {
625  if (other.return_type() != return_types::commutative)
626  throw std::runtime_error("matrix::mul_scalar(): non-commutative scalar");
627 
628  exvector prod(row * col);
629 
630  for (unsigned r=0; r<row; ++r)
631  for (unsigned c=0; c<col; ++c)
632  prod[r*col+c] = m[r*col+c] * other;
633 
634  return matrix(row, col, std::move(prod));
635 }
636 
637 
639 matrix matrix::pow(const ex & expn) const
640 {
641  if (col!=row)
642  throw (std::logic_error("matrix::pow(): matrix not square"));
643 
644  if (is_exactly_a<numeric>(expn)) {
645  // Integer cases are computed by successive multiplication, using the
646  // obvious shortcut of storing temporaries, like A^4 == (A*A)*(A*A).
647  if (expn.info(info_flags::integer)) {
648  numeric b = ex_to<numeric>(expn);
649  matrix A(row,col);
650  if (expn.info(info_flags::negative)) {
651  b *= -1;
652  A = this->inverse();
653  } else {
654  A = *this;
655  }
656  matrix C(row,col);
657  for (unsigned r=0; r<row; ++r)
658  C(r,r) = _ex1;
659  if (b.is_zero())
660  return C;
661  // This loop computes the representation of b in base 2 from right
662  // to left and multiplies the factors whenever needed. Note
663  // that this is not entirely optimal but close to optimal and
664  // "better" algorithms are much harder to implement. (See Knuth,
665  // TAoCP2, section "Evaluation of Powers" for a good discussion.)
666  while (b!=*_num1_p) {
667  if (b.is_odd()) {
668  C = C.mul(A);
669  --b;
670  }
671  b /= *_num2_p; // still integer.
672  A = A.mul(A);
673  }
674  return A.mul(C);
675  }
676  }
677  throw (std::runtime_error("matrix::pow(): don't know how to handle exponent"));
678 }
679 
680 
686 const ex & matrix::operator() (unsigned ro, unsigned co) const
687 {
688  if (ro>=row || co>=col)
689  throw (std::range_error("matrix::operator(): index out of range"));
690 
691  return m[ro*col+co];
692 }
693 
694 
700 ex & matrix::operator() (unsigned ro, unsigned co)
701 {
702  if (ro>=row || co>=col)
703  throw (std::range_error("matrix::operator(): index out of range"));
704 
706  return m[ro*col+co];
707 }
708 
709 
713 {
714  exvector trans(this->cols()*this->rows());
715 
716  for (unsigned r=0; r<this->cols(); ++r)
717  for (unsigned c=0; c<this->rows(); ++c)
718  trans[r*this->rows()+c] = m[c*this->cols()+r];
719 
720  return matrix(this->cols(), this->rows(), std::move(trans));
721 }
722 
737 ex matrix::determinant(unsigned algo) const
738 {
739  if (row!=col)
740  throw (std::logic_error("matrix::determinant(): matrix not square"));
741  GINAC_ASSERT(row*col==m.capacity());
742 
743  // Gather some statistical information about this matrix:
744  bool numeric_flag = true;
745  bool normal_flag = false;
746  unsigned sparse_count = 0; // counts non-zero elements
747  for (auto r : m) {
748  if (!r.info(info_flags::numeric))
749  numeric_flag = false;
750  exmap srl; // symbol replacement list
751  ex rtest = r.to_rational(srl);
752  if (!rtest.is_zero())
753  ++sparse_count;
756  normal_flag = true;
757  }
758 
759  // Here is the heuristics in case this routine has to decide:
760  if (algo == determinant_algo::automatic) {
761  // Minor expansion is generally a good guess:
763  // Does anybody know when a matrix is really sparse?
764  // Maybe <~row/2.236 nonzero elements average in a row?
765  if (row>3 && 5*sparse_count<=row*col)
767  // Purely numeric matrix can be handled by Gauss elimination.
768  // This overrides any prior decisions.
769  if (numeric_flag)
771  }
772 
773  // Trap the trivial case here, since some algorithms don't like it
774  if (this->row==1) {
775  // for consistency with non-trivial determinants...
776  if (normal_flag)
777  return m[0].normal();
778  else
779  return m[0].expand();
780  }
781 
782  // Compute the determinant
783  switch(algo) {
785  ex det = 1;
786  matrix tmp(*this);
787  int sign = tmp.gauss_elimination(true);
788  for (unsigned d=0; d<row; ++d)
789  det *= tmp.m[d*col+d];
790  if (normal_flag)
791  return (sign*det).normal();
792  else
793  return (sign*det).normal().expand();
794  }
796  matrix tmp(*this);
797  int sign;
798  sign = tmp.fraction_free_elimination(true);
799  if (normal_flag)
800  return (sign*tmp.m[row*col-1]).normal();
801  else
802  return (sign*tmp.m[row*col-1]).expand();
803  }
805  matrix tmp(*this);
806  int sign;
807  sign = tmp.division_free_elimination(true);
808  if (sign==0)
809  return _ex0;
810  ex det = tmp.m[row*col-1];
811  // factor out accumulated bogus slag
812  for (unsigned d=0; d<row-2; ++d)
813  for (unsigned j=0; j<row-d-2; ++j)
814  det = (det/tmp.m[d*col+d]).normal();
815  return (sign*det);
816  }
818  default: {
819  // This is the minor expansion scheme. We always develop such
820  // that the smallest minors (i.e, the trivial 1x1 ones) are on the
821  // rightmost column. For this to be efficient, empirical tests
822  // have shown that the emptiest columns (i.e. the ones with most
823  // zeros) should be the ones on the right hand side -- although
824  // this might seem counter-intuitive (and in contradiction to some
825  // literature like the FORM manual). Please go ahead and test it
826  // if you don't believe me! Therefore we presort the columns of
827  // the matrix:
828  typedef std::pair<unsigned,unsigned> uintpair;
829  std::vector<uintpair> c_zeros; // number of zeros in column
830  for (unsigned c=0; c<col; ++c) {
831  unsigned acc = 0;
832  for (unsigned r=0; r<row; ++r)
833  if (m[r*col+c].is_zero())
834  ++acc;
835  c_zeros.push_back(uintpair(acc,c));
836  }
837  std::sort(c_zeros.begin(),c_zeros.end());
838  std::vector<unsigned> pre_sort;
839  for (auto & i : c_zeros)
840  pre_sort.push_back(i.second);
841  std::vector<unsigned> pre_sort_test(pre_sort); // permutation_sign() modifies the vector so we make a copy here
842  int sign = permutation_sign(pre_sort_test.begin(), pre_sort_test.end());
843  exvector result(row*col); // represents sorted matrix
844  unsigned c = 0;
845  for (auto & it : pre_sort) {
846  for (unsigned r=0; r<row; ++r)
847  result[r*col+c] = m[r*col+it];
848  ++c;
849  }
850 
851  if (normal_flag)
852  return (sign*matrix(row, col, std::move(result)).determinant_minor()).normal();
853  else
854  return sign*matrix(row, col, std::move(result)).determinant_minor();
855  }
856  }
857 }
858 
859 
867 {
868  if (row != col)
869  throw (std::logic_error("matrix::trace(): matrix not square"));
870 
871  ex tr;
872  for (unsigned r=0; r<col; ++r)
873  tr += m[r*col+r];
874 
877  return tr.normal();
878  else
879  return tr.expand();
880 }
881 
882 
894 ex matrix::charpoly(const ex & lambda) const
895 {
896  if (row != col)
897  throw (std::logic_error("matrix::charpoly(): matrix not square"));
898 
899  bool numeric_flag = true;
900  for (auto & r : m) {
901  if (!r.info(info_flags::numeric)) {
902  numeric_flag = false;
903  break;
904  }
905  }
906 
907  // The pure numeric case is traditionally rather common. Hence, it is
908  // trapped and we use Leverrier's algorithm which goes as row^3 for
909  // every coefficient. The expensive part is the matrix multiplication.
910  if (numeric_flag) {
911 
912  matrix B(*this);
913  ex c = B.trace();
914  ex poly = power(lambda, row) - c*power(lambda, row-1);
915  for (unsigned i=1; i<row; ++i) {
916  for (unsigned j=0; j<row; ++j)
917  B.m[j*col+j] -= c;
918  B = this->mul(B);
919  c = B.trace() / ex(i+1);
920  poly -= c*power(lambda, row-i-1);
921  }
922  if (row%2)
923  return -poly;
924  else
925  return poly;
926 
927  } else {
928 
929  matrix M(*this);
930  for (unsigned r=0; r<col; ++r)
931  M.m[r*col+r] -= lambda;
932 
933  return M.determinant().collect(lambda);
934  }
935 }
936 
937 
940 {
942 }
943 
950 matrix matrix::inverse(unsigned algo) const
951 {
952  if (row != col)
953  throw (std::logic_error("matrix::inverse(): matrix not square"));
954 
955  // This routine actually doesn't do anything fancy at all. We compute the
956  // inverse of the matrix A by solving the system A * A^{-1} == Id.
957 
958  // First populate the identity matrix supposed to become the right hand side.
959  matrix identity(row,col);
960  for (unsigned i=0; i<row; ++i)
961  identity(i,i) = _ex1;
962 
963  // Populate a dummy matrix of variables, just because of compatibility with
964  // matrix::solve() which wants this (for compatibility with under-determined
965  // systems of equations).
966  matrix vars(row,col);
967  for (unsigned r=0; r<row; ++r)
968  for (unsigned c=0; c<col; ++c)
969  vars(r,c) = symbol();
970 
971  matrix sol(row,col);
972  try {
973  sol = this->solve(vars, identity, algo);
974  } catch (const std::runtime_error & e) {
975  if (e.what()==std::string("matrix::solve(): inconsistent linear system"))
976  throw (std::runtime_error("matrix::inverse(): singular matrix"));
977  else
978  throw;
979  }
980  return sol;
981 }
982 
983 
996  const matrix & rhs,
997  unsigned algo) const
998 {
999  const unsigned m = this->rows();
1000  const unsigned n = this->cols();
1001  const unsigned p = rhs.cols();
1002 
1003  // syntax checks
1004  if ((rhs.rows() != m) || (vars.rows() != n) || (vars.cols() != p))
1005  throw (std::logic_error("matrix::solve(): incompatible matrices"));
1006  for (unsigned ro=0; ro<n; ++ro)
1007  for (unsigned co=0; co<p; ++co)
1008  if (!vars(ro,co).info(info_flags::symbol))
1009  throw (std::invalid_argument("matrix::solve(): 1st argument must be matrix of symbols"));
1010 
1011  // build the augmented matrix of *this with rhs attached to the right
1012  matrix aug(m,n+p);
1013  for (unsigned r=0; r<m; ++r) {
1014  for (unsigned c=0; c<n; ++c)
1015  aug.m[r*(n+p)+c] = this->m[r*n+c];
1016  for (unsigned c=0; c<p; ++c)
1017  aug.m[r*(n+p)+c+n] = rhs.m[r*p+c];
1018  }
1019 
1020  // Eliminate the augmented matrix:
1021  auto colid = aug.echelon_form(algo, n);
1022 
1023  // assemble the solution matrix:
1024  matrix sol(n,p);
1025  for (unsigned co=0; co<p; ++co) {
1026  unsigned last_assigned_sol = n+1;
1027  for (int r=m-1; r>=0; --r) {
1028  unsigned fnz = 1; // first non-zero in row
1029  while ((fnz<=n) && (aug.m[r*(n+p)+(fnz-1)].normal().is_zero()))
1030  ++fnz;
1031  if (fnz>n) {
1032  // row consists only of zeros, corresponding rhs must be 0, too
1033  if (!aug.m[r*(n+p)+n+co].normal().is_zero()) {
1034  throw (std::runtime_error("matrix::solve(): inconsistent linear system"));
1035  }
1036  } else {
1037  // assign solutions for vars between fnz+1 and
1038  // last_assigned_sol-1: free parameters
1039  for (unsigned c=fnz; c<last_assigned_sol-1; ++c)
1040  sol(colid[c],co) = vars.m[colid[c]*p+co];
1041  ex e = aug.m[r*(n+p)+n+co];
1042  for (unsigned c=fnz; c<n; ++c)
1043  e -= aug.m[r*(n+p)+c]*sol.m[colid[c]*p+co];
1044  sol(colid[fnz-1],co) = (e/(aug.m[r*(n+p)+fnz-1])).normal();
1045  last_assigned_sol = fnz;
1046  }
1047  }
1048  // assign solutions for vars between 1 and
1049  // last_assigned_sol-1: free parameters
1050  for (unsigned ro=0; ro<last_assigned_sol-1; ++ro)
1051  sol(colid[ro],co) = vars(colid[ro],co);
1052  }
1053 
1054  return sol;
1055 }
1056 
1058 unsigned matrix::rank() const
1059 {
1060  return rank(solve_algo::automatic);
1061 }
1062 
1065 unsigned matrix::rank(unsigned solve_algo) const
1066 {
1067  // Method:
1068  // Transform this matrix into upper echelon form and then count the
1069  // number of non-zero rows.
1070  GINAC_ASSERT(row*col==m.capacity());
1071 
1072  matrix to_eliminate = *this;
1073  to_eliminate.echelon_form(solve_algo, col);
1074 
1075  unsigned r = row*col; // index of last non-zero element
1076  while (r--) {
1077  if (!to_eliminate.m[r].is_zero())
1078  return 1+r/col;
1079  }
1080  return 0;
1081 }
1082 
1083 
1084 // protected
1085 
1097 {
1098  const unsigned n = this->cols();
1099 
1100  // This algorithm can best be understood by looking at a naive
1101  // implementation of Laplace-expansion, like this one:
1102  // ex det;
1103  // matrix minorM(this->rows()-1,this->cols()-1);
1104  // for (unsigned r1=0; r1<this->rows(); ++r1) {
1105  // // shortcut if element(r1,0) vanishes
1106  // if (m[r1*col].is_zero())
1107  // continue;
1108  // // assemble the minor matrix
1109  // for (unsigned r=0; r<minorM.rows(); ++r) {
1110  // for (unsigned c=0; c<minorM.cols(); ++c) {
1111  // if (r<r1)
1112  // minorM(r,c) = m[r*col+c+1];
1113  // else
1114  // minorM(r,c) = m[(r+1)*col+c+1];
1115  // }
1116  // }
1117  // // recurse down and care for sign:
1118  // if (r1%2)
1119  // det -= m[r1*col] * minorM.determinant_minor();
1120  // else
1121  // det += m[r1*col] * minorM.determinant_minor();
1122  // }
1123  // return det.expand();
1124  // What happens is that while proceeding down many of the minors are
1125  // computed more than once. In particular, there are binomial(n,k)
1126  // kxk minors and each one is computed factorial(n-k) times. Therefore
1127  // it is reasonable to store the results of the minors. We proceed from
1128  // right to left. At each column c we only need to retrieve the minors
1129  // calculated in step c-1. We therefore only have to store at most
1130  // 2*binomial(n,n/2) minors.
1131 
1132  // we store the minors in maps, keyed by the rows they arise from
1133  typedef std::vector<unsigned> keyseq;
1134  typedef std::map<keyseq, ex> Rmap;
1135 
1136  Rmap M, N; // minors used in current and next column, respectively
1137  // populate M with dummy unit, to be used as factor in rightmost column
1138  M[keyseq{}] = _ex1;
1139 
1140  // keys to identify minor of M and N (Mkey is a subsequence of Nkey)
1141  keyseq Mkey, Nkey;
1142  Mkey.reserve(n-1);
1143  Nkey.reserve(n);
1144 
1145  ex det;
1146  // proceed from right to left through matrix
1147  for (int c=n-1; c>=0; --c) {
1148  Nkey.clear();
1149  Mkey.clear();
1150  for (unsigned i=0; i<n-c; ++i)
1151  Nkey.push_back(i);
1152  unsigned fc = 0; // controls logic for minor key generator
1153  do {
1154  det = _ex0;
1155  for (unsigned r=0; r<n-c; ++r) {
1156  // maybe there is nothing to do?
1157  if (m[Nkey[r]*n+c].is_zero())
1158  continue;
1159  // Mkey is same as Nkey, but with element r removed
1160  Mkey.clear();
1161  Mkey.insert(Mkey.begin(), Nkey.begin(), Nkey.begin() + r);
1162  Mkey.insert(Mkey.end(), Nkey.begin() + r + 1, Nkey.end());
1163  // add product of matrix element and minor M to determinant
1164  if (r%2)
1165  det -= m[Nkey[r]*n+c]*M[Mkey];
1166  else
1167  det += m[Nkey[r]*n+c]*M[Mkey];
1168  }
1169  // prevent nested expressions to save time
1170  det = det.expand();
1171  // if the next computed minor is zero, don't store it in N:
1172  // (if key is not found, operator[] will just return a zero ex)
1173  if (!det.is_zero())
1174  N[Nkey] = det;
1175  // compute next minor key
1176  for (fc=n-c; fc>0; --fc) {
1177  ++Nkey[fc-1];
1178  if (Nkey[fc-1]<fc+c)
1179  break;
1180  }
1181  if (fc<n-c && fc>0)
1182  for (unsigned j=fc; j<n-c; ++j)
1183  Nkey[j] = Nkey[j-1]+1;
1184  } while(fc);
1185  // if N contains no minors, then they all vanished
1186  if (N.empty())
1187  return _ex0;
1188 
1189  // proceed to next column: switch roles of M and N, clear N
1190  M = std::move(N);
1191  }
1192 
1193  return det;
1194 }
1195 
1196 std::vector<unsigned>
1197 matrix::echelon_form(unsigned algo, int n)
1198 {
1199  // Here is the heuristics in case this routine has to decide:
1200  if (algo == solve_algo::automatic) {
1201  // Gather some statistical information about the augmented matrix:
1202  bool numeric_flag = true;
1203  for (const auto & r : m) {
1204  if (!r.info(info_flags::numeric)) {
1205  numeric_flag = false;
1206  break;
1207  }
1208  }
1209  unsigned density = 0;
1210  for (const auto & r : m) {
1211  density += !r.is_zero();
1212  }
1213  unsigned ncells = col*row;
1214  if (numeric_flag) {
1215  // For numerical matrices Gauss is good, but Markowitz becomes
1216  // better for large sparse matrices.
1217  if ((ncells > 200) && (density < ncells/2)) {
1218  algo = solve_algo::markowitz;
1219  } else {
1220  algo = solve_algo::gauss;
1221  }
1222  } else {
1223  // For symbolic matrices Markowitz is good, but Bareiss/Divfree
1224  // is better for small and dense matrices.
1225  if ((ncells < 120) && (density*5 > ncells*3)) {
1226  if (ncells <= 12) {
1227  algo = solve_algo::divfree;
1228  } else {
1229  algo = solve_algo::bareiss;
1230  }
1231  } else {
1232  algo = solve_algo::markowitz;
1233  }
1234  }
1235  }
1236  // Eliminate the augmented matrix:
1237  std::vector<unsigned> colid(col);
1238  for (unsigned c = 0; c < col; c++) {
1239  colid[c] = c;
1240  }
1241  switch(algo) {
1242  case solve_algo::gauss:
1244  break;
1245  case solve_algo::divfree:
1247  break;
1248  case solve_algo::bareiss:
1250  break;
1251  case solve_algo::markowitz:
1252  colid = markowitz_elimination(n);
1253  break;
1254  default:
1255  throw std::invalid_argument("matrix::echelon_form(): 'algo' is not one of the solve_algo enum");
1256  }
1257  return colid;
1258 }
1259 
1269 int matrix::gauss_elimination(const bool det)
1270 {
1272  const unsigned m = this->rows();
1273  const unsigned n = this->cols();
1274  GINAC_ASSERT(!det || n==m);
1275  int sign = 1;
1276 
1277  unsigned r0 = 0;
1278  for (unsigned c0=0; c0<n && r0<m-1; ++c0) {
1279  int indx = pivot(r0, c0, true);
1280  if (indx == -1) {
1281  sign = 0;
1282  if (det)
1283  return 0; // leaves *this in a messy state
1284  }
1285  if (indx>=0) {
1286  if (indx > 0)
1287  sign = -sign;
1288  for (unsigned r2=r0+1; r2<m; ++r2) {
1289  if (!this->m[r2*n+c0].is_zero()) {
1290  // yes, there is something to do in this row
1291  ex piv = this->m[r2*n+c0] / this->m[r0*n+c0];
1292  for (unsigned c=c0+1; c<n; ++c) {
1293  this->m[r2*n+c] -= piv * this->m[r0*n+c];
1294  if (!this->m[r2*n+c].info(info_flags::numeric))
1295  this->m[r2*n+c] = this->m[r2*n+c].normal();
1296  }
1297  }
1298  // fill up left hand side with zeros
1299  for (unsigned c=r0; c<=c0; ++c)
1300  this->m[r2*n+c] = _ex0;
1301  }
1302  if (det) {
1303  // save space by deleting no longer needed elements
1304  for (unsigned c=r0+1; c<n; ++c)
1305  this->m[r0*n+c] = _ex0;
1306  }
1307  ++r0;
1308  }
1309  }
1310  // clear remaining rows
1311  for (unsigned r=r0+1; r<m; ++r) {
1312  for (unsigned c=0; c<n; ++c)
1313  this->m[r*n+c] = _ex0;
1314  }
1315 
1316  return sign;
1317 }
1318 
1319 /* Perform Markowitz-ordered Gaussian elimination (with full
1320  * pivoting) on a matrix, constraining the choice of pivots to
1321  * the first n columns (this simplifies handling of augmented
1322  * matrices). Return the column id vector v, such that v[column]
1323  * is the original number of the column before shuffling (v[i]==i
1324  * for i >= n). */
1325 std::vector<unsigned>
1327 {
1328  GINAC_ASSERT(n <= col);
1329  std::vector<int> rowcnt(row, 0);
1330  std::vector<int> colcnt(col, 0);
1331  // Normalize everything before start. We'll keep all the
1332  // cells normalized throughout the algorithm to properly
1333  // handle unnormal zeros.
1334  for (unsigned r = 0; r < row; r++) {
1335  for (unsigned c = 0; c < col; c++) {
1336  if (!m[r*col + c].is_zero()) {
1337  m[r*col + c] = m[r*col + c].normal();
1338  rowcnt[r]++;
1339  colcnt[c]++;
1340  }
1341  }
1342  }
1343  std::vector<unsigned> colid(col);
1344  for (unsigned c = 0; c < col; c++) {
1345  colid[c] = c;
1346  }
1347  exvector ab(row);
1348  for (unsigned k = 0; (k < col) && (k < row - 1); k++) {
1349  // Find the pivot that minimizes (rowcnt[r]-1)*(colcnt[c]-1).
1350  unsigned pivot_r = row + 1;
1351  unsigned pivot_c = col + 1;
1352  int pivot_m = row*col;
1353  for (unsigned r = k; r < row; r++) {
1354  for (unsigned c = k; c < n; c++) {
1355  const ex &mrc = m[r*col + c];
1356  if (mrc.is_zero())
1357  continue;
1358  GINAC_ASSERT(rowcnt[r] > 0);
1359  GINAC_ASSERT(colcnt[c] > 0);
1360  int measure = (rowcnt[r] - 1)*(colcnt[c] - 1);
1361  if (measure < pivot_m) {
1362  pivot_m = measure;
1363  pivot_r = r;
1364  pivot_c = c;
1365  }
1366  }
1367  }
1368  if (pivot_m == row*col) {
1369  // The rest of the matrix is zero.
1370  break;
1371  }
1372  GINAC_ASSERT(k <= pivot_r && pivot_r < row);
1373  GINAC_ASSERT(k <= pivot_c && pivot_c < col);
1374  // Swap the pivot into (k, k).
1375  if (pivot_c != k) {
1376  for (unsigned r = 0; r < row; r++) {
1377  m[r*col + pivot_c].swap(m[r*col + k]);
1378  }
1379  std::swap(colid[pivot_c], colid[k]);
1380  std::swap(colcnt[pivot_c], colcnt[k]);
1381  }
1382  if (pivot_r != k) {
1383  for (unsigned c = k; c < col; c++) {
1384  m[pivot_r*col + c].swap(m[k*col + c]);
1385  }
1386  std::swap(rowcnt[pivot_r], rowcnt[k]);
1387  }
1388  // No normalization before is_zero() here, because
1389  // we maintain the matrix normalized throughout the
1390  // algorithm.
1391  ex a = m[k*col + k];
1392  GINAC_ASSERT(!a.is_zero());
1393  // Subtract the pivot row KJI-style (so: loop by pivot, then
1394  // column, then row) to maximally exploit pivot row zeros (at
1395  // the expense of the pivot column zeros). The speedup compared
1396  // to the usual KIJ order is not really significant though...
1397  for (unsigned r = k + 1; r < row; r++) {
1398  const ex &b = m[r*col + k];
1399  if (!b.is_zero()) {
1400  ab[r] = b/a;
1401  rowcnt[r]--;
1402  }
1403  }
1404  colcnt[k] = rowcnt[k] = 0;
1405  for (unsigned c = k + 1; c < col; c++) {
1406  const ex &mr0c = m[k*col + c];
1407  if (mr0c.is_zero())
1408  continue;
1409  colcnt[c]--;
1410  for (unsigned r = k + 1; r < row; r++) {
1411  if (ab[r].is_zero())
1412  continue;
1413  bool waszero = m[r*col + c].is_zero();
1414  m[r*col + c] = (m[r*col + c] - ab[r]*mr0c).normal();
1415  bool iszero = m[r*col + c].is_zero();
1416  if (waszero && !iszero) {
1417  rowcnt[r]++;
1418  colcnt[c]++;
1419  }
1420  if (!waszero && iszero) {
1421  rowcnt[r]--;
1422  colcnt[c]--;
1423  }
1424  }
1425  }
1426  for (unsigned r = k + 1; r < row; r++) {
1427  ab[r] = m[r*col + k] = _ex0;
1428  }
1429  }
1430  return colid;
1431 }
1432 
1442 {
1444  const unsigned m = this->rows();
1445  const unsigned n = this->cols();
1446  GINAC_ASSERT(!det || n==m);
1447  int sign = 1;
1448 
1449  unsigned r0 = 0;
1450  for (unsigned c0=0; c0<n && r0<m-1; ++c0) {
1451  int indx = pivot(r0, c0, true);
1452  if (indx==-1) {
1453  sign = 0;
1454  if (det)
1455  return 0; // leaves *this in a messy state
1456  }
1457  if (indx>=0) {
1458  if (indx>0)
1459  sign = -sign;
1460  for (unsigned r2=r0+1; r2<m; ++r2) {
1461  for (unsigned c=c0+1; c<n; ++c)
1462  this->m[r2*n+c] = (this->m[r0*n+c0]*this->m[r2*n+c] - this->m[r2*n+c0]*this->m[r0*n+c]).normal();
1463  // fill up left hand side with zeros
1464  for (unsigned c=r0; c<=c0; ++c)
1465  this->m[r2*n+c] = _ex0;
1466  }
1467  if (det) {
1468  // save space by deleting no longer needed elements
1469  for (unsigned c=r0+1; c<n; ++c)
1470  this->m[r0*n+c] = _ex0;
1471  }
1472  ++r0;
1473  }
1474  }
1475  // clear remaining rows
1476  for (unsigned r=r0+1; r<m; ++r) {
1477  for (unsigned c=0; c<n; ++c)
1478  this->m[r*n+c] = _ex0;
1479  }
1480 
1481  return sign;
1482 }
1483 
1484 
1496 {
1497  // Method:
1498  // (single-step fraction free elimination scheme, already known to Jordan)
1499  //
1500  // Usual division-free elimination sets m[0](r,c) = m(r,c) and then sets
1501  // m[k+1](r,c) = m[k](k,k) * m[k](r,c) - m[k](r,k) * m[k](k,c).
1502  //
1503  // Bareiss (fraction-free) elimination in addition divides that element
1504  // by m[k-1](k-1,k-1) for k>1, where it can be shown by means of the
1505  // Sylvester identity that this really divides m[k+1](r,c).
1506  //
1507  // We also allow rational functions where the original prove still holds.
1508  // However, we must care for numerator and denominator separately and
1509  // "manually" work in the integral domains because of subtle cancellations
1510  // (see below). This blows up the bookkeeping a bit and the formula has
1511  // to be modified to expand like this (N{x} stands for numerator of x,
1512  // D{x} for denominator of x):
1513  // N{m[k+1](r,c)} = N{m[k](k,k)}*N{m[k](r,c)}*D{m[k](r,k)}*D{m[k](k,c)}
1514  // -N{m[k](r,k)}*N{m[k](k,c)}*D{m[k](k,k)}*D{m[k](r,c)}
1515  // D{m[k+1](r,c)} = D{m[k](k,k)}*D{m[k](r,c)}*D{m[k](r,k)}*D{m[k](k,c)}
1516  // where for k>1 we now divide N{m[k+1](r,c)} by
1517  // N{m[k-1](k-1,k-1)}
1518  // and D{m[k+1](r,c)} by
1519  // D{m[k-1](k-1,k-1)}.
1520 
1522  const unsigned m = this->rows();
1523  const unsigned n = this->cols();
1524  GINAC_ASSERT(!det || n==m);
1525  int sign = 1;
1526  if (m==1)
1527  return 1;
1528  ex divisor_n = 1;
1529  ex divisor_d = 1;
1530  ex dividend_n;
1531  ex dividend_d;
1532 
1533  // We populate temporary matrices to subsequently operate on. There is
1534  // one holding numerators and another holding denominators of entries.
1535  // This is a must since the evaluator (or even earlier mul's constructor)
1536  // might cancel some trivial element which causes divide() to fail. The
1537  // elements are normalized first (yes, even though this algorithm doesn't
1538  // need GCDs) since the elements of *this might be unnormalized, which
1539  // makes things more complicated than they need to be.
1540  matrix tmp_n(*this);
1541  matrix tmp_d(m,n); // for denominators, if needed
1542  exmap srl; // symbol replacement list
1543  auto tmp_n_it = tmp_n.m.begin(), tmp_d_it = tmp_d.m.begin();
1544  for (auto & it : this->m) {
1545  ex nd = it.normal().to_rational(srl).numer_denom();
1546  *tmp_n_it++ = nd.op(0);
1547  *tmp_d_it++ = nd.op(1);
1548  }
1549 
1550  unsigned r0 = 0;
1551  for (unsigned c0=0; c0<n && r0<m-1; ++c0) {
1552  // When trying to find a pivot, we should try a bit harder than expand().
1553  // Searching the first non-zero element in-place here instead of calling
1554  // pivot() allows us to do no more substitutions and back-substitutions
1555  // than are actually necessary.
1556  unsigned indx = r0;
1557  while ((indx<m) &&
1558  (tmp_n[indx*n+c0].subs(srl, subs_options::no_pattern).expand().is_zero()))
1559  ++indx;
1560  if (indx==m) {
1561  // all elements in column c0 below row r0 vanish
1562  sign = 0;
1563  if (det)
1564  return 0;
1565  } else {
1566  if (indx>r0) {
1567  // Matrix needs pivoting, swap rows r0 and indx of tmp_n and tmp_d.
1568  sign = -sign;
1569  for (unsigned c=c0; c<n; ++c) {
1570  tmp_n.m[n*indx+c].swap(tmp_n.m[n*r0+c]);
1571  tmp_d.m[n*indx+c].swap(tmp_d.m[n*r0+c]);
1572  }
1573  }
1574  for (unsigned r2=r0+1; r2<m; ++r2) {
1575  for (unsigned c=c0+1; c<n; ++c) {
1576  dividend_n = (tmp_n.m[r0*n+c0]*tmp_n.m[r2*n+c]*
1577  tmp_d.m[r2*n+c0]*tmp_d.m[r0*n+c]
1578  -tmp_n.m[r2*n+c0]*tmp_n.m[r0*n+c]*
1579  tmp_d.m[r0*n+c0]*tmp_d.m[r2*n+c]).expand();
1580  dividend_d = (tmp_d.m[r2*n+c0]*tmp_d.m[r0*n+c]*
1581  tmp_d.m[r0*n+c0]*tmp_d.m[r2*n+c]).expand();
1582  bool check = divide(dividend_n, divisor_n,
1583  tmp_n.m[r2*n+c], true);
1584  check &= divide(dividend_d, divisor_d,
1585  tmp_d.m[r2*n+c], true);
1586  GINAC_ASSERT(check);
1587  }
1588  // fill up left hand side with zeros
1589  for (unsigned c=r0; c<=c0; ++c)
1590  tmp_n.m[r2*n+c] = _ex0;
1591  }
1592  if (c0<n && r0<m-1) {
1593  // compute next iteration's divisor
1594  divisor_n = tmp_n.m[r0*n+c0].expand();
1595  divisor_d = tmp_d.m[r0*n+c0].expand();
1596  if (det) {
1597  // save space by deleting no longer needed elements
1598  for (unsigned c=0; c<n; ++c) {
1599  tmp_n.m[r0*n+c] = _ex0;
1600  tmp_d.m[r0*n+c] = _ex1;
1601  }
1602  }
1603  }
1604  ++r0;
1605  }
1606  }
1607  // clear remaining rows
1608  for (unsigned r=r0+1; r<m; ++r) {
1609  for (unsigned c=0; c<n; ++c)
1610  tmp_n.m[r*n+c] = _ex0;
1611  }
1612 
1613  // repopulate *this matrix:
1614  tmp_n_it = tmp_n.m.begin();
1615  tmp_d_it = tmp_d.m.begin();
1616  for (auto & it : this->m)
1617  it = ((*tmp_n_it++)/(*tmp_d_it++)).subs(srl, subs_options::no_pattern);
1618 
1619  return sign;
1620 }
1621 
1622 
1636 int matrix::pivot(unsigned ro, unsigned co, bool symbolic)
1637 {
1638  unsigned k = ro;
1639  if (symbolic) {
1640  // search first non-zero element in column co beginning at row ro
1641  while ((k<row) && (this->m[k*col+co].expand().is_zero()))
1642  ++k;
1643  } else {
1644  // search largest element in column co beginning at row ro
1645  GINAC_ASSERT(is_exactly_a<numeric>(this->m[k*col+co]));
1646  unsigned kmax = k+1;
1647  numeric mmax = abs(ex_to<numeric>(m[kmax*col+co]));
1648  while (kmax<row) {
1649  GINAC_ASSERT(is_exactly_a<numeric>(this->m[kmax*col+co]));
1650  numeric tmp = ex_to<numeric>(this->m[kmax*col+co]);
1651  if (abs(tmp) > mmax) {
1652  mmax = tmp;
1653  k = kmax;
1654  }
1655  ++kmax;
1656  }
1657  if (!mmax.is_zero())
1658  k = kmax;
1659  }
1660  if (k==row)
1661  // all elements in column co below row ro vanish
1662  return -1;
1663  if (k==ro)
1664  // matrix needs no pivoting
1665  return 0;
1666  // matrix needs pivoting, so swap rows k and ro
1668  for (unsigned c=0; c<col; ++c)
1669  this->m[k*col+c].swap(this->m[ro*col+c]);
1670 
1671  return k;
1672 }
1673 
1677 {
1678  for (auto & i : m)
1679  if (!i.is_zero())
1680  return false;
1681  return true;
1682 }
1683 
1685 {
1686  // Find number of rows and columns
1687  size_t rows = l.nops(), cols = 0;
1688  for (auto & itr : l) {
1689  if (!is_a<lst>(itr))
1690  throw (std::invalid_argument("lst_to_matrix: argument must be a list of lists"));
1691  if (itr.nops() > cols)
1692  cols = itr.nops();
1693  }
1694 
1695  // Allocate and fill matrix
1696  matrix & M = dynallocate<matrix>(rows, cols);
1697 
1698  unsigned i = 0;
1699  for (auto & itr : l) {
1700  unsigned j = 0;
1701  for (auto & itc : ex_to<lst>(itr)) {
1702  M(i, j) = itc;
1703  ++j;
1704  }
1705  ++i;
1706  }
1707 
1708  return M;
1709 }
1710 
1711 ex diag_matrix(const lst & l)
1712 {
1713  size_t dim = l.nops();
1714 
1715  // Allocate and fill matrix
1716  matrix & M = dynallocate<matrix>(dim, dim);
1717 
1718  unsigned i = 0;
1719  for (auto & it : l) {
1720  M(i, i) = it;
1721  ++i;
1722  }
1723 
1724  return M;
1725 }
1726 
1727 ex diag_matrix(std::initializer_list<ex> l)
1728 {
1729  size_t dim = l.size();
1730 
1731  // Allocate and fill matrix
1732  matrix & M = dynallocate<matrix>(dim, dim);
1733 
1734  unsigned i = 0;
1735  for (auto & it : l) {
1736  M(i, i) = it;
1737  ++i;
1738  }
1739 
1740  return M;
1741 }
1742 
1743 ex unit_matrix(unsigned r, unsigned c)
1744 {
1745  matrix & Id = dynallocate<matrix>(r, c);
1747  for (unsigned i=0; i<r && i<c; i++)
1748  Id(i,i) = _ex1;
1749 
1750  return Id;
1751 }
1752 
1753 ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name, const std::string & tex_base_name)
1754 {
1755  matrix & M = dynallocate<matrix>(r, c);
1757 
1758  bool long_format = (r > 10 || c > 10);
1759  bool single_row = (r == 1 || c == 1);
1760 
1761  for (unsigned i=0; i<r; i++) {
1762  for (unsigned j=0; j<c; j++) {
1763  std::ostringstream s1, s2;
1764  s1 << base_name;
1765  s2 << tex_base_name << "_{";
1766  if (single_row) {
1767  if (c == 1) {
1768  s1 << i;
1769  s2 << i << '}';
1770  } else {
1771  s1 << j;
1772  s2 << j << '}';
1773  }
1774  } else {
1775  if (long_format) {
1776  s1 << '_' << i << '_' << j;
1777  s2 << i << ';' << j << "}";
1778  } else {
1779  s1 << i << j;
1780  s2 << i << j << '}';
1781  }
1782  }
1783  M(i, j) = symbol(s1.str(), s2.str());
1784  }
1785  }
1786 
1787  return M;
1788 }
1789 
1790 ex reduced_matrix(const matrix& m, unsigned r, unsigned c)
1791 {
1792  if (r+1>m.rows() || c+1>m.cols() || m.cols()<2 || m.rows()<2)
1793  throw std::runtime_error("minor_matrix(): index out of bounds");
1794 
1795  const unsigned rows = m.rows()-1;
1796  const unsigned cols = m.cols()-1;
1797  matrix & M = dynallocate<matrix>(rows, cols);
1799 
1800  unsigned ro = 0;
1801  unsigned ro2 = 0;
1802  while (ro2<rows) {
1803  if (ro==r)
1804  ++ro;
1805  unsigned co = 0;
1806  unsigned co2 = 0;
1807  while (co2<cols) {
1808  if (co==c)
1809  ++co;
1810  M(ro2,co2) = m(ro, co);
1811  ++co;
1812  ++co2;
1813  }
1814  ++ro;
1815  ++ro2;
1816  }
1817 
1818  return M;
1819 }
1820 
1821 ex sub_matrix(const matrix&m, unsigned r, unsigned nr, unsigned c, unsigned nc)
1822 {
1823  if (r+nr>m.rows() || c+nc>m.cols())
1824  throw std::runtime_error("sub_matrix(): index out of bounds");
1825 
1826  matrix & M = dynallocate<matrix>(nr, nc);
1828 
1829  for (unsigned ro=0; ro<nr; ++ro) {
1830  for (unsigned co=0; co<nc; ++co) {
1831  M(ro,co) = m(ro+r,co+c);
1832  }
1833  }
1834 
1835  return M;
1836 }
1837 
1838 } // namespace GiNaC
Interface to GiNaC's sums of expressions.
Archiving of GiNaC expressions.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
This class stores all properties needed to record/retrieve the state of one object of class basic (or...
Definition: archive.h:49
This class is the ABC (abstract base class) of GiNaC's class hierarchy.
Definition: basic.h:105
virtual size_t nops() const
Number of operands/members.
Definition: basic.cpp:229
virtual bool info(unsigned inf) const
Information about the object.
Definition: basic.cpp:222
void ensure_if_modifiable() const
Ensure the object may be modified without hurting others, throws if this is not the case.
Definition: basic.cpp:894
friend class ex
Definition: basic.h:108
virtual ex op(size_t i) const
Return operand/member at position i.
Definition: basic.cpp:238
const basic & hold() const
Stop further evaluation.
Definition: basic.cpp:887
bool is_equal(const basic &other) const
Test for syntactic equality.
Definition: basic.cpp:863
const basic & setflag(unsigned f) const
Set some status_flags.
Definition: basic.h:288
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
Definition: basic.cpp:719
void do_print_tree(const print_tree &c, unsigned level) const
Tree output to stream.
Definition: basic.cpp:175
virtual ex expand(unsigned options=0) const
Expand expression, i.e.
Definition: basic.cpp:796
virtual ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const
Default implementation of ex::normal().
Definition: normal.cpp:2191
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
size_t nops() const override
Number of operands/members.
Definition: container.h:118
@ automatic
Let the system choose.
Definition: flags.h:93
@ divfree
Division-free elimination.
Definition: flags.h:114
@ laplace
Laplace elimination.
Definition: flags.h:120
@ gauss
Gauss elimination.
Definition: flags.h:102
@ bareiss
Bareiss fraction-free elimination.
Definition: flags.h:136
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
ex to_rational(exmap &repl) const
Rationalization of non-rational functions.
Definition: normal.cpp:2600
ex expand(unsigned options=0) const
Definition: ex.cpp:73
ex numer_denom() const
Get numerator and denominator of an expression.
Definition: normal.cpp:2567
bool is_equal(const ex &other) const
Definition: ex.h:345
ex normal() const
Normalization of rational functions.
Definition: normal.cpp:2492
unsigned return_type() const
Definition: ex.h:230
size_t nops() const
Definition: ex.h:135
bool info(unsigned inf) const
Definition: ex.h:132
bool is_zero() const
Definition: ex.h:213
ex collect(const ex &s, bool distributed=false) const
Definition: ex.h:181
ex op(size_t i) const
Definition: ex.h:136
This class holds one index of an indexed object.
Definition: idx.h:36
ex get_dim() const
Get dimension of index space.
Definition: idx.h:81
ex get_value() const
Get value of index.
Definition: idx.h:72
This class holds an indexed expression.
Definition: indexed.h:40
@ crational_polynomial
Definition: flags.h:259
Symbolic matrices.
Definition: matrix.h:38
matrix inverse() const
Inverse of this matrix, with automatic algorithm selection.
Definition: matrix.cpp:939
int gauss_elimination(const bool det=false)
Perform the steps of an ordinary Gaussian elimination to bring the m x n matrix into an upper echelon...
Definition: matrix.cpp:1269
ex scalar_mul_indexed(const ex &self, const numeric &other) const override
Product of an indexed matrix with a number.
Definition: matrix.cpp:433
ex eval_indexed(const basic &i) const override
Automatic symbolic evaluation of an indexed matrix.
Definition: matrix.cpp:320
unsigned cols() const
Get number of columns.
Definition: matrix.h:77
ex charpoly(const ex &lambda) const
Characteristic Polynomial.
Definition: matrix.cpp:894
void do_print_latex(const print_latex &c, unsigned level) const
Definition: matrix.cpp:191
exvector m
representation (cols indexed first)
Definition: matrix.h:117
ex determinant(unsigned algo=determinant_algo::automatic) const
Determinant of square matrix.
Definition: matrix.cpp:737
const ex & operator()(unsigned ro, unsigned co) const
operator() to access elements for reading.
Definition: matrix.cpp:686
void archive(archive_node &n) const override
Save (a.k.a.
Definition: matrix.cpp:152
ex add_indexed(const ex &self, const ex &other) const override
Sum of two indexed matrices.
Definition: matrix.cpp:397
bool is_zero_matrix() const
Function to check that all elements of the matrix are zero.
Definition: matrix.cpp:1676
ex trace() const
Trace of a matrix.
Definition: matrix.cpp:866
bool match_same_type(const basic &other) const override
Returns true if the attributes of two objects are similar enough for a match.
Definition: matrix.cpp:309
matrix add(const matrix &other) const
Sum of matrices.
Definition: matrix.cpp:555
matrix pow(const ex &expn) const
Power of a matrix.
Definition: matrix.cpp:639
ex subs(const exmap &m, unsigned options=0) const override
Substitute a set of objects by arbitrary expressions.
Definition: matrix.cpp:228
matrix(unsigned r, unsigned c)
Very common ctor.
Definition: matrix.cpp:71
std::vector< unsigned > markowitz_elimination(unsigned n)
Definition: matrix.cpp:1326
matrix solve(const matrix &vars, const matrix &rhs, unsigned algo=solve_algo::automatic) const
Solve a linear system consisting of a m x n matrix and a m x p right hand side by applying an elimina...
Definition: matrix.cpp:995
void do_print_python_repr(const print_python_repr &c, unsigned level) const
Definition: matrix.cpp:198
ex imag_part() const override
Definition: matrix.cpp:273
void read_archive(const archive_node &n, lst &syms) override
Read (a.k.a.
Definition: matrix.cpp:134
matrix mul_scalar(const ex &other) const
Product of matrix and scalar expression.
Definition: matrix.cpp:623
void do_print(const print_context &c, unsigned level) const
Definition: matrix.cpp:184
void print_elements(const print_context &c, const char *row_start, const char *row_end, const char *row_sep, const char *col_sep) const
Definition: matrix.cpp:168
size_t nops() const override
nops is defined to be rows x columns.
Definition: matrix.cpp:206
ex real_part() const override
Definition: matrix.cpp:264
unsigned rank() const
Compute the rank of this matrix.
Definition: matrix.cpp:1058
ex determinant_minor() const
Recursive determinant for small matrices having at least one symbolic entry.
Definition: matrix.cpp:1096
matrix mul(const matrix &other) const
Product of matrices.
Definition: matrix.cpp:589
bool contract_with(exvector::iterator self, exvector::iterator other, exvector &v) const override
Contraction of an indexed matrix with something else.
Definition: matrix.cpp:448
unsigned rows() const
Get number of rows.
Definition: matrix.h:75
ex op(size_t i) const override
returns matrix entry at position (i/col, icol).
Definition: matrix.cpp:212
int division_free_elimination(const bool det=false)
Perform the steps of division free elimination to bring the m x n matrix into an upper echelon form.
Definition: matrix.cpp:1441
unsigned col
number of columns
Definition: matrix.h:116
ex & let_op(size_t i) override
returns writable matrix entry at position (i/col, icol).
Definition: matrix.cpp:220
matrix transpose() const
Transposed of an m x n matrix, producing a new n x m matrix object that represents the transposed.
Definition: matrix.cpp:712
matrix sub(const matrix &other) const
Difference of matrices.
Definition: matrix.cpp:572
std::vector< unsigned > echelon_form(unsigned algo, int n)
Definition: matrix.cpp:1197
ex conjugate() const override
Complex conjugate every matrix entry.
Definition: matrix.cpp:239
int pivot(unsigned ro, unsigned co, bool symbolic=true)
Partial pivoting method for matrix elimination schemes.
Definition: matrix.cpp:1636
int fraction_free_elimination(const bool det=false)
Perform the steps of Bareiss' one-step fraction free elimination to bring the matrix into an upper ec...
Definition: matrix.cpp:1495
unsigned row
number of rows
Definition: matrix.h:115
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Definition: numeric.h:82
bool is_odd() const
True if object is an exact odd integer.
Definition: numeric.cpp:1182
bool is_zero() const
True if object is zero.
Definition: numeric.cpp:1129
This class holds a two-component object, a basis and and exponent representing exponentiation.
Definition: power.h:39
Base class for print_contexts.
Definition: print.h:103
Context for latex-parsable output.
Definition: print.h:123
Context for python-parsable output.
Definition: print.h:139
Switch to control algorithm for linear system solving.
Definition: flags.h:141
@ bareiss
Bareiss fraction-free elimination.
Definition: flags.h:185
@ markowitz
Markowitz-ordered Gaussian elimination.
Definition: flags.h:193
@ automatic
Let the system choose.
Definition: flags.h:146
@ divfree
Division-free elimination.
Definition: flags.h:168
@ gauss
Gauss elimination.
Definition: flags.h:156
@ evaluated
.eval() has already done its job
Definition: flags.h:203
@ not_shareable
don't share instances of this object between different expressions unless explicitly asked to (used b...
Definition: flags.h:206
@ no_pattern
disable pattern matching
Definition: flags.h:51
Basic CAS symbol.
Definition: symbol.h:39
upoly poly
Definition: factor.cpp:1474
unsigned options
Definition: factor.cpp:2480
vector< int > k
Definition: factor.cpp:1466
ex x
Definition: factor.cpp:1641
size_t n
Definition: factor.cpp:1463
size_t c
Definition: factor.cpp:770
size_t r
Definition: factor.cpp:770
mvec m
Definition: factor.cpp:771
Interface to GiNaC's indices.
Interface to GiNaC's indexed expressions.
Definition of GiNaC's lst.
Interface to symbolic matrices.
Definition: add.cpp:38
bool is_zero(const ex &thisex)
Definition: ex.h:820
ex symbolic_matrix(unsigned r, unsigned c, const std::string &base_name, const std::string &tex_base_name)
Create an r times c matrix of newly generated symbols consisting of the given base name plus the nume...
Definition: matrix.cpp:1753
ex sub_matrix(const matrix &m, unsigned r, unsigned nr, unsigned c, unsigned nc)
Return the nr times nc submatrix starting at position r, c of matrix m.
Definition: matrix.cpp:1821
std::map< ex, ex, ex_is_less > exmap
Definition: basic.h:50
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
const ex _ex1
Definition: utils.cpp:385
bool are_ex_trivially_equal(const ex &e1, const ex &e2)
Compare two objects of class quickly without doing a deep tree traversal.
Definition: ex.h:684
ex rhs(const ex &thisex)
Definition: ex.h:817
ex diag_matrix(const lst &l)
Convert list of diagonal elements to matrix.
Definition: matrix.cpp:1711
const numeric * _num2_p
Definition: utils.cpp:388
unsigned rows(const matrix &m)
Definition: matrix.h:132
bool is_dummy_pair(const idx &i1, const idx &i2)
Check whether two indices form a dummy pair.
Definition: idx.cpp:502
print_func< print_context >(&varidx::do_print). print_func< print_latex >(&varidx
Definition: idx.cpp:45
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(add, expairseq, print_func< print_context >(&add::do_print). print_func< print_latex >(&add::do_print_latex). print_func< print_csrc >(&add::do_print_csrc). print_func< print_tree >(&add::do_print_tree). print_func< print_python_repr >(&add::do_print_python_repr)) add
Definition: add.cpp:40
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT_T(lst, basic, print_func< print_context >(&lst::do_print). print_func< print_tree >(&lst::do_print_tree)) template<> bool lst GINAC_BIND_UNARCHIVER(lst)
Specialization of container::info() for lst.
Definition: lst.cpp:42
unsigned cols(const matrix &m)
Definition: matrix.h:135
const numeric * _num1_p
Definition: utils.cpp:384
void swap(ex &e1, ex &e2)
Definition: ex.h:823
int permutation_sign(It first, It last)
Definition: utils.h:77
ex lst_to_matrix(const lst &l)
Convert list of lists to matrix.
Definition: matrix.cpp:1684
bool divide(const ex &a, const ex &b, ex &q, bool check_args)
Exact polynomial division of a(X) by b(X) in Q[X].
Definition: normal.cpp:595
const ex _ex0
Definition: utils.cpp:369
std::vector< ex > exvector
Definition: basic.h:46
ex unit_matrix(unsigned r, unsigned c)
Create an r times c unit matrix.
Definition: matrix.cpp:1743
ex reduced_matrix(const matrix &m, unsigned r, unsigned c)
Return the reduced matrix that is formed by deleting the rth row and cth column of matrix m.
Definition: matrix.cpp:1790
Definition: ex.h:972
void swap(GiNaC::ex &a, GiNaC::ex &b)
Specialization of std::swap() for ex objects.
Definition: ex.h:976
This file defines several functions that work on univariate and multivariate polynomials and rational...
Makes the interface to the underlying bignum package available.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
Interface to GiNaC's symbolic objects.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...

This page is part of the GiNaC developer's reference. It was generated automatically by doxygen. For an introduction, see the tutorial.