pdf版本下載㈠ Fibonacci數
剛接觸Fibonacci數的時候，在網上看到“矩陣法”，看到要先實現一個矩陣乘法，感覺太麻煩了。后來仔細觀察Fibonacci數列，發現有下面的規律：
F(n) = F(k)*F(n+1-k) + F(k-1)*F(n-k) =>
F(2*n) = F(n+1) * F(n) + F(n) * F(n - 1)
F(2*n+1) = F(n+1) * F(n+1) + F(n) * F(n)
根據該公式：要計算F(n)，只需先計算出F(n/2)和F(n/2+1)，于是得出一個數的O(log n)解法。（例如：計算F(13) => 計算F(6)、F(7) => 計算F(3)、F(4) => 計算F(1)、F(2)。）
再后來無意間發現，“矩陣法”根本就不必實現一個矩陣，網上廣為流傳的糟糕的做法，掩蓋了“矩陣法”的優美。
先回顧下Fibonacci數列的矩陣法：
上式中，對系數矩陣A求n次方，有O(log n)解法，因而整個算法是O(log n)。
某些介紹矩陣法的文章，會“偷懶”采用上面的第二種寫法，而不是第一種寫法。偷懶的結果，總是要付出代價的。對上面矩陣法的實現，存在兩個盲點，也正由于這兩個盲點，使“矩陣法”的實現代碼看起來很復雜，失去了簡潔之美。
盲點之一：對系數矩陣A求n次方，可以不采用矩陣乘法來實現。
將F(1) = F(2) = 1， F(0) = 0代入上面的公式1，得到：
上式，對任意 n >=1都成立，也就是說A的任意n次方，只要用兩個變量表示，根本沒必要去實現矩陣乘法。
另外，由 A^n = A^k * A^(n-k)，結合上式，很容易就得到前面提到的公式：
F(n) = F(k)*F(n+1-k) + F(k-1)*F(n-k)，
盲點之二： A的n次方計算方法。
計算一個數m的n次方，
若采用迭代法的話，一般是將m^n，拆分成m、m^2、m^4、m^8…中的幾個的乘積。
若采用遞歸的話，則是將m^n拆分成計算m^(n/2)
//迭代法：
int pow1(int m, unsigned n)
{
int result = 1;
int factor = m;
while (n) {
if (n & 1) { result *= factor; }
factor *= factor;
n /= 2u;
}
return result;
}
//遞歸法
int pow2(int m, unsigned n)
{
if (n == 0) return 1;
int square_root = pow2(m, n / 2);
int result = square_root * square_root;
if (n & 1) result *= m;
return result;
}
對于計算一個整數的n次方，顯然第一種解法效率高，但對計算矩陣的n次方，第二種解法（遞歸法）則更簡單。該遞歸算法也可寫成迭代形式：
int pow3(int m, unsigned n)
{
if (n == 0) return 1;
unsigned flag = n; //小等于n的最大的2的k次冪
for (unsigned value = n; value &= (value - 1); ) flag = value;
int result = m;
while (flag >>= 1) {
result *= result;
if (n & flag) result *= m;
}
return result;
}
（求小等于n的最大的2的k次冪（或求二進制表示中的最高/左位1），有兩種不通用的O(1)方法：一種是使用位掃描匯編指令、另外一種是利用浮點數的二進制表示。）
unsigned extract_leftmost_one(unsigned num)
{
union {
unsigned i;
float f;
} u;
u.f = (float)num;
return u.i >> 23;
}
最后可得到如下代碼：
① 采用一般迭代法計算A^n
static inline void matrix_multiply(uint& b1, uint& b2, uint a1, uint a2)
{
const uint r1 = a1 * b1 + a1 * b2 + a2 * b1;
const uint r2 = a1 * b1 + a2 * b2;
b1 = r1;
b2 = r2;
}
uint fib_matrix(uint num)
{
uint b1 = 0, b2 = 1;
uint a1 = 1, a2 = 0;
for (; num != 0; num >>= 1) {
if (num & 1) matrix_multiply(b1, b2, a1, a2);
matrix_multiply(a1, a2, a1, a2);
}
return b1;
}
② 采用新的迭代法計算A^n
typedef unsigned uint;
uint fibonacci(uint num)
{
if (num == 0) return 0;
uint flag = num; //extract_leftmost_one
for (uint value = num; value &= value - 1; ) flag = value;
uint a1 = 1, a2 = 0;
while (flag >>= 1) {
const uint r1 = a1 * a1 + 2 * a1 * a2;
const uint r2 = a1 * a1 + a2 * a2;
a1 = r1;
a2 = r2;
if (num & flag) {
a1 = r1 + r2;
a2 = r1;
}
}
return a1;
}
上面提到的方法，很容易擴展到三階矩陣，下面是《編程之美》書上的一道擴展題的解法：
（具體分析見下一節）
假設：A(0)=1, A(1)=2, A(2)=2，對n>2都有A(n)=A(n-1)+A(n-2)+A(n-3),
1. 對于任何一個給定的n，如何計算出A(n)？
2. 對于n非常大的情況，如n=2^60的時候，如何計算A(n) mod M （M < 100000)呢？
typedef unsigned uint;
typedef unsigned long long uint64;
uint fib_ex(uint64 num, uint M)
{
assert(M != 0);
const uint g0 = 1, g1 = 2, g2 = 2;
if (num == 0) return g0;
uint64 flag = num;
for (uint64 value = num; value &= value - 1; ) flag = value;
uint64 a1 = 0, a2 = 1, a3 = 0;
while (flag >>= 1) {
const uint64 r1 = 2 * (a1 + a2 + a3) * a1 + a2 * a2;
const uint64 r2 = 2 * (a1 + a2) * a1 + 2 * a2 * a3;
const uint64 r3 = (a1 + 2 * a2) * a1 + a3 * a3;
a1 = r1;
a2 = r2;
a3 = r3;
if (num & flag) {
a1 = r1 + r2;
a2 = r1 + r3;
a3 = r1;
}
a1 %= M;
a2 %= M;
a3 %= M;
}
return (a1 * g2 + a2 * g1 + a3 * g0) % M;
}
㈡ 擴展數列的通解：
下面將前面的結果擴展到任意m階數列：
例子：
① m=2： g(n) = f1 * g(n-1) + f2 * g(n-2)， 初始值為：g0 = g(0)， g1=g(1)
設系數矩陣為A，An的最后一行為(a1 a2)，則
倒數第二行為：(f1*a1 + a2 f2*a1)
即：
系數矩陣A An
f1 f2 f1*a1 + a2 f2*a1
1 0 a1 a2
typedef unsigned uint;
uint fib_matrix2(uint num)
{
if (num == 0) return g0;
uint flag = num;
for (uint value = num; value &= value - 1; ) flag = value;
/*
A A^n
f1 f2 f1*a1 + a2 f2*a1
1 0 a1 a2
*/
uint a1 = 1, a2 = 0; // 0 0 ... 1 0
while (flag >>= 1) {
const uint r1 = f1 * a1 * a1 + 2 * a1 * a2;
const uint r2 = f2 * a1 * a1 + a2 * a2;
a1 = r1;
a2 = r2;
if (num & flag) {
a1 = f1 * r1 + r2;
a2 = f2 * r1;
}
}
return a1 * g1 + a2 * g0;
}
② m=3： g(n) = f1 * g(n-1) + f2 * g(n-2) + f3*g(n-3)，初始值為：g0 = g(0)，g1=g(1), g2=g(2)
設系數矩陣為A，An的最后一行為(a1 a2 a3)，則
倒數第二行為：(f1*a1 + a2 f2*a1 + a3 f3*a1)
倒數第三行為：((f1*f1+f2)*a1 + f1*a2 + a3 (f1*f2+f3)*a1 + f2*a2 f1*f3*a1 + f3*a2)
即：
系數矩陣A An
f1 f2 f3 (f1*f1+f2)*a1 + f1*a2 + a3 (f1*f2+f3)*a1 + f2*a2 f1*f3*a1 + f3*a2
1 0 0 f1*a1 + a2 f2*a1 + a3 f3*a1
0 1 0 a1 a2 a3
typedef unsigned uint;
uint fib_matrix3(uint num)
{
if (num == 0) return g0;
uint flag = num;
for (uint value = num; value &= value - 1; ) flag = value;
/*
A A^n
f1 f2 f3 (f1*f1+f2)*a1 + f1*a2 + a3 (f1*f2+f3)*a1 + f2*a2 f1*f3*a1 + f3*a2
1 0 0 f1*a1 + a2 f2*a1 + a3 f3*a1
0 1 0 a1 a2 a3
*/
uint a1 = 0, a2 = 1, a3 = 0; // 0 0 ... 1 0
while (flag >>= 1) {
const uint r1 = (f1 * f1 + f2) * a1 * a1 + 2 * f1 * a1 * a2 + 2 * a1 * a3 + a2 * a2;
const uint r2 = (f1 * f2 + f3) * a1 * a1 + 2 * f2 * a1 * a2 + 2 * a2 * a3;
const uint r3 = (f1 * f3) * a1 * a1 + 2 * f3 * a1 * a2 + a3 * a3;
a1 = r1;
a2 = r2;
a3 = r3;
if (num & flag) {
a1 = f1 * r1 + r2;
a2 = f2 * r1 + r3;
a3 = f3 * r1;
}
}
return a1 * g2 + a2 * g1 + a3 * g0;
}