盒子
盒子
文章目录
  1. 从 Fibonacci 到 矩阵快速幂 ;)
  2. 1000000007
  3. 总结
  4. 参考资料
  5. 经历

矩阵快速幂

ACM校赛,get到新技能–关键词:矩阵快速幂、1000000007。╮(╯_╰)╭,为了方便理解,就使用 Python 代码了,最后附上 C 代码。

吐槽时间:
为啥会写这篇文章? →_→ ACMer 的 blog 不堪入目 -代码凌乱不说,人话也不多说几句- 找了好多篇,不知所以然。而且,网上的实现代码大都是 C++ 的,但我看不懂唉。更何况,C++ 创建一个矩阵类,重载乘法运算,方便多了。用 C 呢,指针乱飞……

从 Fibonacci 到 矩阵快速幂 ;)

Fibonacci 数列的定义:
$$ \begin{equation} \rm A(n) = \left\{\begin{array}{ll} A(n-1) + A(n-2) & x > 1 \\ 1 & x = 0, 1 \end{array} \right. \end{equation} $$
那么,怎么求数列的第 n 项呢?首先想到的是递归:

1
2
3
4
5
def fib_recursion(n):
if (n == 0 or n == 1):
return 1
else:
return fib_recursion(n-1) + fib_recursion(n-2)

考虑到递归一般都可以改成循环:

1
2
3
4
5
6
def fib_loop(n):
a, b = 1, 1
while(n > 1):
a, b = a+b, a
n -= 1
return a

恩,$O(n)$,线性的时间复杂度,现在看起来蛮不错。但当 n 是 $2^{32}$ 这种大数呢,$O(n)$都嫌弃了。OK,开始介绍矩阵快速幂–把$O(n)$变成$log(n)$的神奇方法。

它其实就是一个公式,重点在于矩阵的构造:
$$ \begin{bmatrix} A(n+1) \\ A(n) \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} * \begin{bmatrix} A(n) \\ A(n-1) \end{bmatrix} $$
更一般的:
$$ \begin{bmatrix} A(n+1) \\ A(n) \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} ^n * \begin{bmatrix} A(1) \\ A(0) \end{bmatrix} $$
其中, n 又可以这么表示(就是转换成二进制啦):
$$ n = \sum_{i=0}^k 2^i $$
代码是这样的:

1
2
3
4
5
6
7
8
9
10
11
12
import numpy as np
def fib_fast_matrix_power(n):
frequency = n - 1 # frequency is less one than the index of factor
result = np.matrix([[1], [1]])
factor = np.matrix([[1, 1], [1, 0]])
while(frequency > 0):
if(frequency%2 == 1):
result = factor * result
factor = factor * factor
frequency //= 2 # be divided with no remainder
print("frequency:", frequency)
return np.array(result)[0][0] # return A(n)

tips: 若是求从 0 到 n 项的和呢?有人给他起了个名字——拓展的 Fibonacci 数列。公式是这样的:
$$ \begin{bmatrix} S(n+1) \\ A(n+1) \\ A(n) \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 1 & 0 \end{bmatrix} ^n * \begin{bmatrix} S(1) \\ A(1) \\ A(0) \end{bmatrix} $$
那么,这种递推公式呢?
$$ \begin{equation} \rm A(n) = \left\{\begin{array}{ll} \alpha A(n-1) + \beta A(n-2) + \gamma A(n-3) & x > 2 \\ 1 & x = 0, 1, 2 \end{array} \right. \end{equation} $$
懒得构造矩阵了~,现在运算时间变短了,但若尝试大数的话,会出现溢出哦。

1000000007

这是一个质数,恩,很大的一个质数 - $10^9+7$ - 出现这种数,肯定是关于高精度整数求模啦。这是几个公式:

1
2
3
4
( a + b ) % c = (( a % c ) + ( b % c )) % c
( a * b ) % c = (( a % c ) * ( b % c )) % c
( a – b ) % c = (( a % c ) – ( b % c )) % c
( a / b ) % c != (( a % c ) / ( b % c )) % c

一定要注意,第四个公式并不成立唉!
据说,这些公式在离散数学里面有提及,但是完全没有印象=_=。

总结

  • 满足递归公式的数列,可以利用矩阵快速幂算法减少运算时间
  • 矩阵快速幂的重点在于矩阵的构造
  • $10^9+7$这样大质数的存在往往伴随着溢出
  • 3+1 个公式

参考资料

斐波那契数列:BestCoder Round #29 1002 || hdu 5171
“OUTPUT THE ANSWER MODULO 10^9 + 7”
What is special/different with number 10^9+7? As most of the coding competitive websites asks to output for large number as modulo 10^9+7.

经历

下面娱乐时间:写一段经历——两天速成矩阵快速幂,也当作矩阵快速幂的一个栗子啦。省略不必要的废话,题目是这样的:

求 Tribonacci 数列从 l 到 r 项的和,结果取$10^9+7$的模。Tribonacci 数列定义如下:
$$ \begin{equation} \rm A(n) = \left\{\begin{array}{ll} A(n-1) + A(n-2) + A(n-3) & x > 2 \\ 1 & x = 0, 1, 2 \end{array} \right. \end{equation} $$
关键字:

  • 大质数 - $10^9+7$ -
  • 很多项的和

傻瓜都知道起码得用循环来解,于是第一版本来了:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include <stdio.h>
#define MOD 1000000007

void tribonacci(long long int l, long long int r){
long long int x=1, y=1, z=1, tmp;
long long int result=0;
long long int i;

if (l < 3){
for(i = l; i <= r && i < 3; i++)
result += 1;
}
for(i = 3; i <= r; i++){ // the i number is z
tmp = x;
x = y;
y = z;
z = (tmp + x + y) % MOD;
if(i >= l)
result = (result + z) % MOD;
}
printf("%lld\n", result);
return;
}

恩,开始测试(就别纠结同数量级下 clock 的大小了,我也不知道为啥前三个会递减)。


r l result clock
0 10 423 19
0 100 440169199 12
0 1000 397128969 10
0 10000 749926090 73
0 100000 142414638 619
0 1000000 817893089 5502
0 10000000 628546478 55406
0 100000000 455282242 540773
0 1000000000 38716501 5394326
0 10000000000 906890177 55216246

到 10 个 0 的时候就有点慢了,但是,给的最大的数据量可是会有 18 个 0 ……我要好好鄙视下@石老板,做出来了也不告诉我关键词。还好有@我荇

我不会告诉你使用矩阵快速幂后差距会这么大:


r l result clock
0 10 423 167
0 100 440169199 24
0 1000 397128969 47
0 10000 749926090 72
0 100000 142414638 77
0 1000000 817893089 94
0 10000000 628546478 92
0 100000000 455282242 117
0 1000000000 38716501 119
0 10000000000 906890177 155
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include <stdio.h>
#define MOD 1000000007

void tribonacci(long long int l, long long int r){
long long int x=1, y=1, z=1, tmp;
long long int result=0;
long long int i;


if (l < 3){
for(i = l; i <= r && i < 3; i++)
result += 1;
}
else{
before_start_sum = 0;
if(l < 3){
for(i=0; i < l; i++)
before_start_sum +=1;
}
else{
before_start_sum = fast_matrices(l-3); // frequency == l-3
}
end_sum = fast_matrices(r-2); // frequency == r-2
result = end_sum - before_start_sum;
if(result < 0)
result = (end_sum + MOD) - before_start_sum;
}
printf("%lld\n", result);
return;
}

long long int fast_matrices(long long int frequency){
martix result = init(3);
martix factor = init(2);
long long int ans;

while(frequency){
if(frequency & 1){ // frequency %= 2
martix_mul(factor, result);
}
martix_mul(factor, factor);
frequency >>= 1; // frequency /= 2
}
martix_des(factor);
ans = result->array[0][0];
martix_des(result);
return ans;
}

还有个小坑,体现在if(result < 0) result = (end_sum + MOD) - before_start_sum;。再次鄙视一下@石老板——因我意识到错误了 - debug 时候偶然输入35 35测试 - 后寻问他,他给我了个错误的答案 =_=

能看到这估计也不容易,在这里说说一些思考的过程以及其他的解法:

先是推导出了一个公式:
$$ S(n+2) = 2 * \sum_{i=1}^{n-1} A_i + A_0 + A_2 + A_n $$
但是不能用,因为得用到除法……

@石老板的思路:

从第六项开始,S(n)也是一个 Tribonacci 数列……

附录:

测试代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#include <stdio.h>
#include <time.h>

// clock() 返回从“开启这个程序进程”到“程序中调用 clock() 函数”时之间的CPU时钟计时单元(clock tick)数

int main(void){
int i;
long long int l, r, start, end;
for(l=0, i=r=1; i < 10; i++){
r *= 10;
start = clock()
tribonacci(l, r);
end = clock();
printf("end-start: %lld\n", end-start);
}
return 0;
}

矩阵相关

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
#include <stdio.h>
#include <stdlib.h>

typedef struct {
long long int **array;
int params[2]; // martix's shape
} martix_, *martix;

// initialization
martix init(int flag){
long long int
ones[4][4] = {
{1, 0, 0, 0},
{0, 1, 0, 0},
{0, 0, 1, 0},
{0, 0, 0, 1}
},
factor[4][4] = {
{1, 1, 1, 1},
{0, 1, 1, 1},
{0, 1, 0, 0},
{0, 0, 1, 0}
},
first[4][1] = {
{3},
{1},
{1},
{1}
};

martix array = (martix)malloc(sizeof(martix_));
if(array == NULL)
exit(-1);
array->params[2] = {4};
switch(flag){
case 1:
case 2:
array->params[1] = 4;
break;
case 3:
array->params[1] = 1;
break;
}
array->array = (long long int **)malloc(array->params[0] * sizeof(long long int *));
if (NULL == array->array)
exit(-1);
int i, j;
for(i = 0; i < array->params[0]; i++){
array->array[i] = (long long int *)malloc(array->params[1] * sizeof(long long int));
if(array->array[i] == NULL)
exit(-1);
for(j = 0; j < array->params[1]; j++){
switch(flag){
case 1:
array->array[i][j] = ones[i][j];
break;
case 2:
array->array[i][j] = factor[i][j];
break;
case 3:
array->array[i][j] = first[i][j];
break;
}
}
}
return array;
}

// y = (x + y) % MOD = ((x % MOD) + (y % MOD)) % MOD
void martix_sum(martix x, martix y){
int i, j;
for(i = 0; i < x->params[0]; i++)
for(j = 0; j < x->params[1]; j++)
y->array[i][j] = (x->array[i][j] + y->array[i][j]) % MOD;
return;
}

// y = (x * y) % MOD = ((x % MOD) * (y % MOD)) % MOD
void martix_mul(martix x, martix y){
long long int tmp[x->params[0]][y->params[1]];
int i, j, k;
for(i = 0; i < x->params[0]; i++)
for(j = 0; j < y->params[1]; j++){
tmp[i][j] = 0;
for(k = 0; k < x->params[1]; k++){
tmp[i][j] += (x->array[i][k] * y->array[k][j]) % MOD;
}
tmp[i][j] %= MOD;
}
for(i = 0; i < x->params[0]; i++)
for(j = 0; j < y->params[1]; j++){
y->array[i][j] = tmp[i][j];
}
y->params[0] = x->params[0];
return;
}

// destroy
void martix_des(martix x){
int i;
for(i = 0; i < x->params[0]; i++)
free(x->array[i]);
free(x->array);
free(x);
return;
}