Efficient calculation of Fibonacci series

Question

I'm working on Project Euler problem 2: the one about the sum of the even Fibonacci numbers up to four million.

My code:

def Fibonacci(n):
    if n == 0:
        return 0
    elif n == 1:
        return 1
    else:
        return Fibonacci(n-1) + Fibonacci(n-2)

list1 = [x for x in range(39)]
list2 = [i for i in list1 if Fibonacci(i) % 2 == 0]

The problem's solution can be easily found by printing sum(list2). However, it is taking a lot of time to come up with the list2 I'm guessing.

Is there any way to make this faster?
Or is it okay even this way...

(the problem: By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.)

Answer 1

Since matrix exponentiation can give (n-1)-th, n-th and (n+1)-th items at once we can utilize numpy.

[Image Reference]

import numpy as np

def fib_matrix_exp(n):
    A = np.array([[1, 1], [1, 0]], dtype=np.uint64)
    # return upper triangle, which will give (n-1)-th, n-th and (n+1)-th items
    return np.linalg.matrix_power(A, n)[np.triu_indices(2)][::-1] 

def fib_series(n):
    fib_first_three = np.array([0, 1, 1])
    if n < 4:
        return fib_first_three[:n]

    fib_series = []
    for i in range(1, n+1, 3):        
        fib_nums = fib_matrix_exp(i)
        fib_series.append(fib_nums)
    return np.concatenate(fib_series)[:n]

fib_values = fib_series(39)
sum([i for i in fib_values if i%2==0])

The code can be parallelized like this:

from multiprocessing import Pool

n=1000

pool = Pool(4)
fib_values=np.concatenate(pool.map(fib_matrix_exp, range(1, n+1, 3)))[:n]

sum([i for i in fib_values if i%2==0])

Answer 2

Since Python 3.9, you can use the cache decorator. From the docs:

Returns the same as lru_cache(maxsize=None), creating a thin wrapper around a dictionary lookup for the function arguments. Because it never needs to evict old values, this is smaller and faster than lru_cache() with a size limit.

from functools import cache
@cache
def fibonacci(n):
    if (n==2) or (n==1):
        return 1
    else:
        return fibonacci(n-1) + fibonacci(n-2)

Answer 3

Here's an answer with fully runnable code and benchmark, that provides a objectively the most efficient (so far) way to calculate Fibo(N). The algorithms can be implemented in whichever language. Benchmark results below, after summary of methods

I implement the following methods:

• Binet formula (spoiler: it looks like O(1) but it's very slow)
• Matrix fast exponentiation
• Fast doubling
• Fast doubling in group of 3 trick.
• Code to calculate Fibonacci using GMP (a big math library) that is well tested, to check correctness of my code

It can be downloaded here (can only be run on Linux, oops): https://github.com/lehuyduc/fast-fibonacci/tree/main

Summary of methods and their time complexity:

Space complexity of 3 best methods are O(n). In fact, for N >= 10^7, total memory usage for calculating Fibo(N) is guaranteed to be <= N bytes for those methods

1. F[n] = F[n-1] + F[n-2] => O(n^2)

Binet's formula has time complexity O(n * log(n)^2), while other methods have are O(n * log(n)), so a better way to approximate the cost is by using the number of times FFT (Fast Fourier Transform) is called. FFT is used to multiply big numbers in computer.

2. Matrix exponentiation => 48 * log(n) FFT
3. Binet's formula: idk how to evaluate this but it's a lot slower than the next method
4. Fast doubling formula: 9 * (log(n) - 2) FFT

This formula has a few forms, but let's use:

F[2k] = F[k] * (F[k] + 2F[k - 1])
F[2k + 1] = (2F[k] + F[k-1]) * (2F[k] - F[k-1]) + 2*(-1)^k

5. Fast doubling in group of 3: 6 * (log(n) - 2) FFT

F[2k] = F[k] * (F[k] + 2F[k - 1])
F[2k + 1] = (2F[k] + F[k-1]) * (2F[k] - F[k-1]) + 2*(-1)^k
F[2k - 1] = F[2k + 1] - F[2k]
F[2k + 2] = F[2k] + F[2k + 1]

6. Fast doubling in group of 3, slightly different formula: 4 * (log(n) - 2) FFT

F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k
F[2k-1] =   F[k]^2 + F[k-1]^2
F[2k] = F[2k+1] - F[2k-1]

Summary of empirical measurement:

I measure time to calculate Fibo(N) with N = 10^7, 10^8, 10^9, unit is millisecond

N = 10^7

dp cost = 27.7669 (almost similar method 3, a bit slower)
mpz_fib_ui cost = 16.284 (GMP math library)
dp no-recursion cost = 16.274 (method 6. It has the same speed as mpz_fib_ui, but for larger numbers it uses 2 threads, making it faster)
binet cost = 657.771 (Binet's formula)
matrix cost = 232.092 (Matrix fast exponentiation)
cost to convert number to base10 string = 105.358

N = 10^8

dp cost = 397.618
mpz_fib_ui cost = 268.026
dp no-recursion cost = 239.268
binet cost = 12968.6
matrix cost = 3609.47
cost to convert number to base10 string = 1904.29
Output string cost = 124.072

N = 10^9

dp cost = 5390.61
mpz_fib_ui cost = 3801.59
dp no-recursion cost = 2980.57
binet cost = 183846
matrix cost = 50054.3
cost to convert number to base10 string = 30636.7

For the 3 fast methods, the actual result does not fit the cost estimation super well, that's because:

a) We just count number of FFTs, meaning we treat all FFT to cost the same. But in reality, the last few FFTs take up most of the time. The actual formula is something like:

(N/2)*log(N / 2) + 2 * (N/4) * log(N / 4) + 2 * (N / 8) * log(N / 8) + ... 
which is O(n * log(n))

The first 2 terms (N/2)*log(N / 2) + 2 * (N/4) * log(N / 4) is the same for all methods, and cost the most.

b) We're ignoring the cost of a few other computations

But the main point is to show the reason why those methods are faster, so that estimation is good enough.

---

Actual code (need to install libgmp or just use full_run.sh script in the repo)

#include <iostream>
#include <gmpxx.h>
#include <unordered_map>
#include <string>
#include <sstream>
#include <chrono>
#include <map>
#include <fstream>
#include <thread>
#include <mutex>
#include <memory>
#include <vector>
#include "gmp-impl.h"
#include "longlong.h"
using namespace std;

// Important functions:
// gmp_fibo
// F(int n), dp_fibo
// best_fibo
// binet_fibo
// matrix_fibo

class MyTimer {
    std::chrono::time_point<std::chrono::system_clock> start;

public:
    void startCounter() {
        start = std::chrono::system_clock::now();        
    }

    int64_t getCounterNs() {
        return std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::system_clock::now() - start).count();        
    }

    int64_t getCounterMs() {
        return std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now() - start).count();        
    }

    double getCounterMsPrecise() {
        return std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::system_clock::now() - start).count()
                / 1000000.0;        
    }
};

bool ktt = false;

mpz_class gmp_fibo(int n) {
    mpz_class fib_n, fib_n1;
    mpz_fib_ui(fib_n.get_mpz_t(), n);
    return fib_n;
}

unordered_map<int, mpz_class> dp;
mpz_class& F(int n) {
    if (n <= 2) return dp[n];
    auto it = dp.find(n);
    if (it != dp.end()) return it->second;

    int k = n / 2;
    auto Fk = F(k);
    auto Fk1 = F(k - 1);    

    if (n % 2 == 0) {           
        return dp[n] = Fk * (Fk + 2 * Fk1);
    }

    int sign = (k % 2 == 0) ? 1 : - 1;
    return dp[n] = (2 * Fk + Fk1) * (2 * Fk - Fk1) + 2 * sign;
}

mpz_class dp_fibo(int n) {
    return F(n);
}

void list_dependency(map<int,int>& mp, int n) {    
    if (n <= 1 || mp[n]) return;
    mp[n] = 1;
    list_dependency(mp, n / 2);
    list_dependency(mp, n / 2 - 1);
}

void mpz_class_reserve_bits(mpz_class& x, mp_bitcnt_t bitcount)
{
    // mpz_realloc2 is a C function, so we must call it on x.get_mpz_t().
    mpz_realloc2(x.get_mpz_t(), bitcount);
}

void mpz_rsblsh2(mpz_class& c, mpz_class& a, mpz_class& b) {
    // Access the internal representation of a, b, and c
    auto a_mpz = a.get_mpz_t();
    auto b_mpz = b.get_mpz_t();
    auto c_mpz = c.get_mpz_t();

    // Get the limb data and sizes    
    mp_size_t a_size = abs(a_mpz->_mp_size); // Number of limbs in a
    mp_size_t b_size = abs(b_mpz->_mp_size); // Number of limbs in b

    // Compute the maximum size needed
    // Idk why it has to realloc every time. If I alloc one big block at the start, its size is reset to a small size anyway
    mp_size_t max_size = std::max(a_size, b_size) + 2; // +2 for carry/borrow
    mpz_realloc2(c_mpz, (max_size) * GMP_LIMB_BITS);
    mpz_realloc2(a_mpz, (max_size) * GMP_LIMB_BITS); 
    mpz_realloc2(b_mpz, (max_size) * GMP_LIMB_BITS);

    mp_ptr a_limbs = a_mpz->_mp_d;
    mp_ptr b_limbs = b_mpz->_mp_d;
    for (mp_size_t i = a_size; i < max_size; ++i) {
        a_mpz->_mp_d[i] = 0;
    }

    for (mp_size_t i = b_size; i < max_size; ++i) {
        b_mpz->_mp_d[i] = 0;
    }

    // Ensure c's limb data is clean
    std::fill(c_mpz->_mp_d + a_size, c_mpz->_mp_d + max_size, 0);
    
    mp_limb_t carry = mpn_rsblsh2_n(c_mpz->_mp_d, b_limbs, a_limbs, a_size);

    // Determine the number of significant limbs
    mp_size_t result_size = max_size;
    while (result_size > 0 && c_mpz->_mp_d[result_size - 1] == 0) {
        result_size--; // Trim trailing zeros
    }

    // Handle carry propagation correctly
    if (carry > 0) {
        c_mpz->_mp_d[result_size] = carry; // Add carry as a new highest limb
        result_size++;
    }

    // Update the size of c
    c_mpz->_mp_size = result_size;

    // Final sanity check for size mismatches
    if (result_size == 0 || result_size > max_size) {
        throw std::logic_error("Unexpected result size mismatch");
    }
}


mpz_class best_fibo(int n) {
    if (n <= 100) return gmp_fibo(n);
    map<int,int> mp;
    list_dependency(mp, n);
    
    mpz_class f[3], dummy;
    bool started = false;
    bool flag = 0;
    map<int, mpz_class*> temps;

    for (int i = 0; i < 3; i++) {
        f[i] = 0;
        mpz_class_reserve_bits(f[i], 2 * n + 64);
    }
    dummy = 0;
    mpz_class_reserve_bits(dummy, 2 * n + 64);

    MyTimer timer;
    timer.startCounter();
    // for (auto [key, value] : mp) cout << key << "\n";
    // cout << std::endl;
    vector<thread> threads;

    for (auto &[key, value] : mp)
        if (key >= 20 && mp.count(key - 1) && !mp.count(key - 2))
    {
        int N = key;
        // cout << "key = " << N << std::endl;

        if (!started) {
            f[0] = gmp_fibo(N - 1);
            f[1] = gmp_fibo(N);
            f[2] = f[0] + f[1];
            temps[N - 1] = &f[0];
            temps[N] = &f[1];
            temps[N + 1] = &f[2];
            started = true;
            continue;
        }

        // edge cases: 160, 170
        // 2 3 4 5, 8 9 10, 18 19 20, 38 39 40, 79 80, 160
        // 2 3 4 5, 8 9 10, 19 20 21, 41 42, 84 85, 170
        // 13 14 15, 29 30 31, 61 62 63, 125 126 252
        // 13 14 15, 29 30 31, 61 62 63, 125 126 253
        // 13 14 15, 29 30 31, 60 61 62, 123 124 125, 248 249 250, 499 500, 1000
        // F[2k] = F[k]*(F[k]+2F[k-1])
        // F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k
        // OR
        // F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k
        // F[2k-1] =   F[k]^2 + F[k-1]^2
        // F[2k] = F[2k+1] - F[2k-1]

        int k = N / 2;
        auto &Fk = *temps[k];
        auto &Fk1 = *temps[k - 1];                
        int sign = (k % 2 == 0) ? 1 : -1;

        if (N % 2 == 1) {
            // in this case, previous F[k+1] is unused. We that to store temporary result
            auto& Fkb = *temps[k + 1];            
            // Use f[k + 1] to store F[n - 1], f[k] = F[n], F[k - 1] = F[n + 1]

            if (n >= 50'000'000) {
                threads.clear();
                threads.emplace_back([&]() {
                    Fk *= Fk;
                });
                Fk1 *= Fk1;     
                threads[0].join();
            } else {
                Fk *= Fk;
                Fk1 *= Fk1;
            }
            
            //Fkb = (4 * Fk + 2 * sign) - Fk1; // Fkb = F[2 * k + 1] = F[N]
            mpz_rsblsh2(Fkb, Fk, Fk1);            
            Fkb += 2 * sign;
            Fk1 += Fk; // Fk1 = F[2 * k - 1] = F[n - 2]
            Fk = Fkb - Fk1; // F[k] = F[2 * k] = F[n - 1]            

            if (mp.count(N + 1)) Fk1 = Fkb + Fk;            
            
            temps.clear();
            temps[N - 1] = &Fk;
            temps[N] = &Fkb;
            temps[N + 1] = &Fk1;
        } else {
            // in this case, F[k - 2] is unused. Use it to store F[n - 1]
            auto& Fk2 = *temps[k - 2];            

            if (n >= 50'000'000) {
                threads.clear();
                threads.emplace_back([&]() {
                    Fk *= Fk;
                });
                Fk1 *= Fk1;     
                threads[0].join();
            } else {
                Fk *= Fk;
                Fk1 *= Fk1;
            }
    
            // Fk2 = (4 * Fk + 2 * sign) - Fk1; // Fk2 = F[2k + 1] = F[N + 1] 
            mpz_rsblsh2(Fk2, Fk, Fk1);            
            Fk2 += 2 * sign;
            Fk1 += Fk; // Fk1 = F[2k -1]; F[2k - 1] = F[N - 1]
            Fk = Fk2 - Fk1; // Fk = F[2k] = F[N]


            temps[N - 1] = &Fk1;    
            temps[N] = &Fk;
            temps[N + 1] = &Fk2;
        }
    }

    int k = n / 2;
    auto& Fk = *temps[k];
    auto& Fk1 = *temps[k - 1];
    int sign = (k % 2 == 0) ? 1 : -1;

    if (n % 2 == 0) {
        //return Fk * (Fk + 2 * Fk1);        
        Fk1 *= 2;
        Fk *= (Fk + Fk1);
        return std::move(Fk);
    }

    //return (2 * Fk + Fk1) * (2 * Fk - Fk1) + 2 * sign;
    Fk *= 2;
    Fk1 = (Fk + Fk1) * (Fk - Fk1) + 2 * sign;
    return std::move(Fk1);
}

mpz_class binet_fibo(int n)
{
    // Increase default precision so we don't lose accuracy for large n.    
    
    mpf_set_default_prec(n + 64);

    // Use mpf_class for floating-point operations
    mpf_class sqrt5(5);
    sqrt5 = sqrt(sqrt5); // sqrt(5)

    mpf_class phi = (mpf_class(1) + sqrt5) / 2; // (1 + sqrt(5)) / 2

    // power = phi^n
    mpf_class power;//(0, 2 * n + 32);
    mpf_pow_ui(power.get_mpf_t(), phi.get_mpf_t(), n);

    // result_float = power / sqrt(5) + 0.5
    mpf_class result_float = power / sqrt5;
    result_float += 0.5;

    // Convert the floating-point approximation to an integer (mpz_class)
    mpz_class result = mpz_class(result_float);

    return result;
}

struct Matrix {
    int n;
    vector<mpz_class> data;

    Matrix(int n) {
        this->n = n;
        data.resize(n * n);
    }

    Matrix(int n, const vector<int> &a) {
        this->n = n;
        data.resize(n * n);
        for (int i = 0; i < n * n; i++) data[i] = a[i];        
    }

    mpz_class& get(int row, int col) {
        return data[row * n + col];
    }

    const mpz_class& get(int row, int col) const {
        return data[row * n + col];
    }


    Matrix& operator*=(const Matrix& other) {
        Matrix dummy(n);

        for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++) {
            dummy.get(i, j) = 0;
            for (int k = 0; k < n; k++)
                dummy.get(i, j) += get(i, k) * other.get(k, j);
        }

        for (int i = 0; i < n * n; i++) data[i] = std::move(dummy.data[i]);        
        return *this;
    }
};

mpz_class matrix_fibo(int n) {
    if (n <= 2) return n;
    n--;
    Matrix pow = Matrix(2, {1, 1, 1, 0});
    Matrix res = Matrix(2, {1, 0, 0, 1});
    
    while (n > 0) {
        if (n & 1) res *= pow;
        pow *= pow;
        n >>= 1;
    }
    return res.data[0];
}

//-----------------------

// sum of all even fibonacci number <= Fib(n)
mpz_class sum_even_fib(int n) {
    if (n == 0) return 0;
    n -= n % 3;
    n /= 3;
    return (best_fibo(3 * n + 2) - 1) / 2;
}

//-----------------------

bool test_even(int N)
{
    mpz_class tmp, tmp1 = 1, tmp2 = 0;
    mpz_class sum = 0;
    for (int n = 2; n <= N; n++) {
        tmp = tmp1 + tmp2;
        tmp2 = tmp1;
        tmp1 = tmp;
        if (tmp % 2 == 0) sum += tmp;
        if (sum != sum_even_fib(n)) {
            cout << "test_even wrong at " << n << "\n";
            cout << sum << "\n" << sum_even_fib(n) << "\n";
            exit(0);
        }
    }
    cout << "test_even correct\n";
    return true;
}

bool test(int L, int R)
{
    for (int n = L; n <= R; n++) {
        cout << "n = " << n << "\n";
        auto res1 = gmp_fibo(n);
        auto res2 = best_fibo(n);
        string s1 = res1.get_str();
        string s2 = res2.get_str();
        
        if (s1.length() != s2.length()) cout << "Wrong length\n";

        if (s1 != s2) {
            cout << s1 << " " << s2 << "\n";
            cout << "Fail at n = " << n << "\n";
            return false;
        }        
    }

    cout << "Pass all\n";
    return true;
}

bool test(int n) {
    MyTimer timer;    

    timer.startCounter();
    auto res2 = dp_fibo(n);
    double cost2 = timer.getCounterMsPrecise();
    cout << "dp cost = " << cost2 << std::endl;

    timer.startCounter();
    auto res1 = gmp_fibo(n);
    double cost1 = timer.getCounterMsPrecise();    
    cout << "mpz_fib_ui cost = " << cost1 << std::endl;


    timer.startCounter();
    auto res3 = best_fibo(n);
    double cost3 = timer.getCounterMsPrecise();    
    cout << "dp no-recursion cost = " << cost3 << std::endl;    

    timer.startCounter();
    //auto res4 = res1;
    auto res4 = binet_fibo(n);
    double cost4 = timer.getCounterMsPrecise();
    cout << "binet cost = " << cost4 << std::endl;

    timer.startCounter();
    //auto res5 = res1;
    auto res5 = matrix_fibo(n);
    double cost5 = timer.getCounterMsPrecise();
    cout << "matrix cost = " << cost5 << std::endl;

    timer.startCounter();
    string s1 = res1.get_str();
    cout << "cost to convert number to base10 string = " << timer.getCounterMsPrecise() << std::endl;
    string s2 = res2.get_str();
    string s3 = res3.get_str();
    string s4 = res4.get_str();
    string s5 = res5.get_str();
    
    bool ok = true;
    if (s2 != s1) {cout << "DP wrong answer\n"; ok = false;};
    if (s3 != s1) {cout << "Non-recursive DP wrong answer\n"; ok = false;}
    if (s4 != s1) {cout << "Binet wrong answer\n"; ok = false;}
    if (s5 != s1) {cout << "Matrix wrong answer\n"; ok = false;}
    
    timer.startCounter();
    ofstream fo("output.txt");
    fo << s1 << "\n";
    fo.close();
    cout << "Output string cost = " << timer.getCounterMsPrecise() << "\n";

    return ok;
}

int main(int argc, char* argv[])
{
    std::ios_base::sync_with_stdio(false);
    cin.tie(0);

    int L = (argc > 1) ? atoi(argv[1]) : 10000000;
    int R = (argc > 2) ? atoi(argv[2]) : L;
    R = max(L, R);

    dp[0] = 0;
    dp[1] = 1;
    dp[2] = 1;
    dp[3] = 2;

    // test_even(100'000);    
    if (L < 0) {
        L = abs(L);
        MyTimer timer;
        timer.startCounter();
        auto res = sum_even_fib(L);
        cout << "Time to compute sum of Fibonacci numbers <= Fibo(" << L << ") with even values: " << timer.getCounterMsPrecise() << "ms" << std::endl;

        timer.startCounter();
        ofstream fo("output.txt");
        string s = res.get_str();
        fo << s << "\n";
        fo.close();
        cout << "Time to output result as string to output.txt " << timer.getCounterMsPrecise() << "ms\n";

        if (s.length() <= 100) cout << "Value = " << s << "\n";
        return 0;
    }


    bool result;
    if (L == R) result = test(L);
    else result = test(L, R);    

    if (result) cout << "Correct\n";
    else cout << "Wrong\n";
    return 0;
}

---

METHODS TO CALCULATE FIBONACCI

Some prior information:

• To add 2 numbers with N digits, it takes O(N). Multiplying a big number with a small number (int32 or int64), also takes O(N)
• Use log(n) to mean log2(n)
• To multiply 2 numbers with N digits, it takes N*log(N) using FFT or some other algorithms. It requires FFT 3 times: once for each number in forward direction, and once in reverse direction to get the result.
• The n-th Fibonacci number has around n / 3 digits (apply log10 to Binet's Formula). So the cost to multiply 2 Fibonacci number close together is around O(n / 3 * log(n / 3)) = O(n * log(n))
• We focus on very large n. For smaller numbers, addition/multiplication can use int64 which happens in O(1), instead of O(number_of_digits)
• The slowest part of calculating Fibonacci number is multiplication, and its slowest part is performing FFT. So, we can evaluate the speed of an algorithm by the number of times FFT is called.

1. Naive method: F[n] = F[n-1] + F[n-2]

• This takes n additions, and addition i takes O(i / 3) steps (because Fib(i) has length around i / 3

Total cost: 1 / 3 + 2 / 3 + 3 / 3 + ... + n / 3 = (1 + 2 + 3 ... + n) / 3 = O(n^2 / 3) = O(n^2)

2. Matrix exponentiation: The code to do fast power of a matrix is something like:

func fibo(int n) {
  Matrix pow = [[1, 1], [1, 0]];
  Matrix res = [[1, 0], [0, 1]]; // identity matrix
  while (n > 0) {
    if (n % 2 == 1) res = res * pow;
    pow = pow * pow;
    n /= 2;
  }
  return res;
}

pow = pow * pow is called called log(n) times.

if (n % 2 == 0) res = res * pow; is called number_of_bit_1_in(n) times, can treat it as log(n).

Let's count the actual number of multiplication needed to multiply two 2x2 matrices:

res(1, 1) = A(1,1) * B(1,1) + A(1,2) * B(2,1)
res(1, 2) = A(2,1) * B(1,1) + A(2,2) * B(2,1)
...

In total, it takes 8 multiplication.

So the number of times FFT called is around: 2 * log(n) * 8 * 3 = 48 * log(n)

2 * log(n): number of times res and pow are multiplied

8: number of multiplication per matrix multiply

3: number of FFT per multiplication So, the actual time complexity for this method is around O(n * log(n)).

We'll see later that this is the best possible time complexity to calculate Fibonacci number. But the difference in constant factors matters a lot in the real world.

3. Binet's formula: Though it looks O(1), it requires multiplication of irrational integer with around N bits. Cost per multiplication is still n * log(n)) Number of multiplication to do fast exponentiation: log(n) => O(n * log(n)^2)

In practice, this is much slower than the fast doubling method shown below.

4. Fast doubling formula:

F[2k] = F[k] * (F[k] + 2F[k - 1])
F[2k + 1] = (2F[k] + F[k-1]) * (2F[k] - F[k-1]) + 2*(-1)^k

The code looks something like:

unordered_map<int, mpz_class> dp;
dp[1] = 1;
dp[2] = 1;
mpz_class& F(int n) {
    if (n <= 2) return dp[n];
    auto it = dp.find(n);
    if (it != dp.end()) return it->second;

    int k = n / 2;
    auto Fk = F(k);
    auto Fk1 = F(k - 1);    

    if (n % 2 == 0) {           
        return dp[n] = Fk * (Fk + 2 * Fk1);
    }

    int sign = (k % 2 == 0) ? 1 : - 1;
    return dp[n] = (2 * Fk + Fk1) * (2 * Fk - Fk1) + 2 * sign;
}

=> Cost = number_of_times_F(k)_is_called * 3 FFT.

A good approximation for number_of_times_F(k)_is_called is is 3 * (log(n) - 2)

=> Cost = 9 * (log(n) - 2) FFT, which is a lot better than 48 * log(n) of matrix exponentiation

Let's look at a few examples to see why we can do this approximation:

F(253): 3 5 6 7 13 14 15 29 30 31 61 62 63 125 126 253
F(160): 3 4 5 8 9 10 18 19 20 38 39 40 79 80 160

The sequence can be split into 3 parts:

a) <= 2 leftover number of the beginning (in case of F(253), it's a single 3). They can be pre-computed in O(1) and doesn't need FFT

b) Groups of 3 consecutive numbers

c) (n / 2 - 1), (n / 2), and n

In each iteration, the center number in a group of 3 is doubled. And the center number always start from >= 4, so the number of groups is log(n) - 2 => F(k) is called around 3 * (log(n) - 2) times.

=> Total cost = 9 * (log(n) - 2) FFT

5. Final optimization:

Each group of 3 can be represented by: F(k-1), F(k), F(k+1)

Instead of using the fast doubling formula for all 3 numbers, we just need to calculate F(k-1), F(k), then F(k+1) = F(k) + F(k-1)

=> 2 multiplication = 6 FFT per group instead of 9 FFT per group

=> Total cost = 6 * (log(n) - 2) FFT

6. Final final optimization:

In method (4) and (5), we used the formula:

F[2k] = F[k] * (F[k] + 2F[k - 1])
F[2k + 1] = (2F[k] + F[k-1]) * (2F[k] - F[k-1]) + 2*(-1)^k

However, we can also transform the formula to:

F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k
F[2k-1] =   F[k]^2 + F[k-1]^2
F[2k] = F[2k+1] - F[2k-1]

Notice that, in a group, we need to calculate F[k]^2 and F[k-1]^2.

Multiplying 2 different numbers require 3 FFT, but squaring a number of takes 2 FFT

=> Cost per group = 4 FFT

=> Total cost = 4 * (log(n) - 2) FFT.

Note that for each group of 3, their values only depend on the previous group. So, it's possible to just use 3 bignum variables in total to calculate Fibo(n), instead of storing all the intermediate results. See the code for details

There are more ways to make computation faster, but the real bottleneck is in converting the result to decimal string. So I stop here.

---

Fibonacci numbers with even value are the ones at position F[3k].

Sum of the first K even Fibonacci is equal to: (Fibo(3 * K + 2) - 1) / 2. Proof here: https://math.stackexchange.com/questions/323058/closed-form-for-the-sum-of-even-fibonacci-numbers

So to calculate the sum of Fibonacci numbers with even value that is <= Fibo(n), the code is:

mpz_class sum_even_fib(int n) {
    if (n == 0) return 0;
    n -= n % 3;
    n /= 3;
    return (best_fibo(3 * n + 2) - 1) / 2;
}

Basically, it has same cost as finding Fibo(n)

Answer 4

This might be helpful:

fib_dict = {}

def fib(n): 
    try:
        return fib_dict[n]
    except:
        if n<=1:
            fib_dict[n] = n
            return n
        else:
            fib_dict[n] = fib(n-1) + fib (n-2)
            return fib(n-1) + fib (n-2)
    
print fib(100)

This is much faster than the traditional way.

Answer 5

Fibonacci series; add, the number before the last number to the last number and the resulting sum is the new number in the series.

You defined a recursive function, but it would be best to use a while loop.

i = 0
x = [0,1,2]
y =[]
n = int(input('upper limit of fibonacci number is '))
while i< n:
    z= x[-1]+x[-2]
    x.append(z)
    if z%2 == 0:
        y.append(z)
    i = x[-1]
    print(i)
print('The sum of even fibbunacci num in given fibonacci number is ',sum(y)+2)

Answer 6

I figured I would add my two cents in too. I added the code below and will now explain it.

def fibonacci(n):
    a, b = 1,1
    for char in bin(n)[3:]:
        a, b = (a*a + 5*b*b)//2, a*b
        if(int(char)):
            a, b=(a + 5*b)//2, (a+b)//2
    return(b)

The essential idea is to use Binet's formula (description here) but in a way that avoids radicals. Any number in Binet's formula will be of the form (a + b sqrt(5))/2 = a 1/2 + b sqrt(5)/2 so the pair a and b in the above are multiples of 1/2 and sqrt(5)/2 respectively.

There are a couple of observations. We know from applying the Binomial theorem that only odd powers of sqrt(5) will contribute to the final formula because the even powers will cancel as per Binet's formula. That means we only need to focus on the first term in the formula and in particular the number b will be the final answer once we have performed the exponentiation.

Because I used some of the ideas of Piotr Dabkowski's answer (which Piotr's answer helped me improve it even further), but the method is not really explained, I will explain it here.

We want to know X^n where X is some object with given multiplication properties. In Piotr's case X is a matrix. In this case X is the pair (a,b) which has multiplication properties I will explain below.

First we convert n to its bitwise representation, which is the second line in the function. In python bit strings start with '0b' and, since we know the first bit is 1, we start from the second bit in the number.

If we consider any bitwise representation as

n = a_m 2^m + a_{m-1} 2^{m-1} + ... + a_0 2^0

this can be written as

n = 2(...(2(2(2a_m + a_{m-1}) +a_{m-2}) ... ) + a_0

where the as are 0 or 1, and we start from the inner-most bracket and expand outward. This means we will continually double the previous value and add a_i (1 or 0 depending on its value) as we move across the bits until we reach the last bit. As n is an exponent in X^n, doubling means X^l x X^l (where l is the current value we have in the exponent) and adding 1 means multiplying as X^l x X. However, we don't care about the exponent but the above X^l at any given time. We start with a_m=1.

Now for the multiplication properties of a,b. The first case is multiplying by 1,1 (which is 1/2 + sqrt(5)/2) and this corresponds to ''add 1'' above. We see that:

(a/2 + b sqrt(5)/2)(1/2 + sqrt(5)/2) = 1/2 (a+5*b)/2 + sqrt(5)/2 (a+b)/2  

hence if the current 'char'==1 then we perform this operation, keeping in mind that a and b are in multiples of 1/2 and sqrt(5)/2. In the case of doubling we have the result:

(a/2 + b sqrt(5)/2)(a/2 + b sqrt(5)/2) = 1/2(a^2 + 5b^2)/2 + sqrt(5)/2(ab)

which is reflected in the code once more. One point to note is the division by 2. From the results above we see that if a and b are both even or both odd, then the result for for the first equation (prior to division) is that they are both even. Moreover, after the division the new a and b are both even odd or both even. Thus the division by two is always possible as the initial conditions are a=b=1 (an induction argument).

It probably is not necessary to prove this, but for the second equation you can also show if a and b are both even (odd) then the new a and b are both even (odd).

In testing the code it seems to be comparable speed to the "faster" implementations here.

As for the answer to the question, if you note that

F(2n) = F(2n+1) - F(2n-1)

you can get telescoping cancellations

F(2n-2) = F(2n-1) - F(2n-3)

F(2n-4) = F(2n-3) - F(2n-5)

. 

.

.

F(2) = F(3) - F(1)

so that sum_i^n F(2i) = F(2n+1) - F(1) and therefore the sum since F(32) is the highest even term not exceeding four million, the sum of even Fibonacci numbers not exceeding four million is F(33)-1 = 3524577.

Edit:

It appears I misunderstood the task and it is to find the sum of even Fibonacci terms not indices. That is slightly more complicated but can be done in the above method. First note that even terms have indices divisible by 3, so all even terms have the form F(3n). Using the above approach

F(3n)   = F(3n+1) - F(3n-1) = F(3n-1) - F(3n-2) - F(3n-3)
F(3n-3) = F(3n-2) - F(3n-4) = F(3n-2) - F(3n-5) - F(3n-6)
F(3n-6) = F(3n-5) - F(3n-7) = F(3n-5) - F(3n-8) - F(3n-9)
.
.
.
F(6) = F(7) - F(5) = F(7) - F(4) - F(3)
F(3) = F(4) - F(2) = F(4) - F(2) 

Adding all these up, we see terms on the right cancel so we have

Sum_i^n F(3i) = F(3n+1) - sum_i^{n-1} F(3i) - F(2)

or that

sum_i^{n-1} F(3i) = 1/2(F(3n+1) - F(3n) - F(2)) = 1/2(F(3n-1) - F(2)) 

replacing n-1 with n this gives

sum_i^{n} F(3i) = 1/2(F(3n+2) - F(2))

which is again straightforward. Now it is F(33) which is the highest even term smaller than four million, so that the sum of the even terms up to F(33) is

(F(35) - F(2))/2= 4613732

Edit 2:

So I've made further changes, removing as much redundant multiplication as possible.

def fibonacci(n):
  if(n==0):
    return 0
  if(n==1):
    return 1
  v1, v2 = 1,1
  binary_n = bin(n)[3:]
  n_length = len(binary_n)
  for i in range(0,n_length):
    if(i==n_length-1):
      if(binary_n[i]=='0'):
        return v1*v2
      elif(binary_n[i]=='1'):
        val2 = v2*v2
        return val2 + ((v1*v1+val2)//2 + v1*v2)//2
    if(binary_n[i]=='0'):
      val2 = v2*v2
      v1, v2 =  2*val2 + (v1*v1+val2)//2, v1*v2
    elif(binary_n[i]=='1'):
      val2 = v2*v2
      val1= v1*v2
      val3 = v1*v1
      v1, v2 = val2 +2*val1 + ((val3+val2)//2+val1)//2 ,val2 + ((val3+val2)//2 + val1)//2
  return(v2)

Answer 7

A closed-form solution

It turns out that there is a nice recursive formula for the sum of even Fibonacci numbers. The nth term in the sequence of sums of even Fibonacci numbers is S_{n} = 4*S_{n-1} + S_{n-2} + 2 Proof is left to the reader, but involves proving 1) even Fibo numbers are every third one, 2) proof of the formula above with induction using the definition of Fibo numbers. Using the logic from here, we can derive a closed-form formula for this with a little effort:

S_{n} = -1/2 + (1/4 + 3*sqrt(5)/20)*(2+sqrt(5))**n + (1/4 - 3*sqrt(5)/20)*(2-sqrt(5))**n

Despite the sqrt, this is integral for integral n, so this can be conveniently computed using the handy functions from my previous answer, or using a package such as sympy to handle the roots exactly.

import sympy as sp
one = sp.sympify(1) #to force casting to sympy types
k1 = -one/2
k2 = one/4 + 3*sp.sqrt(5)/20
k3 = one/4 - 3*sp.sqrt(5)/20
r1 = one
r2 = 2 + sp.sqrt(5)
r3 = 2 - sp.sqrt(5)
def even_sum_fibo(n):
  #get the nth number in the sequence of even sums of Fibonacci numbers.  If you want the sum of Fibos up to some number m, use n = m/3 (integer division)
  return sp.simplify(k1*r1**n + k2*r2**n + k3*r3**n)

Answer 8

I implemented the algorithm mentioned in SICP exercise 1.19 as:

def fib(self, n: int) -> int:
    a,b,p,q = 1,0,0,1
    while n>0:
        if n%2 == 0:
            p,q = p*p+q*q, 2*p*q+q*q
            n//=2
        else:
            a,b = b*q + a*q + a*p, b*p + a*q
            n-=1
    return b

which is of the time complexity as O(logn)

Answer 9

I realize this question was asked 8 years ago and it's been thoroughly answered… sorry to bounce it back up to the top. But there is always more to be said. I came across this in a search to improve my own algorithm, which I'd like to share.

I'd like to offer my own take since I see this wasn't really brought up. I think my algorithm is unique amongst the contributors thus far. I make use of well known Fibonacci number equations (wikipedia) in order to scale down the index. One or two others briefly cover a basic version, but I take it a step further.

This is a recursive algorithm, but I'm able to calculate Fib(2 million) in 0.15 seconds, Fib(10 million) in under 2 seconds, and Fib(100 million) in 75 seconds. All without error. I will say this, it isn't the fastest for calculating a whole list of consecutive Fibonacci numbers; this is best for picking out individuals that are very large.

Most algorithms mentioned so far - no matter how fast they may be - struggle to get above Fib(100) without recursion depth issues. A couple of contributors have eluded to parts of my algorithm, though they have some disadvantages that mine doesn't. Not saying mines the best or anything, but I think it's quite fast and can calculate really large fibs. I think it's worth adding to the discussion.

Best of all, I don't make any use of memory. No lists, dictionaries or arrays of any kind. No caches or memoization. Not even a single persistent saved constant. No special packages imported. Just basic, plain, python with basic integer types. Ive also extended the function to compute negative fibs with negligible impact to run time.

I should warn though… I'm a mathematician, not a programmer. I have no doubts this can be improved further. And I have no idea what the Big O is.

def fib(n):
      if n<0: return int(pow(-1, (n&1)+1))*fib(-n)
      if n == 0: return 0
      if n==1 or n==2: return 1
      if n==3: return 2
      
      # n is multiple of 3
      if n%3 == 0:
            third = n//3
            fibthird = fib(third)
            return 5*pow(fibthird,3) + 3*pow(-1, third)*fibthird
            
      # even n
      if n&1==0:
            return pow(fib((n>>1) + 1),2) - pow(fib((n>>1) - 1), 2)

      # for odd n
      return ( pow(fib((n>>1)+1),2) + pow(fib(n>>1),2) )

Run the code, tell me what you think. I'd love to hear from the community. I'm impressed by it, personally, and have been running it for a while. Can't find a way in my limited (programming) knowledge to improve it though. Trying to add lists, memoization, caches, etc., either fails to improve anything, or makes runtime worse. In the rare instance I find something that improves runtime, the benefits to runtime are negligible and the costs to memory are significant, and I don't think it's a fair trade.

---

Prime testing

For added fun, I include a basic probabilistic is_prime test below that relates to Fibonacci numbers:

def is_prime_fib(n):
      # Fibonacci Probabilistic is_prime test. Compositeness deterministic.
      if n==1: return False
      if n==5: return True
      if n%5 in [1,4] and fib(n-1) % n == 0: return True
      if n%5 in [2,3] and fib(n+1) % n == 0: return True
      return False

I include this just for fun even though its off topic. Its a well-known primality test using fibonacci numbers, but unfortunately it goes unused precisely because most fibonacci calculating algorithms are slow, cause recursion error, or otherwise produce inaccuracies, thus making the test unreliable and we naturally resort to other algorithms. I think the game can be changed a bit though.

On its own, the Fibonacci primality test is probabilistic. The n=1 and n=5 cases are oddities that fail to produce correct results, but they are too obvious to worry about. Aside from that, a False is deterministic in compositeness, a True is probabilistic in primeness. A composite that passes as true by this test is a Fibonacci Pseudoprime. In conjunction with other probabilistic tests, we can achieve emergent determinism.

Answer 10

I know this is an old question but I figured I would give it a go anyways.

First, some basics. Every third Fibonacci number is even. Since F(1)+F(2)=F(3), F(4)+F(5)=F(6), etc, all the even Fibonacci numbers make up half the total sum of all Fibonacci numbers up to F(3X). We already have an easy way to find the sum of all Fibonacci numbers up to F(X). The answer is F(X+2)-1. All we have to do is divide that term by two and we have our answer.

Now a little sidetracking to how we solve Fibonacci in O(log2(X)) time. Phi is a very special number. Phi=(sqrt(5)+1)/2. Phi^2=1+Phi. In fact, Phi^X=F(X-1)+F(X)Phi. Bringing back highschool algebra, we know Phi^2X=(Phi^X)^2 = (F(X-1)+F(X)Phi)^2 = F(X-1)^2+2F(X-1)F(X)Phi+(F(X)^2)(Phi^2). We know Phi^2, so substitute and distribute. F(2X-1)+F(2X)Phi=Phi^2X=F(X-1)^2+F(X)^2+Phi(2F(X-1)F(X)+F(X)^2). Since Fibonacci numbers are integers that don't contain Phi, we now know that F(2X-1)=F(X-1)^2+F(X)^2. With the additional fact that F(2X+1)=F(X)+F(X+1), we can find F(2X)=F(X+1)^2-F(X-1)^2. Now lets code!

import math

def fibonacci(x):
    a=1 #start at F(-1)
    b=0 #start at F(0)
    c=1 #start at F(1)
    bits=int((math.log2(x)+1)//1) #number of bits in x
    for i in range(bits,-1,-1):
        #Times 2
        d=b*b+a*a
        e=c*c-a*a
        f=d+e
        a=d
        b=e
        c=f
        bit=(x//(2**i))%2
        #Plus 1
        if bit==1:
            a=b
            b=c
            c=a+b
    return b
    
def fibsum(x):
    y=x-(x%3)
    return (fibonacci(y+2)-1)//2

print(fibsum(600))

Answer 11

I had done a little research and found out about a formula called Binet's formula. This formula can calculate the nth number of the fibonacci sequence easily in O(1) time.

Here is my Java code translated to Python:

def fibonacci(n):
    five_sqrt = 5 ** 0.5

    return int(round((((1 + five_sqrt)/2) ** n)/five_sqrt))

for i in range(1, 21):
    print(fibonacci(i))

Output:

1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765

Answer 12

Just another one fast solution:

def fibonnaci(n):
    a = []  
    while n != 1: 
        a.append(n&1)
        n >>= 1
    f1 = 1
    f2 = 1
    while a:
        t = f1 * (f2 * 2 - f1)
        f2 = f2 * f2 + f1 * f1
        if a.pop() is 1:
            f1 = f2
            f2 += t
        else:
            f1 = t
    return f1

Answer 13

Here are some more formulas, from OEIS:

• F(n) = ((1+sqrt(5))^n - (1-sqrt(5))^n)/(2^n*sqrt(5))
• Alternatively, F(n) = ((1/2+sqrt(5)/2)^n - (1/2-sqrt(5)/2)^n)/sqrt(5)
• F(n) = round(phi^n/sqrt(5)); where phi is (1 + sqrt(5)) / 2
• F(n+1) = Sum_{j=0..floor(n/2)} binomial(n-j, j)

Some of these formulas have implementations in the other comments above.

Answer 14

You can use the equation with square roots to compute this if you don't use floating point arithmetic, but keep track of the coefficients some other way as you go. This gives an O(log n) arithmetic operation (as opposed to O(n log n) operations for memoization) algorithm.

def rootiply(a1,b1,a2,b2,c):
    ''' multipy a1+b1*sqrt(c) and a2+b2*sqrt(c)... return a,b'''
    return a1*a2 + b1*b2*c, a1*b2 + a2*b1

def rootipower(a,b,c,n):
    ''' raise a + b * sqrt(c) to the nth power... returns the new a,b and c of the result in the same format'''
    ar,br = 1,0
    while n != 0:
        if n%2:
            ar,br = rootiply(ar,br,a,b,c)
        a,b = rootiply(a,b,a,b,c)
        n /= 2
    return ar,br

def fib(k):
    ''' the kth fibonacci number'''
    a1,b1 = rootipower(1,1,5,k)
    a2,b2 = rootipower(1,-1,5,k)
    a = a1-a2
    b = b1-b2
    a,b = rootiply(0,1,a,b,5)
    # b should be 0!
    assert b == 0
    return a/2**k/5

if __name__ == "__main__":
    assert rootipower(1,2,3,3) == (37,30) # 1+2sqrt(3) **3 => 13 + 4sqrt(3) => 39 + 30sqrt(3)
    assert fib(10)==55

Answer 15

This is much faster than everything above

from sympy import fibonacci
%timeit fibonacci(10000)

262 ns ± 10.8 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

Answer 16

import time


def calculate_fibonacci_1(n):
    if n == 0:
        return 0
    if n == 1:
        return 1
    return calculate_fibonacci_1(n - 1) + calculate_fibonacci_1(n - 2)


def calculate_fibonacci_2(n):
    fib = [0] * n
    fib[0] = 1
    fib[1] = 1
    for i in range(2, n):
        fib[i] = fib[i - 1] + fib[i - 2]
    return fib[n-1]


def calculate_fibonacci_3(n):
    a, b = 0, 1
    for _ in range(n):
        a, b = b, a + b
    return a


def calculate_fibonacci_4(n):
    v1, v2, v3 = 1, 1, 0
    for rec in bin(n)[3:]:
        calc = v2*v2
        v1, v2, v3 = v1*v1+calc, (v1+v3)*v2, calc+v3*v3
        if rec == '1':
            v1, v2, v3 = v1+v2, v1, v2
    return v2


def calculate_fibonacci_5(n):
    if n == 0:
        return (0, 1)
    else:
        a, b = calculate_fibonacci_5(n // 2)
        c = a * (b * 2 - a)
        d = a * a + b * b
        if n % 2 == 0:
            return (c, d)
        else:
            return (d, c + d)

    n = 30

    start = time.time()
    calculate_fibonacci_1(n)
    end = time.time()
    print(end - start)

    start = time.time()
    calculate_fibonacci_2(n)
    end = time.time()
    print(end - start)

    start = time.time()
    calculate_fibonacci_3(n)
    end = time.time()
    print(end - start)

    start = time.time()
    calculate_fibonacci_4(n)
    end = time.time()
    print(end - start)

    start = time.time()
    calculate_fibonacci_5(n)
    end = time.time()
    print(end - start)

for n=30:

0.264876127243
6.19888305664e-06
8.10623168945e-06
7.15255737305e-06
4.05311584473e-06

for n=300:

>10s
3.19480895996e-05
1.78813934326e-05
7.15255737305e-06
6.19888305664e-06

for n=3000:

>10s
0.000766038894653
0.000277996063232
1.78813934326e-05
1.28746032715e-05

for n=30000:

>10s
0.0550990104675
0.0153529644012
0.000290870666504
0.000216007232666

for n=300000:

>10s
3.35211610794
0.979753017426
0.012097120285
0.00845909118652

for n=3000000:

>10s
>10s
>10s
0.466345071793
0.355515003204

for n=30000000:

>100s
>100s
>100s
16.4943139553
12.6505448818

disclaimer: codes of functions no. 4 and 5 were not written by me

Answer 17

O(1) SOLUTION

The formula is also called Binet's Formula (read more)

Basically, we can write it in python like this:

def fib(n):
    a = ((1 + (5 ** 0.5)) / 2)**int(n)
    b = ((1 - (5 ** 0.5)) / 2)**int(n)
    return round((a - b) / (5 ** 0.5))

However, Because of the relatively low value of b, we can ignore it and the function can be as simple as

def fib(n):
    return round((((1+(5**0.5))/2)**int(n))/(5**0.5))

Answer 18

Here is an Optimized Solution with the Dictionary

def Fibonacci(n):
    if n<2 : return n
    elif not n in fib_dict :
            fib_dict[n]= Fibonacci(n-1) + Fibonacci(n-2)
    return fib_dict[n]

#dictionary which store Fibonacci values with the Key
fib_dict = {}
print(Fibonacci(440))

Answer 19

Spoiler alert: don't read this if you're doing Project Euler Question 2 until you've had a crack at it yourself.

Closed-form series-summation-based approaches aside, this seems more efficient than most/all of what I've seen posted, as it only needs one rather cheap loop iteration per even Fibonacci number, so only 12 iterations to get to 4,000,000.

def sumOfEvenFibonacciNumbersUpTo(inclusiveLimit):
    even = 0
    next = 1
    sum  = 0
    while even<=inclusiveLimit:
        sum  += even
        even += next<<1
        next  = (even<<1)-next
    return sum

Answer 20

Here's a simple one without recursion and O(n)

def fibonacci(n):
    a, b = 0, 1
    for _ in range(n):
        a, b = b, a + b
    return a

Answer 21

Given the starting number and the maximum number; I think the following solution for fibonacci would be interesting. The good thing is that it doesn't include recursion - thus reducing memory burden.

# starting number is a
# largest number in the fibonacci sequence is b

def fibonacci(a,b):
    fib_series = [a, a]

    while sum(fib_series[-2:]) <=b:
        next_fib = sum(fib_series[-2:])
        fib_series.append(next_fib)

    return fib_series

print('the fibonacci series for the range %s is %s'
      %([3, 27], fibonacci(3, 27)))

the fibonacci series for the range [1, 12] is [3, 3, 6, 9, 15, 24]

Answer 22

There is an O(1) solution: https://en.wikipedia.org/wiki/Fibonacci_number#Computation_by_rounding

import math

PHI = (1 + math.sqrt(5)) / 2
SQRT5 = math.sqrt(5)


def fast_fib(n):
    if n < 0:
        raise ValueError('Fibs for negative values are not defined.')
    return round(math.pow(PHI, n) / SQRT5)

Answer 23

This is some improved version of Fibonacci where we compute Fibonacci of number only once:

dicFib = { 0:0 ,1 :1 }
iterations = 0
def fibonacci(a):
    if  (a in dicFib):      
        return dicFib[a]    
    else :
        global iterations               
        fib = fibonacci(a-2)+fibonacci(a-1)
        dicFib[a] = fib
        iterations += 1
        return fib

print ("Fibonacci of 10 is:" , fibonacci(10))
print ("Fibonacci of all numbers:" ,dicFib)
print ("iterations:" ,iterations)

# ('Fibonacci of 10 is:', 55)
# ('Fibonacci of all numbers:', {0: 0, 1: 1, 2: 1, 3: 2, 4: 3, 5: 5, 6: 8, 7: 13, 8: 21, 9: 34, 10: 55})
# ('iterations:', 9)

Here we are storing Fibonacci of each number in dictionary. So you can see it calculates only once for each iteration and for Fibonacci(10) it is only 9 times.

Answer 24

I based this on an article on Fibonacci numbers on Wikipedia. The idea is to avoid looping and recursion and simply calculate the value as needed.

Not being a math wiz, selected one of the formulas and rendered it to code and tweaked it until the values came out right.

import cmath

def getFib(n):
    #Given which fibonacci number we want, calculate its value
    lsa = (1 / cmath.sqrt(5)) * pow(((1 + cmath.sqrt(5)) / 2), n)
    rsa = (1 / cmath.sqrt(5)) * pow(((1 - cmath.sqrt(5)) / 2), n)
    fib = lsa-rsa
    #coerce to real so we can round the complex result
    fn = round(fib.real) 
    return fn 

#Demo using the function
s = ''
for m in range(0,30):
    s = s + '(' + str(m) + ')' + str(getFib(m)) + ' '

print(s)

Answer 25

Solution in R, benchmark calculates 1 to 1000th Fibonacci number series in 1.9 seconds. Would be much faster in C++ or Fortran, in fact, since writing the initial post, I wrote an equivalent function in C++ which completed in an impressive 0.0033 seconds, even python completed in 0.3 seconds.

#Calculate Fibonnaci Sequence
fib <- function(n){
  if(n <= 2)
    return(as.integer(as.logical(n)))
  k = as.integer(n/2)
  a = fib(k + 1)
  b = fib(k)
  if(n %% 2 == 1)
    return(a*a + b*b)
  return(b*(2*a - b))
}

#Function to do every fibonacci number up to nmax
doFib <- function(nmax = 25,doPrint=FALSE){
    res = sapply(0:abs(nmax),fib)
    if(doPrint)
        print(paste(res,collapse=","))
    return(res)
}

#Benchmark
system.time(doFib(1000))

#user  system elapsed 
#  1.874   0.007   1.892 

Answer 26

Based on the fact that fib(n) = fib(n-1)+fib(n-2), the straightforward solution is

def fib(n):
    if (n <=1):
        return(1)
    else:
        return(fib(n-1)+fib(n-2))

however, the problem here is that some values are calculated multiple times, and therefore it is very inefficient. The reason can be seen in this sketch:

Essentially, each recursive call to fib function has to compute all the previous fibonacci numbers for its own use. So, the most computed value will be fib(1) since it has to appear in all the leaf nodes of the tree shown by answer of @kqr. The complexity of this algorithm is the number of nodes of the tree, which is $O(2^n)$.

Now a better way is to keep track of two numbers, the current value and the previous value, so each call does not have to compute all the previous values. This is the second algorithm in the sketch, and can be implemented as follows

def fib(n):
   if (n==0):
       return(0,1)
   elif (n==1):
       return(1,1)
   else:
       a,b = fib(n-1)
       return(b,a+b)

The complexity of this algorithm is linear $O(n)$, and some examples will be

>>> fib(1)
(1, 1)
>>> fib(2)
(1, 2)
>>> fib(4)
(3, 5)
>>> fib(6)
(8, 13)

Answer 27

Python doesn't optimize tail recursion, thus most solutions presented here will fail with Error: maximum recursion depth exceeded in comparison if n is too big (and by big, I mean 1000).

The recursion limit can be increased, but it will make Python crash on stack overflow in the operating system.

Note the difference in performance between fib_memo / fib_local and fib_lru / fib_local_exc: LRU cache is a lot slower and didn't even complete, because it produces a runtime error already for n = ~500:

import functools
from time import clock
#import sys
#sys.setrecursionlimit()

@functools.lru_cache(None)
def fib_lru(n):
    if n < 2:
        return n
    return fib_lru(n-1) + fib_lru(n-2)

def fib_memo(n, computed = {0: 0, 1: 1}):
    if n not in computed:
        computed[n] = fib_memo(n-1, computed) + fib_memo(n-2, computed)
    return computed[n]

def fib_local(n):
    computed = {0: 0, 1: 1}
    def fib_inner(n):
        if n not in computed:
            computed[n] = fib_inner(n-1) + fib_inner(n-2)
        return computed[n]
    return fib_inner(n)

def fib_local_exc(n):
    computed = {0: 0, 1: 1}
    def fib_inner_x(n):
        try:
            computed[n]
        except KeyError:
            computed[n] = fib_inner_x(n-1) + fib_inner_x(n-2)
        return computed[n]

    return fib_inner_x(n)

def fib_iter(n):
    a, b = 0, 1
    for i in range(n):
        a, b = b, a + b
    return a

def benchmark(n, *args):
    print("-" * 80)
    for func in args:
        print(func.__name__)
        start = clock()
        try:
            ret = func(n)
            #print("Result:", ret)
        except RuntimeError as e:
            print("Error:", e)
        print("Time:", "{:.8f}".format(clock() - start))
        print()

benchmark(500, fib_iter, fib_memo, fib_local, fib_local_exc, fib_lru)

Results:

fib_iter
Time: 0.00008168

fib_memo
Time: 0.00048622

fib_local
Time: 0.00044645

fib_local_exc
Time: 0.00146036

fib_lru
Error: maximum recursion depth exceeded in comparison
Time: 0.00112552

The iterative solution is by far the fastest and does not corrupt the stack even for n=100k (0.162 seconds). It does not return the intermediate Fibonacci numbers indeed.

If you want to compute the nth even Fibonacci number, you could adapt the iterative approach like this:

def fib_even_iter(n):
    a, b = 0, 1
    c = 1
    while c < n:
        a, b = b, a + b
        if a % 2 == 0:
            c += 1
    return a

Or if you are interested in every even number on the way, use a generator:

def fib_even_gen(n):
    a, b = 0, 1
    c = 1
    yield a
    while c < n:
        a, b = b, a + b
        if a % 2 == 0:
            yield a
            c += 1
    return a

for i, f in enumerate(fib_even_gen(100), 1):
    print("{:3d}.  {:d}".format(i, f))

Result:

  1.  0
  2.  2
  3.  8
  4.  34
  5.  144
  6.  610
  7.  2584
  8.  10946
  9.  46368
 10.  196418
 11.  832040
 12.  3524578
 13.  14930352
 14.  63245986
 15.  267914296
 16.  1134903170
 17.  4807526976
 18.  20365011074
 19.  86267571272
 20.  365435296162
 21.  1548008755920
 22.  6557470319842
 23.  27777890035288
 24.  117669030460994
 25.  498454011879264
 26.  2111485077978050
 27.  8944394323791464
 28.  37889062373143906
 29.  160500643816367088
 30.  679891637638612258
 31.  2880067194370816120
 32.  12200160415121876738
 33.  51680708854858323072
 34.  218922995834555169026
 35.  927372692193078999176
 36.  3928413764606871165730
 37.  16641027750620563662096
 38.  70492524767089125814114
 39.  298611126818977066918552
 40.  1264937032042997393488322
 41.  5358359254990966640871840
 42.  22698374052006863956975682
 43.  96151855463018422468774568
 44.  407305795904080553832073954
 45.  1725375039079340637797070384
 46.  7308805952221443105020355490
 47.  30960598847965113057878492344
 48.  131151201344081895336534324866
 49.  555565404224292694404015791808
 50.  2353412818241252672952597492098
 51.  9969216677189303386214405760200
 52.  42230279526998466217810220532898
 53.  178890334785183168257455287891792
 54.  757791618667731139247631372100066
 55.  3210056809456107725247980776292056
 56.  13598018856492162040239554477268290
 57.  57602132235424755886206198685365216
 58.  244006547798191185585064349218729154
 59.  1033628323428189498226463595560281832
 60.  4378519841510949178490918731459856482
 61.  18547707689471986212190138521399707760
 62.  78569350599398894027251472817058687522
 63.  332825110087067562321196029789634457848
 64.  1409869790947669143312035591975596518914
 65.  5972304273877744135569338397692020533504
 66.  25299086886458645685589389182743678652930
 67.  107168651819712326877926895128666735145224
 68.  453973694165307953197296969697410619233826
 69.  1923063428480944139667114773918309212080528
 70.  8146227408089084511865756065370647467555938
 71.  34507973060837282187130139035400899082304280
 72.  146178119651438213260386312206974243796773058
 73.  619220451666590135228675387863297874269396512
 74.  2623059926317798754175087863660165740874359106
 75.  11111460156937785151929026842503960837766832936
 76.  47068900554068939361891195233676009091941690850
 77.  199387062373213542599493807777207997205533596336
 78.  844617150046923109759866426342507997914076076194
 79.  3577855662560905981638959513147239988861837901112
 80.  15156039800290547036315704478931467953361427680642
 81.  64202014863723094126901777428873111802307548623680
 82.  271964099255182923543922814194423915162591622175362
 83.  1152058411884454788302593034206568772452674037325128
 84.  4880197746793002076754294951020699004973287771475874
 85.  20672849399056463095319772838289364792345825123228624
 86.  87571595343018854458033386304178158174356588264390370
 87.  370959230771131880927453318055001997489772178180790104
 88.  1571408518427546378167846658524186148133445300987550786
 89.  6656593304481317393598839952151746590023553382130993248
 90.  28197781736352815952563206467131172508227658829511523778
 91.  119447720249892581203851665820676436622934188700177088360
 92.  505988662735923140767969869749836918999964413630219877218
 93.  2143402371193585144275731144820024112622791843221056597232
 94.  9079598147510263717870894449029933369491131786514446266146
 95.  38461794961234640015759308940939757590587318989278841661816
 96.  162926777992448823780908130212788963731840407743629812913410
 97.  690168906931029935139391829792095612517948949963798093315456
 98.  2923602405716568564338475449381171413803636207598822186175234
 99.  12384578529797304192493293627316781267732493780359086838016392
100.  52461916524905785334311649958648296484733611329035169538240802

Time: 0.00698620

That's the first 100 even Fibonacci numbers in ~7ms and includes the overhead of printing to terminal (easy to underestimate on Windows).

Answer 28

This is a very fast algorithm and it can find n-th Fibonacci number much faster than simple iterative approach presented in other answers, it is quite advanced though:

def fib(n):
    v1, v2, v3 = 1, 1, 0    # initialise a matrix [[1,1],[1,0]]
    for rec in bin(n)[3:]:  # perform fast exponentiation of the matrix (quickly raise it to the nth power)
        calc = v2*v2
        v1, v2, v3 = v1*v1+calc, (v1+v3)*v2, calc+v3*v3
        if rec=='1':    v1, v2, v3 = v1+v2, v1, v2
    return v2

You can read some more about involved math here.

---

Answer 29

To find the sum of the first n even-valued fibonacci numbers directly, put 3n + 2 in your favourite method to efficiently compute a single fibonacci number, decrement by one and divide by two (fib((3*n+2) - 1)/2)). How did math dummies survive before OEIS?

Answer 30

Haskell 1 liner :-

fibs = 0 : (f 1 1) where f a b = a : f b (a+b)

This code is extremely efficient and calculates Fibonacci numbers up to (10^1000) in less than a second ! This code will also be useful for this problem in Project Euler.