问题：

斐波那契数列由如下递推关系式定义：F(0) = 0，F(1)=1，F(n)=F(n-1)+F(n-2) if n>1。

解法：

斐波那契数列是二阶递推数列，所以存在一个2*2的矩阵A，使得：

(Fn, Fn-1) = (Fn-1, Fn-2)*A

求得A=(1   1)

            (1   0)

那么求数列的第n项就是等于求矩阵A的第n-1次幂，计算的速度非常快，时间复杂度为O(logn)。

首先我们用long long 型表示数列中的元素，它只能表示20位的整数，能表示的范围太小，最多第92个元素。

[cpp] 
view plain
copy

1. #include <iostream>  
2. #include <cassert>  
3. using namespace std;  
4.   
5. const int MAXLENGTH = 10;  
6.   
7. struct Matrix   
8. {  
9.     unsigned side;  
10.     long long dat[MAXLENGTH*MAXLENGTH];  
11. };  
12.   
13. // 方阵的乘法   
14. void MatrixMult(const Matrix a, const Matrix b, Matrix &m)  
15. {  
16.     assert(a.side == b.side);  
17.     m.side = a.side;  
18.     for (int i=0; i<m.side; ++i)  
19.         for (int j=0; j<m.side; ++j)  
20.         {  
21.             m.dat[i*m.side+j] = 0;  
22.             for (int k=0; k<m.side; ++k)  
23.                 m.dat[i*m.side+j] += a.dat[i*a.side+k]*b.dat[k*b.side+j];  
24.         }  
25. }  
26.   
27. long long Fibonaci(unsigned n)  
28. {  
29.     if (n==0) return 0;  
30.     –n;    // 计算矩阵prod的n-1次幂  
31.     Matrix res; res.side = 2;  
32.     res.dat[0] = 1; res.dat[1] = 0;   
33.     res.dat[2] = 0; res.dat[3] = 1;  
34.     Matrix prod;        prod.side = 2;  
35.     prod.dat[0] = 1;  prod.dat[1] = 1;  
36.     prod.dat[2] = 1;  prod.dat[3] = 0;  
37.     // 只需要O(logn)的复杂度就能算出x的n次幂    
38.     while (n)  
39.     {  
40.         // 如果n的最低二进制位为1，则乘上对应的幂次prod  
41.         if (n&1) MatrixMult(res, prod, res);  
42.         MatrixMult(prod, prod, prod);  
43.         n >>= 1;  
44.     }  
45.     return res.dat[0];  
46. }  
47.   
48. int main()  
49. {  
50.     long long res = Fibonaci(93);  
51.     cout << res << endl;  
52. }  

如果想输出更高项的值，就要用大整数来表示矩阵的元素。这里我们可以表示1000位的整数，表达能力是上面算法的50倍，能计算出第4000多个元素，而且计算速度非常快。

[cpp] 
view plain
copy

1. #include <iostream>  
2. #include <cstring>  
3. #include <cassert>  
4. using namespace std;  
5.   
6. // 大整数类型  
7. const int MAXLENGTH = 1000;  
8. struct HP{int len, s[MAXLENGTH];};  
9.   
10. // 字符串转大整数  
11. void Str2HP(const char* s, HP& x)  
12. {  
13.     x.len = strlen(s);  
14.     for (int i=0; i<x.len; ++i)  
15.     {  
16.         assert(s[i]>=‘0’ && s[i]<=‘9’);  
17.         x.s[x.len-1-i] = s[i]-‘0’;  
18.     }  
19.     if (x.len == 0)  
20.     {  
21.         x.len = 1;  
22.         x.s[0] = 0;  
23.     }  
24. }  
25.   
26. void PrintHP(const HP& x)  
27. {  
28.     for (int i=x.len-1; i>=0; –i)  
29.         cout << x.s[i];  
30. }  
31.   
32. // 大整数的加法  
33. void PlusHP(const HP x, const HP y, HP &z)  
34. {  
35.     int i;  z.s[0] = 0;  
36.     // 大整数x,y的加法操作和输出大整数z的进位操作  
37.     for (i=0; i<x.len || i<y.len; ++i)  
38.     {  
39.         if (i<x.len) z.s[i] += x.s[i];  
40.         if (i<y.len) z.s[i] += y.s[i];  
41.         z.s[i+1] = z.s[i]/10;  
42.         z.s[i] %= 10;  
43.     }  
44.     // 第i位不会再进位，这里可省去  
45.     while (z.s[i])  
46.     {  
47.         z.s[i+1] = z.s[i]/10;  
48.         z.s[i] %= 10;  
49.         ++i;  
50.     }  
51.     z.len = i;  
52. }  
53.   
54. // 大整数的减法   
55. void SubtractHP(const HP x, const HP y, HP &z)  
56. {  
57.     int i, j;  
58.     // j表示是否要对高位进行借位，1表示借1位，0表示不借位  
59.     for (i=0, j=0; i<x.len; ++i)  
60.     {  
61.         z.s[i] = x.s[i] – j;  
62.         if (i<y.len) z.s[i] -= y.s[i];  
63.         if (z.s[i]<0)   
64.         {  
65.             // 向高位借位，该位补10  
66.             j=1;   
67.             z.s[i] += 10;  
68.         }  
69.         else  
70.             j=0;  
71.     }  
72.     do –i;  
73.     while (i>0 && !z.s[i]);  
74.     z.len = i+1;  
75. }  
76.   
77. // 大整数的比较  
78. int CompareHP(const HP &x, const HP &y)  
79. {  
80.     if (x.len > y.len) return 1;  
81.     if (x.len < y.len) return -1;  
82.     int i = x.len-1;  
83.     while (i>0 && x.s[i]==y.s[i]) –i;  
84.     return x.s[i]-y.s[i];  
85. }  
86.   
87. // 大整数的乘法  
88. void MultiHP(const HP x, const HP y, HP &z)  
89. {  
90.     int i, j;  
91.     // 对乘法结果大整数z初值化为0，以方便之后的+=运算  
92.     z.len = x.len + y.len;  
93.     for (i=0; i<z.len; ++i) z.s[i] = 0;  
94.     for (i=0; i<x.len; ++i)  
95.         for (j=0; j<y.len; ++j)  
96.             z.s[i+j] += x.s[i]*y.s[j];  
97.     // 大整数z进位  
98.     for (i=0; i<z.len-1; ++i) {z.s[i+1] += z.s[i]/10; z.s[i] %= 10;}  
99.     // 最高位继续进位，这步不会执行可省去  
100.     while (z.s[i]) {z.s[i+1] = z.s[i]/10; z.s[i] %= 10; i++;}  
101.     // 直到最高位不为0 ,以确定大整数的长度  
102.     while (i>0 && !z.s[i]) –i;  
103.     z.len = i+1;  
104. }  
105.   
106.   
107. void DivideHP(const HP x, const HP y, HP &u, HP &v)  
108. {  
109.     int i, j;  
110.     v.len = 1; v.s[0] = 0;  
111.     // u表示x的前i位除y的商  
112.     // v表示x的前i位除y的余数  
113.     for (i=x.len-1; i>=0; –i)  
114.     {  
115.         // 余数v先向左移一位，再加上x.s[i]  
116.         if (!(v.len==1 && v.s[0]==0))  
117.         {  
118.             for (j=v.len-1; j>=0; –j)  
119.                 v.s[j+1] = v.s[j];  
120.             ++v.len;  
121.         }  
122.         v.s[0] = x.s[i];  
123.         // 每次循环都计算出商的第i位u.s[i]  
124.         u.s[i] = 0;  
125.         while ((j=CompareHP(v, y))>=0)  
126.         {  
127.             // 余数v大于等于除数y时，就进行减操作  
128.             SubtractHP(v, y, v);  
129.             ++u.s[i];  
130.             if (j==0) break;  
131.         }  
132.     }  
133.     i = x.len – 1;  
134.     while (i>0 && !u.s[i]) –i;  
135.     u.len = i+1;  
136. }  
137.   
138. // 大整数置0  
139. inline void ZeroHP(HP& x)  
140. {  
141.     x.len = 1; x.s[0] = 0;  
142. }  
143.   
144. // 大整数置1  
145. inline void OneHP(HP& x)  
146. {  
147.     x.len = 1; x.s[0] = 1;  
148. }  
149.   
150. // 二维数组（方阵）  
151. const int MAXSIDE = 2;  
152. struct Matrix{int side; HP a[MAXSIDE*MAXSIDE];};  
153.   
154. // 方阵的乘法  
155. void MultiMatrix(const Matrix x, const Matrix y, Matrix &z)  
156. {  
157.     assert(x.side == y.side);  
158.     z.side = x.side;  
159.     HP tmp;  
160.     for (int i=0; i<z.side; ++i)  
161.         for (int j=0; j<z.side; ++j)  
162.         {  
163.             ZeroHP(z.a[i*z.side+j]);  
164.             for (int k=0; k<z.side; ++k)  
165.             {  
166.                 MultiHP(x.a[i*x.side+k], y.a[k*y.side+j], tmp);  
167.                 PlusHP(z.a[i*z.side+j], tmp, z.a[i*z.side+j]);  
168.             }  
169.         }  
170. }  
171.   
172. const HP& Fibonaci(int n)  
173. {  
174.     HP zero;  
175.     ZeroHP(zero);  
176.     if (n==0) return zero;  
177.     –n;  
178.     Matrix tmp; tmp.side = 2;  
179.     OneHP(tmp.a[0]);    OneHP(tmp.a[1]);  
180.     OneHP(tmp.a[2]);    ZeroHP(tmp.a[3]);  
181.     Matrix res; res.side = 2;  
182.     OneHP(res.a[0]);    ZeroHP(res.a[1]);  
183.     ZeroHP(res.a[2]);   OneHP(res.a[3]);  
184.     while (n)  
185.     {  
186.         if (n&1) MultiMatrix(res, tmp, res);  
187.         MultiMatrix(tmp, tmp, tmp);  
188.         n >>= 1;  
189.     }  
190.     return res.a[0];  
191. }  
192.   
193. int main()  
194. {  
195.     const HP& res = Fibonaci(4000);  
196.     PrintHP(res);  
197.     cout << endl;  
198. }  

转载自：http://blog.csdn.net/linyunzju/article/details/7706896