nth fibonacci number in sublinear time
Following from Pillsy's reference to matrix exponentiation, such that for the matrix
M = [1 1]
[1 0]
then
fib(n) = Mn1,2
Raising matrices to powers using repeated multiplication is not very efficient.
Two approaches to matrix exponentiation are divide and conquer which yields Mn in O(ln n) steps, or eigenvalue decomposition which is constant time, but may introduce errors due to limited floating point precision.
If you want an exact value greater than the precision of your floating point implementation, you have to use the O ( ln n ) approach based on this relation:
Mn = (Mn/2)2 if n even = M·Mn-1 if n is odd
The eigenvalue decomposition on M finds two matrices U and Λ such that Λ is diagonal and
M = U Λ U-1Mn = ( U Λ U-1) n
= U Λ U-1U Λ U-1U Λ U-1 ... n times
= U Λ Λ Λ ... U-1
= U Λ nU-1
Raising a the diagonal matrix Λ to the nth power is a simple matter of raising each element in Λ to the nth, so this gives an O(1) method of raising M to the nth power. However, the values in Λ are not likely to be integers, so some error will occur.
Defining Λ for our 2x2 matrix as
Λ = [ λ1 0 ] = [ 0 λ2 ]
To find each λ, we solve
|M - λI| = 0
which gives
|M - λI| = -λ ( 1 - λ ) - 1 λ² - λ - 1 = 0
using the quadratic formula
λ = ( -b ± √ ( b² - 4ac ) ) / 2a
= ( 1 ± √5 ) / 2
{ λ1, λ2 } = { Φ, 1-Φ } where Φ = ( 1 + √5 ) / 2
If you've read Jason's answer, you can see where this is going to go.
Solving for the eigenvectors X1 and X2:
if X1 = [ X1,1, X1,2 ] M.X1 1 = λ1X1X1,1 + X1,2 = λ1X1,1X1,1 = λ1X1,2 => X1 = [ Φ, 1 ] X2 = [ 1-Φ, 1 ]
These vectors give U:
U = [ X1,1, X2,2 ]
[ X1,1, X2,2 ]
= [ Φ, 1-Φ ]
[ 1, 1 ]
Inverting U using
A = [ a b ]
[ c d ]
=>
A-1 = ( 1 / |A| ) [ d -b ]
[ -c a ]
so U-1 is given by
U-1 = ( 1 / ( Φ - ( 1 - Φ ) ) [ 1 Φ-1 ]
[ -1 Φ ]
U-1 = ( √5 )-1 [ 1 Φ-1 ]
[ -1 Φ ]
Sanity check:
UΛU-1 = ( √5 )-1 [ Φ 1-Φ ] . [ Φ 0 ] . [ 1 Φ-1 ]
[ 1 1 ] [ 0 1-Φ ] [ -1 Φ ]
let Ψ = 1-Φ, the other eigenvalue
as Φ is a root of λ²-λ-1=0
so -ΨΦ = Φ²-Φ = 1
and Ψ+Φ = 1
UΛU-1 = ( √5 )-1 [ Φ Ψ ] . [ Φ 0 ] . [ 1 -Ψ ]
[ 1 1 ] [ 0 Ψ ] [ -1 Φ ]
= ( √5 )-1 [ Φ Ψ ] . [ Φ -ΨΦ ]
[ 1 1 ] [ -Ψ ΨΦ ]
= ( √5 )-1 [ Φ Ψ ] . [ Φ 1 ]
[ 1 1 ] [ -Ψ -1 ]
= ( √5 )-1 [ Φ²-Ψ² Φ-Ψ ]
[ Φ-Ψ 0 ]
= [ Φ+Ψ 1 ]
[ 1 0 ]
= [ 1 1 ]
[ 1 0 ]
= M
So the sanity check holds.
Now we have everything we need to calculate Mn1,2:
Mn = UΛnU-1
= ( √5 )-1 [ Φ Ψ ] . [ Φn 0 ] . [ 1 -Ψ ]
[ 1 1 ] [ 0 Ψn ] [ -1 Φ ]
= ( √5 )-1 [ Φ Ψ ] . [ Φn -ΨΦn ]
[ 1 1 ] [ -Ψn ΨnΦ ]
= ( √5 )-1 [ Φ Ψ ] . [ Φn Φn-1 ]
[ 1 1 ] [ -Ψn -Ψn-1 ] as ΨΦ = -1
= ( √5 )-1 [ Φn+1-Ψn+1 Φn-Ψn ]
[ Φn-Ψn Φn-1-Ψn-1 ]
so
fib(n) = Mn1,2
= ( Φn - (1-Φ)n ) / √5
Which agrees with the formula given elsewhere.
You can derive it from a recurrance relation, but in engineering computing and simulation calculating the eigenvalues and eigenvectors of large matrices is an important activity, as it gives stability and harmonics of systems of equations, as well as allowing raising matrices to high powers efficiently.
The nth Fibonacci number is given by
f(n) = Floor(phi^n / sqrt(5) + 1/2)
where
phi = (1 + sqrt(5)) / 2
Assuming that the primitive mathematical operations (+, -, * and /) are O(1) you can use this result to compute the nth Fibonacci number in O(log n) time (O(log n) because of the exponentiation in the formula).
In C#:
static double inverseSqrt5 = 1 / Math.Sqrt(5);
static double phi = (1 + Math.Sqrt(5)) / 2;
/* should use
const double inverseSqrt5 = 0.44721359549995793928183473374626
const double phi = 1.6180339887498948482045868343656
*/
static int Fibonacci(int n) {
return (int)Math.Floor(Math.Pow(phi, n) * inverseSqrt5 + 0.5);
}
If you want the exact number (which is a "bignum", rather than an int/float), then I'm afraid that
It's impossible!
As stated above, the formula for Fibonacci numbers is:
fib n = floor (phin/√5 + 1/2)
fib n ~= phin/√5
How many digits is fib n?
numDigits (fib n) = log (fib n) = log (phin/√5) = log phin - log √5 = n * log phi - log √5
numDigits (fib n) = n * const + const
it's O(n)
Since the requested result is of O(n), it can't be calculated in less than O(n) time.
If you only want the lower digits of the answer, then it is possible to calculate in sub-linear time using the matrix exponentiation method.
One of the exercises in SICP is about this, which has the answer described here.
In the imperative style, the program would look something like
Function Fib(count)
a ← 1
b ← 0
p ← 0
q ← 1
While count > 0 Do
If Even(count) Then
p ← p² + q²
q ← 2pq + q²
count ← count ÷ 2
Else
a ← bq + aq + ap
b ← bp + aq
count ← count - 1
End If
End While
Return b
End Function