Matrix power in base R using exponentiation by squaring

matrix_pow.R

matrix_pow <- function(x, n) {
  if (abs(n - round(n)) < .Machine$double.eps^0.5) {
    stop(sprintf("Power %.1f is not an integer", n))
  }
  if (n < 0) {
    return(matrix_exp(solve(x), -n))
  } else if (n == 0) {
    return(1) 
  } else if (n == 1) {
    return(x) 
  } else if (n %% 2 == 0) {
    return(matrix_exp(x %*% x, n / 2))
  } else {
    return(x %*% matrix_exp(x %*% x, (n - 1) / 2))
  }
}