存档

文章标签 ‘数列’

差分序列

2011年9月29日 8 条评论

不小心遇到一题4次求和，推公式推了好久。。想想实在是没有必要，每次遇到还得推一边实在是累。。所以，回来就写了一个差分的版本，单次计算的效率不比直接公式差，不容易出错，而且，再也不用蛋去推多项式公式了

至于什么是差分序列，这个我就不解释了，wiki上应该可以找到。我的写法是借鉴了牛顿插值的，使用了一个递推的写法（一般求多项式都会用这个， 比如4n^3 + n^2 + 2 * n + 3，一般写成 ((4n + 1) * n + 2) * n + 3）

初始化细节说明：order是生成的多项式的阶，请确保它大于或等于实际的多项式阶。order阶的多项式需要用order+1项进行init，请把这前order+1项存入数组cof中。最后调用init函数，就可以用calc来计算任意项了。

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

/******************************************************************************
  * Differential Sequence
  * init: O(order^2)
  * you should init the cof[i] with the value of i-th value of sequence
  * the sequene start from startfrom, default is 0, so a[0] is the first one
  * 	of sequence a
  * calc: O(order)
  ****************************************************************************/
 
const int ORDERLIM = 100;
 
double cof[ORDERLIM + 1];
int order;
const int startfrom = 0;
 
void init()
{
	int i, j;
	for(i = 1; i <= order; i++) {
		for(j = order; j >= i; j--) {
			cof[j] -= cof[j - 1];
		}
	}
}
 
double calc(double n)
{
	int i;
	double res = 0;
	for(i = order; i > 0; i--) {
		res += cof[i];
		res *= (n - i + 1 - startfrom) / i;	
	}
	res += cof[i];
	return res;
}
 
/******************************************************************************
  * MOD version of Differential Sequence
  * cof should not overflow MOD
  * MOD should should not less than order
  * MOD should be prime
  * need inverse in mod MOD system
  ****************************************************************************/
typec MOD = 1000000007;
typec cof[ORDERLIM + 1];
typec inv[ORDERLIM + 1];
int order;
const int startfrom = 0;
 
void init()
{
	int i, j;
	for(i = 1; i <= order; i++) {
		inv[i] = inverse(i, MOD);
		for(j = order; j >= i; j--) {
			cof[j] -= cof[j - 1];
		}
	}
}
 
typec calc(typec n)
{
	int i;
	typec res = 0;
	for(i = order; i > 0; i--) {
		res += cof[i], res %= MOD;
		res *= (n - i + 1 - startfrom) % MOD * inv[i] % MOD;
		res %= MOD;
	}
	res += cof[i], res %= MOD;
	while(res < 0) res += MOD;
	return res;
}
分类: ACM, C, Math 标签: ,

Fibonacci序列 及其前n项和

2010年5月7日 3 条评论

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;
}
分类: Cpp, Math 标签: ,