dense_matrix.cpp
1 #include <symengine/matrix.h>
2 #include <symengine/number.h>
3 #include <symengine/add.h>
4 #include <symengine/functions.h>
5 #include <symengine/pow.h>
6 #include <symengine/subs.h>
7 #include <symengine/symengine_exception.h>
8 #include <symengine/polys/uexprpoly.h>
9 #include <symengine/solve.h>
10 #include <symengine/test_visitors.h>
11 
12 namespace SymEngine
13 {
14 
15 // Constructors
16 DenseMatrix::DenseMatrix() {}
17 
18 DenseMatrix::DenseMatrix(unsigned row, unsigned col) : row_(row), col_(col)
19 {
20  m_ = std::vector<RCP<const Basic>>(row * col);
21 }
22 
23 DenseMatrix::DenseMatrix(unsigned row, unsigned col, const vec_basic &l)
24  : m_{l}, row_(row), col_(col){SYMENGINE_ASSERT(m_.size() == row * col)}
25 
26  DenseMatrix::DenseMatrix(const vec_basic &column_elements)
27  : m_(column_elements), row_(static_cast<unsigned>(column_elements.size())),
28  col_(1)
29 {
30 }
31 
32 // Resize the Matrix
33 void DenseMatrix::resize(unsigned row, unsigned col)
34 {
35  row_ = row;
36  col_ = col;
37  m_.resize(row * col);
38 }
39 
40 // Get and set elements
41 RCP<const Basic> DenseMatrix::get(unsigned i, unsigned j) const
42 {
43  SYMENGINE_ASSERT(i < row_ and j < col_);
44  return m_[i * col_ + j];
45 }
46 
47 void DenseMatrix::set(unsigned i, unsigned j, const RCP<const Basic> &e)
48 {
49  SYMENGINE_ASSERT(i < row_ and j < col_);
50  m_[i * col_ + j] = e;
51 }
52 
53 // Flat (row major) view of Matrix
54 vec_basic DenseMatrix::as_vec_basic() const
55 {
56  return m_;
57 }
58 
59 RCP<const Basic> DenseMatrix::trace() const
60 {
61  SYMENGINE_ASSERT(row_ == col_);
62  unsigned offset = 0;
63  vec_basic diag;
64  for (unsigned i = 0; i < row_; i++) {
65  diag.push_back(m_[offset]);
66  offset += row_ + 1;
67  }
68  auto sum = add(diag);
69  return sum;
70 }
71 
72 unsigned DenseMatrix::rank() const
73 {
74  throw NotImplementedError("Not Implemented");
75 }
76 
77 bool DenseMatrix::is_lower() const
78 {
79  auto A = *this;
80  unsigned n = A.nrows();
81  for (unsigned i = 1; i < n; ++i) {
82  for (unsigned j = 0; j < i; ++j) {
83  if (not is_number_and_zero(*A.get(i, j))) {
84  return false;
85  }
86  }
87  }
88  return true;
89 }
90 
91 bool DenseMatrix::is_upper() const
92 {
93  auto A = *this;
94  unsigned n = A.nrows();
95  for (unsigned i = 0; i < n - 1; ++i) {
96  for (unsigned j = i + 1; j < n; ++j) {
97  if (not is_number_and_zero(*A.get(i, j))) {
98  return false;
99  }
100  }
101  }
102  return true;
103 }
104 
105 tribool DenseMatrix::is_zero() const
106 {
107  auto A = *this;
108  tribool cur = tribool::tritrue;
109  for (auto &e : m_) {
110  cur = and_tribool(cur, SymEngine::is_zero(*e));
111  if (is_false(cur)) {
112  return cur;
113  }
114  }
115  return cur;
116 }
117 
118 tribool DenseMatrix::is_diagonal() const
119 {
120  auto A = *this;
121  if (not A.is_square()) {
122  return tribool::trifalse;
123  }
124  unsigned ncols = A.ncols();
125  unsigned offset;
126  tribool cur = tribool::tritrue;
127  for (unsigned i = 0; i < ncols; i++) {
128  offset = i * ncols;
129  for (unsigned j = 0; j < ncols; j++) {
130  if (j != i) {
131  auto &e = m_[offset];
132  cur = and_tribool(cur, SymEngine::is_zero(*e));
133  if (is_false(cur)) {
134  return cur;
135  }
136  }
137  offset++;
138  }
139  }
140  return cur;
141 }
142 
143 tribool DenseMatrix::is_real(const Assumptions *assumptions) const
144 {
145  RealVisitor visitor(assumptions);
146  tribool cur = tribool::tritrue;
147  for (auto &e : m_) {
148  cur = and_tribool(cur, visitor.apply(*e));
149  if (is_false(cur)) {
150  return cur;
151  }
152  }
153  return cur;
154 }
155 
156 tribool DenseMatrix::is_symmetric() const
157 {
158  auto A = *this;
159  if (not A.is_square()) {
160  return tribool::trifalse;
161  }
162  unsigned ncols = A.ncols();
163  tribool cur = tribool::tritrue;
164  for (unsigned i = 0; i < ncols; i++) {
165  for (unsigned j = 0; j <= i; j++) {
166  if (j != i) {
167  auto &e = m_[i * ncols + j];
168  auto &e2 = m_[j * ncols + i];
169  cur = and_tribool(cur, SymEngine::is_zero(*sub(e, e2)));
170  }
171  if (is_false(cur)) {
172  return cur;
173  }
174  }
175  }
176  return cur;
177 }
178 
179 tribool DenseMatrix::is_hermitian() const
180 {
181  auto A = *this;
182  if (not A.is_square()) {
183  return tribool::trifalse;
184  }
185  unsigned ncols = A.ncols();
186  tribool cur = tribool::tritrue;
187  for (unsigned i = 0; i < ncols; i++) {
188  for (unsigned j = 0; j <= i; j++) {
189  auto &e = m_[i * ncols + j];
190  if (j != i) {
191  auto &e2 = m_[j * ncols + i];
192  cur = and_tribool(
194  } else {
195  cur = and_tribool(cur, SymEngine::is_real(*e));
196  }
197  if (is_false(cur)) {
198  return cur;
199  }
200  }
201  }
202  return cur;
203 }
204 
205 tribool DenseMatrix::is_weakly_diagonally_dominant() const
206 {
207  auto A = *this;
208  if (not A.is_square()) {
209  return tribool::trifalse;
210  }
211 
212  unsigned ncols = A.ncols();
213  RCP<const Basic> diag;
214  RCP<const Basic> sum;
215  tribool diagdom = tribool::tritrue;
216  for (unsigned i = 0; i < ncols; i++) {
217  sum = zero;
218  for (unsigned j = 0; j < ncols; j++) {
219  auto &e = m_[i * ncols + j];
220  if (i == j) {
221  diag = abs(e);
222  } else {
223  sum = add(sum, abs(e));
224  }
225  }
226  diagdom = and_tribool(diagdom, is_nonnegative(*sub(diag, sum)));
227  if (is_false(diagdom)) {
228  return diagdom;
229  }
230  }
231  return diagdom;
232 }
233 
234 tribool DenseMatrix::is_strictly_diagonally_dominant() const
235 {
236  auto A = *this;
237  if (not A.is_square()) {
238  return tribool::trifalse;
239  }
240 
241  unsigned ncols = A.ncols();
242  RCP<const Basic> diag;
243  RCP<const Basic> sum;
244  tribool diagdom = tribool::tritrue;
245  for (unsigned i = 0; i < ncols; i++) {
246  sum = zero;
247  for (unsigned j = 0; j < ncols; j++) {
248  auto &e = m_[i * ncols + j];
249  if (i == j) {
250  diag = abs(e);
251  } else {
252  sum = add(sum, abs(e));
253  }
254  }
255  diagdom = and_tribool(diagdom, is_positive(*sub(diag, sum)));
256  if (is_false(diagdom)) {
257  return diagdom;
258  }
259  }
260  return diagdom;
261 }
262 
263 tribool DenseMatrix::shortcut_to_posdef() const
264 {
265  tribool is_diagonal_positive = tribool::tritrue;
266  unsigned offset = 0;
267  for (unsigned i = 0; i < row_; i++) {
268  is_diagonal_positive
269  = and_tribool(is_diagonal_positive, is_positive(*m_[offset]));
270  if (is_false(is_diagonal_positive))
271  return is_diagonal_positive;
272  offset += row_ + 1;
273  }
274  if (is_true(and_tribool(is_diagonal_positive,
275  this->is_strictly_diagonally_dominant())))
276  return tribool::tritrue;
277  return tribool::indeterminate;
278 }
279 
280 tribool DenseMatrix::is_positive_definite_GE()
281 {
282  auto size = row_;
283  for (unsigned i = 0; i < size; i++) {
284  auto ispos = is_positive(*m_[i * size + i]);
285  if (!is_true(ispos))
286  return ispos;
287  for (unsigned j = i + 1; j < size; j++) {
288  for (unsigned k = i + 1; k < size; k++) {
289  m_[j * size + k] = sub(mul(m_[i * size + i], m_[j * size + k]),
290  mul(m_[j * size + i], m_[i * size + k]));
291  }
292  }
293  }
294  return tribool::tritrue;
295 }
296 
297 tribool DenseMatrix::is_positive_definite() const
298 {
299 
300  auto A = *this;
302  const DenseMatrix *H;
303  if (!is_true(A.is_hermitian())) {
304  if (!A.is_square())
305  return tribool::trifalse;
306  DenseMatrix tmp1 = DenseMatrix(A.row_, A.col_);
307  B = std::unique_ptr<DenseMatrix>(new DenseMatrix(A.row_, A.col_));
308  A.conjugate_transpose(tmp1);
309  add_dense_dense(A, tmp1, *B.get());
310  H = B.get();
311  } else {
312  H = this;
313  }
314 
315  tribool shortcut = H->shortcut_to_posdef();
316  if (!is_indeterminate(shortcut))
317  return shortcut;
318 
319  if (!B) {
320  B = std::unique_ptr<DenseMatrix>(new DenseMatrix(A));
321  }
322  return B->is_positive_definite_GE();
323 }
324 
325 tribool DenseMatrix::is_negative_definite() const
326 {
327  auto res = DenseMatrix(row_, col_);
328  mul_dense_scalar(*this, integer(-1), res);
329  return res.is_positive_definite();
330 }
331 
332 RCP<const Basic> DenseMatrix::det() const
333 {
334  return det_bareis(*this);
335 }
336 
337 void DenseMatrix::inv(MatrixBase &result) const
338 {
339  if (is_a<DenseMatrix>(result)) {
340  DenseMatrix &r = down_cast<DenseMatrix &>(result);
341  inverse_pivoted_LU(*this, r);
342  }
343 }
344 
345 void DenseMatrix::add_matrix(const MatrixBase &other, MatrixBase &result) const
346 {
347  SYMENGINE_ASSERT(row_ == result.nrows() and col_ == result.ncols());
348 
349  if (is_a<DenseMatrix>(other) and is_a<DenseMatrix>(result)) {
350  const DenseMatrix &o = down_cast<const DenseMatrix &>(other);
351  DenseMatrix &r = down_cast<DenseMatrix &>(result);
352  add_dense_dense(*this, o, r);
353  }
354 }
355 
356 void DenseMatrix::mul_matrix(const MatrixBase &other, MatrixBase &result) const
357 {
358  SYMENGINE_ASSERT(row_ == result.nrows()
359  and other.ncols() == result.ncols());
360 
361  if (is_a<DenseMatrix>(other) and is_a<DenseMatrix>(result)) {
362  const DenseMatrix &o = down_cast<const DenseMatrix &>(other);
363  DenseMatrix &r = down_cast<DenseMatrix &>(result);
364  mul_dense_dense(*this, o, r);
365  }
366 }
367 
368 void DenseMatrix::elementwise_mul_matrix(const MatrixBase &other,
369  MatrixBase &result) const
370 {
371  SYMENGINE_ASSERT(row_ == result.nrows() and col_ == result.ncols()
372  and row_ == other.nrows() and col_ == other.ncols());
373 
374  if (is_a<DenseMatrix>(other) and is_a<DenseMatrix>(result)) {
375  const DenseMatrix &o = down_cast<const DenseMatrix &>(other);
376  DenseMatrix &r = down_cast<DenseMatrix &>(result);
377  elementwise_mul_dense_dense(*this, o, r);
378  }
379 }
380 
381 // Add a scalar
382 void DenseMatrix::add_scalar(const RCP<const Basic> &k,
383  MatrixBase &result) const
384 {
385  if (is_a<DenseMatrix>(result)) {
386  DenseMatrix &r = down_cast<DenseMatrix &>(result);
387  add_dense_scalar(*this, k, r);
388  }
389 }
390 
391 // Multiply by a scalar
392 void DenseMatrix::mul_scalar(const RCP<const Basic> &k,
393  MatrixBase &result) const
394 {
395  if (is_a<DenseMatrix>(result)) {
396  DenseMatrix &r = down_cast<DenseMatrix &>(result);
397  mul_dense_scalar(*this, k, r);
398  }
399 }
400 
401 // Matrix conjugate
402 void DenseMatrix::conjugate(MatrixBase &result) const
403 {
404  if (is_a<DenseMatrix>(result)) {
405  DenseMatrix &r = down_cast<DenseMatrix &>(result);
406  conjugate_dense(*this, r);
407  }
408 }
409 
410 // Matrix transpose
411 void DenseMatrix::transpose(MatrixBase &result) const
412 {
413  if (is_a<DenseMatrix>(result)) {
414  DenseMatrix &r = down_cast<DenseMatrix &>(result);
415  transpose_dense(*this, r);
416  }
417 }
418 
419 // Matrix conjugate transpose
420 void DenseMatrix::conjugate_transpose(MatrixBase &result) const
421 {
422  if (is_a<DenseMatrix>(result)) {
423  DenseMatrix &r = down_cast<DenseMatrix &>(result);
424  conjugate_transpose_dense(*this, r);
425  }
426 }
427 
428 // Extract out a submatrix
429 void DenseMatrix::submatrix(MatrixBase &result, unsigned row_start,
430  unsigned col_start, unsigned row_end,
431  unsigned col_end, unsigned row_step,
432  unsigned col_step) const
433 {
434  if (is_a<DenseMatrix>(result)) {
435  DenseMatrix &r = down_cast<DenseMatrix &>(result);
436  submatrix_dense(*this, r, row_start, col_start, row_end, col_end,
437  row_step, col_step);
438  }
439 }
440 
441 // LU factorization
442 void DenseMatrix::LU(MatrixBase &L, MatrixBase &U) const
443 {
444  if (is_a<DenseMatrix>(L) and is_a<DenseMatrix>(U)) {
445  DenseMatrix &L_ = down_cast<DenseMatrix &>(L);
446  DenseMatrix &U_ = down_cast<DenseMatrix &>(U);
447  SymEngine::LU(*this, L_, U_);
448  }
449 }
450 
451 // LDL factorization
452 void DenseMatrix::LDL(MatrixBase &L, MatrixBase &D) const
453 {
454  if (is_a<DenseMatrix>(L) and is_a<DenseMatrix>(D)) {
455  DenseMatrix &L_ = down_cast<DenseMatrix &>(L);
456  DenseMatrix &D_ = down_cast<DenseMatrix &>(D);
457  SymEngine::LDL(*this, L_, D_);
458  }
459 }
460 
461 // Solve Ax = b using LU factorization
462 void DenseMatrix::LU_solve(const MatrixBase &b, MatrixBase &x) const
463 {
464  if (is_a<DenseMatrix>(b) and is_a<DenseMatrix>(x)) {
465  const DenseMatrix &b_ = down_cast<const DenseMatrix &>(b);
466  DenseMatrix &x_ = down_cast<DenseMatrix &>(x);
467  SymEngine::LU_solve(*this, b_, x_);
468  }
469 }
470 
471 // Fraction free LU factorization
472 void DenseMatrix::FFLU(MatrixBase &LU) const
473 {
474  if (is_a<DenseMatrix>(LU)) {
475  DenseMatrix &LU_ = down_cast<DenseMatrix &>(LU);
476  fraction_free_LU(*this, LU_);
477  }
478 }
479 
480 // Fraction free LDU factorization
481 void DenseMatrix::FFLDU(MatrixBase &L, MatrixBase &D, MatrixBase &U) const
482 {
483  if (is_a<DenseMatrix>(L) and is_a<DenseMatrix>(D)
484  and is_a<DenseMatrix>(U)) {
485  DenseMatrix &L_ = down_cast<DenseMatrix &>(L);
486  DenseMatrix &D_ = down_cast<DenseMatrix &>(D);
487  DenseMatrix &U_ = down_cast<DenseMatrix &>(U);
488  fraction_free_LDU(*this, L_, D_, U_);
489  }
490 }
491 
492 // QR factorization
493 void DenseMatrix::QR(MatrixBase &Q, MatrixBase &R) const
494 {
495  if (is_a<DenseMatrix>(Q) and is_a<DenseMatrix>(R)) {
496  DenseMatrix &Q_ = down_cast<DenseMatrix &>(Q);
497  DenseMatrix &R_ = down_cast<DenseMatrix &>(R);
498  SymEngine::QR(*this, Q_, R_);
499  }
500 }
501 
502 // Cholesky decomposition
503 void DenseMatrix::cholesky(MatrixBase &L) const
504 {
505  if (is_a<DenseMatrix>(L)) {
506  DenseMatrix &L_ = down_cast<DenseMatrix &>(L);
507  SymEngine::cholesky(*this, L_);
508  }
509 }
510 
511 // ---------------------------- Jacobian -------------------------------------//
512 
513 void jacobian(const DenseMatrix &A, const DenseMatrix &x, DenseMatrix &result,
514  bool diff_cache)
515 {
516  SYMENGINE_ASSERT(A.col_ == 1);
517  SYMENGINE_ASSERT(x.col_ == 1);
518  SYMENGINE_ASSERT(A.row_ == result.nrows() and x.row_ == result.ncols());
519  bool error = false;
520 #pragma omp parallel for
521  for (unsigned i = 0; i < result.row_; i++) {
522  for (unsigned j = 0; j < result.col_; j++) {
523  if (is_a<Symbol>(*(x.m_[j]))) {
524  const RCP<const Symbol> x_
525  = rcp_static_cast<const Symbol>(x.m_[j]);
526  result.m_[i * result.col_ + j] = A.m_[i]->diff(x_, diff_cache);
527  } else {
528  error = true;
529  break;
530  }
531  }
532  }
533  if (error) {
534  throw SymEngineException(
535  "'x' must contain Symbols only. "
536  "Use sjacobian for SymPy style differentiation");
537  }
538 }
539 
540 void sjacobian(const DenseMatrix &A, const DenseMatrix &x, DenseMatrix &result,
541  bool diff_cache)
542 {
543  SYMENGINE_ASSERT(A.col_ == 1);
544  SYMENGINE_ASSERT(x.col_ == 1);
545  SYMENGINE_ASSERT(A.row_ == result.nrows() and x.row_ == result.ncols());
546 #pragma omp parallel for
547  for (unsigned i = 0; i < result.row_; i++) {
548  for (unsigned j = 0; j < result.col_; j++) {
549  if (is_a<Symbol>(*(x.m_[j]))) {
550  const RCP<const Symbol> x_
551  = rcp_static_cast<const Symbol>(x.m_[j]);
552  result.m_[i * result.col_ + j] = A.m_[i]->diff(x_, diff_cache);
553  } else {
554  // TODO: Use a dummy symbol
555  const RCP<const Symbol> x_ = symbol("x_");
556  result.m_[i * result.col_ + j] = ssubs(
557  ssubs(A.m_[i], {{x.m_[j], x_}})->diff(x_, diff_cache),
558  {{x_, x.m_[j]}});
559  }
560  }
561  }
562 }
563 
564 // ---------------------------- Diff -------------------------------------//
565 
566 void diff(const DenseMatrix &A, const RCP<const Symbol> &x, DenseMatrix &result,
567  bool diff_cache)
568 {
569  SYMENGINE_ASSERT(A.row_ == result.nrows() and A.col_ == result.ncols());
570 #pragma omp parallel for
571  for (unsigned i = 0; i < result.row_; i++) {
572  for (unsigned j = 0; j < result.col_; j++) {
573  result.m_[i * result.col_ + j]
574  = A.m_[i * result.col_ + j]->diff(x, diff_cache);
575  }
576  }
577 }
578 
579 void sdiff(const DenseMatrix &A, const RCP<const Basic> &x, DenseMatrix &result,
580  bool diff_cache)
581 {
582  SYMENGINE_ASSERT(A.row_ == result.nrows() and A.col_ == result.ncols());
583 #pragma omp parallel for
584  for (unsigned i = 0; i < result.row_; i++) {
585  for (unsigned j = 0; j < result.col_; j++) {
586  if (is_a<Symbol>(*x)) {
587  const RCP<const Symbol> x_ = rcp_static_cast<const Symbol>(x);
588  result.m_[i * result.col_ + j]
589  = A.m_[i * result.col_ + j]->diff(x_, diff_cache);
590  } else {
591  // TODO: Use a dummy symbol
592  const RCP<const Symbol> x_ = symbol("_x");
593  result.m_[i * result.col_ + j]
594  = ssubs(ssubs(A.m_[i * result.col_ + j], {{x, x_}})
595  ->diff(x_, diff_cache),
596  {{x_, x}});
597  }
598  }
599  }
600 }
601 
602 // ----------------------------- Matrix Conjugate ----------------------------//
603 void conjugate_dense(const DenseMatrix &A, DenseMatrix &B)
604 {
605  SYMENGINE_ASSERT(B.col_ == A.col_ and B.row_ == A.row_);
606 
607  for (unsigned i = 0; i < A.row_; i++)
608  for (unsigned j = 0; j < A.col_; j++)
609  B.m_[i * B.col_ + j] = conjugate(A.m_[i * A.col_ + j]);
610 }
611 
612 // ----------------------------- Matrix Transpose ----------------------------//
613 void transpose_dense(const DenseMatrix &A, DenseMatrix &B)
614 {
615  SYMENGINE_ASSERT(B.row_ == A.col_ and B.col_ == A.row_);
616 
617  for (unsigned i = 0; i < A.row_; i++)
618  for (unsigned j = 0; j < A.col_; j++)
619  B.m_[j * B.col_ + i] = A.m_[i * A.col_ + j];
620 }
621 
622 // ----------------------------- Matrix Conjugate Transpose -----------------//
623 void conjugate_transpose_dense(const DenseMatrix &A, DenseMatrix &B)
624 {
625  SYMENGINE_ASSERT(B.row_ == A.col_ and B.col_ == A.row_);
626 
627  for (unsigned i = 0; i < A.row_; i++)
628  for (unsigned j = 0; j < A.col_; j++)
629  B.m_[j * B.col_ + i] = conjugate(A.m_[i * A.col_ + j]);
630 }
631 
632 // ------------------------------- Submatrix ---------------------------------//
633 void submatrix_dense(const DenseMatrix &A, DenseMatrix &B, unsigned row_start,
634  unsigned col_start, unsigned row_end, unsigned col_end,
635  unsigned row_step, unsigned col_step)
636 {
637  SYMENGINE_ASSERT(row_end >= row_start and col_end >= col_start);
638  SYMENGINE_ASSERT(row_end < A.row_);
639  SYMENGINE_ASSERT(col_end < A.col_);
640  SYMENGINE_ASSERT(B.row_ == row_end - row_start + 1
641  and B.col_ == col_end - col_start + 1);
642 
643  unsigned row = B.row_, col = B.col_;
644 
645  for (unsigned i = 0; i < row; i += row_step)
646  for (unsigned j = 0; j < col; j += col_step)
647  B.m_[i * col + j] = A.m_[(row_start + i) * A.col_ + col_start + j];
648 }
649 
650 // ------------------------------- Matrix Addition ---------------------------//
651 void add_dense_dense(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C)
652 {
653  SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_ and A.row_ == C.row_
654  and A.col_ == C.col_);
655 
656  unsigned row = A.row_, col = A.col_;
657 
658  for (unsigned i = 0; i < row; i++) {
659  for (unsigned j = 0; j < col; j++) {
660  C.m_[i * col + j] = add(A.m_[i * col + j], B.m_[i * col + j]);
661  }
662  }
663 }
664 
665 void add_dense_scalar(const DenseMatrix &A, const RCP<const Basic> &k,
666  DenseMatrix &B)
667 {
668  SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
669 
670  unsigned row = A.row_, col = A.col_;
671 
672  for (unsigned i = 0; i < row; i++) {
673  for (unsigned j = 0; j < col; j++) {
674  B.m_[i * col + j] = add(A.m_[i * col + j], k);
675  }
676  }
677 }
678 
679 // ------------------------------- Matrix Multiplication ---------------------//
680 void mul_dense_dense(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C)
681 {
682  SYMENGINE_ASSERT(A.col_ == B.row_ and C.row_ == A.row_
683  and C.col_ == B.col_);
684 
685  unsigned row = A.row_, col = B.col_;
686 
687  if (&A != &C and &B != &C) {
688  for (unsigned r = 0; r < row; r++) {
689  for (unsigned c = 0; c < col; c++) {
690  C.m_[r * col + c] = zero; // Integer Zero
691  for (unsigned k = 0; k < A.col_; k++)
692  C.m_[r * col + c]
693  = add(C.m_[r * col + c],
694  mul(A.m_[r * A.col_ + k], B.m_[k * col + c]));
695  }
696  }
697  } else {
698  DenseMatrix tmp = DenseMatrix(A.row_, B.col_);
699  mul_dense_dense(A, B, tmp);
700  C = tmp;
701  }
702 }
703 
704 void elementwise_mul_dense_dense(const DenseMatrix &A, const DenseMatrix &B,
705  DenseMatrix &C)
706 {
707  SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_ and A.row_ == C.row_
708  and A.col_ == C.col_);
709 
710  unsigned row = A.row_, col = A.col_;
711 
712  for (unsigned i = 0; i < row; i++) {
713  for (unsigned j = 0; j < col; j++) {
714  C.m_[i * col + j] = mul(A.m_[i * col + j], B.m_[i * col + j]);
715  }
716  }
717 }
718 
719 void mul_dense_scalar(const DenseMatrix &A, const RCP<const Basic> &k,
720  DenseMatrix &B)
721 {
722  SYMENGINE_ASSERT(A.col_ == B.col_ and A.row_ == B.row_);
723 
724  unsigned row = A.row_, col = A.col_;
725 
726  for (unsigned i = 0; i < row; i++) {
727  for (unsigned j = 0; j < col; j++) {
728  B.m_[i * col + j] = mul(A.m_[i * col + j], k);
729  }
730  }
731 }
732 
733 // ---------------------------- Joining Operations ---------------------------//
734 void DenseMatrix::row_join(const DenseMatrix &B)
735 {
736  this->col_insert(B, col_);
737 }
738 
739 void DenseMatrix::col_join(const DenseMatrix &B)
740 {
741  this->row_insert(B, row_);
742 }
743 
744 void DenseMatrix::row_insert(const DenseMatrix &B, unsigned pos)
745 {
746  SYMENGINE_ASSERT(col_ == B.col_ and pos <= row_)
747 
748  unsigned row = row_, col = col_;
749  this->resize(row_ + B.row_, col_);
750 
751  for (unsigned i = row; i-- > pos;) {
752  for (unsigned j = col; j-- > 0;) {
753  this->m_[(i + B.row_) * col + j] = this->m_[i * col + j];
754  }
755  }
756 
757  for (unsigned i = 0; i < B.row_; i++) {
758  for (unsigned j = 0; j < col; j++) {
759  this->m_[(i + pos) * col + j] = B.m_[i * col + j];
760  }
761  }
762 }
763 
764 void DenseMatrix::col_insert(const DenseMatrix &B, unsigned pos)
765 {
766  SYMENGINE_ASSERT(row_ == B.row_ and pos <= col_)
767 
768  unsigned row = row_, col = col_;
769  this->resize(row_, col_ + B.col_);
770 
771  for (unsigned i = row; i-- > 0;) {
772  for (unsigned j = col; j-- > 0;) {
773  if (j >= pos) {
774  this->m_[i * (col + B.col_) + j + B.col_]
775  = this->m_[i * col + j];
776  } else {
777  this->m_[i * (col + B.col_) + j] = this->m_[i * col + j];
778  }
779  }
780  }
781 
782  for (unsigned i = 0; i < row; i++) {
783  for (unsigned j = 0; j < B.col_; j++) {
784  this->m_[i * (col + B.col_) + j + pos] = B.m_[i * B.col_ + j];
785  }
786  }
787 }
788 
789 void DenseMatrix::row_del(unsigned k)
790 {
791  SYMENGINE_ASSERT(k < row_)
792 
793  if (row_ == 1)
794  this->resize(0, 0);
795  else {
796  for (unsigned i = k; i < row_ - 1; i++) {
797  row_exchange_dense(*this, i, i + 1);
798  }
799  this->resize(row_ - 1, col_);
800  }
801 }
802 
803 void DenseMatrix::col_del(unsigned k)
804 {
805  SYMENGINE_ASSERT(k < col_)
806 
807  if (col_ == 1)
808  this->resize(0, 0);
809  else {
810  unsigned row = row_, col = col_, m = 0;
811 
812  for (unsigned i = 0; i < row; i++) {
813  for (unsigned j = 0; j < col; j++) {
814  if (j != k) {
815  this->m_[m] = this->m_[i * col + j];
816  m++;
817  }
818  }
819  }
820  this->resize(row_, col_ - 1);
821  }
822 }
823 
824 // -------------------------------- Row Operations ---------------------------//
825 void row_exchange_dense(DenseMatrix &A, unsigned i, unsigned j)
826 {
827  SYMENGINE_ASSERT(i != j and i < A.row_ and j < A.row_);
828 
829  unsigned col = A.col_;
830 
831  for (unsigned k = 0; k < A.col_; k++)
832  std::swap(A.m_[i * col + k], A.m_[j * col + k]);
833 }
834 
835 void row_mul_scalar_dense(DenseMatrix &A, unsigned i, RCP<const Basic> &c)
836 {
837  SYMENGINE_ASSERT(i < A.row_);
838 
839  unsigned col = A.col_;
840 
841  for (unsigned j = 0; j < A.col_; j++)
842  A.m_[i * col + j] = mul(c, A.m_[i * col + j]);
843 }
844 
845 void row_add_row_dense(DenseMatrix &A, unsigned i, unsigned j,
846  RCP<const Basic> &c)
847 {
848  SYMENGINE_ASSERT(i != j and i < A.row_ and j < A.row_);
849 
850  unsigned col = A.col_;
851 
852  for (unsigned k = 0; k < A.col_; k++)
853  A.m_[i * col + k] = add(A.m_[i * col + k], mul(c, A.m_[j * col + k]));
854 }
855 
856 void permuteFwd(DenseMatrix &A, permutelist &pl)
857 {
858  for (auto &p : pl) {
859  row_exchange_dense(A, p.first, p.second);
860  }
861 }
862 
863 // ----------------------------- Column Operations ---------------------------//
864 void column_exchange_dense(DenseMatrix &A, unsigned i, unsigned j)
865 {
866  SYMENGINE_ASSERT(i != j and i < A.col_ and j < A.col_);
867 
868  unsigned col = A.col_;
869 
870  for (unsigned k = 0; k < A.row_; k++)
871  std::swap(A.m_[k * col + i], A.m_[k * col + j]);
872 }
873 
874 // ------------------------------ Gaussian Elimination -----------------------//
875 void pivoted_gaussian_elimination(const DenseMatrix &A, DenseMatrix &B,
876  permutelist &pl)
877 {
878  SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
879 
880  unsigned row = A.row_, col = A.col_;
881  unsigned index = 0, i, j, k;
882  B.m_ = A.m_;
883 
884  RCP<const Basic> scale;
885 
886  for (i = 0; i < col - 1; i++) {
887  if (index == row)
888  break;
889 
890  k = pivot(B, index, i);
891  if (k == row)
892  continue;
893  if (k != index) {
894  row_exchange_dense(B, k, index);
895  pl.push_back({k, index});
896  }
897 
898  scale = div(one, B.m_[index * col + i]);
899  row_mul_scalar_dense(B, index, scale);
900 
901  for (j = i + 1; j < row; j++) {
902  for (k = i + 1; k < col; k++)
903  B.m_[j * col + k]
904  = sub(B.m_[j * col + k],
905  mul(B.m_[j * col + i], B.m_[i * col + k]));
906  B.m_[j * col + i] = zero;
907  }
908 
909  index++;
910  }
911 }
912 
913 // Algorithm 1, page 12, Nakos, G. C., Turner, P. R., Williams, R. M. (1997).
914 // Fraction-free algorithms for linear and polynomial equations.
915 // ACM SIGSAM Bulletin, 31(3), 11–19. doi:10.1145/271130.271133.
916 void fraction_free_gaussian_elimination(const DenseMatrix &A, DenseMatrix &B)
917 {
918  SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
919 
920  unsigned col = A.col_;
921  B.m_ = A.m_;
922 
923  for (unsigned i = 0; i < col - 1; i++)
924  for (unsigned j = i + 1; j < A.row_; j++) {
925  for (unsigned k = i + 1; k < col; k++) {
926  B.m_[j * col + k]
927  = sub(mul(B.m_[i * col + i], B.m_[j * col + k]),
928  mul(B.m_[j * col + i], B.m_[i * col + k]));
929  if (i > 0)
930  B.m_[j * col + k]
931  = div(B.m_[j * col + k], B.m_[i * col - col + i - 1]);
932  }
933  B.m_[j * col + i] = zero;
934  }
935 }
936 
937 // Pivoted version of `fraction_free_gaussian_elimination`
938 void pivoted_fraction_free_gaussian_elimination(const DenseMatrix &A,
939  DenseMatrix &B, permutelist &pl)
940 {
941  SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
942 
943  unsigned col = A.col_, row = A.row_;
944  unsigned index = 0, i, k, j;
945  B.m_ = A.m_;
946 
947  for (i = 0; i < col - 1; i++) {
948  if (index == row)
949  break;
950 
951  k = pivot(B, index, i);
952  if (k == row)
953  continue;
954  if (k != index) {
955  row_exchange_dense(B, k, index);
956  pl.push_back({k, index});
957  }
958 
959  for (j = i + 1; j < row; j++) {
960  for (k = i + 1; k < col; k++) {
961  B.m_[j * col + k]
962  = sub(mul(B.m_[i * col + i], B.m_[j * col + k]),
963  mul(B.m_[j * col + i], B.m_[i * col + k]));
964  if (i > 0)
965  B.m_[j * col + k]
966  = div(B.m_[j * col + k], B.m_[i * col - col + i - 1]);
967  }
968  B.m_[j * col + i] = zero;
969  }
970 
971  index++;
972  }
973 }
974 
975 void pivoted_gauss_jordan_elimination(const DenseMatrix &A, DenseMatrix &B,
976  permutelist &pl)
977 {
978  SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
979 
980  unsigned row = A.row_, col = A.col_;
981  unsigned index = 0, i, j, k;
982  RCP<const Basic> scale;
983  B.m_ = A.m_;
984 
985  for (i = 0; i < col; i++) {
986  if (index == row)
987  break;
988 
989  k = pivot(B, index, i);
990  if (k == row)
991  continue;
992  if (k != index) {
993  row_exchange_dense(B, k, index);
994  pl.push_back({k, index});
995  }
996 
997  scale = div(one, B.m_[index * col + i]);
998  row_mul_scalar_dense(B, index, scale);
999 
1000  for (j = 0; j < row; j++) {
1001  if (j == index)
1002  continue;
1003 
1004  scale = mul(minus_one, B.m_[j * col + i]);
1005  row_add_row_dense(B, j, index, scale);
1006  }
1007 
1008  index++;
1009  }
1010 }
1011 
1012 // Algorithm 2, page 13, Nakos, G. C., Turner, P. R., Williams, R. M. (1997).
1013 // Fraction-free algorithms for linear and polynomial equations.
1014 // ACM SIGSAM Bulletin, 31(3), 11–19. doi:10.1145/271130.271133.
1015 void fraction_free_gauss_jordan_elimination(const DenseMatrix &A,
1016  DenseMatrix &B)
1017 {
1018  SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
1019 
1020  unsigned row = A.row_, col = A.col_;
1021  unsigned i, j, k;
1022  RCP<const Basic> d;
1023 
1024  B.m_ = A.m_;
1025 
1026  for (i = 0; i < col; i++) {
1027  if (i > 0)
1028  d = B.m_[i * col - col + i - 1];
1029  for (j = 0; j < row; j++)
1030  if (j != i)
1031  for (k = 0; k < col; k++) {
1032  if (k != i) {
1033  B.m_[j * col + k]
1034  = sub(mul(B.m_[i * col + i], B.m_[j * col + k]),
1035  mul(B.m_[j * col + i], B.m_[i * col + k]));
1036  if (i > 0)
1037  B.m_[j * col + k] = div(B.m_[j * col + k], d);
1038  }
1039  }
1040  for (j = 0; j < row; j++)
1041  if (j != i)
1042  B.m_[j * col + i] = zero;
1043  }
1044 }
1045 
1046 void pivoted_fraction_free_gauss_jordan_elimination(const DenseMatrix &A,
1047  DenseMatrix &B,
1048  permutelist &pl)
1049 {
1050  SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
1051 
1052  unsigned row = A.row_, col = A.col_;
1053  unsigned index = 0, i, k, j;
1054  RCP<const Basic> d;
1055 
1056  B.m_ = A.m_;
1057 
1058  for (i = 0; i < col; i++) {
1059  if (index == row)
1060  break;
1061 
1062  k = pivot(B, index, i);
1063  if (k == row)
1064  continue;
1065  if (k != index) {
1066  row_exchange_dense(B, k, index);
1067  pl.push_back({k, index});
1068  }
1069 
1070  for (j = 0; j < row; j++) {
1071  if (j == index)
1072  continue;
1073  for (k = 0; k < col; k++) {
1074  if (k == i)
1075  continue;
1076  B.m_[j * col + k]
1077  = sub(mul(B.m_[index * col + i], B.m_[j * col + k]),
1078  mul(B.m_[j * col + i], B.m_[index * col + k]));
1079  if (index == 0)
1080  continue;
1081  B.m_[j * col + k] = div(B.m_[j * col + k], d);
1082  }
1083  }
1084 
1085  d = B.m_[index * col + i];
1086 
1087  for (j = 0; j < row; j++) {
1088  if (j == index)
1089  continue;
1090  B.m_[j * col + i] = zero;
1091  }
1092  index++;
1093  }
1094 }
1095 
1096 unsigned pivot(DenseMatrix &B, unsigned r, unsigned c)
1097 {
1098  for (unsigned k = r; k < B.row_; k++) {
1099  if (!is_true(is_zero(*(B.m_[k * B.col_ + c])))) {
1100  return k;
1101  }
1102  }
1103  return B.row_;
1104 }
1105 
1106 void reduced_row_echelon_form(const DenseMatrix &A, DenseMatrix &b,
1107  vec_uint &pivot_cols, bool normalize_last)
1108 {
1109  permutelist pl;
1110  if (normalize_last) {
1111  pivoted_fraction_free_gauss_jordan_elimination(A, b, pl);
1112  } else {
1113  pivoted_gauss_jordan_elimination(A, b, pl);
1114  }
1115  unsigned row = 0;
1116  for (unsigned col = 0; col < b.col_ && row < b.row_; col++) {
1117  if (is_true(is_zero(*b.get(row, col))))
1118  continue;
1119  pivot_cols.push_back(col);
1120  if (row == 0 and normalize_last) {
1121  RCP<const Basic> m = div(one, b.get(row, col));
1122  b.mul_scalar(m, b);
1123  }
1124  row++;
1125  }
1126 }
1127 
1128 // --------------------------- Solve Ax = b ---------------------------------//
1129 // Assuming A is a diagonal square matrix
1130 void diagonal_solve(const DenseMatrix &A, const DenseMatrix &b, DenseMatrix &x)
1131 {
1132  SYMENGINE_ASSERT(A.row_ == A.col_);
1133  SYMENGINE_ASSERT(x.row_ == A.col_ and x.col_ == b.col_);
1134 
1135  const unsigned sys = b.col_;
1136 
1137  // No checks are done to see if the diagonal entries are zero
1138  for (unsigned k = 0; k < sys; k++) {
1139  for (unsigned i = 0; i < A.col_; i++) {
1140  x.m_[i * sys + k] = div(b.m_[i * sys + k], A.m_[i * A.col_ + i]);
1141  }
1142  }
1143 }
1144 
1145 void back_substitution(const DenseMatrix &U, const DenseMatrix &b,
1146  DenseMatrix &x)
1147 {
1148  SYMENGINE_ASSERT(U.row_ == U.col_);
1149  SYMENGINE_ASSERT(b.row_ == U.row_);
1150  SYMENGINE_ASSERT(x.row_ == U.col_ and x.col_ == b.col_);
1151 
1152  const unsigned col = U.col_;
1153  const unsigned sys = b.col_;
1154  x.m_ = b.m_;
1155 
1156  for (unsigned k = 0; k < sys; k++) {
1157  for (int i = col - 1; i >= 0; i--) {
1158  for (unsigned j = i + 1; j < col; j++)
1159  x.m_[i * sys + k]
1160  = sub(x.m_[i * sys + k],
1161  mul(U.m_[i * col + j], x.m_[j * sys + k]));
1162  x.m_[i * sys + k] = div(x.m_[i * sys + k], U.m_[i * col + i]);
1163  }
1164  }
1165 }
1166 
1167 void forward_substitution(const DenseMatrix &A, const DenseMatrix &b,
1168  DenseMatrix &x)
1169 {
1170  SYMENGINE_ASSERT(A.row_ == A.col_);
1171  SYMENGINE_ASSERT(b.row_ == A.row_);
1172  SYMENGINE_ASSERT(x.row_ == A.col_ and x.col_ == b.col_);
1173 
1174  unsigned col = A.col_;
1175  const unsigned sys = b.col_;
1176  x.m_ = b.m_;
1177 
1178  for (unsigned k = 0; k < b.col_; k++) {
1179  for (unsigned i = 0; i < col - 1; i++) {
1180  for (unsigned j = i + 1; j < col; j++) {
1181  x.m_[j * sys + k]
1182  = sub(mul(A.m_[i * col + i], x.m_[j * sys + k]),
1183  mul(A.m_[j * col + i], x.m_[i * sys + k]));
1184  if (i > 0)
1185  x.m_[j * sys + k]
1186  = div(x.m_[j * sys + k], A.m_[(i - 1) * col + i - 1]);
1187  }
1188  }
1189  }
1190 }
1191 
1192 void fraction_free_gaussian_elimination_solve(const DenseMatrix &A,
1193  const DenseMatrix &b,
1194  DenseMatrix &x)
1195 {
1196  SYMENGINE_ASSERT(A.row_ == A.col_);
1197  SYMENGINE_ASSERT(b.row_ == A.row_ and x.row_ == A.row_);
1198  SYMENGINE_ASSERT(x.col_ == b.col_);
1199 
1200  int i, j, k, col = A.col_, bcol = b.col_;
1201  DenseMatrix A_ = DenseMatrix(A.row_, A.col_, A.m_);
1202  DenseMatrix b_ = DenseMatrix(b.row_, b.col_, b.m_);
1203 
1204  for (i = 0; i < col - 1; i++)
1205  for (j = i + 1; j < col; j++) {
1206  for (k = 0; k < bcol; k++) {
1207  b_.m_[j * bcol + k]
1208  = sub(mul(A_.m_[i * col + i], b_.m_[j * bcol + k]),
1209  mul(A_.m_[j * col + i], b_.m_[i * bcol + k]));
1210  if (i > 0)
1211  b_.m_[j * bcol + k] = div(b_.m_[j * bcol + k],
1212  A_.m_[i * col - col + i - 1]);
1213  }
1214 
1215  for (k = i + 1; k < col; k++) {
1216  A_.m_[j * col + k]
1217  = sub(mul(A_.m_[i * col + i], A_.m_[j * col + k]),
1218  mul(A_.m_[j * col + i], A_.m_[i * col + k]));
1219  if (i > 0)
1220  A_.m_[j * col + k]
1221  = div(A_.m_[j * col + k], A_.m_[i * col - col + i - 1]);
1222  }
1223  A_.m_[j * col + i] = zero;
1224  }
1225 
1226  for (i = 0; i < col * bcol; i++)
1227  x.m_[i] = zero; // Integer zero;
1228 
1229  for (k = 0; k < bcol; k++) {
1230  for (i = col - 1; i >= 0; i--) {
1231  for (j = i + 1; j < col; j++)
1232  b_.m_[i * bcol + k]
1233  = sub(b_.m_[i * bcol + k],
1234  mul(A_.m_[i * col + j], x.m_[j * bcol + k]));
1235  x.m_[i * bcol + k] = div(b_.m_[i * bcol + k], A_.m_[i * col + i]);
1236  }
1237  }
1238 }
1239 
1240 void fraction_free_gauss_jordan_solve(const DenseMatrix &A,
1241  const DenseMatrix &b, DenseMatrix &x,
1242  bool pivot)
1243 {
1244  SYMENGINE_ASSERT(A.row_ == A.col_);
1245  SYMENGINE_ASSERT(b.row_ == A.row_ and x.row_ == A.row_);
1246  SYMENGINE_ASSERT(x.col_ == b.col_);
1247 
1248  unsigned i, j, k, p, col = A.col_, bcol = b.col_;
1249  RCP<const Basic> d;
1250  DenseMatrix A_ = DenseMatrix(A.row_, A.col_, A.m_);
1251  DenseMatrix b_ = DenseMatrix(b.row_, b.col_, b.m_);
1252 
1253  for (i = 0; i < col; i++) {
1254  if (i > 0)
1255  d = A_.m_[i * col - col + i - 1];
1256  if (pivot) {
1257  p = i;
1258  while (p < col and eq(*A_.m_[p * col + i], *zero)) {
1259  p++;
1260  }
1261  SYMENGINE_ASSERT(p != col);
1262  if (p != i) {
1263  // pivot A
1264  for (k = i; k < col; k++) {
1265  std::swap(A_.m_[p * col + k], A_.m_[i * col + k]);
1266  }
1267  for (k = 0; k < bcol; k++) {
1268  std::swap(b_.m_[p * bcol + k], b_.m_[i * bcol + k]);
1269  }
1270  }
1271  }
1272  for (j = 0; j < col; j++)
1273  if (j != i) {
1274  for (k = 0; k < bcol; k++) {
1275  b_.m_[j * bcol + k]
1276  = sub(mul(A_.m_[i * col + i], b_.m_[j * bcol + k]),
1277  mul(A_.m_[j * col + i], b_.m_[i * bcol + k]));
1278  if (i > 0)
1279  b_.m_[j * bcol + k] = div(b_.m_[j * bcol + k], d);
1280  }
1281 
1282  for (k = 0; k < col; k++) {
1283  if (k != i) {
1284  A_.m_[j * col + k]
1285  = sub(mul(A_.m_[i * col + i], A_.m_[j * col + k]),
1286  mul(A_.m_[j * col + i], A_.m_[i * col + k]));
1287  if (i > 0)
1288  A_.m_[j * col + k] = div(A_.m_[j * col + k], d);
1289  }
1290  }
1291  }
1292 
1293  for (j = 0; j < col; j++)
1294  if (j != i)
1295  A_.m_[j * col + i] = zero;
1296  }
1297 
1298  // No checks are done to see if the diagonal entries are zero
1299  for (k = 0; k < bcol; k++)
1300  for (i = 0; i < col; i++)
1301  x.m_[i * bcol + k] = div(b_.m_[i * bcol + k], A_.m_[i * col + i]);
1302 }
1303 
1304 void fraction_free_LU_solve(const DenseMatrix &A, const DenseMatrix &b,
1305  DenseMatrix &x)
1306 {
1307  DenseMatrix LU = DenseMatrix(A.nrows(), A.ncols());
1308  DenseMatrix x_ = DenseMatrix(b.nrows(), b.ncols());
1309 
1310  fraction_free_LU(A, LU);
1311  forward_substitution(LU, b, x_);
1312  back_substitution(LU, x_, x);
1313 }
1314 
1315 void LU_solve(const DenseMatrix &A, const DenseMatrix &b, DenseMatrix &x)
1316 {
1317  DenseMatrix L = DenseMatrix(A.nrows(), A.ncols());
1318  DenseMatrix U = DenseMatrix(A.nrows(), A.ncols());
1319  DenseMatrix x_ = DenseMatrix(b.nrows(), b.ncols());
1320 
1321  LU(A, L, U);
1322  forward_substitution(L, b, x_);
1323  back_substitution(U, x_, x);
1324 }
1325 
1326 void pivoted_LU_solve(const DenseMatrix &A, const DenseMatrix &b,
1327  DenseMatrix &x)
1328 {
1329  DenseMatrix L = DenseMatrix(A.nrows(), A.ncols());
1330  DenseMatrix U = DenseMatrix(A.nrows(), A.ncols());
1331  DenseMatrix x_ = DenseMatrix(b);
1332  permutelist pl;
1333 
1334  pivoted_LU(A, L, U, pl);
1335  permuteFwd(x_, pl);
1336  forward_substitution(L, x_, x_);
1337  back_substitution(U, x_, x);
1338 }
1339 
1340 void LDL_solve(const DenseMatrix &A, const DenseMatrix &b, DenseMatrix &x)
1341 {
1342  DenseMatrix L = DenseMatrix(A.nrows(), A.ncols());
1343  DenseMatrix D = DenseMatrix(A.nrows(), A.ncols());
1344  DenseMatrix x_ = DenseMatrix(b.nrows(), b.ncols());
1345 
1346  if (not is_symmetric_dense(A))
1347  throw SymEngineException("Matrix must be symmetric");
1348 
1349  LDL(A, L, D);
1350  forward_substitution(L, b, x);
1351  diagonal_solve(D, x, x_);
1352  transpose_dense(L, D);
1353  back_substitution(D, x_, x);
1354 }
1355 
1356 // --------------------------- Matrix Decomposition --------------------------//
1357 
1358 // Algorithm 3, page 14, Nakos, G. C., Turner, P. R., Williams, R. M. (1997).
1359 // Fraction-free algorithms for linear and polynomial equations.
1360 // ACM SIGSAM Bulletin, 31(3), 11–19. doi:10.1145/271130.271133.
1361 // This algorithms is not a true factorization of the matrix A(i.e. A != LU))
1362 // but can be used to solve linear systems by applying forward and backward
1363 // substitutions respectively.
1364 void fraction_free_LU(const DenseMatrix &A, DenseMatrix &LU)
1365 {
1366  SYMENGINE_ASSERT(A.row_ == A.col_ and LU.row_ == LU.col_
1367  and A.row_ == LU.row_);
1368 
1369  unsigned n = A.row_;
1370  unsigned i, j, k;
1371 
1372  LU.m_ = A.m_;
1373 
1374  for (i = 0; i < n - 1; i++)
1375  for (j = i + 1; j < n; j++)
1376  for (k = i + 1; k < n; k++) {
1377  LU.m_[j * n + k] = sub(mul(LU.m_[i * n + i], LU.m_[j * n + k]),
1378  mul(LU.m_[j * n + i], LU.m_[i * n + k]));
1379  if (i)
1380  LU.m_[j * n + k]
1381  = div(LU.m_[j * n + k], LU.m_[i * n - n + i - 1]);
1382  }
1383 }
1384 
1385 // SymPy LUDecomposition algorithm, in
1386 // sympy.matrices.matrices.Matrix.LUdecomposition
1387 // with no pivoting
1388 void LU(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &U)
1389 {
1390  SYMENGINE_ASSERT(A.row_ == A.col_ and L.row_ == L.col_
1391  and U.row_ == U.col_);
1392  SYMENGINE_ASSERT(A.row_ == L.row_ and A.row_ == U.row_);
1393 
1394  unsigned n = A.row_;
1395  unsigned i, j, k;
1396  RCP<const Basic> scale;
1397 
1398  U.m_ = A.m_;
1399 
1400  for (j = 0; j < n; j++) {
1401  for (i = 0; i < j; i++)
1402  for (k = 0; k < i; k++)
1403  U.m_[i * n + j] = sub(U.m_[i * n + j],
1404  mul(U.m_[i * n + k], U.m_[k * n + j]));
1405 
1406  for (i = j; i < n; i++) {
1407  for (k = 0; k < j; k++)
1408  U.m_[i * n + j] = sub(U.m_[i * n + j],
1409  mul(U.m_[i * n + k], U.m_[k * n + j]));
1410  }
1411 
1412  scale = div(one, U.m_[j * n + j]);
1413 
1414  for (i = j + 1; i < n; i++)
1415  U.m_[i * n + j] = mul(U.m_[i * n + j], scale);
1416  }
1417 
1418  for (i = 0; i < n; i++) {
1419  for (j = 0; j < i; j++) {
1420  L.m_[i * n + j] = U.m_[i * n + j];
1421  U.m_[i * n + j] = zero; // Integer zero
1422  }
1423  L.m_[i * n + i] = one; // Integer one
1424  for (j = i + 1; j < n; j++)
1425  L.m_[i * n + j] = zero; // Integer zero
1426  }
1427 }
1428 
1429 // SymPy LUDecomposition_Simple algorithm, in
1430 // sympy.matrices.matrices.Matrix.LUdecomposition_Simple with pivoting.
1431 // P must be an initialized matrix and will be permuted.
1432 void pivoted_LU(const DenseMatrix &A, DenseMatrix &LU, permutelist &pl)
1433 {
1434  SYMENGINE_ASSERT(A.row_ == A.col_ and LU.row_ == LU.col_);
1435  SYMENGINE_ASSERT(A.row_ == LU.row_);
1436 
1437  unsigned n = A.row_;
1438  unsigned i, j, k;
1439  RCP<const Basic> scale;
1440  int pivot;
1441 
1442  LU.m_ = A.m_;
1443 
1444  for (j = 0; j < n; j++) {
1445  for (i = 0; i < j; i++)
1446  for (k = 0; k < i; k++)
1447  LU.m_[i * n + j] = sub(LU.m_[i * n + j],
1448  mul(LU.m_[i * n + k], LU.m_[k * n + j]));
1449  pivot = -1;
1450  for (i = j; i < n; i++) {
1451  for (k = 0; k < j; k++) {
1452  LU.m_[i * n + j] = sub(LU.m_[i * n + j],
1453  mul(LU.m_[i * n + k], LU.m_[k * n + j]));
1454  }
1455  if (pivot == -1 and !is_true(is_zero(*LU.m_[i * n + j])))
1456  pivot = i;
1457  }
1458  if (pivot == -1)
1459  throw SymEngineException("Matrix is rank deficient");
1460  if (pivot - j != 0) { // row must be swapped
1461  row_exchange_dense(LU, pivot, j);
1462  pl.push_back({pivot, j});
1463  }
1464  scale = div(one, LU.m_[j * n + j]);
1465 
1466  for (i = j + 1; i < n; i++)
1467  LU.m_[i * n + j] = mul(LU.m_[i * n + j], scale);
1468  }
1469 }
1470 
1471 // SymPy LUDecomposition algorithm, in
1472 // sympy.matrices.matrices.Matrix.LUdecomposition_Simple with pivoting.
1473 // P must be an initialized matrix and will be permuted.
1474 void pivoted_LU(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &U,
1475  permutelist &pl)
1476 {
1477  SYMENGINE_ASSERT(A.row_ == A.col_ and L.row_ == L.col_
1478  and U.row_ == U.col_);
1479  SYMENGINE_ASSERT(A.row_ == L.row_ and A.row_ == U.row_);
1480 
1481  pivoted_LU(A, U, pl);
1482  unsigned n = A.col_;
1483  for (unsigned i = 0; i < n; i++) {
1484  for (unsigned j = 0; j < i; j++) {
1485  L.m_[i * n + j] = U.m_[i * n + j];
1486  U.m_[i * n + j] = zero; // Integer zero
1487  }
1488  L.m_[i * n + i] = one; // Integer one
1489  for (unsigned j = i + 1; j < n; j++)
1490  L.m_[i * n + j] = zero; // Integer zero
1491  }
1492 }
1493 
1494 // SymPy's fraction free LU decomposition, without pivoting
1495 // sympy.matrices.matrices.MatrixBase.LUdecompositionFF
1496 // W. Zhou & D.J. Jeffrey, "Fraction-free matrix factors: new forms for LU and
1497 // QR factors".
1498 // Frontiers in Computer Science in China, Vol 2, no. 1, pp. 67-80, 2008.
1499 void fraction_free_LDU(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &D,
1500  DenseMatrix &U)
1501 {
1502  SYMENGINE_ASSERT(A.row_ == L.row_ and A.row_ == U.row_);
1503  SYMENGINE_ASSERT(A.col_ == L.col_ and A.col_ == U.col_);
1504 
1505  unsigned row = A.row_, col = A.col_;
1506  unsigned i, j, k;
1507  RCP<const Basic> old = integer(1);
1508 
1509  U.m_ = A.m_;
1510 
1511  // Initialize L
1512  for (i = 0; i < row; i++)
1513  for (j = 0; j < row; j++)
1514  if (i != j)
1515  L.m_[i * col + j] = zero;
1516  else
1517  L.m_[i * col + i] = one;
1518 
1519  // Initialize D
1520  for (i = 0; i < row * row; i++)
1521  D.m_[i] = zero; // Integer zero
1522 
1523  for (k = 0; k < row - 1; k++) {
1524  L.m_[k * col + k] = U.m_[k * col + k];
1525  D.m_[k * col + k] = mul(old, U.m_[k * col + k]);
1526 
1527  for (i = k + 1; i < row; i++) {
1528  L.m_[i * col + k] = U.m_[i * col + k];
1529  for (j = k + 1; j < col; j++)
1530  U.m_[i * col + j]
1531  = div(sub(mul(U.m_[k * col + k], U.m_[i * col + j]),
1532  mul(U.m_[k * col + j], U.m_[i * col + k])),
1533  old);
1534  U.m_[i * col + k] = zero; // Integer zero
1535  }
1536 
1537  old = U.m_[k * col + k];
1538  }
1539 
1540  D.m_[row * col - col + row - 1] = old;
1541 }
1542 
1543 // SymPy's QRecomposition in sympy.matrices.matrices.MatrixBase.QRdecomposition
1544 // Rank check is not performed
1545 void QR(const DenseMatrix &A, DenseMatrix &Q, DenseMatrix &R)
1546 {
1547  unsigned row = A.row_;
1548  unsigned col = A.col_;
1549 
1550  SYMENGINE_ASSERT(Q.row_ == row and Q.col_ == col and R.row_ == col
1551  and R.col_ == col);
1552 
1553  unsigned i, j, k;
1554  RCP<const Basic> t;
1555  std::vector<RCP<const Basic>> tmp(row);
1556 
1557  // Initialize Q
1558  for (i = 0; i < row * col; i++)
1559  Q.m_[i] = zero;
1560 
1561  // Initialize R
1562  for (i = 0; i < col * col; i++)
1563  R.m_[i] = zero;
1564 
1565  for (j = 0; j < col; j++) {
1566  // Use submatrix for this
1567  for (k = 0; k < row; k++)
1568  tmp[k] = A.m_[k * col + j];
1569 
1570  for (i = 0; i < j; i++) {
1571  t = zero;
1572  for (k = 0; k < row; k++)
1573  t = add(t, mul(A.m_[k * col + j], Q.m_[k * col + i]));
1574  for (k = 0; k < row; k++)
1575  tmp[k] = expand(sub(tmp[k], mul(Q.m_[k * col + i], t)));
1576  }
1577 
1578  // calculate norm
1579  t = zero;
1580  for (k = 0; k < row; k++)
1581  t = add(t, pow(tmp[k], integer(2)));
1582 
1583  t = pow(t, div(one, integer(2)));
1584 
1585  R.m_[j * col + j] = t;
1586  for (k = 0; k < row; k++)
1587  Q.m_[k * col + j] = div(tmp[k], t);
1588 
1589  for (i = 0; i < j; i++) {
1590  t = zero;
1591  for (k = 0; k < row; k++)
1592  t = add(t, mul(Q.m_[k * col + i], A.m_[k * col + j]));
1593  R.m_[i * col + j] = t;
1594  }
1595  }
1596 }
1597 
1598 // SymPy's LDL decomposition, Assuming A is a symmetric, square, positive
1599 // definite non singular matrix
1600 void LDL(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &D)
1601 {
1602  SYMENGINE_ASSERT(A.row_ == A.col_);
1603  SYMENGINE_ASSERT(L.row_ == A.row_ and L.col_ == A.row_);
1604  SYMENGINE_ASSERT(D.row_ == A.row_ and D.col_ == A.row_);
1605 
1606  unsigned col = A.col_;
1607  unsigned i, k, j;
1608  RCP<const Basic> sum;
1609  RCP<const Basic> i2 = integer(2);
1610 
1611  // Initialize D
1612  for (i = 0; i < col; i++)
1613  for (j = 0; j < col; j++)
1614  D.m_[i * col + j] = zero; // Integer zero
1615 
1616  // Initialize L
1617  for (i = 0; i < col; i++)
1618  for (j = 0; j < col; j++)
1619  L.m_[i * col + j] = (i != j) ? zero : one;
1620 
1621  for (i = 0; i < col; i++) {
1622  for (j = 0; j < i; j++) {
1623  sum = zero;
1624  for (k = 0; k < j; k++)
1625  sum = add(sum, mul(mul(L.m_[i * col + k], L.m_[j * col + k]),
1626  D.m_[k * col + k]));
1627  L.m_[i * col + j]
1628  = mul(div(one, D.m_[j * col + j]), sub(A.m_[i * col + j], sum));
1629  }
1630  sum = zero;
1631  for (k = 0; k < i; k++)
1632  sum = add(sum, mul(pow(L.m_[i * col + k], i2), D.m_[k * col + k]));
1633  D.m_[i * col + i] = sub(A.m_[i * col + i], sum);
1634  }
1635 }
1636 
1637 // SymPy's cholesky decomposition
1638 void cholesky(const DenseMatrix &A, DenseMatrix &L)
1639 {
1640  SYMENGINE_ASSERT(A.row_ == A.col_);
1641  SYMENGINE_ASSERT(L.row_ == A.row_ and L.col_ == A.row_);
1642 
1643  unsigned col = A.col_;
1644  unsigned i, j, k;
1645  RCP<const Basic> sum;
1646  RCP<const Basic> i2 = integer(2);
1647  RCP<const Basic> half = div(one, i2);
1648 
1649  // Initialize L
1650  for (i = 0; i < col; i++)
1651  for (j = 0; j < col; j++)
1652  L.m_[i * col + j] = zero;
1653 
1654  for (i = 0; i < col; i++) {
1655  for (j = 0; j < i; j++) {
1656  sum = zero;
1657  for (k = 0; k < j; k++)
1658  sum = add(sum, mul(L.m_[i * col + k], L.m_[j * col + k]));
1659 
1660  L.m_[i * col + j]
1661  = mul(div(one, L.m_[j * col + j]), sub(A.m_[i * col + j], sum));
1662  }
1663  sum = zero;
1664  for (k = 0; k < i; k++)
1665  sum = add(sum, pow(L.m_[i * col + k], i2));
1666  L.m_[i * col + i] = pow(sub(A.m_[i * col + i], sum), half);
1667  }
1668 }
1669 
1670 // Matrix Queries
1671 bool is_symmetric_dense(const DenseMatrix &A)
1672 {
1673  if (A.col_ != A.row_)
1674  return false;
1675 
1676  unsigned col = A.col_;
1677  bool sym = true;
1678 
1679  for (unsigned i = 0; i < col; i++)
1680  for (unsigned j = i + 1; j < col; j++)
1681  if (not eq(*(A.m_[j * col + i]), *(A.m_[i * col + j]))) {
1682  sym = false;
1683  break;
1684  }
1685 
1686  return sym;
1687 }
1688 
1689 // ----------------------------- Determinant ---------------------------------//
1690 RCP<const Basic> det_bareis(const DenseMatrix &A)
1691 {
1692  SYMENGINE_ASSERT(A.row_ == A.col_);
1693 
1694  unsigned n = A.row_;
1695 
1696  if (n == 1) {
1697  return A.m_[0];
1698  } else if (n == 2) {
1699  // If A = [[a, b], [c, d]] then det(A) = ad - bc
1700  return sub(mul(A.m_[0], A.m_[3]), mul(A.m_[1], A.m_[2]));
1701  } else if (n == 3) {
1702  // if A = [[a, b, c], [d, e, f], [g, h, i]] then
1703  // det(A) = (aei + bfg + cdh) - (ceg + bdi + afh)
1704  return sub(add(add(mul(mul(A.m_[0], A.m_[4]), A.m_[8]),
1705  mul(mul(A.m_[1], A.m_[5]), A.m_[6])),
1706  mul(mul(A.m_[2], A.m_[3]), A.m_[7])),
1707  add(add(mul(mul(A.m_[2], A.m_[4]), A.m_[6]),
1708  mul(mul(A.m_[1], A.m_[3]), A.m_[8])),
1709  mul(mul(A.m_[0], A.m_[5]), A.m_[7])));
1710  } else {
1711 
1712  if (A.is_lower() or A.is_upper()) {
1713  RCP<const Basic> det = A.m_[0];
1714  for (unsigned i = 1; i < n; ++i) {
1715  det = mul(det, A.m_[i * n + i]);
1716  }
1717  return det;
1718  }
1719 
1720  DenseMatrix B = DenseMatrix(n, n, A.m_);
1721  unsigned i, sign = 1;
1722  RCP<const Basic> d;
1723 
1724  for (unsigned k = 0; k < n - 1; k++) {
1725  if (is_true(is_zero(*B.m_[k * n + k]))) {
1726  for (i = k + 1; i < n; i++)
1727  if (!is_true(is_zero(*B.m_[i * n + k]))) {
1728  row_exchange_dense(B, i, k);
1729  sign *= -1;
1730  break;
1731  }
1732  if (i == n)
1733  return zero;
1734  }
1735 
1736  for (i = k + 1; i < n; i++) {
1737  for (unsigned j = k + 1; j < n; j++) {
1738  d = sub(mul(B.m_[k * n + k], B.m_[i * n + j]),
1739  mul(B.m_[i * n + k], B.m_[k * n + j]));
1740  if (k > 0)
1741  d = div(d, B.m_[(k - 1) * n + k - 1]);
1742  B.m_[i * n + j] = d;
1743  }
1744  }
1745  }
1746 
1747  return (sign == 1) ? B.m_[n * n - 1] : mul(minus_one, B.m_[n * n - 1]);
1748  }
1749 }
1750 
1751 // Returns the coefficients of characterisitc polynomials of leading principal
1752 // minors of Matrix A as elements in `polys`. Principal leading minor of kth
1753 // order is the submatrix of A obtained by deleting last n-k rows and columns
1754 // from A. Here `n` is the dimension of the square matrix A.
1755 void berkowitz(const DenseMatrix &A, std::vector<DenseMatrix> &polys)
1756 {
1757  SYMENGINE_ASSERT(A.row_ == A.col_);
1758 
1759  unsigned col = A.col_;
1760  unsigned i, k, l, m;
1761 
1763  std::vector<DenseMatrix> transforms;
1765 
1766  for (unsigned n = col; n > 1; n--) {
1767  items.clear();
1768  k = n - 1;
1769  DenseMatrix T = DenseMatrix(n + 1, n);
1770  DenseMatrix C = DenseMatrix(k, 1);
1771 
1772  // Initialize T and C
1773  for (i = 0; i < n * (n + 1); i++)
1774  T.m_[i] = zero;
1775  for (i = 0; i < k; i++)
1776  C.m_[i] = A.m_[i * col + k];
1777  items.push_back(C);
1778 
1779  for (i = 0; i < n - 2; i++) {
1780  DenseMatrix B = DenseMatrix(k, 1);
1781  for (l = 0; l < k; l++) {
1782  B.m_[l] = zero;
1783  for (m = 0; m < k; m++)
1784  B.m_[l]
1785  = add(B.m_[l], mul(A.m_[l * col + m], items[i].m_[m]));
1786  }
1787  items.push_back(B);
1788  }
1789 
1790  items_.clear();
1791  for (i = 0; i < n - 1; i++) {
1792  RCP<const Basic> element = zero;
1793  for (l = 0; l < k; l++)
1794  element = add(element, mul(A.m_[k * col + l], items[i].m_[l]));
1795  items_.push_back(mul(minus_one, element));
1796  }
1797  items_.insert(items_.begin(), mul(minus_one, A.m_[k * col + k]));
1798  items_.insert(items_.begin(), one);
1799 
1800  for (i = 0; i < n; i++) {
1801  for (l = 0; l < n - i + 1; l++)
1802  T.m_[(i + l) * n + i] = items_[l];
1803  }
1804 
1805  transforms.push_back(T);
1806  }
1807 
1808  polys.push_back(DenseMatrix(2, 1, {one, mul(A.m_[0], minus_one)}));
1809 
1810  for (i = 0; i < col - 1; i++) {
1811  unsigned t_row = transforms[col - 2 - i].nrows();
1812  unsigned t_col = transforms[col - 2 - i].ncols();
1813  DenseMatrix B = DenseMatrix(t_row, 1);
1814 
1815  for (l = 0; l < t_row; l++) {
1816  B.m_[l] = zero;
1817  for (m = 0; m < t_col; m++) {
1818  B.m_[l] = add(B.m_[l],
1819  mul(transforms[col - 2 - i].m_[l * t_col + m],
1820  polys[i].m_[m]));
1821  B.m_[l] = expand(B.m_[l]);
1822  }
1823  }
1824  polys.push_back(B);
1825  }
1826 }
1827 
1828 RCP<const Basic> det_berkowitz(const DenseMatrix &A)
1829 {
1831 
1832  berkowitz(A, polys);
1833  DenseMatrix poly = polys[polys.size() - 1];
1834 
1835  if (polys.size() % 2 == 1)
1836  return mul(minus_one, poly.get(poly.nrows() - 1, 0));
1837 
1838  return poly.get(poly.nrows() - 1, 0);
1839 }
1840 
1841 void char_poly(const DenseMatrix &A, DenseMatrix &B)
1842 {
1843  SYMENGINE_ASSERT(B.ncols() == 1 and B.nrows() == A.nrows() + 1);
1844  SYMENGINE_ASSERT(A.nrows() == A.ncols());
1845 
1847 
1848  berkowitz(A, polys);
1849  B = polys[polys.size() - 1];
1850 }
1851 
1852 void inverse_fraction_free_LU(const DenseMatrix &A, DenseMatrix &B)
1853 {
1854  SYMENGINE_ASSERT(A.row_ == A.col_ and B.row_ == B.col_
1855  and B.row_ == A.row_);
1856 
1857  unsigned n = A.row_, i;
1858  DenseMatrix LU = DenseMatrix(n, n);
1859  DenseMatrix e = DenseMatrix(n, 1);
1860  DenseMatrix x = DenseMatrix(n, 1);
1861  DenseMatrix x_ = DenseMatrix(n, 1);
1862 
1863  // Initialize matrices
1864  for (i = 0; i < n * n; i++) {
1865  LU.m_[i] = zero;
1866  B.m_[i] = zero;
1867  }
1868  for (i = 0; i < n; i++) {
1869  e.m_[i] = zero;
1870  x.m_[i] = zero;
1871  x_.m_[i] = zero;
1872  }
1873 
1874  fraction_free_LU(A, LU);
1875 
1876  // We solve AX_{i} = e_{i} for i = 1, 2, .. n and combine the row vectors
1877  // X_{1}, X_{2}, ... X_{n} to form the inverse of A. Here, e_{i}'s are the
1878  // elements of the standard basis.
1879  for (unsigned j = 0; j < n; j++) {
1880  e.m_[j] = one;
1881 
1882  forward_substitution(LU, e, x_);
1883  back_substitution(LU, x_, x);
1884 
1885  for (i = 0; i < n; i++)
1886  B.m_[i * n + j] = x.m_[i];
1887 
1888  e.m_[j] = zero;
1889  }
1890 }
1891 
1892 void inverse_LU(const DenseMatrix &A, DenseMatrix &B)
1893 {
1894  SYMENGINE_ASSERT(A.row_ == A.col_ and B.row_ == B.col_
1895  and B.row_ == A.row_);
1896  DenseMatrix e = DenseMatrix(A.row_, A.col_);
1897  eye(e);
1898  LU_solve(A, e, B);
1899 }
1900 
1901 void inverse_pivoted_LU(const DenseMatrix &A, DenseMatrix &B)
1902 {
1903  SYMENGINE_ASSERT(A.row_ == A.col_ and B.row_ == B.col_
1904  and B.row_ == A.row_);
1905  DenseMatrix e = DenseMatrix(A.row_, A.col_);
1906  eye(e);
1907  pivoted_LU_solve(A, e, B);
1908 }
1909 
1910 void inverse_gauss_jordan(const DenseMatrix &A, DenseMatrix &B)
1911 {
1912  SYMENGINE_ASSERT(A.row_ == A.col_ and B.row_ == B.col_
1913  and B.row_ == A.row_);
1914 
1915  unsigned n = A.row_;
1916  DenseMatrix e = DenseMatrix(n, n);
1917 
1918  // Initialize matrices
1919  for (unsigned i = 0; i < n; i++)
1920  for (unsigned j = 0; j < n; j++) {
1921  if (i != j) {
1922  e.m_[i * n + j] = zero;
1923  } else {
1924  e.m_[i * n + i] = one;
1925  }
1926  B.m_[i * n + j] = zero;
1927  }
1928 
1929  fraction_free_gauss_jordan_solve(A, e, B);
1930 }
1931 
1932 // ----------------------- Vector-specific Methods --------------------------//
1933 
1934 void dot(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C)
1935 {
1936  if (A.col_ == B.row_) {
1937  if (B.col_ != 1) {
1938  DenseMatrix tmp1 = DenseMatrix(A.col_, A.row_);
1939  A.transpose(tmp1);
1940  DenseMatrix tmp2 = DenseMatrix(B.col_, B.row_);
1941  B.transpose(tmp2);
1942  C.resize(tmp1.row_, tmp2.col_);
1943  mul_dense_dense(tmp1, tmp2, C);
1944  } else {
1945  C.resize(A.row_, B.col_);
1946  mul_dense_dense(A, B, C);
1947  }
1948  C.resize(1, C.row_ * C.col_);
1949  } else if (A.col_ == B.col_) {
1950  DenseMatrix tmp2 = DenseMatrix(B.col_, B.row_);
1951  B.transpose(tmp2);
1952  dot(A, tmp2, C);
1953  } else if (A.row_ == B.row_) {
1954  DenseMatrix tmp1 = DenseMatrix(A.col_, A.row_);
1955  A.transpose(tmp1);
1956  dot(tmp1, B, C);
1957  } else {
1958  throw SymEngineException("Dimensions incorrect for dot product");
1959  }
1960 }
1961 
1962 void cross(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C)
1963 {
1964  SYMENGINE_ASSERT((A.row_ * A.col_ == 3 and B.row_ * B.col_ == 3)
1965  and (A.row_ == C.row_ and A.col_ == C.col_));
1966 
1967  C.m_[0] = sub(mul(A.m_[1], B.m_[2]), mul(A.m_[2], B.m_[1]));
1968  C.m_[1] = sub(mul(A.m_[2], B.m_[0]), mul(A.m_[0], B.m_[2]));
1969  C.m_[2] = sub(mul(A.m_[0], B.m_[1]), mul(A.m_[1], B.m_[0]));
1970 }
1971 
1972 RCP<const Set> eigen_values(const DenseMatrix &A)
1973 {
1974  unsigned n = A.nrows();
1975  if (A.is_lower() or A.is_upper()) {
1976  RCP<const Set> eigenvals = emptyset();
1977  set_basic x;
1978  for (unsigned i = 0; i < n; ++i) {
1979  x.insert(A.get(i, i));
1980  }
1981  eigenvals = finiteset(x);
1982  return eigenvals;
1983  }
1984 
1985  DenseMatrix B = DenseMatrix(A.nrows() + 1, 1);
1986  char_poly(A, B);
1987  map_int_Expr coeffs;
1988  auto nr = A.nrows();
1989  for (unsigned i = 0; i <= nr; i++)
1990  insert(coeffs, nr - i, B.get(i, 0));
1991  auto lambda = symbol("lambda");
1992  return solve_poly(uexpr_poly(lambda, coeffs), lambda);
1993 }
1994 
1995 // ------------------------- NumPy-like functions ----------------------------//
1996 
1997 // Mimic `eye` function in NumPy
1998 void eye(DenseMatrix &A, int k)
1999 {
2000  if ((k >= 0 and (unsigned) k >= A.col_) or k + A.row_ <= 0) {
2001  zeros(A);
2002  }
2003 
2004  vec_basic v = vec_basic(k > 0 ? A.col_ - k : A.row_ + k, one);
2005 
2006  diag(A, v, k);
2007 }
2008 
2009 // Create diagonal matrices directly
2010 void diag(DenseMatrix &A, vec_basic &v, int k)
2011 {
2012  SYMENGINE_ASSERT(v.size() > 0);
2013 
2014  unsigned k_ = std::abs(k);
2015 
2016  if (k >= 0) {
2017  for (unsigned i = 0; i < A.row_; i++) {
2018  for (unsigned j = 0; j < A.col_; j++) {
2019  if (j != (unsigned)k) {
2020  A.m_[i * A.col_ + j] = zero;
2021  } else {
2022  A.m_[i * A.col_ + j] = v[k - k_];
2023  }
2024  }
2025  k++;
2026  }
2027  } else {
2028  k = -k;
2029 
2030  for (unsigned j = 0; j < A.col_; j++) {
2031  for (unsigned i = 0; i < A.row_; i++) {
2032  if (i != (unsigned)k) {
2033  A.m_[i * A.col_ + j] = zero;
2034  } else {
2035  A.m_[i * A.col_ + j] = v[k - k_];
2036  }
2037  }
2038  k++;
2039  }
2040  }
2041 }
2042 
2043 // Create a matrix filled with ones
2044 void ones(DenseMatrix &A)
2045 {
2046  for (unsigned i = 0; i < A.row_ * A.col_; i++) {
2047  A.m_[i] = one;
2048  }
2049 }
2050 
2051 // Create a matrix filled with zeros
2052 void zeros(DenseMatrix &A)
2053 {
2054  for (unsigned i = 0; i < A.row_ * A.col_; i++) {
2055  A.m_[i] = zero;
2056  }
2057 }
2058 
2059 } // namespace SymEngine
Classes and functions relating to the binary operation of addition.
T begin(T... args)
RCP< const Basic > sub(const RCP< const Basic > &a, const RCP< const Basic > &b)
Substracts b from a.
Definition: add.cpp:495
RCP< const Basic > add(const RCP< const Basic > &a, const RCP< const Basic > &b)
Adds two objects (safely).
Definition: add.cpp:425
T clear(T... args)
T div(T... args)
T get(T... args)
T insert(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
std::enable_if< std::is_integral< T >::value, RCP< const Integer > >::type integer(T i)
Definition: integer.h:200
bool is_number_and_zero(const Basic &b)
Definition: number.h:139
RCP< const Symbol > symbol(const std::string &name)
inline version to return Symbol
Definition: symbol.h:82
bool eq(const Basic &a, const Basic &b)
Checks equality for a and b
Definition: basic-inl.h:21
RCP< const EmptySet > emptyset()
Definition: sets.h:454
RCP< const Basic > sign(const RCP< const Basic > &arg)
Canonicalize Sign.
Definition: functions.cpp:484
RCP< const Basic > abs(const RCP< const Basic > &arg)
Canonicalize Abs:
Definition: functions.cpp:3449
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:347
void insert(T1 &m, const T2 &first, const T3 &second)
Definition: dict.h:83
tribool is_zero(const Basic &b, const Assumptions *assumptions=nullptr)
Check if a number is zero.
RCP< const Set > finiteset(const set_basic &container)
Definition: sets.h:466
RCP< const Basic > expand(const RCP< const Basic > &self, bool deep=true)
Expands self
Definition: expand.cpp:369
RCP< const Basic > ssubs(const RCP< const Basic > &x, const map_basic_basic &subs_dict, bool cache=true)
SymPy compatible subs.
Definition: subs.h:543
RCP< const Basic > conjugate(const RCP< const Basic > &arg)
Canonicalize Conjugate.
Definition: functions.cpp:106
T pow(T... args)
T push_back(T... args)
T size(T... args)
T swap(T... args)