Matrix Expo.cpp

/*
Uses powers of two to exponentiate numbers and matrices. Calculates
n^k in O(log(k)) time when n is a number. If A is an n x n matrix,
calculates A^k in O(n^3*log(k)) time.
*/
#include <bits/stdc++.h>
using namespace std;
typedef double T;
typedef vector<T> VT;
typedef vector<VT> VVT;
T power(T x, int k) {
  T ret = 1;
  while(k) {
    if(k & 1) ret *= x;
    k >>= 1; x *= x;
  }
  return ret;
}
VVT multiply(VVT& A, VVT& B) {
  int n = A.size(), m = A[0].size(), k = B[0].size();
  VVT C(n, VT(k, 0));
  
  for(int i = 0; i < n; i++)
    for(int j = 0; j < k; j++)
      for(int l = 0; l < m; l++)
        C[i][j] += A[i][l] * B[l][j];

  return C;
}
VVT power(VVT& A, int k) {
  int n = A.size();
  VVT ret(n, VT(n)), B = A;
  for(int i = 0; i < n; i++) ret[i][i]=1;
  while(k) {
    if(k & 1) ret = multiply(ret, B);
    k >>= 1; B = multiply(B, B);
  }
  return ret;
}
int main() {
  /* Expected Output:
     2.37^48 = 9.72569e+17
     376 264 285 220 265 
     550 376 529 285 484 
     484 265 376 264 285 
     285 220 265 156 264 
     529 285 484 265 376 */
  double n = 2.37; int k = 48;
  cout << n << "^" << k << " = " << power(n, k) << endl;
  double At[5][5] = {
    { 0, 0, 1, 0, 0 },
    { 1, 0, 0, 1, 0 },
    { 0, 0, 0, 0, 1 },
    { 1, 0, 0, 0, 0 },
    { 0, 1, 0, 0, 0 } };
    
  vector <vector <double> > A(5, vector <double>(5));    
  for(int i = 0; i < 5; i++)
    for(int j = 0; j < 5; j++)
      A[i][j] = At[i][j];
    
  vector <vector <double> > Ap = power(A, k);
  
  cout << endl;
  for(int i = 0; i < 5; i++) {
    for(int j = 0; j < 5; j++)
      cout << Ap[i][j] << " ";
    cout << endl;
  }
}