计算斐波那契数列

Question

计算斐波那契数列

c#-4.0fibonacci

12

我收到了一个很好的非递归函数用来计算斐波那契数列。

alt text

所以我编写了一些c#代码，并能够验证所有小于1474的数字都是正确的。

问题在于尝试计算1475及以上的数字时出现了问题。我的c#数学技能还不足以找到另一种解决方法。那么，有人能用c#表达这个特定的数学函数的更好方式吗？除了传统的递归函数之外？

顺便说一句，我开始使用BigInteger作为返回类型。但当尝试将(1+Math.Sqrt(5) /2)的1475次方时，问题真正出现了。我只看到需要什么数据类型（或者机制）才能使其返回其他内容而不是Infinity。

这是一个起点。

private Double FibSequence(Int32 input) {
    Double part1 = (1 / Math.Sqrt(5));
    Double part2 = Math.Pow(((1 + Math.Sqrt(5)) / 2), input);
    Double part3 = Math.Pow(((1 - Math.Sqrt(5)) / 2), input);

    return (part1 * part2) - (part1 * part3);
}

不，这不是作业。只是一个“简单”的问题，适合慢日子。

- NotMe

3

不要被它的复杂性迷惑了。使用“无限”精度计算一个无理数的第N个幂不是O(1)。这个公式只有在需要近似值时才有意义。 - ruslik

2

嗯...现在我正在遇到困难，无法在没有递归或迭代的情况下计算根。也许我已经疯了。 - 3Dave

@David，即使是加两个整数也不是O(1)。虽然硬件（ALU）能够在固定大小的整数上以O(1)的速度执行它，但通过将其扩展到任意精度，复杂性将变为O(log(N))，其中机器字长是对数的底数。 - ruslik

@David 考虑到输出的大小将是O(n)，不可能在O(1)时间内生成它。 - ruslik

什么？输出位数与达到答案所需的迭代次数无关。直接公式具有固定数量的操作，无论请求的项索引如何。迭代解决方案将需要N个循环才能达到解决方案。固定周期= O（1）。线性依赖于输入值的周期= O（n）。除非我误解了你的陈述？ - 3Dave

显示剩余9条评论

9个回答

5

using System;
using Nat = System.Numerics.BigInteger; // needs a reference to System.Numerics

class Program
{
    static void Main()
    {
        Console.WriteLine(Fibonacci(1000));
    }

    static Nat Fibonacci(Nat n)
    {
        if (n == 0) return 0;
        Nat _, fibonacci = MatrixPower(1, 1, 1, 0, Nat.Abs(n) - 1, out _, out _, out _);
        return n < 0 && n.IsEven ? -fibonacci : fibonacci;
    }

    /// <summary>Calculates matrix power B = A^n of a 2x2 matrix.</summary>
    /// <returns>b11</returns>
    static Nat MatrixPower(
        Nat a11, Nat a12, Nat a21, Nat a22, Nat n,
        out Nat b12, out Nat b21, out Nat b22)
    {
        if (n == 0)
        {
            b12 = b21 = 0; return b22 = 1;
        }

        Nat c12, c21, c22, c11 = MatrixPower(
            a11, a12, a21, a22,
            n.IsEven ? n / 2 : n - 1,
            out c12, out c21, out c22);

        if (n.IsEven)
        {
            a11 = c11; a12 = c12; a21 = c21; a22 = c22;
        }

        b12 = c11 * a12 + c12 * a22;
        b21 = c21 * a11 + c22 * a21;
        b22 = c21 * a12 + c22 * a22;
        return c11 * a11 + c12 * a21;
    }
}

- Vladimir Reshetnikov

3

Double数据类型的上限为1.7 x 10^308。在计算1474时，其中一个步骤包含了约1.1 x 10^308的值。因此，在1475时，您肯定超出了Double可以表示的范围。不幸的是，C#中唯一更大的原始数据类型Decimal（一个128位数字）被设计为具有非常高的精度，但其范围相对较小（只能达到约10^28）。除非设计一种自定义数据类型来处理大于10^308的数字并带有一定程度的小数精度，否则我看不到完成这项任务的方法。话虽如此，可能已经有人制作了这样的类，因为我可以想象它可能会派上用场。请参见double：http://msdn.microsoft.com/en-us/library/678hzkk9(v=VS.80).aspx和decimal：http://msdn.microsoft.com/en-us/library/364x0z75(v=VS.80).aspx。

- DGH

2

'Solver Foundation'库似乎包含一些“大”数类型。它的Rational类型可能为您提供所需的精度和范围。它将有理数表示为两个BigInteger值的比率。（它自带BigInteger - 我猜它是在.NET 4发布之前编写的。）

理论上，这使它能够表示非常大的数字，但也能够表示很高的精度。（显然，您的公式不涉及有理数，但浮点数在这里也是一种近似。）

它提供了一个将Rational提高到其他指数的方法：http://msdn.microsoft.com/en-us/library/microsoft.solverfoundation.common.rational.power(v=VS.93).aspx

- Ian Griffiths

1

正如您所指出的那样，BigInteger 没有实现 Sqrt 方法。不过您可以自己实现：

计算 BigInteger 的平方根 (System.Numerics.BigInteger)

然而，精度仍然可能是您代码中的问题。

- Robert Jeppesen

1

这里有许多答案表明，复杂度可以最小化到O(log(n))。为什么不尝试一下log(n)方法的整数实现呢？

首先，考虑给定斐波那契序列中的两个项：F(n)和F(n+1)。可以说较大的项F(n+k)可以写成F(n)和F(n+1)的线性函数作为：

 F(n+k) = Ck1*F(n) + Ck2*F(n+1)

你可以计算这些系数（仅依赖于 k），（有趣的是，它们也是 Fibonacci 序列！）并使用它们来更快地前进，然后为更大的 k 再次计算它们以便能够更快地前进，以此类推。

- ruslik

1

这正是我所提出的矩阵指数解决方案所做的，但更为简洁 :) - ephemient

2

@ephemient 我讨厌数学。它甚至能使简单的思想变得晦涩难懂。 - ruslik

最终我们达成了一致。 - 3Dave

1

最快（也最骯髒）的方法？ :D

private Double dirty_math_function(Int32 input){
       Double part1 = (1 / Math.Sqrt(5));
       Double part2 = Math.Pow(((1 + Math.Sqrt(5)) / 2), input);
       Double part3 = Math.Pow(((1 - Math.Sqrt(5)) / 2), input);
       return (part1 * part2) - (part1 * part3);
 }

private Double FibSequence(Int32 input) {
  if(input < 1475)
       return dirty_math_function(input);
  else{
       return (FibSequence(input -1) + FibSequence(intput -2));
  }
}

- ykatchou

有趣的解决方案。你测试过了吗？ - NotMe

如果你尝试在大数字上操作，可能会遇到堆栈溢出错误，但不要担心，可以在Stack Overflow上找到解决方案。 - ykatchou

斐波那契数列增长呈指数级别（这一点由闭合形式解决方案很明显），因此在达到最大堆栈深度之前，该函数返回的所有内容都将是Double.PositiveInfinity。 - ephemient

0

这不是一个精确的公式，它只会给你一个估计值。而且，由于浮点运算每个数字仅限于6-8字节，随着数字变得更大，偏差也会增加。
为什么不在循环中使用大整数加法，这应该可以正常工作。比浮点数好得多。

- Nickolay Olshevsky

符号上来说，这个公式是精确的。它是斐波那契递推关系的闭合形式。 - ephemient

嗯，实际上是这样的。谢谢你指出这一点，我懒得去读维基百科。 - Nickolay Olshevsky

0

问题在于(5^(1/2)^1475)很容易溢出int。你需要编写一个“大数学”库来处理从内存（位对位）中进行数学运算，而不是使用硬类型数据类型。我知道这很麻烦。可以查找平方乘法方法。

- JERiv

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- ephemient · Accepted Answer

我认为C#没有足够的浮点精度和范围来单纯地处理这个。

如果你确实想要走这条路，你可以注意到共轭 $\Phi=\phi^{-1}=\phi-1=\frac{-1+\sqrt 5}{2}$ 小于1，因此 $-\frac{(-\Phi)^n}{\sqrt 5}$ 与四舍五入到最近的整数相同，因此你可以简化你的解决方案，找到 $\left\lfloor\frac{\phi^n}{\sqrt 5}+\frac12\right\rfloor$ 。然后使用二项式展开，这样你只需要计算 $\left\lfloor a+b\sqrt 5\right\rfloor$ ，使用适当的a和b（它们是有理数，并且可以通过BigInteger精确计算）。如果你仍然回到Double，你仍然不会比1475更进一步，但你应该能够弄清楚如何仅使用精确的整数运算来完成这部分☺。

$\frac{\phi^n}{\sqrt 5}=\frac{(1+\sqrt 5)^n}{2^n\sqrt 5}=\frac{\sum_{k=0}^n{n\choose k}\sqrt 5^k}{2^n\sqrt 5}$
$=\left(\sum_{k=0}^{\left\lfloor\frac{n-1}2\right\rfloor}\frac{{n\choose 2k+1}5^k}{2^n}\right)+\left(\sum_{k=0}^{\left\lfloor\frac n2\right\rfloor}\frac{{n\choose 2k}5^{k-1}}{2^n}\right)\sqrt 5$

还有一种计算斐波那契数列的巧妙方法，使用矩阵乘法：

$\left(\begin{matrix}1&1\1&0\end{matrix}\right)^n=\left(\begin{matrix}F_n&F_{n-1}\F_{n-1}&F_{n-2}\end{matrix}\right)$

如果你足够聪明，这可以在O(log n)的时间内完成。

我最终使用Haskell实现了这些代码。 fib1 是矩阵幂的算法，而fib2则是上面描述的闭式公式的确切整数转换。当使用GHC 7.0.3编译时，它们各自的运行时间如下，由Criterion测量：
矩阵幂运行时间闭式公式运行时间

import Control.Arrow
import Data.List
import Data.Ratio

newtype Matrix2 a = Matrix2 (a, a, a, a) deriving (Show, Eq)
instance (Num a) => Num (Matrix2 a) where
    Matrix2 (a, b, c, d) * Matrix2 (e, f, g, h) =
        Matrix2 (a*e+b*g, a*f+b*h, c*e+d*g, c*f+d*h)
    fromInteger x = let y = fromInteger x in Matrix2 (y, 0, 0, y)
fib1 n = let Matrix2 (_, x, _, _) = Matrix2 (1, 1, 1, 0) ^ n in x

binom n =
    scanl (\a (b, c)-> a*b `div` c) 1 $
    takeWhile ((/=) 0 . fst) $ iterate (pred *** succ) (n, 1)
evens (x:_:xs) = x : evens xs
evens xs = xs
odds (_:x:xs) = x : odds xs
odds _ = []
iterate' f x = x : (iterate' f $! f x)
powers b = iterate' (b *) 1
esqrt e n = x where
    (_, x):_ = dropWhile ((<=) e . abs . uncurry (-)) $ zip trials (tail trials)
    trials = iterate (\x -> (x + n / x) / 2) n
fib' n = (a, b) where
    d = 2 ^ n
    a = sum (zipWith (*) (odds $ binom n) (powers 5)) % d
    b = sum (zipWith (*) (evens $ binom n) (powers 5)) % d
fib2 n = numerator r `div` denominator r where
    (a, b) = fib' n
    l = lcm (denominator a) (denominator a)
    r = a + esqrt (1 % max 3 l) (b * b / 5) + 1 % 2