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

[iceheart]

问题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项?

[蒋晟]

1 不是有公式么？算乘方的时候可以用二进制法降低复杂度到log(2) N
2 对于5次及以上的一元高次方程没有通用的代数解法，所以K>4的时候是没有公式可以套的

[iceheart]

[Quote=引用 1 楼 jiangsheng 的回复:]
1 不是有公式么？算乘方的时候可以用二进制法降低复杂度到log(2) N
2 对于5次及以上的一元高次方程没有通用的代数解法，所以K>4的时候是没有公式可以套的
[/Quote]
还有些疑问
1、是指矩阵的乘法么？(线性代数毕业时就还给老师了，偶去找回来)
2、请教k小于等于4的时候的公式

非常感谢

[iceheart]

刚看了一下方阵的乘法运算规则.

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在上面那个帖子里那样归纳出来的方法呢?

[medie2005]

是的,求菲波拉契数列的第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变换.

[medie2005]

请教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的四根。

四次以上的一元高次方程无公式解.有的可以 ，有的不可以，伟大的数学家阿贝尔和伽罗瓦已经讨论得很清楚了.

[mathe]

如果精度要求不是很高，那么上面问题总可以通过求特征方程的数值解在O(1)的时间内算出任意一项的值。
但是如果要得到精确的精度，那么乘法本身的时间复杂度就要增加。
但是模范3楼中矩阵方法，可以转化为一个计算k阶方阵n次幂问题，这个可以通过O(log(n))次乘法解决

[medie2005]

当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)的通项表达式。

[njwangchuan]

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表示取最接近整数，公式须证明

[mathe]

证明|((1-sqrt(5))/2)^(n+1)|/sqrt(5)<1/2就可以了
而由于|(1-sqrt(5))/2|<1,所以只要证明对于n=0成立就可以了，而这个结论是显然的。

[medie2005]

那是因为(1-√5)/2)^(n+1)这个数很小。
证明的话，只需要对开始几个小的n值验证一下就可以了。

[zgg___]

这个问题的答案是肯定的，而且过程也不大复杂的。我来把大家的想法总结一下吧，呵呵。

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

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、上面两个一综合，就可以了吧。

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

[mathe]

对于这个问题，分析计算复杂度的一个大问题是计算需要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))

[rover___]

楼上各位分析得透彻。赞一个。

[iceheart]

谢谢楼上各位的指教，这个问题我基本了解了。
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库中的用的是协处理器指令。

[iceheart]

这个方法所依赖的就是最快的求出 f(n)/f(n-1)，我用的是取两项中值的办法。
不知道各位还有没有更好的方法？

[medie2005]

我觉得lz方向搞错了：
你应该求平方根的快速算法，用以快速计算Fibonacci数列；而不是由Fibonacci数列来开发快速平方根算法。
有点本末倒置了。

[medie2005]

当然，要是lz是弄着玩的，就另说了。

另外：这几句：“形如 f(n) = a * f(n-1) + b * f(n-2) + c * f(n-3)...的数列都是等比数列。 n趋于无穷大时 f(n) / f(n-1)趋于固定值。 ”
貌似不对哦。因为有可能存在若干个大于1的特征根。

[iceheart]

多谢楼上关注。
我只是做了一些测试，通过观察到的结果产生这样的想法的，只是猜想，呵呵。
不过对于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的方法取估计值，不过效果并不理想，反而效率更低了。我想可能每次循环取中值就已经达到了比这个方法更好的效果。也许这个结果还可以用到其他优秀的开平方算法上。