Relação de Recorrência com exponenciação de matriz

Relações de recorrência

Relações de recorrência aparecem recorretemente :) em computação. Podemos encontrar o n-ésimo termo de uma relação de recorrência em O(n) utilizando programação dinâmica:
*Construindo uma solução bottom-up com resultados pré-calculados.
*Utilizando memorização.

Podemos também resolver relações de recorrência através de  álgebra linear. Uma relação recorrência pode ser vista como uma sequência de vetores definida por uma equação matricial x(i+1) = M x(i) para uma matriz constante M.

Considere a relação de recorrência que define a sequência de Fibonacci.

F(i+1) = F(i) + F(i-1)

F(i+1) =[a b][F(i)  ]
F(i)       [c d][F(i-1)]

aF(i) + bF(i-1) = F(i+1)
a=b=1
F(i) + F(i-1) = F(i+1)

cF(i) + dF(i-1) = F(i)
c = 1 , d=0
F(i)=F(i)

f(0) = 0
f(1) = 1
[1 1][f(1)=1]=[f(2)=1]
[1 0][f(0)=0]  [f(1)=1]

[1 1][f(2)=1]=[f(3)=2]
[1 0][f(1)=1]  [f(2)=1]

[1 1][f(3)=2]=[f(3)=3]
[1 0][f(2)=1]  [f(2)=2]

([1 1] )^n [F(1)] = [F(n+1)]
([1 0] )     [F(0)]    [F(n)    ]

Podemos calcular o n-ésimo termo da sequência de fibonacci utilizando exponenciação de matriz com complexidade $O(log_2 n)$ da seguinte maneira:

Considere a seguinte propriedade para n >= 1:

M^n = [F(n+1) F(n)]
           [F(n)   F(n-1)]

n=1
[1 1]
[1 0]
é válido
Passo de indução:
M^(n+1) = M*M^n
[1 1][F(n+1) F(n)  ]=[F(n+2) F(n+1)]
[1 0][F(n)   F(n-1)]   [F(n+1) F(n)  ]

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <map>
#include <vector>
#include <iostream>
using namespace std;
template <typename T>
T** multiplica(T  **a, T  **b, int n){
 int i,j,k;
 T** c; 
 c = (T **)malloc(n*sizeof(T *));
 for(i=0;i<2;i++)
  c[i] = (T *)malloc(n*sizeof(T));
 for(i=0;i<n;i++){
  for(j=0;j<n;j++){
    c[i][j] = T();     
    for(k=0;k<n;k++){
 c[i][j] += a[i][k]*b[k][j];
     }
   }
 }
 return c; 
}

template <typename T>
T** exp(T **a,int n, int m){
 int i,j;
 T ** c;
 T ** b;
 c = (T **)malloc(n*sizeof(T *));
 for(i=0;i<2;i++)
  c[i] = (T *)malloc(n*sizeof(T));
 b = (T **)malloc(n*sizeof(T *));
 for(i=0;i<2;i++)
  b[i] = (T *)malloc(n*sizeof(T));
 if(m==0){
   for(i=0;i<n;i++)
     for(j=0;j<n;j++)
       if(i==j) b[i][j] = (T)1;
       else b[i][j] = (T)0;
 }else if(m==1){
   for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      b[i][j] = a[i][j];
  }else if(m%2==0){
    c = exp(a,n,(int)m/2);
    b = multiplica(c,c,n);
  }else if(m%2==1){
    c = exp(a,n,m-1);
    b = multiplica(c,a,n);
  }
  free(c);
  return b;
}

template <typename T>
void show(T** M, int n){
 int i,j;
 
 for(i=0;i<n;i++){
  for(j=0;j<n;j++){
   cout << M[i][j] << " ";
  }
  cout << "\n";
 }
}

int main(){
 int i;
 int ** M;
 
 M = (int **)malloc(2*sizeof(int *));
 for(i=0;i<2;i++)
  M[i] = (int *)malloc(2*sizeof(int));
 
 M[0][0] = 1;
 M[0][1] = 1;
 M[1][0] = 1;
 M[1][1] = 0;
  
 show(M,2);
 
 int **Mn; 
 int n;
 
 n=1;
 Mn = exp(M,2,n);
 printf("n=1\n");
 printf("F[%d]=%d\n",n+1,Mn[0][0]);
 printf("F[%d]=%d\n",n,Mn[0][1]);
 printf("F[%d]=%d\n\n",n-1,Mn[1][1]);
 
 n = 2;
 Mn = exp(M,2,n);
 printf("n=2\n");
 printf("F[%d]=%d\n",n+1,Mn[0][0]);
 printf("F[%d]=%d\n",n,Mn[0][1]);
 printf("F[%d]=%d\n\n",n-1,Mn[1][1]);
 
 n=10;
 Mn = exp(M,2,n);
 printf("n=10\n");
 printf("F[%d]=%d\n",n+1,Mn[0][0]);
 printf("F[%d]=%d\n",n,Mn[0][1]);
 printf("F[%d]=%d\n\n",n-1,Mn[1][1]);
}

1 1
1 0

n=1
F[2]=1
F[1]=1
F[0]=0

n=2
F[3]=2
F[2]=1
F[1]=1

n=10
F[11]=89
F[10]=55
F[9]=34


Referências:
http://en.wikipedia.org/wiki/Recurrence_relation
http://comeoncodeon.wordpress.com/2011/05/08/recurrence-relation-and-matrix-exponentiation/