GiNaC  1.8.0
matrix.cpp
Go to the documentation of this file.
1 
5 /*
6  * GiNaC Copyright (C) 1999-2020 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
GiNaC::print_python_repr
Context for python-parsable output.
Definition: print.h:139
GiNaC::ex::expand
ex expand(unsigned options=0) const
Definition: ex.cpp:73
GiNaC::rhs
ex rhs(const ex &thisex)
Definition: ex.h:817
GiNaC::status_flags::not_shareable
@ not_shareable
don't share instances of this object between different expressions unless explicitly asked to (used b...
Definition: flags.h:206
GiNaC::GINAC_IMPLEMENT_REGISTERED_CLASS_OPT
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::info_flags::integer
@ integer
Definition: flags.h:223
GiNaC::matrix::col
unsigned col
number of columns
Definition: matrix.h:116
r
size_t r
Definition: factor.cpp:770
GiNaC::info_flags::negative
@ negative
Definition: flags.h:227
GiNaC::ex::to_rational
ex to_rational(exmap &repl) const
Rationalization of non-rational functions.
Definition: normal.cpp:2597
GiNaC::ex::normal
ex normal() const
Normalization of rational functions.
Definition: normal.cpp:2489
GiNaC::determinant_algo::automatic
@ automatic
Let the system choose.
Definition: flags.h:93
GiNaC::matrix::division_free_elimination
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
GiNaC::_ex0
const ex _ex0
Definition: utils.cpp:177
GiNaC::matrix::nops
size_t nops() const override
nops is defined to be rows x columns.
Definition: matrix.cpp:206
GiNaC::reduced_matrix
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
GiNaC::swap
void swap(ex &e1, ex &e2)
Definition: ex.h:823
numeric.h
Makes the interface to the underlying bignum package available.
GiNaC::matrix::determinant
ex determinant(unsigned algo=determinant_algo::automatic) const
Determinant of square matrix.
Definition: matrix.cpp:737
GiNaC::print_context
Base class for print_contexts.
Definition: print.h:103
GiNaC::determinant_algo::bareiss
@ bareiss
Bareiss fraction-free elimination.
Definition: flags.h:136
GiNaC::matrix::mul
matrix mul(const matrix &other) const
Product of matrices.
Definition: matrix.cpp:589
GiNaC::matrix::do_print
void do_print(const print_context &c, unsigned level) const
Definition: matrix.cpp:184
add.h
Interface to GiNaC's sums of expressions.
k
vector< int > k
Definition: factor.cpp:1466
GiNaC::matrix::read_archive
void read_archive(const archive_node &n, lst &syms) override
Read (a.k.a.
Definition: matrix.cpp:134
GiNaC::basic::expand
virtual ex expand(unsigned options=0) const
Expand expression, i.e.
Definition: basic.cpp:796
power.h
Interface to GiNaC's symbolic exponentiation (basis^exponent).
GiNaC::ex::collect
ex collect(const ex &s, bool distributed=false) const
Definition: ex.h:181
GiNaC::exvector
std::vector< ex > exvector
Definition: basic.h:46
GiNaC::ex::nops
size_t nops() const
Definition: ex.h:135
GiNaC::rows
unsigned rows(const matrix &m)
Definition: matrix.h:132
GiNaC::power
This class holds a two-component object, a basis and and exponent representing exponentiation.
Definition: power.h:39
GiNaC::status_flags::evaluated
@ evaluated
.eval() has already done its job
Definition: flags.h:203
GiNaC::matrix::do_print_python_repr
void do_print_python_repr(const print_python_repr &c, unsigned level) const
Definition: matrix.cpp:198
poly
upoly poly
Definition: factor.cpp:1474
GiNaC::archive_node
This class stores all properties needed to record/retrieve the state of one object of class basic (or...
Definition: archive.h:49
GiNaC::numeric::is_odd
bool is_odd() const
True if object is an exact odd integer.
Definition: numeric.cpp:1182
GiNaC::matrix::matrix
matrix(unsigned r, unsigned c)
Very common ctor.
Definition: matrix.cpp:71
GiNaC::basic::normal
virtual ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const
Default implementation of ex::normal().
Definition: normal.cpp:2192
idx.h
Interface to GiNaC's indices.
GiNaC::determinant_algo::divfree
@ divfree
Division-free elimination.
Definition: flags.h:114
GiNaC::_ex1
const ex _ex1
Definition: utils.cpp:193
GiNaC::idx::get_value
ex get_value() const
Get value of index.
Definition: idx.h:72
GiNaC::basic::nops
virtual size_t nops() const
Number of operands/members.
Definition: basic.cpp:229
GiNaC::diag_matrix
ex diag_matrix(const lst &l)
Convert list of diagonal elements to matrix.
Definition: matrix.cpp:1711
GiNaC::matrix::pivot
int pivot(unsigned ro, unsigned co, bool symbolic=true)
Partial pivoting method for matrix elimination schemes.
Definition: matrix.cpp:1636
options
unsigned options
Definition: factor.cpp:2480
GiNaC::ex::is_equal
bool is_equal(const ex &other) const
Definition: ex.h:345
GiNaC::matrix::trace
ex trace() const
Trace of a matrix.
Definition: matrix.cpp:866
m
mvec m
Definition: factor.cpp:771
GiNaC::GINAC_BIND_UNARCHIVER
GINAC_BIND_UNARCHIVER(add)
GiNaC::matrix::solve
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
GiNaC::matrix::cols
unsigned cols() const
Get number of columns.
Definition: matrix.h:77
GiNaC
Definition: add.cpp:38
GiNaC::matrix::add_indexed
ex add_indexed(const ex &self, const ex &other) const override
Sum of two indexed matrices.
Definition: matrix.cpp:397
GiNaC::matrix::print_elements
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
GiNaC::idx::get_dim
ex get_dim() const
Get dimension of index space.
Definition: idx.h:81
GiNaC::basic::is_equal
bool is_equal(const basic &other) const
Test for syntactic equality.
Definition: basic.cpp:863
GiNaC::determinant_algo::laplace
@ laplace
Laplace elimination.
Definition: flags.h:120
GiNaC::info_flags::symbol
@ symbol
Definition: flags.h:246
GiNaC::info_flags::nonnegint
@ nonnegint
Definition: flags.h:231
GiNaC::ex::numer_denom
ex numer_denom() const
Get numerator and denominator of an expression.
Definition: normal.cpp:2564
GiNaC::are_ex_trivially_equal
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
matrix.h
Interface to symbolic matrices.
x
ex x
Definition: factor.cpp:1641
GiNaC::ex::op
ex op(size_t i) const
Definition: ex.h:136
GiNaC::matrix::scalar_mul_indexed
ex scalar_mul_indexed(const ex &self, const numeric &other) const override
Product of an indexed matrix with a number.
Definition: matrix.cpp:433
GiNaC::ex::info
bool info(unsigned inf) const
Definition: ex.h:132
utils.h
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...
GiNaC::matrix::do_print_latex
void do_print_latex(const print_latex &c, unsigned level) const
Definition: matrix.cpp:191
GiNaC::info_flags::crational_polynomial
@ crational_polynomial
Definition: flags.h:259
GiNaC::matrix::charpoly
ex charpoly(const ex &lambda) const
Characteristic Polynomial.
Definition: matrix.cpp:894
GiNaC::ex
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
GiNaC::basic::info
virtual bool info(unsigned inf) const
Information about the object.
Definition: basic.cpp:222
lst.h
Definition of GiNaC's lst.
GiNaC::numeric::is_zero
bool is_zero() const
True if object is zero.
Definition: numeric.cpp:1129
GiNaC::matrix::determinant_minor
ex determinant_minor() const
Recursive determinant for small matrices having at least one symbolic entry.
Definition: matrix.cpp:1096
GiNaC::basic::compare_same_type
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
Definition: basic.cpp:719
GiNaC::matrix::markowitz_elimination
std::vector< unsigned > markowitz_elimination(unsigned n)
Definition: matrix.cpp:1326
GiNaC::matrix::gauss_elimination
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
GiNaC::matrix::echelon_form
std::vector< unsigned > echelon_form(unsigned algo, int n)
Definition: matrix.cpp:1197
GiNaC::solve_algo::bareiss
@ bareiss
Bareiss fraction-free elimination.
Definition: flags.h:185
GiNaC::solve_algo::automatic
@ automatic
Let the system choose.
Definition: flags.h:146
GiNaC::matrix::fraction_free_elimination
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
GiNaC::matrix::inverse
matrix inverse() const
Inverse of this matrix, with automatic algorithm selection.
Definition: matrix.cpp:939
GiNaC::sub_matrix
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
GiNaC::matrix::transpose
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
GiNaC::basic::hold
const basic & hold() const
Stop further evaluation.
Definition: basic.cpp:887
GiNaC::matrix::is_zero_matrix
bool is_zero_matrix() const
Function to check that all elements of the matrix are zero.
Definition: matrix.cpp:1676
GiNaC::print_latex
Context for latex-parsable output.
Definition: print.h:123
symbol.h
Interface to GiNaC's symbolic objects.
GiNaC::matrix::eval_indexed
ex eval_indexed(const basic &i) const override
Automatic symbolic evaluation of an indexed matrix.
Definition: matrix.cpp:320
GiNaC::symbolic_matrix
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
GiNaC::info_flags::numeric
@ numeric
Definition: flags.h:220
GiNaC::matrix::pow
matrix pow(const ex &expn) const
Power of a matrix.
Definition: matrix.cpp:639
GiNaC::matrix::real_part
ex real_part() const override
Definition: matrix.cpp:264
GiNaC::return_types::commutative
@ commutative
Definition: flags.h:280
GiNaC::matrix::operator()
const ex & operator()(unsigned ro, unsigned co) const
operator() to access elements for reading.
Definition: matrix.cpp:686
GiNaC::determinant_algo::gauss
@ gauss
Gauss elimination.
Definition: flags.h:102
GiNaC::basic::setflag
const basic & setflag(unsigned f) const
Set some status_flags.
Definition: basic.h:288
GiNaC::container
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
GiNaC::subs_options::no_pattern
@ no_pattern
disable pattern matching
Definition: flags.h:51
GiNaC::cols
unsigned cols(const matrix &m)
Definition: matrix.h:135
GiNaC::matrix::contract_with
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
GiNaC::permutation_sign
int permutation_sign(It first, It last)
Definition: utils.h:77
GiNaC::basic::ensure_if_modifiable
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
GiNaC::matrix::rank
unsigned rank() const
Compute the rank of this matrix.
Definition: matrix.cpp:1058
GiNaC::solve_algo::gauss
@ gauss
Gauss elimination.
Definition: flags.h:156
c
size_t c
Definition: factor.cpp:770
GiNaC::solve_algo::divfree
@ divfree
Division-free elimination.
Definition: flags.h:168
GiNaC::exmap
std::map< ex, ex, ex_is_less > exmap
Definition: basic.h:50
GiNaC::matrix::imag_part
ex imag_part() const override
Definition: matrix.cpp:273
archive.h
Archiving of GiNaC expressions.
std
Definition: ex.h:972
GiNaC::symbol
Basic CAS symbol.
Definition: symbol.h:39
GiNaC::matrix::archive
void archive(archive_node &n) const override
Save (a.k.a.
Definition: matrix.cpp:152
n
size_t n
Definition: factor.cpp:1463
GiNaC::info_flags::rational_function
@ rational_function
Definition: flags.h:260
GiNaC::matrix::m
exvector m
representation (cols indexed first)
Definition: matrix.h:117
GiNaC::lst_to_matrix
ex lst_to_matrix(const lst &l)
Convert list of lists to matrix.
Definition: matrix.cpp:1684
GiNaC::basic::do_print_tree
void do_print_tree(const print_tree &c, unsigned level) const
Tree output to stream.
Definition: basic.cpp:175
GiNaC::matrix::match_same_type
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
GiNaC::matrix
Symbolic matrices.
Definition: matrix.h:38
GiNaC::matrix::sub
matrix sub(const matrix &other) const
Difference of matrices.
Definition: matrix.cpp:572
GiNaC::is_zero
bool is_zero(const ex &thisex)
Definition: ex.h:820
GiNaC::basic::op
virtual ex op(size_t i) const
Return operand/member at position i.
Definition: basic.cpp:238
GiNaC::basic
This class is the ABC (abstract base class) of GiNaC's class hierarchy.
Definition: basic.h:105
GiNaC::basic::ex
friend class ex
Definition: basic.h:108
GiNaC::matrix::conjugate
ex conjugate() const override
Complex conjugate every matrix entry.
Definition: matrix.cpp:239
GiNaC::matrix::rows
unsigned rows() const
Get number of rows.
Definition: matrix.h:75
GiNaC::matrix::row
unsigned row
number of rows
Definition: matrix.h:115
GiNaC::ex::is_zero
bool is_zero() const
Definition: ex.h:213
GiNaC::unit_matrix
ex unit_matrix(unsigned r, unsigned c)
Create an r times c unit matrix.
Definition: matrix.cpp:1743
indexed.h
Interface to GiNaC's indexed expressions.
GiNaC::matrix::subs
ex subs(const exmap &m, unsigned options=0) const override
Substitute a set of objects by arbitrary expressions.
Definition: matrix.cpp:228
std::swap
void swap(GiNaC::ex &a, GiNaC::ex &b)
Specialization of std::swap() for ex objects.
Definition: ex.h:976
GiNaC::ex::return_type
unsigned return_type() const
Definition: ex.h:230
GiNaC::matrix::op
ex op(size_t i) const override
returns matrix entry at position (i/col, icol).
Definition: matrix.cpp:212
GiNaC::abs
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
GiNaC::idx
This class holds one index of an indexed object.
Definition: idx.h:36
GiNaC::container::nops
size_t nops() const override
Number of operands/members.
Definition: container.h:118
GiNaC::_num1_p
const numeric * _num1_p
Definition: utils.cpp:192
GiNaC::is_dummy_pair
bool is_dummy_pair(const idx &i1, const idx &i2)
Check whether two indices form a dummy pair.
Definition: idx.cpp:502
GiNaC::print_func< print_context >
print_func< print_context >(&varidx::do_print). print_func< print_latex >(&varidx
Definition: idx.cpp:45
GiNaC::matrix::let_op
ex & let_op(size_t i) override
returns writable matrix entry at position (i/col, icol).
Definition: matrix.cpp:220
GiNaC::divide
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
operators.h
Interface to GiNaC's overloaded operators.
GiNaC::solve_algo::markowitz
@ markowitz
Markowitz-ordered Gaussian elimination.
Definition: flags.h:193
GiNaC::numeric
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Definition: numeric.h:82
normal.h
This file defines several functions that work on univariate and multivariate polynomials and rational...
GINAC_ASSERT
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
GiNaC::matrix::add
matrix add(const matrix &other) const
Sum of matrices.
Definition: matrix.cpp:555
GiNaC::_num2_p
const numeric * _num2_p
Definition: utils.cpp:196
GiNaC::solve_algo
Switch to control algorithm for linear system solving.
Definition: flags.h:141
GiNaC::indexed
This class holds an indexed expression.
Definition: indexed.h:40
GiNaC::matrix::mul_scalar
matrix mul_scalar(const ex &other) const
Product of matrix and scalar expression.
Definition: matrix.cpp:623

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