[C Programming] 快速Fibonacci數算法

續Fibonacci數列問題。在非遞迴的Fibonacci程式中，
我們在算fn時最多不超過n-2個加法，
請寫作一個速度更快的程式，或許您可以用乘法。
如果每一個乘法用m單位時間，每一個加法用p單位時間，
於是非遞迴的寫法因為最多有n-2個加法，因此最多用(n-2)p時間。
請問，您的改良程式要用多少時間？
假設只考慮加法與乘法而已。

答：

令上面的矩陣為M，在算M^(n-2)時，
我們可以依照n的奇偶數來像之前算m^n一樣來化簡。
           單位矩陣   r=0
M^r =  (M^k)^2     r=2k
           M*(M^2k)  r=2k+1

因此使用遞迴，我們用2 log n次矩陣相乘就可以得到M^n了。(log底數為2)
2x2的矩陣相乘一次需要8個乘法，4個加法，
所以最多總時間不超過(2 log n)*(8m+4p)。
當n夠大的時候，n > log n，其他的m和p和n無關所以可以視為常數。
所以當n非常大了，此方法所用的時間就小於(n-2)p，因此會比較快。

Code：

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

/**                         (n-2)
 +--------+    +--------+
 | aa  bb |    | a    b |
 |        | =  |        |
 | cc  dd |    | c    d |
 +--------+    +--------+

*/
void matrix_power(unsigned long a, unsigned long b,
      unsigned long c, unsigned long d, int n,
      unsigned long *aa, unsigned long *bb,
      unsigned long *cc, unsigned long *dd    )
{
 unsigned long xa, xb, xc, xd;
 if(n==1)
  *aa = a, *bb = b, *cc = c, *dd = d;
 else if(n & 0x01 == 1) // n = 2k+1, break into m*m^2k
 {
  // compute m^2k, then multiply m
  matrix_power(a, b, c, d, n-1, &xa, &xb, &xc, &xd);
  *aa = a*xa + b*xc;
  *bb = a*xb + b*xd;
  *cc = c*xa + d*xc;
  *dd = c*xb + d*xd;
 }
 else // n = 2k, break into (m^k) * (m^k)
 {
  matrix_power(a, b, c, d, n>>1, &xa, &xb, &xc, &xd);
  *aa = xa*xa + xb*xc;
  *bb = xa*xb + xb*xd;
  *cc = xc*xa + xd*xc;
  *dd = xc*xb + xd*xd;
 }

}

unsigned long fibonacci(int n)
{
 unsigned long a, b, c, d;
 if( n <= 2 )
  return 1;
 else
 {
  matrix_power(1UL, 1UL, 1UL, 0UL, n-2, &a, &b, &c, &d);
  return a + b;
 }

}

int main(int argc, char *argv[])
{ 
 int n;
 for(n = 1; n <= 20; n++)
  printf("f(%d) = %ld\n", n, fibonacci(n));
 system("PAUSE");
 return 0;
}