logN的时间内算fibonacci number

 

三、矩阵法

算Fibonacci数精确值的最快的方法应该就是矩阵法，看过的人都觉得这个方法很好。如果你跟我一样，曾经为记住这个方法中的矩阵而烦恼，那今天就来看看怎么进行推导。其实方法非常简单，想清楚了也就自然而然地记住了。

我们把Fibonacci数列中相邻的两项：F(n)和F(n - 1)写成一个2x1的矩阵，然后对其进行变形，看能得到什么：

[FnFn−1]=[Fn−1+Fn−2Fn−1]=[1×Fn−1+1×Fn−21×Fn−1+0×Fn−2]=[1110]×[Fn−1Fn−2]

是不是非常自然呢？把等式最右边继续算下去，最后得到：

[FnFn−1]=[1110]n−1×[F1F0]=[1110]n−1×[10]

因此要求F(n)，只要对这个二阶方阵求n - 1次方，最后取结果方阵第一行第一列的数字就可以了。

看起来有点儿化简为繁的感觉，但关键点在于，幂运算是可以二分加速的。设有一个方阵a，利用分治法求a的n次方，有：

an={an/2×an/2a(n−1)/2×a(n−1)/2×a,if x iseven, if x isodd

可见复杂度满足T(n) = T(n / 2) +O(1)，根据Master定理可得：T(n) = O(log n)。

在实现的时候，可以用循环代替递归实现这里的二分分治，好处是降低了空间复杂度（用递归的话，空间复杂度为O(log n)）。下面的Python程序直接利用的numpy库中的矩阵乘法（当然这个库也实现了矩阵的幂运算，我把它单独写出来是为了强调这里的分治算法）。

 

斐波那契数列Fibonacci

http://blog.csdn.net/ju136/article/details/8062078

• 最简单的处理方式

当n<=50的时候，是可以利用64位int长整型。long long来进行处理。

[cpp] view plaincopy

1. #include <stdio.h> 
2. #include <stdlib.h> 
3.   
4. long long a[71]; 
5. void init(void) { 
6.     a[0] = 0; 
7.     a[1] = 1; a[2] = 2; 
8.     for (int i = 3; i < 71; ++i) { 
9.         a[i] = a[i-1] + a[i-2]; 
10.     } 
11. } 
12.   
13. int main(void) { 
14.     int n;init(); 
15.     while (scanf("%d", &n) != EOF) { 
16.         printf("%lld\n", a[n]); 
17.     } 
18.     return 0; 
19. } 
20. //注意打表及长整型。 

很简单，就把上面那个题给秒掉了。

 

高效的处理

如果你还对前面讲到的斐波那契数有印象的话，可以轻易地得到。

[ FnFn-1] = [F2 F1] pow(A, n - 2);

A =[1 1

       1 0];

 

值得注意的是，针对于原始的斐波那契数而言，其递推之后是到[F1 F0] 截止，所以其指数是n-1。而在这里是到[F2F1] 截止。所以指数为n-2。

此外，假设T = pow(A, n-2);

与前面不同的是，Fn = F2*T[0][0] + F1*T[1][0];

从而可以写出如下代码：

 

[cpp] view plaincopy

 

21. #include <stdio.h> 
22. #include <stdlib.h> 
23.   
24. typedef struct _node 
25. { 
26.     long long a, b; 
27.     long long c, d; 
28. } node; 
29.   
30.   
31. void _multiple(node *x, node *y) { 
32.     node temp; 
33.     temp.a = x->a*y->a + x->b*y->c; 
34.     temp.b = x->a*y->b + x->b*y->d; 
35.     temp.c = x->c*y->a + x->d*y->c; 
36.     temp.d = x->c*y->b + x->d*y->d; 
37.     *x = temp; 
38. } 
39.   
40. long long fib(int n) 
41. { 
42.     if (0 == n) return 0; 
43.     if (1 == n) return 1; 
44.     if (2 == n) return 2; 
45.   
46.     node odd; odd.a = odd.d = 1; odd.c = odd.b = 0; //单位矩阵 
47.     node temp; temp.a = temp.b = temp.c = 1; temp.d = 0; // A矩阵 
48.   
49.     n -= 2; 
50.     while (n) { 
51.         if (n&1) _multiple(&odd, &temp); 
52.         _multiple(&temp, &temp); 
53.         n >>= 1; 
54.     } 
55.   
56.     return odd.a * 2 + odd.c; 
57. } 
58.   
59. int main() 
60. { 
61.     int n; 
62.     while (scanf("%d", &n) != EOF) { 
63.         printf("%lld\n", fib(n)); 
64.     } 
65.     return 0; 
66. } 

大整数的处理

这里，还是可以利用前面提到的大整数的加法模板来进行运算。

[cpp] view plaincopy

 

67. #include <iostream> 
68. #include <string> 
69. #include <stdio.h> 
70. #include <stdlib.h> 
71. using namespace std; 
72.   
73. string &_del_zeros_before_dot(string &a) 
74. { 
75.     if (a.length() <= 0 || a[0] != '0') return a; 
76.     int i = 0; 
77.     while (i < a.length() && a[i] == '0') ++i; 
78.     a = a.substr(i, a.length() - i); 
79.     return a; 
80. } 
81.   
82. string &_string_add_string(const string &a, const string &b, string &res) 
83. { 
84.     int sum_value = 0, add_bit = 0; 
85.     const int alen = a.length(), blen = b.length(); 
86.     res = "0" + (alen > blen ? a : b); 
87.     for (int i = alen-1, j = blen-1, k = res.length() - 1;  
88.          i >= 0 || j >= 0 || add_bit > 0;  
89.          --i, --j, --k){ 
90.         sum_value = (i>=0 ? a[i]-48: 0) + (j>=0 ? b[j]-48: 0) + add_bit; 
91.         add_bit = sum_value / 10;  
92.         res[k] = sum_value%10 + '0'; 
93.     } 
94.     if (res[0] == '0') res = res.substr(1, res.length() - 1); 
95.     return res; 
96. } 
97.   
98. string fib(int n) { 
99.     if (0 == n) return "0"; 
100.     if (1 == n) return "1"; 
101.     if (2 == n) return "2"; 
102.     string a_2 = "1", b_1 = "2", ret = "0"; 
103.     for (int i = 3; i <= n; ++i) { 
104.         _string_add_string(a_2, b_1, ret); // fn = fn-2 + fn-1; 
105.         a_2 = b_1; 
106.         b_1 = ret; 
107.     } 
108.     return ret; 
109. } 
110.   
111. int main(void) { 
112.     int n; 
113.     while (scanf("%d", &n) != EOF) { 
114.         printf("%s\n", fib(n).c_str()); 
115.     } 
116.     return 0; 
117. } 

扩展-1

首先题意变成：如果一下子可以走1，2，3步呢。也就是多了一种走法。根据惯性思维

Fn = Fn-1 + Fn-2 + Fn-3

且

F1 = 1

F2 = 2

F3 = F1 + F2 + 1 = 4

由于有了递推公式。可以轻易地写出如下代码：

[cpp] view plaincopy

 

118. #include <stdio.h> 
119.   
120. long long step(int n) { 
121.     int i = 0; 
122.     long long a = 1, b = 2, c = 4, ret; 
123.     if (0 == n) return 0; 
124.     if (1 == n) return 1; 
125.     if (2 == n) return 2; 
126.     if (3 == n) return 4; 
127.     for (i = 4; i <= n; ++i) { 
128.         ret = a + b + c; 
129.         a = b; 
130.         b = c; 
131.         c = ret; 
132.     } 
133.     return ret; 
134. } 
135.   
136. int main(void) 
137. { 
138.     int n; 
139.     while (scanf("%d", &n) != EOF) { 
140.         printf("%lld\n", step(n)); 
141.     } 
142.     return 0; 
143. } 

同样也可以写成大数模板形式

[cpp] view plaincopy

 

144. #include <iostream> 
145. #include <string> 
146. #include <stdio.h> 
147. #include <stdlib.h> 
148. using namespace std; 
149.   
150. string &_del_zeros_before_dot(string &a) 
151. { 
152.     if (a.length() <= 0 || a[0] != '0') return a; 
153.     int i = 0; 
154.     while (i < a.length() && a[i] == '0') ++i; 
155.     a = a.substr(i, a.length() - i); 
156.     return a; 
157. } 
158.   
159. string &_string_add_string(const string &a, const string &b, string &res) 
160. { 
161.     int sum_value = 0, add_bit = 0; 
162.     const int alen = a.length(), blen = b.length(); 
163.     res = "0" + (alen > blen ? a : b); 
164.     for (int i = alen-1, j = blen-1, k = res.length() - 1;  
165.          i >= 0 || j >= 0 || add_bit > 0;  
166.          --i, --j, --k){ 
167.         sum_value = (i>=0 ? a[i]-48: 0) + (j>=0 ? b[j]-48: 0) + add_bit; 
168.         add_bit = sum_value / 10;  
169.         res[k] = sum_value%10 + '0'; 
170.     } 
171.     if (res[0] == '0') res = res.substr(1, res.length() - 1); 
172.     return res; 
173. } 
174.   
175. string fib(int n) { 
176.     if (0 == n) return "0"; 
177.     if (1 == n) return "1"; 
178.     if (2 == n) return "2"; 
179.     if (3 == n) return "4"; 
180.     string a = "1", b = "2", c = "4", ret = "0", temp="0"; 
181.     for (int i = 4; i <= n; ++i) { 
182.         _string_add_string(a, b, temp); // temp = fn-2 + fn-3; 
183.         _string_add_string(temp, c, ret); // ret = temp + fn-1; 
184.         a = b; 
185.         b = c; 
186.         c = ret; 
187.     } 
188.     return ret; 
189. } 
190.   
191. int main(void) { 
192.     int n; 
193.     while (scanf("%d", &n) != EOF) { 
194.         printf("%s\n", fib(n).c_str()); 
195.     } 
196.     return 0; 
197. } 

那么我们考虑，如何写成O(lgN)的效率呢。

同样，我们也可以写出矩阵的表达式。

[Fn Fn-1 Fn-2] = [Fn-1Fn-2 Fn3] * A;

可以得出A的值为

1  1   0

1  0   1

1  0   0

那么，在这种情况下，如果层次递推，我们可以得到

[Fn Fn-1 Fn-2] = [F3 F2F1] * pow(A, n-3);

如果设T = pow(A, n-3);

那么Fn = F3*T[0,0] + F2*T[1,0] + F1*T[2,0];

 

因此，也可以写出如下的代码。

[cpp] view plaincopy

 

198. #include <stdio.h> 
199. #include <stdlib.h> 
200. #include <string.h> 
201. #define MX 3 
202.   
203. typedef struct _node 
204. { 
205.     long long a[MX][MX]; 
206. } node; 
207.   
208.   
209. void _multiple(node *x, node *y) { 
210.     node temp; 
211.     long long (*a)[MX] = x->a; 
212.     long long (*b)[MX] = y->a; 
213.     long long (*c)[MX] = temp.a; 
214.   
215.     for (int row = 0; row < MX; ++row) {      
216.         for (int col = 0; col < MX; ++col) { 
217.             long long s = 0; 
218.             for (int i = 0; i < MX; ++i) { 
219.                 s += a[row][i] * b[i][col]; 
220.             } 
221.             c[row][col] = s; 
222.         } 
223.     } 
224.     *x = temp; 
225. } 
226.   
227. long long fib(int n) 
228. { 
229.     if (0 == n) return 0; 
230.     if (1 == n) return 1; 
231.     if (2 == n) return 2; 
232.     if (3 == n) return 4; 
233.   
234.     node odd; 
235.     //设置odd为单位矩阵 
236.      memset(odd.a, 0, sizeof(odd.a)); 
237.     for (int i = 0; i < MX; ++i) { 
238.         odd.a[i][i] = 1; 
239.     } 
240.     node temp; // A矩阵 
241.     memset(temp.a, 0, sizeof(temp.a)); 
242.     for (int i = 0; i < MX; ++i) { 
243.         temp.a[i][0] = 1; 
244.         if (i + 1 < MX) temp.a[i][i+1] = 1; 
245.     } 
246.   
247.     n -= 3; 
248.     while (n) { 
249.         if (n&1) _multiple(&odd, &temp); 
250.         _multiple(&temp, &temp); 
251.         n >>= 1; 
252.     } 
253.     return 4*odd.a[0][0] + 2*odd.a[1][0] + odd.a[2][0]; 
254. } 
255.   
256. int main() 
257. { 
258.     int n; 
259.     while (scanf("%d", &n) != EOF) { 
260.         printf("%lld\n", fib(n)); 
261.     } 
262.     return 0; 
263. } 

扩展-2

那么接下的问题是

假设A上台阶，一次可以跨1层，2层，3层..或m层，问A上n层台阶，有多少种走法？

其中，m和n都是正整数，并且 m<= n, m <= 10, n <= 50

由于m, n并不清楚。

 

不过对于1<= x <= m的走法，应该是很清楚的。也就是Fx = Fx-1 + Fx-2 + ...... + F1 + 1

而对于大于m的而言，则是Fy= Fy-1 + Fy-2 + Fy-3 + .... + Fy-m;

 

公式都出来了，就不用去想了。

[cpp] view plaincopy

 

264. #include <stdio.h> 
265. #include <iostream> 
266. #include <list> 
267.   
268. using namespace std; 
269.   
270. long long fib(unsigned int m, int n) { 
271.     list<long long> l; 
272.     long long a = 0, s = 0, ret; 
273.     for (int i = 1; i <= m && i <= n; ++i) { 
274.         a = s + 1; 
275.         l.push_back(a); 
276.         s += a; 
277.     } 
278.     if (n <= m) return l.back(); 
279.   
280.     for (int i = m + 1; i <= n; ++i) { 
281.         s = 0; 
282.         for (list<long long>::iterator iter = l.begin();  
283.             iter != l.end(); ++iter) { 
284.             s += *iter; 
285.         } 
286.         ret = s; 
287.   
288.         l.pop_front(); 
289.         l.push_back(ret); 
290.     } 
291.     return ret; 
292. } 
293.   
294. int main(void) { 
295.     int m, n; 
296.     while (scanf("%d%d", &m, &n) != EOF){ 
297.         printf("%lld\n", fib(m, n)); 
298.     } 
299.     return 0; 
300. } 

如果觉得上面的方法不好理解。没关系。直接采用打表法。

[cpp] view plaincopy

 

301. #include <stdio.h> 
302. int m, n; 
303. long long a[51]; 
304. long long fib(int m, int n) { 
305.     long long s = 0; 
306.     a[0] = 0; 
307.     for (int i = 1; i <=m && i <= n; ++i) { 
308.         a[i] = s + 1; 
309.         s += a[i]; 
310.     } 
311.     if (n <= m) return a[n]; 
312.     for (int i = m + 1; i <= n; ++i) { 
313.         s = 0; 
314.         for (int j = 1; j <= m; ++j) { 
315.             s += a[i-j]; 
316.         } 
317.         a[i] = s; 
318.     } 
319.     return a[n]; 
320. } 
321. int main(void) 
322. { 
323.     while (scanf("%d%d", &m, &n) != EOF &&0 <= m && m <=n && n <= 50) { 
324.         printf("%lld\n", fib(m, n)); 
325.     } 
326.     return 0; 
327. } 

这里也顺便给出大整数求解的方案，也就是当n超过50的时候的处理方式。

[cpp] view plaincopy

 

1. #include <iostream> 
2. #include <string> 
3. #include <stdio.h> 
4. #include <stdlib.h> 
5. #include <list> 
6. using namespace std; 
7.   
8. string &_del_zeros_before_dot(string &a) 
9. { 
10.     if (a.length() <= 0 || a[0] != '0') return a; 
11.     int i = 0; 
12.     while (i < a.length() && a[i] == '0') ++i; 
13.     a = a.substr(i, a.length() - i); 
14.     return a; 
15. } 
16.   
17. string &_string_add_string(const string &a, const string &b, string &res) 
18. { 
19.     int sum_value = 0, add_bit = 0; 
20.     const int alen = a.length(), blen = b.length(); 
21.     res = "0" + (alen > blen ? a : b); 
22.     for (int i = alen-1, j = blen-1, k = res.length() - 1;  
23.          i >= 0 || j >= 0 || add_bit > 0;  
24.          --i, --j, --k){ 
25.         sum_value = (i>=0 ? a[i]-48: 0) + (j>=0 ? b[j]-48: 0) + add_bit; 
26.         add_bit = sum_value / 10;  
27.         res[k] = sum_value%10 + '0'; 
28.     } 
29.     if (res[0] == '0') res = res.substr(1, res.length() - 1); 
30.     return res; 
31. } 
32.   
33. string fib(unsigned int m, int n) { 
34.     list<string> l; 
35.     string a = "0", s = "0", ret, temp; 
36.     for (int i = 1; i <= m && i <= n; ++i) { 
37.         _string_add_string(s, "1", a); 
38.         l.push_back(a); 
39.         _string_add_string(s, a, temp); 
40.         s = temp; 
41.     } 
42.     if (n <= m) return l.back(); 
43.   
44.     for (int i = m + 1; i <= n; ++i) { 
45.         s = "0"; 
46.         for (list<string>::iterator iter = l.begin();  
47.             iter != l.end(); ++iter) { 
48.             _string_add_string(s, *iter, temp); 
49.             s = temp; 
50.         } 
51.         ret = s; 
52.   
53.         l.pop_front(); 
54.         l.push_back(ret); 
55.     } 
56.     return ret; 
57. } 
58.   
59. int main(void) { 
60.     int m, n; 
61.     while (scanf("%d%d", &m, &n) != EOF) { 
62.         printf("%s\n", fib(m, n).c_str()); 
63.     } 
64.     return 0; 
65. } 

可能你还希望找出一个针对于存在m的情况下的高效率的解法。那好吧。这里给个示例。原理与前面O(lgN)的原理几乎一样。

只不过你需要注意一下A矩阵的表达形式。其他的没什么不一样的。还有就是最后有个相乘的情况。

[cpp] view plaincopy

 

66. #include <stdio.h> 
67. #include <stdlib.h> 
68. #include <string.h> 
69.   
70. typedef struct _node { 
71.     long long **a; 
72. }node; 
73. long long **_init(int m) { 
74.     long long **t = NULL; 
75.     t = (long long **)malloc(sizeof(long long*)*m); 
76.     for (int i = 0; i < m; ++i){ 
77.         t[i] = (long long *)malloc(sizeof(long long) * m); 
78.         for (int j = 0; j < m; ++j) { 
79.             t[i][j] = 0; 
80.         } 
81.     } 
82.     return t; 
83. } 
84.   
85. void _destroy(long long **t, int m) { 
86.     for (int i = 0; i < m; ++i) { 
87.         free(t[i]); 
88.     } 
89.     free(t); 
90. } 
91.   
92. void _multiple(node *x, node *y, int m) { 
93.     node temp; temp.a = _init(m); 
94.     long long **a = x->a; 
95.     long long **b = y->a; 
96.     long long **c = temp.a; 
97.   
98.     for (int row = 0; row < m; ++row) {       
99.         for (int col = 0; col < m; ++col) { 
100.             long long s = 0; 
101.             for (int i = 0; i < m; ++i) { 
102.                 s += a[row][i] * b[i][col]; 
103.             } 
104.             c[row][col] = s; 
105.         } 
106.     } 
107.     _destroy(x->a, m); 
108.     x->a = temp.a; 
109. } 
110.   
111. long long fib(int m, int n) 
112. { 
113.     long long ret; 
114.     long long *front = (long long *) malloc(sizeof(long long)*(m + 1)); 
115.   
116.     long long s = 0;  
117.     front[0] = 0; 
118.     for (int i = 1; i <= m && i <= n; ++i) { 
119.         front[i] = s + 1; 
120.         s += front[i]; 
121.     } 
122.   
123.     if (n <= m) { 
124.         ret = front[n]; 
125.         free(front); 
126.         return ret; 
127.     } 
128.   
129.     node odd; odd.a = _init(m); 
130.     //设置odd为单位矩阵 
131.     for (int i = 0; i < m; ++i) { 
132.         odd.a[i][i] = 1; 
133.     } 
134.     //A矩阵. 
135.     node temp; temp.a = _init(m);  
136.     for (int i = 0; i < m; ++i) { 
137.         temp.a[i][0] = 1; 
138.         if (i + 1 < m) temp.a[i][i+1] = 1; 
139.     } 
140.   
141.     n -= m; 
142.     while (n) { 
143.         if (n&1) _multiple(&odd, &temp, m); 
144.         _multiple(&temp, &temp, m); 
145.         n >>= 1; 
146.     } 
147.   
148.     ret = 0; 
149.     for (int i = 1; i <= m; ++i) { 
150.         ret += front[i] * odd.a[m-i][0]; 
151.     } 
152.   
153.     _destroy(odd.a, m); 
154.     _destroy(temp.a, m); 
155.     return ret; 
156. } 
157.   
158. int main() 
159. { 
160.     int m, n; 
161.     while (scanf("%d%d", &m, &n) != EOF) { 
162.         printf("%lld\n", fib(m, n)); 
163.     } 
164.     return 0; 
165. } 

扩展-3

一只青蛙一次可以跳上1级台阶，也可以跳上2级……它也可以跳上n级。求该青蛙跳上一个n级的台阶总共有多少种跳法。

测试链接:http://ac.jobdu.com/problem.php?cid=1039&pid=5
 

这个只是第二个扩展的特殊情况。也就是当m==n的时候的情况。

[cpp] view plaincopy

 

166. #include <stdio.h> 
167. #include <stdlib.h> 
168. #include <string.h> 
169.   
170. long long fib(int m) 
171. { 
172.     long long ret; 
173.     long long *front = (long long *) malloc(sizeof(long long)*(m + 1)); 
174.   
175.     long long s = 0;  
176.     front[0] = 0; 
177.     for (int i = 1; i <= m; ++i) { 
178.         front[i] = s + 1; 
179.         s += front[i]; 
180.     } 
181.     ret = front[m]; 
182.     free(front); 
183.     return ret; 
184. } 
185.   
186. int main() 
187. { 
188.     int m, n; 
189.     while (scanf("%d", &n) != EOF) { 
190.         printf("%lld\n", fib(n)); 
191.     } 
192.     return 0; 
193. } 

这个就没有更高效的办法了。除了再加入大整数之外。~~大整数就不用再加了。都差不多的。

 

 当然，如果再深入地看一下，会发现F(n) = 2^(n-1);

那这个代码就更好写了。

[cpp] view plaincopy

 

194. #include <stdio.h> 
195.   
196. long long pow(long long x, long long n) { 
197.     long odd = 1; 
198.     while (n) { 
199.         if (n&1) odd *= x; 
200.         x *= x; 
201.         n >>= 1; 
202.     } 
203.     return odd; 
204. } 
205.   
206. int main(void) { 
207.     int n; 
208.     while (scanf("%d", &n) != EOF) { 
209.         printf("%lld\n", pow(2, n-1));  
210.     } 
211.     return 0; 
212. }