Calculating fibonacci

Question

I was sent this nice non-recursive function for computing a fibonacci sequence.

So I coded up a bit of c# and was able to verify all numbers up to 1474 were correct.

The problem comes in when attempting to calculate it for 1475 and above. My c# math skills just aren't up to the task of figuring out a different way. So, does someone have a better way of expressing this particular math function in c#? other than the traditional way of doing a recursive function?

Incidentally, I started to use BigInteger as the return type. But the problem really comes in when trying to raise (1+Math.Sqrt(5) /2) to the 1475th power. I just don't see what data type I need (nor mechanism for that matter) to get this to come back with something other than Infinity.

Here's a starting point.

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);
}

And, no, it's not homework. Just a "simple" problem for a slow day.

Answer 1

I don't think C# has a data type with enough floating precision and range to handle this naïvely.

If you really want to go down this path, you can notice that the conjugate is less than one, so does the same as rounding to the nearest integer, thus you can simplify your solution to finding . Then use binomial expansion so that you only have to calculate , with appropriate a and b (which are rational and can be computed exactly with BigInteger). If you still go back to Double for this, you still won't get much further than 1475, but you should be able to figure out how to do even this part with exact integer math only ☺

---

There's another neat method for computing Fibonacci numbers, using matrix exponentiation:

This can be done in O(log n) if you're clever.

---

I ended up implementing these in Haskell. fib1 is matrix exponentiation and fib2 is exact integer translation of the closed-form formula, as described above. Their respective runtimes look like this, as measured by Criterion when compiled by GHC 7.0.3:

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

Answer 2

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;
    }
}