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

Answer 3

The Double data type has an upper value limit of 1.7 x 10^308

The calculation for 1474 includes a value of ~ 1.1 x 10^308 at one step.

So, by 1475 you're definitely exceeding what Double can represent. Unfortunately, C#'s only larger primitive, Decimal (a 128-bit number) is designed to have very high precession but with relatively small range (only up to around 10^28).

Without designing a custom data type that can handle numbers larger than 10^308 with some degree of decimal precision, I don't see a way to do this. That said, someone out there may have already made such a class, as I can imagine scenarios where it could come in handy.

See double: http://msdn.microsoft.com/en-us/library/678hzkk9(v=VS.80).aspx

and decimal: http://msdn.microsoft.com/en-us/library/364x0z75(v=VS.80).aspx

Answer 4

The 'Solver Foundation' library appears to include some 'big' number types. It's possible that its Rational type might provide you with the precision and range that you need. It represents a rational as a ratio of two BigInteger values. (It brings its own BigInteger - I guess it was written before .NET 4 shipped.)

In theory, that enables it to represent very large numbers, but also to represent a great deal of precision. (Obviously your formula doesn't deal with rationals, but then floating point is also an approximation here.)

It offers a method for raising a Rational to the power of something else: http://msdn.microsoft.com/en-us/library/microsoft.solverfoundation.common.rational.power(v=VS.93).aspx

Answer 5

As you correctly pointed out, there is no Sqrt method implemented for BigInteger. You can implement it yourself though:

Calculate square root of a BigInteger (System.Numerics.BigInteger)

The precision should still be a problem with your code though.

Answer 6

Many answers here suggest that the complexity could be minimized to O(log(n)). Why not try an integer implementation of the log(n) approach?

First, consider that you are given two terms from Fibonacchi's sequence: F(n) and F(n+1). It's logic that the larger terms F(n+k) can be written as a linear function of F(n) and F(n+1) as

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

You could just compute these coefficients (that will depend only on k), (interestingly, they are a Fibonacci sequence too!) and use them to advance faster, then compute them again for larger values of k to be able to advance even faster, and so on.

Answer 7

The quickest (and dirtiest) ? :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));
  }
}

Answer 8

the problem is that (5^(1/2)^1475) easily overflows int. What you have to do is write a 'big math' library to handle doing math from memory (bit for bit) rather than using the hard typed data types. its a pain in the but, i know. Look up the square and multiply method.

Answer 9

1. That's not a precise formula, it will give you only estimated value. And, since floating-point arithmetic is limited to 6-8 bytes per number, deviation will raise with bigger numbers.
2. Why not to use big integer addition in cycle, it should work ok. Far better than floating-point.