Polymetic  1.1
A c++ library for polynomial and matrix arithmetic, focused on applications in Kinematics.
 All Classes Namespaces Files Functions Variables Typedefs Friends Macros Pages
MatrixAlgorithms_impl.ipp
Go to the documentation of this file.
1 // Copyright 2018 Dhruvesh Nikhilkumar Patel
2 //
3 // Licensed under the Apache License, Version 2.0 (the "License");
4 // you may not use this file except in compliance with the License.
5 // You may obtain a copy of the License at
6 //
7 // http://www.apache.org/licenses/LICENSE-2.0
8 //
9 // Unless required by applicable law or agreed to in writing, software
10 // distributed under the License is distributed on an "AS IS" BASIS,
11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 // See the License for the specific language governing permissions and
13 // limitations under the License.
14 //
18 #ifndef _MATRIXALGORITHMS_IMPL_IPP
19 #define _MATRIXALGORITHMS_IMPL_IPP
20 
21 #include "../include/MatrixAlgorithms.hpp"
22 #include <memory>
23 #include <stdexcept>
24 #include <vector>
25 
26 namespace {
27  template <typename MatrixT>
28  MatrixT exp_by_squaring(const MatrixT& M,size_t n)
29  {
30  if (!isSquare(M)) {
31  throw (std::invalid_argument("Input to exp_by_squaring is not a square matrix"));
32  }
33  switch(n) {
34  case 0: {
35  throw (std::invalid_argument("exp_by_squaring: 0 is not a valid argument for exponent"));
36  //return Identity_Matrix<MatrixT::value_type>(n);
37  }
38  case 1: {
39  return MatrixT(M);
40  }
41  default:{
42  if (n%2 == 0) {//even
43  return exp_by_squaring<MatrixT>(prod(M,M), n/2);
44  }
45  else {// odd
46  return prod(M,exp_by_squaring<MatrixT>(prod(M,M),(n-1)/2));
47  }
48  }
49  }
50  }
51 
52  template <typename MatrixT>
53  typename MatrixT::value_type trace(const MatrixT& M, size_t n)
54  {
55  if (!isSquare(M)) {
56  throw (std::invalid_argument("Input to exp_by_squaring is not a square matrix"));
57  }
58  std::unique_ptr<MatrixT> pB;
59  std::unique_ptr<MatrixT> pA;
60 
61  /* divide M^n into two almost equal parts */
62  if(n%2 ==0) {// even
63  pB = std::make_unique<MatrixT>(exp(M,n/2));
64  pA = std::make_unique<MatrixT>(trans(*pB));
65  // the output object of trans(B) is passed as rvalue ref
66  // to the move constructor of MatrixT and then move
67  // assignment of unique_ptr is used to move the unique prt
68  // to pA
69 
70  }
71  else {
72  pB = std::make_unique<MatrixT>(exp(M,(n-1)/2));
73  pA = std::make_unique<MatrixT>(trans(prod(*pB,M))); // M^((n+1)/2)
74  }
75  MatrixT R = element_prod(*pA,*pB);
76  typename MatrixT::value_type tr=0;
78  for(typename MatrixT::iterator1 it1 = R.begin1(); it1 != R.end1(); ++it1) {
79  for(typename MatrixT::iterator2 it2 = it1.begin(); it2 != it1.end(); ++it2){
80  tr += *it2;
81  }
82  }
83  return tr;
84  }
85 }
90 template <typename MatrixT>
91 MatrixT exp(const MatrixT& M, size_t n)
92 {
93  return exp_by_squaring(M,n);
94 }
95 
99 template <typename MatrixT>
100 typename MatrixT::value_type det(const MatrixT& M)
101 {
102  if(!isSquare(M)) {
103  throw (std::invalid_argument("Input to det is not a square matrix"));
104  }
105  /* base */
106  using valueT = typename MatrixT::value_type;
107  std::vector<valueT> invariants;
108  std::vector<valueT> traces;
109  size_t size = M.size1();
110  invariants.reserve(size);
111  traces.reserve(size);
112  valueT temp = trace(M);
113  invariants.push_back(temp);
114  traces.push_back(temp);
115  for(size_t k = 2; k <= size ; ++k) { // remember : k is 1 based but index of vector is 0 based
117  valueT outer_multiplier = (k%2==0) ? (-1/double(k)) : (1/double(k));
118  /* first inner term */
119  traces.push_back(trace(M,k));
120  /* second inner term */
121  valueT sum = 0;
122  for(size_t i =1 ;i<k;++i) {
123  sum += (i%2==0) ? invariants[i-1]*traces[k-i-1] : invariants[i-1]*traces[k-i-1]*(valueT(-1));
124  }
125  invariants.push_back((traces[k-1] + sum)*(outer_multiplier));
126 
127  }
128  return (valueT(invariants[size-1]));
129 }
130 
131 
132 
133 
134 
135 #endif // _MATRIXALGORITHMS_H
136 
137 
138 
139 
140 
MatrixT::value_type trace(const MatrixT &M, size_t n)
MatrixT exp_by_squaring(const MatrixT &M, size_t n)
MatrixT exp(const MatrixT &M, size_t n)
Finds the M^n for a square matrix M.
MatrixT::value_type det(const MatrixT &M)
Finds the determinant of a square matrix.
bool isSquare(const MatrixT &M)
Convinience function to check if a matrix is square or not.
Definition: Matrix_impl.ipp:44