Fibonacci序列：1，1，2，3，5，8，13，21，34，55...
递推式：fib(n)=fib(n-1)+fib(n-2),(n>2)
至于通项就不说了，找起来很容易
现在我们要计算Fibonacci的第n项，用递推是容易的，这是一个O(n)的算法，但是当n很大的时候显然不行（当然这个时候fib(n)本身也非常大，暂且不管它）。我刚学会了一种新方法——矩阵
其中的fib(n)=Fn.以次类推.这个公式很好证明.当然只是证明是不够的，我还没吃透它，先不管这点，看看怎么使用它.
令左边的矩阵为A，右边的矩阵为F(n+1), 则F(n)=A^(n-1)
定义矩阵的乘法之后，用逐次平方计算即可得到F(n),也就能得到fib(n),我在下面给出代码
那么我们怎么计算fib序列的前n项和呢？
考虑到S(n)=F(1)+F(2)+...+F(n-1)+F(n)=A^0+A^1+A^2+..+A^(n-1)
注意到中学的等比数列公式没?
(A-I)S(n)=A^n-I
其中I是二阶单位矩阵，左乘A-I的逆(A-I)^-1,就能得到S(n)了.
你算下就知道，逆就是A.于是S(n)=A*(A^n-I)逐次平方搞定，类似的还可以计算连续n项的和，提出第一项对应的矩阵就明白了，下面给出相关矩阵的代码，以及逐次平方实现
****推荐下英文的wiki，上面对于fib序列的描述是令人满意的
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
| struct matrix2
{
long long num[2][2];
matrix2(){num[0][0]=1;num[1][1]=1;num[0][1]=0;num[1][0]=0;}// 初始化为单位矩阵
};
int p=2010;//以为fib序列很大，所以结果全部mod p
matrix2 operator*(matrix2 x,matrix2 y)
{
matrix2 solution;
solution.num[0][0]=((x.num[0][0]*y.num[0][0])%p+x.num[0][1]*y.num[1][0])%p;
solution.num[0][1]=((x.num[0][0]*y.num[0][1])%p+x.num[0][1]*y.num[1][1])%p;
solution.num[1][0]=((x.num[1][0]*y.num[0][0])%p+x.num[1][1]*y.num[1][0])%p;
solution.num[1][1]=((x.num[1][0]*y.num[0][1])%p+x.num[1][1]*y.num[1][1])%p;
return solution;
}
matrix2 operator-(matrix2 x,matrix2 y)
{
matrix2 solution;
solution.num[0][0]=x.num[0][0]-y.num[0][0];
solution.num[0][1]=x.num[0][1]-y.num[0][1];
solution.num[1][0]=x.num[1][0]-y.num[1][0];
solution.num[1][1]=x.num[1][1]-y.num[1][1];
return solution;
}
matrix2 operator^(matrix2 h,int n) //计算h^n,逐次平方法O(lgn).令h=A，计算F(n+1)
{
matrix2 solution;
while(n)
{
if(n%2)
{
solution=solution*h;
}
h=h*h;
n>>=1;
}
return solution;
} |
首先给出一个特殊情况的分数m/n循环节输出(代码,基于小学就学过的除法笔算算法，p是进制，当然你可以全部用10替换)，我们先假设它的循环节是从小数点第一位开始（在下面我会用进制变换证明）
void printDecimalRec(int m,int n,int p) //这里假设m<n,且(n,p)=1,且输出省略“0.”
{
int flag=m;
do
{
m*=p;
printf("%d",m/n);
m%=n;
}while(m!=flag);
}
给指针分配的一块内存，都应该是连续的。这必然会导致一个问题：当我分配的空间过大，系统将需要整理出一块有效的地方，这显然是不太合适的（事实上，系统不会支持过大的内存请求的，在我这里不能超过一百万）
我的策略是通过分配较小的单元用另一个指针数组对这些单元索引，类似下图
如果单元的大小为k，指针数组的大小为t，那么就可以处理k*t大小的空间。比如k=100,t=100，就是大小10000空间，只要我建立合适的函数operator[]，就可以对这部分空间索引，实现一个模拟的数组
当然这只是双层的索引（基本的数组就是单层的），我在下面给出了一个三层的索引，能够处理lengthunit^3的空间，这其实是我为大数类设计的，用初始3*lengthunit的空间，实现一个开放数组（虽然是有上限的，但它是使用时才分配的）
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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
| //############numarray.h#################################
//###初始化注意：确保ptr1和ptr2指向的位置不是NULL，并且已经分配了空间######
#ifndef NUMARRAY_H_INCLUDED
#define NUMARRAY_H_INCLUDED
typedef unsigned int numunit;
const int lengthunit=1000;
class numArray
{
private:
numunit ***num;
int ptr1; //顶层的索引，已分配空间和未分配空间的分界线
int ptr2; //第二层索引，类似ptr1，该位置相对于ptr1
public:
numArray();
numArray(const numArray&);
~numArray();
numArray& operator=(const numArray&);
//###########
numunit& operator[](int n);//若空间不存在，会自动创建，并初始化为0
const numunit& operator[](int n) const;//该函数不会分配空间
void resize(int n); //可以额外维护一个长度，因为长度的定义不属于该类
private:
void resize(int,int);//resize(int)就是调用它实现的
};
#endif // NUMARRAY_H_INCLUDED
//#################numarray.cpp###################
#include "numarray.h"
#include <stdio.h>
numArray::numArray() //构造函数，负责最基本的初始化。
{ //注意这里的初始化原则，ptr1和ptr2指向的位置必须分配空间
num=new numunit**[lengthunit];//为顶层的索引数组分配空间
num[0]=new numunit*[lengthunit];//为顶层第一个位置指定空间
num[0][0]=new numunit[lengthunit];//为第二层索引的第一个位置分配空间
num[0][0][0]=0; //初始化归0，下同
for(int i=1;i<lengthunit;i++)
{
num[i]=NULL;
num[0][i]=NULL;
num[0][0][i]=0;
}
ptr1=0;
ptr2=0;
}
numArray::numArray(const numArray &in)
{
num=new numunit**[lengthunit];
for(int i=0;i<lengthunit;i++)
{
if(in.num[i]==NULL)
{
num[i]=NULL;
}
else
{
num[i]=new numunit*[lengthunit];
for(int j=0;j<lengthunit;j++)
{
if(in.num[i][j]==NULL)
{
num[i][j]=NULL;
}
else
{
num[i][j]=new numunit[lengthunit];
for(int t=0;t<lengthunit;t++)
{
num[i][j][t]=in.num[i][j][t];
}
}
}
}
}
ptr1=in.ptr1;
ptr2=in.ptr2;
}
numArray::~numArray()
{
for(int i=0;i<lengthunit;i++)
{
if(num[i]!=NULL)
{
for(int j=0;j<lengthunit;j++)
{
if(num[i][j]!=NULL)
{
delete [] num[i][j];
}
}
delete [] num[i];
}
}
delete [] num;
}
numArray& numArray::operator=(const numArray &in)
{
resize(in.ptr1,in.ptr2);
for(int i=0;i<=ptr1;i++)
{
for(int j=0;j<=ptr2;j++)
{
for(int t=0;t<lengthunit;t++)
{
num[i][j][t]=in.num[i][j][t];
}
}
}
return *this;
}
numunit& numArray::operator[](int n)
{
int p1=(n/lengthunit)/lengthunit;
int p2=(n/lengthunit)%lengthunit;
if(ptr1<p1||(ptr1==p1&&ptr2<p2))
{
resize(n);
}
return num[p1][p2][n%lengthunit];
}
const numunit& numArray::operator[](int n) const
{
int p1=(n/lengthunit)/lengthunit;
int p2=(n/lengthunit)%lengthunit;
if(num[p1][p2]==NULL)
{
printf("error:out of space!\n");
system("pause");
}
return num[p1][p2][n%lengthunit];
}
void numArray::resize(int n)
{
int p1=(n/lengthunit)/lengthunit;
int p2=(n/lengthunit)%lengthunit;
resize(p1,p2);
}
void numArray::resize(int p1,int p2)
{
if(ptr1>p1)
{
for(int i=p1+1;i<=ptr1;i++)
{
for(int j=0;j<=lengthunit;j++)
{
if(num[i][j]!=NULL)
{
delete [] num[i][j];
}
}
delete [] num[i];
num[i]=NULL;
}
}
else if(ptr1<p1)
{
for(int j=1;j<lengthunit;j++)//分配num[ptr1]的空间，注意该空间被初始化过
{
if(num[ptr1][j]==NULL)
{
num[ptr1][j]=new numunit[lengthunit];
for(int t=0;t<lengthunit;t++)
{
num[ptr1][j][t]=0;
}
}
}
for(int i=ptr1+1;i<p1;i++)//这部分空间是完全用0填充的
{
num[i]=new numunit*[lengthunit];
for(int j=0;j<lengthunit;j++)
{
num[i][j]=new numunit[lengthunit];
for(int t=0;t<lengthunit;t++)
{
num[i][j][t]=0;
}
}
}
num[p1]=new numunit*[lengthunit];//num[p2]是要初始化的，不能随便填0
ptr2=0;
num[p1][0]=new numunit[lengthunit];//确保了ptr2指向的位置分配了空间
for(int t=0;t<lengthunit;t++)
{
num[p1][0][t]=0;
}
for(int j=1;j<lengthunit;j++)
{
num[p1][j]=NULL;
}
} //到此为止，p1和ptr1应该是相等的
ptr1=p1;
if(ptr2<p2) //开始处理ptr2
{
for(int j=ptr2+1;j<=p2;j++)
{
if(num[ptr1][j]==NULL)
{
num[ptr1][j]=new numunit[lengthunit];
for(int t=0;t<lengthunit;t++)
{
num[ptr1][j][t]=0;
}
}
}
}
else if(ptr2>p2)
{
for(int j=p2+1;j<=ptr2;j++)
{
if(num[ptr1][j]!=NULL)
{
delete [] num[ptr1][j];
num[ptr1][j]=NULL;
}
}
}
ptr2=p2;
} |
这个计算的目的是什么？当组合数较大的时候，显然不能很好保存，因此把它转化为p的模（p是任意的）。
方案很简单，仍然利用公式的除法优先（我的另一篇文章组合数计算），不同的是，这次需要保证除法的绝对优先，这将导致不那么高的效率
首先明确一点a*b/c mod t和(a mod t)*b/c mod t可能是不同的，因为如果a能被c整除，a mod t却不一定能被c整除，(a mod t)*b/c这个式子甚至的结果甚至不是整数（如果那个除号不是整除的话。）比如：4*5/2 mod 3=1，而 4 mod 3=1 ，1*5/2不是一个整数。如果你看了我的那篇文章就知道，在那个算法中我对这一点的依赖
因此，我选择用数组保存了公式
上边的因子。现在，我的目标是，把下面的因子全部约掉，我使用gcd（最大公约数）找上面元素和下面元素的因子，并且约掉该因子，扫描……显然最后下面的因子全部归1
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
| int gcd(int a,int b)
{
int swap;
while(b!=0)
{
swap=b;
b=a%b;
a=swap;
}
return a;
}
int combinationNumMod_p(int n,int k,int p)
{
int solution=1;
int *dividend;
k=(n-k>k)?k:n-k;
dividend=new int[k];
for(int i=0;i<k;i++)
{
dividend[i]=n-i;
}
for(int i=2;i<=k;i++) //i就是生成的下因子
{
int s,g;
s=i; //把i的值转给s，因i的值不能随便改（i负责生成下一个s）
for(int j=0;j<k;j++) //对上因子扫描
{
if(s==1) //如果当前下因子已经是1了，就跳出循环
{
break;
}
g=gcd(s,dividend[j]); //这个算法的还依赖于gcd，g是保存的最大公约数
if(g!=1) //公约数不等于1就把公约数约掉
{
s/=g;
dividend[j]/=g;
}
}
}
for(int i=0;i<k;i++) //把处理好的上因子乘起来
{
solution*=dividend[i]%p; //公式a*b mod p=(a mod p)*(b mod p) mod p
solution%=p;
}
return solution;
} |
这是一个困扰我很久的问题，刚上高中就开始找它的计数方案（证明），直到高三才搞定。我的把这种计数的情况称为潜在空位和潜在禁止位，它能解决一类组合计数问题。
当然，先看一下，可重复的组合。
有n种元素的一个集合，每种的数目假设无限（至少比k大），现在我们要选取其中的k个元素，令该组合所有情况的集合为T，求组合数。在书上我们很容易查到，它等于C(n+k-1,k)
先考虑一种一般的策略（在非可重复组合中）：有n个座位，我要选择其中的k个（我把这个集合称为C(n,k)），为了方便，同时我把座位标号（1，2，3，……），我可以这样做，首先选择一个座位，然后在该座的后边（标号更大）选择第二个，……依次选择，直到最后一个。很容易证明它没有重复，它对C(n,k)能建立一一映射
对于可重复组合，我只是加入了一个潜在位而已，比如，当我选择1号位，那么在1号位的后面将产生一个新的1号位，可供后面选择，选择k位的情况将会产生一个n+k个位子的现象，但是请注意，选择的第k个位子后面产生的潜在位，不能供任何其它位选择（或者说，那个位子在任何情况下都不会被选择到），所以实际上该位不是一个供选择的位。这样产生了k-1个“虚位”，很容易建立它与C(n+k-1,k)的一一映射。而且，它对T也是一一映射的（在下图中，最后一句话的意思是，在n+k-1位的等价图中，不存在两种不同情况，他们解释到T的时候，是一样的）
现在来看另一种情况，即我所谓的潜在禁止位：当我选择位子的时候，有这样的要求，相邻位子的必须有一个最小间隔，甚至不同位子间最小间隔可以不同。这其实并不困难，我们可以这样理解：每个选择了的位子后面有一个给定长度s的禁止位子，当我选择某个位子的时候，后面s个位子就等于消失了。（或者你可以理解为我选择的位子长度不是1，而是s+1……）这与重复组合的情况是类似的，同样注意最后一个被选择的位子。容易得到，对于给定的最小间隔s，有公式C(n-s(k-1),k),其中的s(k-1)是对禁止位的计数（很奇怪吧，这些禁止位的位置是不确定的，但他们确实存在，而且就是有那么多位，对于不同间隔的情况，把它的禁止位数替换上去就可以了）