首页 > ACM, C, Math > 差分序列

差分序列

2011年9月29日 发表评论 阅读评论

不小心遇到一题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 标签: ,
  1. bnuzhanyu
    2011年10月9日17:02 | #1
    回复 | 引用

    cof怎么初始化的呢？这里cof[0]程序并没给出。

  2. bnuzhanyu
    2011年10月9日18:46 | #2
    回复 | 引用

    好像懂了，不过您的calc函数貌似有点问题。
    我改了一下：
    int function(int k)
    {
    return k*k*k + k*k + 1;
    }

    void init()
    {
    int i, j;
    for(int i=0; i<=order; ++i)
    cof[i] = function(i);//e.g. f(i)=i^4

    inv[0] = 1;
    for(i = 1; i = 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;
    }
    while(res < 0) res += MOD;
    return res;
    }
    欢迎指正~

  3. bnuzhanyu
    2011年10月9日18:47 | #3
    回复 | 引用

    inv[0] = 1;
    for(i = 1; i = i; j--) {
    cof[j] -= cof[j - 1];
    }
    }

  4. bnuzhanyu
    2011年10月9日18:49 | #4
    回复 | 引用

    不知道怎么回事~ 那一段是inv[0]的值是1，inv[i]的值是i+1对于MOD的逆

  5. Earthson
    2011年10月21日12:25 | #5
    回复 | 引用

    @bnuzhanyu
    cof直接用前n+1项初始化~

  6. Earthson
    2011年10月21日12:30 | #6
    回复 | 引用

    @bnuzhanyu
    我是直接用的inv[i]是i对MOD的逆，浪费了inv[0]的^^

  7. Earthson
    2011年10月21日12:31 | #7
    回复 | 引用

    @bnuzhanyu
    你的代码里面，没有看到inv的初始化哎，只初始化了inv[0]~

  8. Earthson
    2011年10月21日12:35 | #8
    回复 | 引用

    @bnuzhanyu
    不好意思啊，好久没上blog了~ 你那个代码，在init的时候貌似写错了比较多的地方，包括循环参数什么的。
    我的代码，好多地方我用参数调的（懒得把循环参数想太清楚，所以就看运行结果，再调一下，看对了就不管了^^），所以有些循环参数看起来会比较别扭~

  1. 本文目前尚无任何 trackbacks 和 pingbacks.