回复 Devil_W
呵呵,看来不解释一下算法都对不起仁兄的置顶卡了。不知道仁兄的置顶卡哪来的...
进入正题
首先要更正一下Fibonacci数列的起始项,应该是
1, 1, 2, 3, 5, 8, 13, 21,...
现在讨论每一项的奇偶性质。
根据数列的定义,每一项是前两项的和。那么该项的奇偶性与前两项有关。
对于两个数的和有下面的性质(别嫌简单)
奇 + 奇 = 偶
奇 + 偶 = 奇
偶 + 奇 = 奇
奇 + 奇 = 偶
由此数列的奇偶序列是
奇 奇 偶 奇 奇 偶 奇 奇 偶 ...
(这也是一个有限状态自动机)
由此也证明Fibonacci数列每三项出现一个偶数。
因为每个偶数都是它前面两个奇数的和,所以某偶数项之前的偶数和就是全部和的一半。
为了讨论方便,定义几个常数
A = (1 + √5) / 2
B = (1 - √5) / 2
C = 1 / √5
这样,Fibonacci的通项公式就是
F(n) = C * (A^n - B^n)
n 是项的索引,也可以叫序号,从1开始。
换个符号更顺眼一些。
F(x) = C * (A^x - B^x)
对于前 N 项的和就是求F(x)在[1,N]的积分。。。。
开个玩笑,这个问题用不着做积分。当然,其本质就是积分。
前 x 项的和为
sum(F(x)) = sum(C * (A^x - B^x))
= C * (sum(A^x) - sum(B^x))
而sum(A^x)是一个初始值为A,比率为A的等比数列,它前x项的和为
sum(A^x) = A * (A^x - 1) / (A - 1)
同理
sum(B^x) = B * (B^x - 1) / (B - 1)
将它们以及A B C代入F(x)整理一下就得到我算法中最后那个求和公式
现在的问题就是如何根据已知的项值 A 求出它对应的索引 x
也就是求F(x)的反函数F^-1(x)
这本来是个理想的方法,但不得不承认毕业太久了,我的数学水平退步了。
我用了两个小时也没有得到一个显式的反函数,只能得到隐函数。
虽然用迭代法也可以从隐函数中求得 x,但消耗的计算量是我不能接受的,也不想为这么简单的问题动用这么复杂的算法。
于是我着手简化原函数,用一个近似函数代替原函数。
注意F(x) 中 B^x 这个幂函数。
B 是一个小于 0 大于 -1 的负数, B^x 是一个振荡函数。
为了便于讨论,换个写法
B^x = (-1)^x * (-B)^x
这样,(-1)^x 是一个符号函数(x的定义域是自然数集),0 < -B < 1
当 x 增大时 (-B)^x 减小趋于 0,而且速度很快。这一点可以有很严格的数学证明,这里我就不废话了。
基于此,我略去了原函数中这一项得到一个近似函数
Q(x) = C * A^x
它的反函数就很简单了
Q^-1(x) = ln(x / C) / ln(A)
用它我可以代入项值求出一个近似的索引值,就是算法中的第一句
由于我们由于我们已知的值也不是一个真正的项值,真正的项值可能比这个值小,但下一项的值一定比这个值大。
得到的索引值也是一个近似值。
所以算法的中间部分就是对这个近似索引值的修正,以得到精确的索引值。
其实之前发的代码是错误的,只是很偶然得到了正确的结果。现在发修正后的代码。
仁兄可以对比一下,有兴趣的话仁兄来解释一下为什么
程序代码:
int Test2(int a)
{
double x, y, sf;
int n, t;
x = log(SQRT_5 * a) / log(A);
x += 0.5;
n = (int)x;
y = (pow(A, n) - pow(B, n)) / SQRT_5;
y += 0.5;
t = (int)y;
if(t > a) n--;
n -= n % 3;
sf = ((SQRT_5 + 1) / (SQRT_5 - 1) * (pow(A, n) - 1) - (1 - SQRT_5) / (1 + SQRT_5) * (1 - pow(B, n))) / SQRT_5;
sf += 0.5;
return ((int)sf) / 2;
}
不知道分析明白没,有质疑我们还可以再讨论