1 #ifndef STAN_MATH_PRIM_FUN_MATRIX_POWER_HPP
2 #define STAN_MATH_PRIM_FUN_MATRIX_POWER_HPP
3
4 #include <stan/math/prim/err.hpp>
5 #include <stan/math/prim/fun/Eigen.hpp>
6
7 namespace stan {
8 namespace math {
9
10 /**
11 * Returns the nth power of the specific matrix. M^n = M * M * ... * M.
12 *
13 * @tparam T type of the matrix
14 *
15 * @param[in] M a square matrix
16 * @param[in] n exponent
17 * @return nth power of M
18 * @throw std::domain_error if the matrix contains NaNs or infinities.
19 * @throw std::invalid_argument if the exponent is negative or the matrix is not
20 * square.
21 */
22 template <typename EigMat, require_eigen_t<EigMat>* = nullptr,
23 require_not_vt_var<EigMat>* = nullptr>
24 inline Eigen::Matrix<value_type_t<EigMat>, EigMat::RowsAtCompileTime,
25 EigMat::ColsAtCompileTime>
matrix_power(const EigMat & M,const int n)26 matrix_power(const EigMat& M, const int n) {
27 using T = value_type_t<EigMat>;
28 constexpr int R = EigMat::RowsAtCompileTime;
29 constexpr int C = EigMat::ColsAtCompileTime;
30
31 check_square("matrix_power", "M", M);
32 check_nonnegative("matrix_power", "n", n);
33 Eigen::Matrix<T, R, C> MM = M;
34 check_finite("matrix_power", "M", MM);
35 if (n == 0)
36 return Eigen::Matrix<T, R, C>::Identity(M.rows(), M.cols());
37 Eigen::Matrix<T, R, C> result = MM;
38 for (int nn = n - 1; nn > 0; nn /= 2) {
39 if (nn % 2 == 1) {
40 result = result * MM;
41 --nn;
42 }
43 MM = MM * MM;
44 }
45 return result;
46 }
47
48 template <typename EigMat, require_eigen_t<EigMat>* = nullptr>
49 inline Eigen::Matrix<value_type_t<EigMat>, EigMat::RowsAtCompileTime,
50 EigMat::ColsAtCompileTime>
operator ^(const EigMat & M,const int n)51 operator^(const EigMat& M, const int n) {
52 return matrix_power(M, n);
53 }
54
55 } // namespace math
56 } // namespace stan
57
58 #endif
59