由菲波拉契数列想到的问题

iceheart 2008-03-29 01:06:15
问题1:
对于这样一个数列
f(1) = 1; f(2) = 1; f(3) = (1); f(4) = 2;... f(n) = f(n - 1) + f(n - 3);
能否找出一种时间复杂度为O(logN)的方法求数列的第n项?

问题2:
对如下更一般的数列
(1)
f(1) = f(2) = ... = f(k) = 1;
f(n) = f(n - 1) + f(n - 2) + ... + f(n - k);
(2)
f(1) = f(2) = ... = f(k) = 1;
f(n) = f(n - 1) + f(n - k);
能否找出时间复杂度为O(logN)的方法求数列的第n项?

...全文
815 点赞 收藏 18
写回复
18 条回复
切换为时间正序
当前发帖距今超过3年,不再开放新的回复
发表回复
iceheart 2008-04-11
多谢楼上关注。
我只是做了一些测试,通过观察到的结果产生这样的想法的,只是猜想,呵呵。
不过对于f(n) = f(n-1) + a * f(n-2) 中 f(n)/f(n-1) == f(n-1)/f(n-2),当a>0时
我觉得这个猜想应该是成立的,其他情况还没仔细考虑,我也是猜测,没证明过,如果错了就再想,呵呵。
对于这种方法求平方根,我知道肯定根标准做法大相径庭,不过我觉得这是个思路。
我的意思是不只把菲波拉契数列当个数列来算,那样也没太大意思。有可能的话就从这个数列的特性中挖掘些有用的东西,如果能在其他领域有所突破那岂不是更好。希望能起到抛砖引玉的作用:)

平方根的快速算法,从网上搜集到的就是卡马克在quake3中用的方法,据说比C库函数的快百分之二十,精度不高,
不过足够实用了。
上面写的菲波拉契的方法中用到了取中值的思想,目的是让数列快速收敛,目前看到的效果是求一个平方根大概要平均循环15到20次,后来试了下指数除2的方法取估计值,不过效果并不理想,反而效率更低了。我想可能每次循环取中值就已经达到了比这个方法更好的效果。也许这个结果还可以用到其他优秀的开平方算法上。
回复
medie2005 2008-04-09
当然,要是lz是弄着玩的,就另说了。

另外:这几句:“形如 f(n) = a * f(n-1) + b * f(n-2) + c * f(n-3)...的数列都是等比数列。 n趋于无穷大时 f(n) / f(n-1)趋于固定值。 ”
貌似不对哦。因为有可能存在若干个大于1的特征根。
回复
medie2005 2008-04-09
我觉得lz方向搞错了:
你应该求平方根的快速算法,用以快速计算Fibonacci数列;而不是由Fibonacci数列来开发快速平方根算法。
有点本末倒置了。
回复
iceheart 2008-04-08
这个方法所依赖的就是最快的求出 f(n)/f(n-1),我用的是取两项中值的办法。
不知道各位还有没有更好的方法?
回复
iceheart 2008-04-08
谢谢楼上各位的指教,这个问题我基本了解了。
mathe和medie2005 说的非常清楚了,尤其是mathe后面的复杂度分析,谢谢!
不过我觉得能写出充分好的复杂度为O(n*log(n))的大整数乘法难度也很高,
考虑整数的复制等等因素,可能会导致有更多的时间浪费。

不好意思,这贴稍等等再结,分都是上面几位的,呵呵。
偶然发现了个通过菲波拉契数列求浮点数平方根的方法,不开新贴了,直接在这里说了,呵呵
根据观察到的现象,当n趋近于无穷大时,
形如 f(n) = a * f(n-1) + b * f(n-2) + c * f(n-3)...的数列都是等比数列。
n趋于无穷大时 f(n) / f(n-1)趋于固定值。
而对于 f(n) = f(n-1) + a * f(n-2)这样的数列,置f(0) = 1, f(1) = f(n) / f(n-1),
则此数列第n项就是f(1)^n;

对于数列 a * f(n-2) + f(n-1) = f(n) 有:
f(n) / f(n - 1) = (1 + sqrt(4 * a + 1) ) / 2;
于是sqrt(4 * a + 1) = 2 * (f(n)/f(n-1)) - 1;
如果能快速求出f(n)/f(n-1), 就能快速求出 4*a+1的平方根

已知a,浮点数求数列f(n)=f(n-1) + a*f(n-2)中的 f(n)/f(n-1)
double f(int n, double x) {
double buf[2] = {1, a};
if ( n < 2 ) return 1;
for ( int i = 2; i <= n; i++ ) {
buf[i & 1] = buf[i & 1] * x + buf[(i + 1) & 1];
}
return buf[n & 1]/buf[(n+1) & 1];
}

进而演化出以下方法

double sqrt2( double x )
{

if ( x <= 0 ) return 0;
int n = 2;
double a = 1;

double y = (x - 1) / 4;
double buf[2];
double b = 0;
for ( int i = 0; a > b || a < b; i++ ) {
buf[0] = 1; buf[1] = a;
double last = 1;
//这个循环求数列n、n-1、n-2项
for ( int j = 2; j <= n; j++ ) {
last = buf[j & 1];
buf[j & 1] = last * y + buf[(j + 1) & 1];
}
b = a;
//(f(n)/(f(n-1) + f(n-1)/f(n-2)) / 2,结果作为f(1)
//这是因为f(n)/f(n-1) 和f(n-1)/f(n-2) 总是在极限值左右震荡,取中值逼近更快
a = (buf[(n-1) & 1] / last + buf[n & 1] / buf[(n-1) & 1]) / 2;
}
return 2 * a - 1;
}

对这个方法再稍微优化以下

double sqrt3( double x )
{
if ( x <= 0 ) return 0;
if ( x < 1 ) return 1 / sqrt3(1 / x);
if ( x < 100 ) return sqrt3( 100 * x ) * 0.1;
double a = 1;
double y = (x - 1) / 4;
double last = 0;
while (a > last || a < last) {
last = a;
a = (a + y / a + 1) / 2;
}
return (2 * a - 1);
}
统计了一下,算1~0x100000 花的时间大概是 c库函数的15倍左右。
c库中的用的是协处理器指令。


回复
rover___ 2008-04-01
楼上各位分析得透彻。赞一个。
回复
mathe 2008-04-01
对于这个问题,分析计算复杂度的一个大问题是计算需要log(n)次大数乘法。
如果假设这个递推公式特征方程中范数(或者说绝对值)最大的根的范数为M
那么我们可以有个估计公式
|f(n)|<=b*M^n (其中b为一个常数)
也就是说,f(n)的长度同n成正比。通常我们认为对于大数乘法,充分好的算法做一次乘法的时间复杂度近似为O(nlog(n))
所以上面算法过程的时间复杂度可以认为大概为O(n*log(n)^2).
而如果直接用递推公式累加,那么时间复杂为O(n^2) (一次加法时间为O(n))
回复
zgg___ 2008-04-01
这个问题的答案是肯定的,而且过程也不大复杂的。我来把大家的想法总结一下吧,呵呵。

我们可以注意到下面的事实:

1、对于某递推关系,都可以求得矩阵m,使得关系X(i+1)=mX(i)成立。

其中X(i)是n维矢量,有X(i)={x(i+n-1),...,x(i)},m是n阶方阵。
拿lz的数列为例:对于数列1,1,1,2,3,4,6……一定存在3乘3的方阵m,
使得m满足{2,1,1}=m{1,1,1},{3,2,1}=m{2,1,1},{4,3,2}=m{3,2,1},{6,4,3}=m{4,3,2}……
我们可以看到m可以是{{1, 0, 1}, {1, 0, 0}, {0, 1, 0}},这个已经在6层及之前都说清了。
因此,对于某个线性递推关系,问题都可以转化成为求给定n阶方阵m的N次幂的问题。

2、求某矩阵m的N次幂,存在log(N)的算法。

这个问题在论坛中讨论过了。不光对于矩阵,对于一个数,对于一个多项式,对于一个向量,只要是定义了乘法的对象,我们都可以用共同的方法求这个对象的N次幂。(如果有觉得由有理数构成的矩阵比实数更象“数”的同学,我们可以扩展讨论,呵呵。)具体的方法就是:选择性的进行平方、乘m的操作就可以了。另一方面,矩阵的大小会影响到每一次计算的速度,但是这个和算法是log(N)的并不矛盾,这里假设为N远大于n的情况。还有就是矩阵的每一项或许会变得很大的问题,因此而涉及的存贮、效率和精度的问题,或许也不是这个问题所要解决的了。

3、上面两个一综合,就可以了吧。

对于求矩阵的幂,我们可以用求特征值的方法。我们知道求多项式根的近似值也有很好的算法的。

回复
medie2005 2008-03-31
那是因为(1-√5)/2)^(n+1)这个数很小。
证明的话,只需要对开始几个小的n值验证一下就可以了。
回复
mathe 2008-03-31
证明|((1-sqrt(5))/2)^(n+1)|/sqrt(5)<1/2就可以了
而由于|(1-sqrt(5))/2|<1,所以只要证明对于n=0成立就可以了,而这个结论是显然的。
回复
njwangchuan 2008-03-30
f(n) = 1 /√5 * [((1 + √5) / 2)^(n+1) - ((1 - √5) / 2)^(n+1)]

不过这个还不是最快的

f(n)=math.round(1/√5*((1+√5)/2)^n)),math.round表示取最接近整数,公式须证明
回复
medie2005 2008-03-30
当k<5时,的确可以通过解特征方程来获得数列的通项表达式。
但当k>=5时,特征方程是最高次数大于4的一元高次方程了,如我在4楼说的:
“四次以上的一元高次方程无公式解.有的可以 ,有的不可以,伟大的数学家阿贝尔和伽罗瓦已经讨论得很清楚了.”
因此,用手工求特征方程的解是不太可能了,也就是说,得不到数列的通项表达式。

当然,如果你对解的精度要求不是很高,可以用计算机求得特征方程的渐进数值解。

当k<5时,以菲波拉契数列为例:
递归式:
f(1) = f(2) = 1; (1)
f(n) = f(n - 1) + f(n - 2); (2)
递归关系(2)确定了以下的特征方程(原理是Z变换):
x^2-x-1=0
求上面的方程的解,得:
x1=(1+sqrt(5))/2;
x2=(1-sqrt(5))/2;
于是f(n)=C1*x1^n+C2*x2^n (C1、C2由初始值确定,即由(1)式确定)。
再令n=1和n=2,得到方程组:
C1*x1+C2*x2=1
C1*x1^2+C2*x2^2=1
解方程组,就得到了C1和C2。于是得到了f(n)的通项表达式。
回复
mathe 2008-03-29
如果精度要求不是很高,那么上面问题总可以通过求特征方程的数值解在O(1)的时间内算出任意一项的值。
但是如果要得到精确的精度,那么乘法本身的时间复杂度就要增加。
但是模范3楼中矩阵方法,可以转化为一个计算k阶方阵n次幂问题,这个可以通过O(log(n))次乘法解决
回复
medie2005 2008-03-29
请教k小于等于4的时候的公式
========================
2次方程就不说了,韦达定理.

一元三次方程实根的解法:
这是意大利数学家达达利得出的,并没有公开,后来意大利数学家卡当向他讨教,并发誓永不外传,1544年才从达达利那里获得了这个解法。谁知,1545年,卡当就把这个解法写进他的著作《大法》里面了。所以,从此这个方法叫做卡当法。
对任意
ax^3+3bx^2+3cx+d=0..................(1)

x=y-b/a
p=c/a-(b^2)/(a^2)
q=(b^3)/(x^3)-(3bc)/(2a^2)+d/a ...........(2)
带入(1)有,
y^3+3py+2q=0
引入方程
y=u-p/u
=>
(u^3)^2+2qu^3-p^3=0...................(3)
解关于(u^3)的一元二次方程:
u^3=-q+sqr(q^2+p^3)
or
u^3=-q-sqr(q^2+p^3)
=>
y=(-q+sqr(q^2+p^3))^(1/3)+(-q-sqr(q^2+p^3))^(1/3)
x=y-b/a

一元四次方程实根的解法:
这个方法是卡当的学生费拉里给出的。
对任意
x^4+bx^3+cx^2+dx+e=0
变化为:
x^2*(x+b/2)^2-(b^2)*(x^2)/4+c*x^2+dx+e=0
待定变量y
((x^2+bx/2)+y/2)^2=((b^2)*(x^2)/4)*cx^2-dx-e+(x^2+bx/2)y+(1/4)y^2=(b^2/4-c+y)x^2+(by/2-d)x+y^2/4-e(*)
取适当的y,使右边为一个完全平方式
即delta=0
(by/2-d)^2-4(b^2/4-c+y)*(y^2/4-e)=0
这是一个三次方程,解法上面已经给出,得出y值
带入(*) 可以化为关于x的两个一元二次方程。
既可以得到x的四根。

四次以上的一元高次方程无公式解.有的可以 ,有的不可以,伟大的数学家阿贝尔和伽罗瓦已经讨论得很清楚了.
回复
medie2005 2008-03-29
是的,求菲波拉契数列的第n项可以用矩阵乘法在O(log(n))内求解(不考虑大数乘法的复杂度).

以菲波拉契数列来说.
fib(1) = fib(2) = 1; (1)
fib(n) = f(n - 1) + f(n - 2); (2)
把(2)写成矩阵相乘的形式:

| fib(n) | |1 1| |fib(n-1)|
| |=| |*| | (3)
|fib(n-1)| |1 0| |fib(n-2)|

我们将(3)展开就得到:

| fib(n) | |1 1|^(n-2) |fib(2)|
| |=| | * | | (4)
|fib(n-1)| |1 0| |fib(1)|


|1 1|^(n-2)
| |
|1 0|
类似与a^n.我们用计算a^n的方法来计算这个式子就可以得到O(log(n))[前提是不超过内置类型的表示范围].

lz可以仿照上面的方法自己解决问题1和问题2.
实际上就是求解差分方程,可以学习一下Z变换.
回复
iceheart 2008-03-29
刚看了一下方阵的乘法运算规则.

http://topic.csdn.net/u/20071216/22/ffd33e30-64d3-40c3-92ca-8b31f7743a04.html
两位在上面这个帖子里的精彩回答我也看过了.
非常感谢楼上几位的帮助, 总算明白点了,呵呵.

当k为3时
对于 数列f(n) = f(n-1) + f(n - 3),对应的矩阵乘法应该是下面的形式
|f(n) | | 1 0 1 |^n-3 | f(3) |
|f(n-1)| = | 1 0 0 | * | f(2) |
|f(n-2)| | 0 1 0 | | f(1) |

对于 数列f(n) = f(n-1) + f(n-2) + f(n-3),对应的矩阵乘法如下
|f(n) | | 1 1 1 |^n-3 | f(3) |
|f(n-1)| = | 1 0 0 | * | f(2) |
|f(n-2)| | 0 1 0 | | f(1) |


方阵的n次方的算法有点迷茫, 每次两个k阶方阵相乘,都要计算k^3次乘法和(k - 1) * k^2 次加法
所以计算k阶数列第n项的算法的时间复杂度变成了O((k^3)*log2(n));

当n很小时,还不如递推的速度快些.
有没有像mathe在上面那个帖子里那样归纳出来的方法呢?
回复
iceheart 2008-03-29
[Quote=引用 1 楼 jiangsheng 的回复:]
1 不是有公式么?算乘方的时候可以用二进制法降低复杂度到log(2) N
2 对于5次及以上的一元高次方程没有通用的代数解法,所以K>4的时候是没有公式可以套的
[/Quote]
还有些疑问
1、是指矩阵的乘法么?(线性代数毕业时就还给老师了,偶去找回来)
2、请教k小于等于4的时候的公式

非常感谢
回复
蒋晟 2008-03-29
1 不是有公式么?算乘方的时候可以用二进制法降低复杂度到log(2) N
2 对于5次及以上的一元高次方程没有通用的代数解法,所以K>4的时候是没有公式可以套的
回复
相关推荐
发帖
数据结构与算法
创建于2007-08-27

3.2w+

社区成员

数据结构与算法相关内容讨论专区
申请成为版主
帖子事件
创建了帖子
2008-03-29 01:06
社区公告
暂无公告