[C Programming] 快速Fibonacci數算法-延伸(上)

請參見快速Fibonacci數算法。
延伸問題：
1. 請擴充這個方法，讓它可以算出 k 階 Fibonacci 數。
2. (續上題) 請評估一下上一個習題程式的時間值，
    當k愈來愈大時，計算時間會以什麼方式增加；
    這個方法與上一個問題中習題3的做法相比又如何呢？
3. 在上一題與本題的解答中分別有兩個計算Fibonacci數的程式，
    一個很短，很簡單，但另一個卻很長、很複雜，請用比較大的n，
    並且透過反覆計算1000次的方式，看看哪一個程式比較快。

我們今天先解答第一題，剩下兩小題留待明天。

個人解答：
1.
我們在k階的Fibonacci中，
如同原本的方法一樣，將fn....f(n-k+1)和f(n-1)...f(n-k)的關係寫出來。
我們寫出來的矩陣應該會像：
1 1 1 1 ... 1 1
1 0 0 0 ... 0 0
0 1 0 0 ... 0 0
0 0 1 0 ... 0 0
.....................
0 0 0 0 ... 1 0

所以最後化簡到fk ... f1的時候，
我們需要將這個矩陣的(n-k)次方，
而且加總是把第一行的所有值共k個加起來。
這裡我們的m矩陣就是最開始的矩陣，
xm矩陣則是用來承載中間過度的矩陣，
跟原本2階的xa~xd是一樣的概念；
mm則是乘法完成時的矩陣。
由於全部使用陣列，相當於傳遞指標，
所以就不用擔心reference的問題了XD!
原則上下面在做的就是將原本的四行乘法，
擴充到k階的矩陣乘法而已。

別忘了要初始化m矩陣(unik的二維陣列)，
當然如果一次執行只想求一個f(n,k)值，
可以直接令為{0}後，將1的部分填入就好，
但為了檢驗是否中間都有算對，
我們還是得乖一點每次都確實初始化才行~

Code：

#include <stdio.h>
#include <stdlib.h>


#define MAX_K 20
unsigned long unik[MAX_K+1][MAX_K+1];

/**                         (n-2)
 +--------+    +--------+
 |        |    |        |
 |   mm   | =  |   m    |   (m == unik)
 |        |    |        |
 +--------+    +--------+

*/


void kmatrix_power(unsigned long m[MAX_K+1][MAX_K+1], int n, int k,
       unsigned long mm[MAX_K+1][MAX_K+1]        )
{
 unsigned long xm[MAX_K+1][MAX_K+1];
 int i, j, l;
 if(n==1)
 {
  for(i=1;i<=k;i++)
   for(j=1;j<=k;j++)
    mm[i][j] = m[i][j];
 }
 else if(n & 0x01 == 1) // n = 2k+1, break into m*m^2k
 {
  // compute m^2k, then multiply m
  kmatrix_power(m, n-1, k, xm);
  for(i=1;i<=k;i++)
   for(j=1;j<=k;j++)
   {
    mm[i][j] = m[i][1] * xm[1][j];
    for(l=2;l<=k;l++)
     mm[i][j] += m[i][l] * xm[l][j];
   }
 }
 else // n = 2k, break into (m^k) * (m^k)
 {
  kmatrix_power(m, n>>1, k, xm);
  for(i=1;i<=k;i++)
   for(j=1;j<=k;j++)
   {
    mm[i][j] = xm[i][1] * xm[1][j];
    for(l=2;l<=k;l++)
     mm[i][j] += xm[i][l] * xm[l][j];
   }
 }
}

unsigned long fibonacci(int n, int k)
{
 int i;
 unsigned long mm[MAX_K+1][MAX_K+1];
 if( n <= k )
  return 1;
 else
 {
  kmatrix_power(unik, n-k, k, mm);
  unsigned long tmp;
  for(i = 2, tmp = mm[1][1]; i <= k; i++)
   tmp += mm[1][i];
  return tmp;
 }
}

int main(int argc, char *argv[])
{ 
 int n, i, j, k = 5;
 for(n = 1; n <= 30; n++)
 {
  for(i = 1; i <= k; i++)   // first row is 1 1 1 1 ... 1 
   unik[1][i] = 1UL;
  for(i = 2; i <= k; i++)      // then 1 0 0 0 0 ... 0,
  {        // 0 1 0 0 0 ... 0, ...
   for(j = 1; j <= k; j++){
    unik[i][j] = 0UL;
   }
   unik[i][i-1] = 1UL;
  } 
  printf("f(%d) = %ld\n", n, fibonacci(n, k));
 }
 system("PAUSE");
 return 0;
}