71 matrix::matrix(
unsigned r,
unsigned c) : row(
r), col(
c),
m(
r*
c,
_ex0)
100 : row(l.size()), col(l.begin()->size())
105 for (
const auto &
r : l) {
107 for (
const auto & e :
r) {
112 throw std::invalid_argument(
"matrix::matrix{{}}: wrong dimension");
120 : row(
r), col(
c),
m(m2)
125 : row(
r), col(
c),
m(
std::move(m2))
136 inherited::read_archive(
n, sym_lst);
138 if (!(
n.find_unsigned(
"row",
row)) || !(
n.find_unsigned(
"col",
col)))
139 throw (std::runtime_error(
"unknown matrix dimensions in archive"));
143 auto range =
n.find_property_range(
"m",
"m");
144 for (
auto i=range.begin; i != range.end; ++i) {
146 n.find_ex_by_loc(i, e, sym_lst);
154 inherited::archive(
n);
155 n.add_unsigned(
"row",
row);
156 n.add_unsigned(
"col",
col);
170 for (
unsigned ro=0; ro<
row; ++ro) {
172 for (
unsigned co=0; co<
col; ++co) {
173 m[ro*
col+co].print(
c);
193 c.s <<
"\\left(\\begin{array}{" << std::string(
col,
'c') <<
"}";
195 c.s <<
"\\end{array}\\right)";
200 c.s << class_name() <<
'(';
208 return static_cast<size_t>(
row) * static_cast<size_t>(
col);
231 for (
unsigned r=0;
r<
row; ++
r)
232 for (
unsigned c=0;
c<
col; ++
c)
241 std::unique_ptr<exvector> ev(
nullptr);
242 for (
auto i=
m.begin(); i!=
m.end(); ++i) {
243 ex x = i->conjugate();
252 ev->reserve(
m.size());
253 for (
auto j=
m.begin(); j!=i; ++j) {
269 v.push_back(i.real_part());
278 v.push_back(i.imag_part());
291 return row < o.
rows() ? -1 : 1;
295 return col < o.
cols() ? -1 : 1;
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;
332 throw (std::runtime_error(
"matrix::eval_indexed(): vector must have exactly 1 index"));
334 const idx & i1 = ex_to<idx>(i.
op(1));
340 throw (std::runtime_error(
"matrix::eval_indexed(): dimension of index must match number of vector elements"));
343 if (all_indices_unsigned) {
346 throw (std::runtime_error(
"matrix::eval_indexed(): value of index exceeds number of vector elements"));
347 return (*
this)(n1, 0);
354 throw (std::runtime_error(
"matrix::eval_indexed(): dimension of index must match number of vector elements"));
357 if (all_indices_unsigned) {
360 throw (std::runtime_error(
"matrix::eval_indexed(): value of index exceeds number of vector elements"));
361 return (*
this)(0, n1);
365 }
else if (i.
nops() == 3) {
368 const idx & i1 = ex_to<idx>(i.
op(1));
369 const idx & i2 = ex_to<idx>(i.
op(2));
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"));
381 if (all_indices_unsigned) {
384 throw (std::runtime_error(
"matrix::eval_indexed(): value of first index exceeds number of rows"));
386 throw (std::runtime_error(
"matrix::eval_indexed(): value of second index exceeds number of columns"));
387 return (*
this)(n1, n2);
391 throw (std::runtime_error(
"matrix::eval_indexed(): matrix must have exactly 2 indices"));
405 if (is_a<matrix>(other.
op(0))) {
408 const matrix &self_matrix = ex_to<matrix>(
self.op(0));
409 const matrix &other_matrix = ex_to<matrix>(other.
op(0));
411 if (
self.
nops() == 2 && other.
nops() == 2) {
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));
418 }
else if (
self.
nops() == 3 && other.
nops() == 3) {
421 return indexed(self_matrix.add(other_matrix),
self.op(1),
self.op(2));
423 return indexed(self_matrix.add(other_matrix.transpose()),
self.
op(1),
self.
op(2));
439 const matrix &self_matrix = ex_to<matrix>(
self.op(0));
441 if (
self.
nops() == 2)
456 if (!is_a<matrix>(other->op(0)))
461 const matrix &self_matrix = ex_to<matrix>(
self->op(0));
462 const matrix &other_matrix = ex_to<matrix>(other->op(0));
464 if (self->nops() == 2) {
466 if (other->nops() == 2) {
468 if (self_matrix.col == 1) {
469 if (other_matrix.col == 1) {
471 *
self = self_matrix.
transpose().
mul(other_matrix)(0, 0);
474 *
self = other_matrix.
mul(self_matrix)(0, 0);
477 if (other_matrix.col == 1) {
479 *
self = self_matrix.
mul(other_matrix)(0, 0);
482 *
self = self_matrix.
mul(other_matrix.transpose())(0, 0);
492 if (self_matrix.row == 1)
493 *
self =
indexed(self_matrix.mul(other_matrix), other->op(2));
495 *
self =
indexed(self_matrix.transpose().mul(other_matrix), other->op(2));
502 if (self_matrix.col == 1)
503 *
self =
indexed(other_matrix.mul(self_matrix), other->op(1));
505 *
self =
indexed(other_matrix.mul(self_matrix.transpose()), other->op(1));
511 }
else if (other->nops() == 3) {
515 *
self =
indexed(self_matrix.mul(other_matrix),
self->op(1), other->op(2));
522 *
self =
indexed(self_matrix.mul(other_matrix.transpose()), self->op(1), other->op(1));
529 *
self =
indexed(self_matrix.transpose().mul(other_matrix),
self->op(2), other->op(2));
536 *
self =
indexed(other_matrix.mul(self_matrix), other->op(1),
self->op(2));
558 throw std::logic_error(
"matrix::add(): incompatible matrices");
561 auto ci = other.
m.begin();
575 throw std::logic_error(
"matrix::sub(): incompatible matrices");
578 auto ci = other.
m.begin();
592 throw std::logic_error(
"matrix::mul(): incompatible matrices");
596 for (
unsigned r1=0; r1<this->
rows(); ++r1) {
597 for (
unsigned c=0;
c<this->
cols(); ++
c) {
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]);
614 for (
unsigned r=0;
r<
row; ++
r)
615 for (
unsigned c=0;
c<
col; ++
c)
626 throw std::runtime_error(
"matrix::mul_scalar(): non-commutative scalar");
630 for (
unsigned r=0;
r<
row; ++
r)
631 for (
unsigned c=0;
c<
col; ++
c)
642 throw (std::logic_error(
"matrix::pow(): matrix not square"));
644 if (is_exactly_a<numeric>(expn)) {
648 numeric b = ex_to<numeric>(expn);
657 for (
unsigned r=0;
r<
row; ++
r)
677 throw (std::runtime_error(
"matrix::pow(): don't know how to handle exponent"));
689 throw (std::range_error(
"matrix::operator(): index out of range"));
703 throw (std::range_error(
"matrix::operator(): index out of range"));
716 for (
unsigned r=0;
r<this->
cols(); ++
r)
717 for (
unsigned c=0;
c<this->
rows(); ++
c)
740 throw (std::logic_error(
"matrix::determinant(): matrix not square"));
744 bool numeric_flag =
true;
745 bool normal_flag =
false;
746 unsigned sparse_count = 0;
749 numeric_flag =
false;
751 ex rtest =
r.to_rational(srl);
777 return m[0].normal();
779 return m[0].expand();
788 for (
unsigned d=0; d<
row; ++d)
789 det *= tmp.
m[d*
col+d];
791 return (sign*det).normal();
793 return (sign*det).normal().expand();
800 return (sign*tmp.
m[
row*
col-1]).normal();
802 return (sign*tmp.
m[
row*
col-1]).expand();
812 for (
unsigned d=0; d<
row-2; ++d)
813 for (
unsigned j=0; j<
row-d-2; ++j)
828 typedef std::pair<unsigned,unsigned> uintpair;
829 std::vector<uintpair> c_zeros;
830 for (
unsigned c=0;
c<
col; ++
c) {
832 for (
unsigned r=0;
r<
row; ++
r)
835 c_zeros.push_back(uintpair(acc,
c));
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);
845 for (
auto & it : pre_sort) {
846 for (
unsigned r=0;
r<
row; ++
r)
854 return sign*
matrix(
row,
col, std::move(result)).determinant_minor();
869 throw (std::logic_error(
"matrix::trace(): matrix not square"));
872 for (
unsigned r=0;
r<
col; ++
r)
897 throw (std::logic_error(
"matrix::charpoly(): matrix not square"));
899 bool numeric_flag =
true;
902 numeric_flag =
false;
915 for (
unsigned i=1; i<
row; ++i) {
916 for (
unsigned j=0; j<
row; ++j)
930 for (
unsigned r=0;
r<
col; ++
r)
953 throw (std::logic_error(
"matrix::inverse(): matrix not square"));
960 for (
unsigned i=0; i<
row; ++i)
961 identity(i,i) =
_ex1;
967 for (
unsigned r=0;
r<
row; ++
r)
968 for (
unsigned c=0;
c<
col; ++
c)
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"));
999 const unsigned m = this->
rows();
1000 const unsigned n = this->
cols();
1001 const unsigned p =
rhs.cols();
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)
1009 throw (std::invalid_argument(
"matrix::solve(): 1st argument must be matrix of symbols"));
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)
1025 for (
unsigned co=0; co<p; ++co) {
1026 unsigned last_assigned_sol =
n+1;
1027 for (
int r=
m-1;
r>=0; --
r) {
1029 while ((fnz<=
n) && (aug.
m[
r*(
n+p)+(fnz-1)].normal().is_zero()))
1033 if (!aug.
m[
r*(
n+p)+
n+co].normal().is_zero()) {
1034 throw (std::runtime_error(
"matrix::solve(): inconsistent linear system"));
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;
1050 for (
unsigned ro=0; ro<last_assigned_sol-1; ++ro)
1051 sol(colid[ro],co) = vars(colid[ro],co);
1072 matrix to_eliminate = *
this;
1077 if (!to_eliminate.m[
r].is_zero())
1098 const unsigned n = this->
cols();
1133 typedef std::vector<unsigned> keyseq;
1134 typedef std::map<keyseq, ex> Rmap;
1147 for (
int c=
n-1;
c>=0; --
c) {
1150 for (
unsigned i=0; i<
n-
c; ++i)
1155 for (
unsigned r=0;
r<
n-
c; ++
r) {
1161 Mkey.insert(Mkey.begin(), Nkey.begin(), Nkey.begin() +
r);
1162 Mkey.insert(Mkey.end(), Nkey.begin() +
r + 1, Nkey.end());
1165 det -=
m[Nkey[
r]*
n+
c]*M[Mkey];
1167 det +=
m[Nkey[
r]*
n+
c]*M[Mkey];
1176 for (fc=
n-
c; fc>0; --fc) {
1178 if (Nkey[fc-1]<fc+
c)
1182 for (
unsigned j=fc; j<
n-
c; ++j)
1183 Nkey[j] = Nkey[j-1]+1;
1196 std::vector<unsigned>
1202 bool numeric_flag =
true;
1203 for (
const auto &
r :
m) {
1205 numeric_flag =
false;
1209 unsigned density = 0;
1210 for (
const auto &
r :
m) {
1211 density += !
r.is_zero();
1213 unsigned ncells =
col*
row;
1217 if ((ncells > 200) && (density < ncells/2)) {
1225 if ((ncells < 120) && (density*5 > ncells*3)) {
1237 std::vector<unsigned> colid(
col);
1238 for (
unsigned c = 0;
c <
col;
c++) {
1255 throw std::invalid_argument(
"matrix::echelon_form(): 'algo' is not one of the solve_algo enum");
1272 const unsigned m = this->
rows();
1273 const unsigned n = this->
cols();
1278 for (
unsigned c0=0; c0<
n && r0<
m-1; ++c0) {
1279 int indx =
pivot(r0, c0,
true);
1288 for (
unsigned r2=r0+1; r2<
m; ++r2) {
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];
1299 for (
unsigned c=r0;
c<=c0; ++
c)
1304 for (
unsigned c=r0+1;
c<
n; ++
c)
1311 for (
unsigned r=r0+1;
r<
m; ++
r) {
1312 for (
unsigned c=0;
c<
n; ++
c)
1325 std::vector<unsigned>
1329 std::vector<int> rowcnt(
row, 0);
1330 std::vector<int> colcnt(
col, 0);
1334 for (
unsigned r = 0;
r <
row;
r++) {
1335 for (
unsigned c = 0;
c <
col;
c++) {
1343 std::vector<unsigned> colid(
col);
1344 for (
unsigned c = 0;
c <
col;
c++) {
1348 for (
unsigned k = 0; (
k <
col) && (
k <
row - 1);
k++) {
1350 unsigned pivot_r =
row + 1;
1351 unsigned pivot_c =
col + 1;
1353 for (
unsigned r =
k;
r <
row;
r++) {
1354 for (
unsigned c =
k;
c <
n;
c++) {
1360 int measure = (rowcnt[
r] - 1)*(colcnt[
c] - 1);
1361 if (measure < pivot_m) {
1368 if (pivot_m ==
row*
col) {
1376 for (
unsigned r = 0;
r <
row;
r++) {
1383 for (
unsigned c =
k;
c <
col;
c++) {
1397 for (
unsigned r =
k + 1;
r <
row;
r++) {
1404 colcnt[
k] = rowcnt[
k] = 0;
1405 for (
unsigned c =
k + 1;
c <
col;
c++) {
1410 for (
unsigned r =
k + 1;
r <
row;
r++) {
1413 bool waszero =
m[
r*
col +
c].is_zero();
1415 bool iszero =
m[
r*
col +
c].is_zero();
1416 if (waszero && !iszero) {
1420 if (!waszero && iszero) {
1426 for (
unsigned r =
k + 1;
r <
row;
r++) {
1444 const unsigned m = this->
rows();
1445 const unsigned n = this->
cols();
1450 for (
unsigned c0=0; c0<
n && r0<
m-1; ++c0) {
1451 int indx =
pivot(r0, c0,
true);
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();
1464 for (
unsigned c=r0;
c<=c0; ++
c)
1469 for (
unsigned c=r0+1;
c<
n; ++
c)
1476 for (
unsigned r=r0+1;
r<
m; ++
r) {
1477 for (
unsigned c=0;
c<
n; ++
c)
1522 const unsigned m = this->
rows();
1523 const unsigned n = this->
cols();
1543 auto tmp_n_it = tmp_n.
m.begin(), tmp_d_it = tmp_d.
m.begin();
1544 for (
auto & it : this->m) {
1546 *tmp_n_it++ = nd.
op(0);
1547 *tmp_d_it++ = nd.
op(1);
1551 for (
unsigned c0=0; c0<
n && r0<
m-1; ++c0) {
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]);
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]*
1580 dividend_d = (tmp_d.
m[r2*
n+c0]*tmp_d.
m[r0*
n+
c]*
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);
1589 for (
unsigned c=r0;
c<=c0; ++
c)
1592 if (c0<
n && r0<
m-1) {
1594 divisor_n = tmp_n.
m[r0*
n+c0].expand();
1595 divisor_d = tmp_d.
m[r0*
n+c0].expand();
1598 for (
unsigned c=0;
c<
n; ++
c) {
1608 for (
unsigned r=r0+1;
r<
m; ++
r) {
1609 for (
unsigned c=0;
c<
n; ++
c)
1614 tmp_n_it = tmp_n.
m.begin();
1615 tmp_d_it = tmp_d.
m.begin();
1616 for (
auto & it : this->m)
1646 unsigned kmax =
k+1;
1650 numeric tmp = ex_to<numeric>(this->
m[kmax*
col+co]);
1651 if (
abs(tmp) > mmax) {
1668 for (
unsigned c=0;
c<
col; ++
c)
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)
1699 for (
auto & itr : l) {
1701 for (
auto & itc : ex_to<lst>(itr)) {
1713 size_t dim = l.
nops();
1716 matrix & M = dynallocate<matrix>(dim, dim);
1719 for (
auto & it : l) {
1729 size_t dim = l.size();
1732 matrix & M = dynallocate<matrix>(dim, dim);
1735 for (
auto & it : l) {
1745 matrix & Id = dynallocate<matrix>(
r,
c);
1747 for (
unsigned i=0; i<
r && i<
c; i++)
1755 matrix & M = dynallocate<matrix>(
r,
c);
1758 bool long_format = (
r > 10 ||
c > 10);
1759 bool single_row = (
r == 1 ||
c == 1);
1761 for (
unsigned i=0; i<
r; i++) {
1762 for (
unsigned j=0; j<
c; j++) {
1763 std::ostringstream s1, s2;
1765 s2 << tex_base_name <<
"_{";
1776 s1 <<
'_' << i <<
'_' << j;
1777 s2 << i <<
';' << j <<
"}";
1780 s2 << i << j <<
'}';
1783 M(i, j) =
symbol(s1.str(), s2.str());
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");
1795 const unsigned rows =
m.rows()-1;
1796 const unsigned cols =
m.cols()-1;
1810 M(ro2,co2) =
m(ro, co);
1823 if (
r+nr>
m.rows() ||
c+nc>
m.cols())
1824 throw std::runtime_error(
"sub_matrix(): index out of bounds");
1826 matrix & M = dynallocate<matrix>(nr, nc);
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);
virtual ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const
Default implementation of ex::normal().
matrix transpose() const
Transposed of an m x n matrix, producing a new n x m matrix object that represents the transposed...
ex & let_op(size_t i) override
returns writable matrix entry at position (i/col, icol).
matrix sub(const matrix &other) const
Difference of matrices.
ex eval_indexed(const basic &i) const override
Automatic symbolic evaluation of an indexed matrix.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
ex collect(const ex &s, bool distributed=false) const
ex to_rational(exmap &repl) const
Rationalization of non-rational functions.
ex normal() const
Normalization of rational functions.
Interface to GiNaC's indexed expressions.
Interface to GiNaC's symbolic objects.
ex scalar_mul_indexed(const ex &self, const numeric &other) const override
Product of an indexed matrix with a number.
bool match_same_type(const basic &other) const override
Returns true if the attributes of two objects are similar enough for a match.
virtual ex expand(unsigned options=0) const
Expand expression, i.e.
int to_int(const numeric &x)
void do_print(const print_context &c, unsigned level) const
unsigned rows(const matrix &m)
ex diag_matrix(const lst &l)
Convert list of diagonal elements to matrix.
unsigned rank() const
Compute the rank of this matrix.
don't share instances of this object between different expressions unless explicitly asked to (used b...
const basic & hold() const
Stop further evaluation.
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.
ex subs(const exmap &m, unsigned options=0) const override
Substitute a set of objects by arbitrary expressions.
bool is_equal(const ex &other) const
const basic & setflag(unsigned f) const
Set some status_flags.
ex expand(unsigned options=0) const
int permutation_sign(It first, It last)
ex get_value() const
Get value of index.
ex unit_matrix(unsigned r, unsigned c)
Create an r times c unit matrix.
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...
Bareiss 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...
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...
Archiving of GiNaC expressions.
Interface to GiNaC's sums of expressions.
This class is the ABC (abstract base class) of GiNaC's class hierarchy.
Context for python-parsable output.
void print_elements(const print_context &c, const char *row_start, const char *row_end, const char *row_sep, const char *col_sep) const
bool are_ex_trivially_equal(const ex &e1, const ex &e2)
Compare two objects of class quickly without doing a deep tree traversal.
void swap(ex &e1, ex &e2)
matrix mul_scalar(const ex &other) const
Product of matrix and scalar expression.
Division-free elimination.
matrix add(const matrix &other) const
Sum of matrices.
virtual size_t nops() const
Number of operands/members.
bool contract_with(exvector::iterator self, exvector::iterator other, exvector &v) const override
Contraction of an indexed matrix with something else.
unsigned row
number of rows
Interface to GiNaC's indices.
const ex & operator()(unsigned ro, unsigned co) const
operator() to access elements for reading.
bool info(unsigned inf) const
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...
int pivot(unsigned ro, unsigned co, bool symbolic=true)
Partial pivoting method for matrix elimination schemes.
bool is_zero(const ex &thisex)
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].
Bareiss fraction-free elimination.
unsigned cols(const matrix &m)
unsigned cols() const
Get number of columns.
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Context for latex-parsable output.
const numeric mul(const numeric &other) const
Numerical multiplication method.
print_func< print_context >(&varidx::do_print). print_func< print_latex >(&varidx
This class holds an indexed expression.
virtual ex op(size_t i) const
Return operand/member at position i.
size_t nops() const override
Number of operands/members.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
bool is_zero_matrix() const
Function to check that all elements of the matrix are zero.
matrix mul(const matrix &other) const
Product of matrices.
ex conjugate() const override
Complex conjugate every matrix entry.
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...
exvector m
representation (cols indexed first)
matrix(unsigned r, unsigned c)
Very common ctor.
bool is_odd() const
True if object is an exact odd integer.
ex numer_denom() const
Get numerator and denominator of an expression.
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
bool is_dummy_pair(const idx &i1, const idx &i2)
Check whether two indices form a dummy pair.
matrix pow(const ex &expn) const
Power of a matrix.
Interface to symbolic matrices.
const numeric abs(const numeric &x)
Absolute value.
std::vector< ex > exvector
std::map< ex, ex, ex_is_less > exmap
Interface to GiNaC's overloaded operators.
Division-free elimination.
ex add_indexed(const ex &self, const ex &other) const override
Sum of two indexed matrices.
unsigned col
number of columns
ex determinant(unsigned algo=determinant_algo::automatic) const
Determinant of square matrix.
Base class for print_contexts.
Definition of GiNaC's lst.
void ensure_if_modifiable() const
Ensure the object may be modified without hurting others, throws if this is not the case...
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...
bool is_zero() const
True if object is zero.
void do_print_tree(const print_tree &c, unsigned level) const
Tree output to stream.
void do_print_latex(const print_latex &c, unsigned level) const
This class stores all properties needed to record/retrieve the state of one object of class basic (or...
Lightweight wrapper for GiNaC's symbolic objects.
This class holds a two-component object, a basis and and exponent representing exponentiation.
ex determinant_minor() const
Recursive determinant for small matrices having at least one symbolic entry.
ex get_dim() const
Get dimension of index space.
void do_print_python_repr(const print_python_repr &c, unsigned level) const
ex imag_part() const override
ex op(size_t i) const override
returns matrix entry at position (i/col, icol).
ex trace() const
Trace of a matrix.
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
unsigned return_type() const
matrix inverse() const
Inverse of this matrix, with automatic algorithm selection.
Wrapper template for making GiNaC classes out of STL containers.
void read_archive(const archive_node &n, lst &syms) override
Read (a.k.a.
This file defines several functions that work on univariate and multivariate polynomials and rational...
std::vector< unsigned > echelon_form(unsigned algo, int n)
This class holds one index of an indexed object.
size_t nops() const override
nops is defined to be rows x columns.
Makes the interface to the underlying bignum package available.
Switch to control algorithm for linear system solving.
ex lst_to_matrix(const lst &l)
Convert list of lists to 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...
virtual bool info(unsigned inf) const
Information about the object.
void archive(archive_node &n) const override
Save (a.k.a.
unsigned rows() const
Get number of rows.
ex real_part() const override
ex charpoly(const ex &lambda) const
Characteristic Polynomial.
.eval() has already done its job
std::vector< unsigned > markowitz_elimination(unsigned n)
void swap(GiNaC::ex &a, GiNaC::ex &b)
Specialization of std::swap() for ex objects.
GINAC_BIND_UNARCHIVER(add)
Markowitz-ordered Gaussian elimination.