Project Euler | Problem 2 - Even Fibonacci numbers
A tutorial about how to solve problem 2 from Project Euler.
Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with $1$ and $2$, the first $10$ terms will be:
$1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...$
By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.
There is a similar exercise on HackerRank:
Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with $1$ and $2$, the first $10$ terms will be:
$1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...$
By considering the terms in the Fibonacci sequence whose values do not exceed $N$, find the sum of the even-valued terms.
Input Format
First line contains $T$ that denotes the number of test cases. This is followed by $T$ lines, each containing an integer, $N$.
Constraints
$1 \leq T \leq 10^5$ and $10 \leq N \leq 4 \cdot 10^6$
The HackerRank exercise is a generic version of the Project Euler exercise, since our maximum value is variable here (through the variable $N$). In the Project Euler version, in contrast, the maximum value is fixed at $4000000$. In the following I will consider a maximum number of $4000000$ (Project Euler exercise) and $4 \cdot 10^6$ (upper bound of the HankerRank exercise). To avoid any confusion, I will call the maximum number $max \_ num$ in the following and not $N$ as HackerRank does. I will use $n$ for the $n$th fibonacci number in the fibonacci sequence and calling the maximum number $N$ at the same time might be confusing.
The fibonacci sequence is defined as follows (for $n > 1$):
$$F_0 = 0, \quad F_1 = 1, \quad F_n = F_{n-1} + F_{n-2}$$
As a result, the first fibonacci numbers are (as stated in the exercise):
$0, \quad 1, \quad 1, \quad 2, \quad 3, \quad 5, \quad 8, \quad 13, \quad 21, \quad 34, \quad 55, \quad 89, \quad 144, \quad ...$
We should only consider the even fibonacci numbers that are smaller than $max \_ num$:
$2, \quad 8, \quad 34, \quad 144, \quad ...$
Finally, we should take the sum of the even numbers smaller than $max \_ num$.
My first approach was to simply use brute force. We create all the fibonacci numbers smaller than $max \_ num$ using the formula $F_n = F_{n-1}+F_{n-2}$, look for the even ones among them and add them together. This could be implemented using the following steps:
- We create a $sum$ variable with an intial value of $0$ and set $n = 3$. We set $n$ to $3$, since $F_3$ is the first even fibonacci number ($F_3 = 2$).
- We create the $n$th fibonacci number using the formula $F_n = F_{n-1} + F_{n-2}$.
- If $F_n \leq max \_ num$, we will go to step 4. Otherwise we stop and return $sum$.
- We check if the number is even using the modulo operator (I explained how to do this in my previous post about Project Euler problem 1).
- If the number is even, we add it to $sum$.
- We increase $n$ by $1$ and jump back to step 2.
In Python we can implement the approach as follows:
def fibonacci_sum_bf(max_num):
sum = 0
fn_1 = 1
fn_2 = 1
fn = fn_1 + fn_2
while fn <= max_num:
if fn % 2 == 0:
sum += fn
fn_1, fn_2 = fn, fn_1
fn = fn_1 + fn_2
return sum
Let's calculate the sum for $max \_ num = 4000000$.
fibonacci_sum_bf(4000000)
The sum is 4613732. Now, let's also check how long the brute force version takes. I use the %timeit magic for this.
%timeit -n10000 -r100 fibonacci_sum_bf(4000000)
The parameter -n10000 means it runs the function 10000 times in a loop. The parameter -r100 means it repeats the loop of 10000 runs 100 times. Finally, it returns the mean and standard deviation of how long the function fibonacci_sum_bf needed to run. As we can see, the function took 2.83 microseconds on average. Let's also calculate the sum for $max \_ num = 4 \cdot 10^{16}$.
fibonacci_sum_bf(4 * 10**16)
The sum is 49597426547377748. Let's also check how long it takes.
%timeit -n10000 -r100 fibonacci_sum_bf(4 * 10**16)
The function took 8.02 microseconds.
I was wondering if there is a better approach and I found a blog post from Oussama Zaki describing a nice idea. Our brute force approach creates all fibonacci numbers smaller than $max \_ num$ and checks after each created number if it is even. However, it would be much better, if we could directly create only the even fibonacci numbers and add them together. This way we don't create any odd numbers that we don't need anyway, and we also don't need to check if a number is even. But how can we create only the even fibonacci numbers? Let's take a look at the fibonacci sequence again and focus on where the even numbers are:
$0, \quad 1, \quad 1, \quad $ $\bf2$ $, \quad 3, \quad 5, \quad $ $\bf8$ $, \quad 13, \quad 21, \quad $ $\bf34$ $, \quad 55, \quad 89, \quad $ $\bf144$ $, \quad ...$
We can see a pattern here:
$odd, \quad odd, \quad $ $\bf even$ $, \quad odd, \quad odd, \quad $ $\bf even$ $, \quad odd, \quad odd, \quad $ $\bf even$ $, \quad ...$
This follows from the addition laws on even and odd numbers:
$$odd + odd = even, \quad even + odd = odd$$
Now, let's assume $F_n$ is even. Then, $F_{n-1}$ and $F_{n-2}$ must be odd, $F_{n-3}$ must be even, $F_{n-4}$ and $F_{n-5}$ must be odd, and $F_{n-6}$ must be even. Let's try to express $F_n$ only using the even numbers $F_{n-3}$ and $F_{n-6}$. Originally, $F_n$ is defined as
$F_n = F_{n-1} + F_{n-2}$
We know that $F_{n-1} = F_{n-2} + F_{n-3}$ and $F_{n-2} = F_{n-3} + F_{n-4}$. So, let's rewrite the equation.
$F_n = F_{n-2} + F_{n-3} + F_{n-3} + F_{n-4} = 2 \cdot F_{n-3} + F_{n-2} + F_{n-4}$
We know that $F_{n-2} = F_{n-3} + F_{n-4}$ and $F_{n-4} = F_{n-5} + F_{n-6}$. So, let's rewrite the equation again.
$F_n = 2 \cdot F_{n-3} + F_{n-3} + F_{n-4} + F_{n-5} + F_{n-6} = 3 \cdot F_{n-3} + F_{n-4} + F_{n-5} + F_{n-6}$
We know that $F_{n-4} + F_{n-5} = F_{n-3}$. So, let's rewrite the equation one more time.
$F_n = 3 \cdot F_{n-3} + F_{n-3} + F_{n-6} = 4 \cdot F_{n-3} + F_{n-6}$
Now, we have a formula to create even fibonacci numbers. This results in a sequence that Oussama Zaki calls EvenFibonacci sequence ${EF}_n$. It is defined as follows (for $n > 1$):
$${EF}_0 = 2, \quad {EF}_1 = 8, \quad {EF}_n = 4 \cdot {EF}_{n-1} + {EF}_{n-2}$$
The fibonacci numbers $2$ and $8$ are the first two even fibonacci numbers. All other even fibonacci numbers are calculated using our derived formula from above. We can implemented the sum of the even fibonacci numbers using the following steps:
- We create a $sum$ variable with an intial value of $10$ (the sum of the first 2 even fibonacci numbers $2$ and $8$) and set $n = 3$.
- We create the $n$th even fibonacci number using the formula ${EF}_n = 4 \cdot {EF}_{n-1} + {EF}_{n-2}$.
- If ${EF}_n \leq max \_ num$, we will go to step 4. Otherwise we stop and return $sum$.
- We add ${EF}_n$ to $sum$.
- We increase $n$ by $1$ and jump back to step 2.
In Python we can implement the approach as follows:
def fibonacci_sum_impr(max_num):
sum = 10
fn_1 = 8
fn_2 = 2
fn = 4 * fn_1 + fn_2
while fn <= max_num:
sum += fn
fn_1, fn_2 = fn, fn_1
fn = 4 * fn_1 + fn_2
return sum
Let's calculate the sum for $max \_ num = 4000000$.
fibonacci_sum_impr(4000000)
The sum is 4613732. It is the same solution that we also received using our brute force approach. Now, let's check how long the improved version fibonacci_sum_impr takes. Could we speed up?
%timeit -n10000 -r100 fibonacci_sum_impr(4000000)
The function took 973 nanoseconds on average. This is faster than our brute force approach, which took 2.83 microseconds. So, we could speed up a little. Now, let's also calculate the sum for $max \_ num = 4 \cdot 10^{16}$.
fibonacci_sum_impr(4 * 10**16)
The sum is 49597426547377748. Not surprisingly, this is the same result that we already received using our brute force approach. Let's also check how long the function takes for $max \_ num = 4 \cdot 10^{16}$.
%timeit -n10000 -r100 fibonacci_sum_impr(4 * 10**16)
The function took 2.89 microseconds. This is much better than our brute force approach, which took 8.02 microseconds.
It would be nice to improve even further. Especially, a pure mathematical approach, that runs in constant time (i.e., the run time is independent of $max \_ num$), would be nice. Oussama Zaki describes such an approach in his blog post as well. I will explain it in the following.
Fibonacci numbers can also be calculated using the following closed format known as Binet's formula:
$$F_n = \frac{\varphi^n - \psi^n}{\sqrt{5}}$$
$\varphi$ is called the golden ratio. $\psi$ is called the silver ratio or golden ratio conjugate.
$$\varphi = \frac{1 + \sqrt{5}}{2}, \quad \psi = \frac{1 - \sqrt{5}}{2}$$
We want to add up all even fibonacci numbers. As we have seen before, that is every third number:
$F_3 \quad + \quad F_6 \quad + \quad F_9 \quad + \quad ... $
So, we can write the sum $S$ as follows:
$$S = \sum^k_{i=1} F_{3i}$$
Now, let's use the Binet's formula:
$$S \quad = \quad \sum^k_{i=1} F_{3i} \quad = \quad \sum^k_{i=1} \frac{\varphi^{3i} - \psi^{3i}}{\sqrt{5}}$$
Let's change the expression as follows:
$$S \quad = \quad \sum^k_{i=1} \frac{\varphi^{3i} - \psi^{3i}}{\sqrt{5}} \quad = \quad \frac{1}{\sqrt{5}} \cdot \bigg(\sum^k_{i=1} (\varphi^3)^i - \sum^k_{i=1} (\psi^3)^i\bigg)$$
Now, we can apply the formula for the geometric series to adjust our expression further. The formula for the geometric series is the following one:
$$\sum^k_{i=1} r^i \quad = \quad \frac{r^1 - r^{k+1}}{1 - r}$$
Now, we can replace $r$ with $\varphi^3$:
$$\sum^k_{i=1} (\varphi^3)^i \quad = \quad \frac{(\varphi^3)^1 - (\varphi^3)^{k+1}}{1 - \varphi^3} \quad = \quad \varphi^3 \cdot \frac{1 - (\varphi^3)^k}{1 - \varphi^3}$$
Let's do the same for $\psi^3$:
$$\sum^k_{i=1} (\psi^3)^i \quad = \quad \frac{(\psi^3)^1 - (\psi^3)^{k+1}}{1 - \psi^3} \quad = \quad \psi^3 \cdot \frac{1 - (\psi^3)^k}{1 - \psi^3}$$
Finally, we can adjust our expression above as follows:
$$S \quad = \quad \frac{1}{\sqrt{5}} \cdot \bigg(\sum^k_{i=1} (\varphi^3)^i - \sum^k_{i=1} (\psi^3)^i\bigg) \quad = \quad \frac{1}{\sqrt{5}} \cdot \bigg( \varphi^3 \cdot \frac{1 - (\varphi^3)^k}{1 - \varphi^3} - \psi^3 \cdot \frac{1 - (\psi^3)^k}{1 - \psi^3} \bigg)$$
We are almost done. There is only one piece missing. What is $k$? Well, we know that $3 \cdot k \leq n_{max}$, where $n_{max}$ is the $n$th fibonacci number that is smaller than $max \_ num$. Why is $3 \cdot k$ smaller than $n_{max}$ and not equal? Well, $n_{max}$ is not necessarily an even fibonacci number. How can we get $k$ then? We can get $k$ as follows:
$$3 \cdot k \leq n_{max} \quad \Rightarrow \quad k = \biggl\lfloor \frac{n_{max}}{3} \biggr\rfloor$$
The brackets $\lfloor ... \rfloor$ mean that we take the floor. This makes sure that $3 \cdot k$ is the next smaller even fibonacci number to $n_{max}$ (I have used the same approach in my blog post about solving Project Euler problem 1). But what is $n_{max}$? Well, we can calculate $n_{max}$ via the inverse of the closed form fibonacci function. I haven't really understood how Oussama Zaki did this, but then I found a post on Stackoverflow from PengOne that describes it quite well. First, we need to look at the fibonacci function in closed form again for $n_{max}$:
$$F_{n_{max}} = \frac{\varphi^{n_{max}} - \psi^{n_{max}}}{\sqrt{5}}$$
As already mentioned above, $\psi$ is defined as follows:
$$\psi = \frac{1 - \sqrt{5}}{2}$$
Thus, $\psi^{n_{max}}$ is:
$$\psi^{n_{max}} = \bigg(\frac{1 - \sqrt{5}}{2}\bigg)^{n_{max}} = \frac{(1 - \sqrt{5})^{n_{max}}}{2^{n_{max}}}$$
The absolute value of $1 - \sqrt{5}$ is smaller than $2$. If we take the exponential of both expressions (i.e., $(1 - \sqrt{5})^{n_{max}}$ and $2^{n_{max}}$), $2^{n_{max}}$ will always be higher than $(1 - \sqrt{5})^{n_{max}}$. For large ${n_{max}}$, the expression $2^{n_{max}}$ will be even much higher and in our case here $n_{max}$ is indeed large. As a result, $\psi^{n_{max}}$ will be close to $0$ for large ${n_{max}}$, since the denominator $2^{n_{max}}$ is much higher than the numerator $(1 - \sqrt{5})^{n_{max}}$. Hence, we can ignore $\psi^{n_{max}}$ in $F_{n_{max}}$ for large $n$ resulting in:
$$F_{n_{max}} = \frac{\varphi^{n_{max}}}{\sqrt{5}} \quad \Rightarrow \quad \varphi^{n_{max}} = \sqrt{5} \cdot F_{n_{max}}$$
Now, we can transform the formula to the inverse function:
$${n_{max}} = \log_{\varphi} \big( \sqrt{5} \cdot F_{n_{max}} \big)$$
But there is one more problem! We actually don't have $F_{n_{max}}$! However, according to the post on Stackoverflow, also any number close to a large fibonacci number $F_{n_{max}}$ gives the correct result. Now, we have all pieces together to calculate the sum $S$:
$$S \quad = \quad \frac{1}{\sqrt{5}} \cdot \bigg( \varphi^3 \cdot \frac{1 - (\varphi^3)^k}{1 - \varphi^3} - \psi^3 \cdot \frac{1 - (\psi^3)^k}{1 - \psi^3} \bigg)$$
$$k = \biggl\lfloor \frac{n_{max}}{3} \biggr\rfloor, \quad {n_{max}} = \log_{\varphi} \big( \sqrt{5} \cdot F_{n_{max}} \big)$$
Let's see how this looks in Python code:
from math import sqrt, floor, log
def fibonacci_sum_math(max_num):
# phi and psi
phi = (1 + sqrt(5)) / 2
psi = (1 - sqrt(5)) / 2
# get k through the inverse fibonacci function
def inverse_fib(fn):
return int(round(log(fn * sqrt(5), phi)))
k = inverse_fib(max_num) // 3
# compute sum
phi3 = phi ** 3
psi3 = psi ** 3
sum = int(round((1 / sqrt(5)) * (
phi3 * ((1 - phi3 ** k) / (1 - phi3)) -
psi3 * ((1 - psi3 ** k) / (1 - psi3))
)))
return sum
As you can see, I round the sum and the result of the inverse fibonacci function. This is necessary, since we probably don't get exact integer values due to floating point rounding errors.
Let's calculate the sum for $max \_ num = 4000000$.
fibonacci_sum_math(4000000)
The sum is 4613732. This is the same result that we also received using our brute force and our improved approach. Now, let's check how long the pure mathematical version takes to see if we could speed up.
%timeit -n10000 -r100 fibonacci_sum_math(4000000)
It took 1.54 microseconds. Oddly, this is faster than the brute force approach (2.83 microseconds) but a little slower than the improved approach (973 nanoseconds). My guess is that operations like log or sqrt can take up some time, and for $max \_ num = 4000000$ we actually don't need to add up so many numbers. So, let's calculate the solution for $max \_ num = 4 \cdot 10^{16}$. For $max \_ num = 4 \cdot 10^{16}$ we need to add up more numbers. Let's see if there is a difference in the run times.
%timeit -n10000 -r100 fibonacci_sum_math(4 * 10**16)
It took 1.62 microseconds. This is faster than our brute force approach (8.02 microseconds) and our improved approach (2.89 microseconds). Moreover, it is almost as fast as calculating the sum for $max \_ num = 4000000$ (1.54 microseconds). This makes sense, since our pure mathematical approach has a complexity of O(1) (we compute the sum directly, i.e., we are independent of $max \_ num$), whereas our brute force and our improved approach has the complexity of O(n) (we need to loop through the fibonacci sequence, i.e., we are dependent on $max \_ num$). Let's see if we also get the correct result.
O(1) means constant time. O(n) means linear time. This notation is called Big-O-Notation and it can be used to express the complexity of an algorithm. You can get more information about it here.
fibonacci_sum_math(4 * 10**16)
We get 49597426547377776. However, this is not the same result that we got before! Before, we got 49597426547377748. Let's try a few different values for $max \_ num$ to obtain a few more insights about the problem.
import pandas as pd
max_nums = [10, 1000, 1000000, 4000000, 10**15, 10**16, 4 * 10**16]
result = {}
for i in max_nums:
result[i] = {
'naive': fibonacci_sum_bf(i),
'improved': fibonacci_sum_impr(i),
'math': fibonacci_sum_math(i)
}
pd.DataFrame(result)
Okay, the sums differ already for $max \_ num = 10^{16}$. For the smaller values of $max \_ num$ the sum is correct. Why is that? Well, our pure mathematical approach suffers from floating point issues. Float variables in standard Python can be of IEEE 754 double precision at maximum. However, we need higher precision here. What can we do? Well, there is a Python module called decimal, which allows a higher precision (for additional information about decimal you can also check out this blog post). So, let's adjust our code by integrating the decimal module.
from math import sqrt, floor, log
from decimal import *
def fibonacci_sum_math2(max_num):
context = Context(prec=3000, rounding=ROUND_05UP)
setcontext(context)
# phi and psi
phi = Decimal(1 + Decimal(5).sqrt()) / Decimal(2)
psi = Decimal(1 - Decimal(5).sqrt()) / Decimal(2)
# get k through the inverse fibonacci function
def inverse_fib(fn):
return int(round(log(fn * sqrt(5), phi)))
k = inverse_fib(max_num) // 3
# compute sum
phi3 = context.power(phi, 3)
psi3 = context.power(psi, 3)
sum = int(round((Decimal(1) / Decimal(5).sqrt()) * (
phi3 * ((1 - context.power(phi3, k)) / (1 - phi3)) -
psi3 * ((1 - context.power(psi3, k)) / (1 - psi3))
)))
return sum
There are actually not too many changes. We just need to set the precison over a context object, and wrap float values with the Decimal class. Mathematical operations are executed either over the Decimal or the context objects. Okay, let's calculate the sum for $max \_ num = 4 \cdot 10^{16}$.
fibonacci_sum_math2(4 * 10**16)
The sum is 49597426547377748. This is the correct result finally. Let's also check the other values for $max \_ num$.
import pandas as pd
max_nums = [10, 1000, 1000000, 4000000, 10**15, 10**16, 4 * 10**16]
result = {}
for i in max_nums:
result[i] = {
'naive': fibonacci_sum_bf(i),
'improved': fibonacci_sum_impr(i),
'math': fibonacci_sum_math2(i)
}
pd.DataFrame(result)
All the numbers look good now. Next, let's also check the run time.
%timeit -n1000 -r10 fibonacci_sum_math2(4 * 10**16)
The function took 15.3 milliseconds! This is a lot more than before! This is actually the disadvantage of the decimal module. It is slow! To have a faster execution time, Oussama Zaki suggests to use a programming language that supports a high floating point precision like the Wolfram Language. However, the Wolfram Language is proprietary. In my understanding, it can only be used in the Wolfram Mathematica software system. So, I was looking for an alternative and I found out that the programming language Fortran offers a higher floating point precision as well. Well, Fortran is typically used for numerical computing. The original BLAS library is written in Fortran. BLAS is a linear algebra library, that is included in e.g., NumPy or Mathematica! So, let's transform our pure mathematical approach from Python into Fortran.
function reverse_fib(fn) result(n)
! function that calculates the reverse fibonacci number, i.e., given the fibonacci number
! it calculates which fibonacci number it is in the fibonacci sequence (e.g., the 3rd, 10th)
implicit none
! input (fibonacci number fn) and output (place of fn in fibonacci sequence)
real*16, intent(in) :: fn
real*16 :: n
! constants as long double (necessary because otherwise we get rounding errors)
real*16, parameter :: c_1 = 1
real*16, parameter :: c_2 = 2
real*16, parameter :: c_4 = 4
real*16, parameter :: c_5 = 5
! calculating phi
real*16 :: phi
phi = (c_1 + sqrt(c_5)) / c_2
! calcluate output
n = anint(log((fn * sqrt(c_5) + sqrt(c_5 * (fn**c_2) - c_4)) / c_2) / log(phi))
end function
program fib_sum
! program that calculates the sum of the even fibonacci numbers which do not exceed
! the number 4 * 10**16
implicit none
! the maximum number the fibonacci numbers are not allowed to exceed
real*16 :: max_num = 4.0 * 10.0**16
! define some variables that are needed for calculation
real*16 :: phi
real*16 :: psi
real*16 :: phi3
real*16 :: psi3
real*16 :: k
real*16 :: reverse_fib
real*16 :: temp_1
real*16 :: temp_2
real*16 :: temp_3
! declare our sum as long integer (otherwise with a normal integer we always get the
! maximum normal integer when we have big numbers, since they do not fit in)
integer, parameter :: LargeInt_K = selected_int_kind (18)
integer (kind=LargeInt_K) :: sum
! constants as long double (necessary because otherwise we get rounding errors)
real*16 :: c_1 = 1
real*16 :: c_2 = 2
real*16 :: c_3 = 3
real*16 :: c_5 = 5
! calculate the sum of the even fibonacci numbers
phi = (c_1 + sqrt(c_5)) / c_2
psi = (c_1 - sqrt(c_5)) / c_2
phi3 = phi ** c_3
psi3 = psi ** c_3
k = floor(reverse_fib(max_num) / c_3)
temp_1 = c_1 / sqrt(c_5)
temp_2 = phi3 * ((c_1 - phi3**k) / (c_1 - phi3))
temp_3 = psi3 * ((c_1 - psi3**k) / (c_1 - psi3))
sum = anint(temp_1 * (temp_2 - temp_3))
print *,"result: ", sum
end program
The code looks quite long, but basically it is the same as our Python version. I just needed to take care of the following things:
- To initialize variables with a higher precision we need to use the data type
real*16. However, depending on your Fortran compiler, you might have problems with that data type. If so, then you can check out this post on Stackoverflow. - I have not found a way to calculate the log of a number with an arbitrary base. However, this is not a problem. You can easily calculate it as shown in the code. For more information, you can check out this Wikipedia article or this post on Stackoverflow.
- I had to wrap constant numbers as
real*16as well (c_1,c_2,c_3,c_4,c_5). Before, I just put them in the code, but I got a small error. I assume constant numbers are of a lower precision datatype by default and when they are used in an expression the whole expression seems to be converted to this lower precision datatype as well. Probably, there is a better way to implement this, but I have not found out how to do this, yet. - When we cast the floating value to integer at the end, we also need a bigger integer datatype. How to initialize bigger integers in Fortran is described here.
It's actually the first time I wrote Fortran code. Thus, the code probably looks quite bad to a Fortran programmer (you can find a Fortran tutorial here), but I simply wanted to show if we can use Fortran to obtain the correct result. So, let's try it out. Save the above code in a file called fib_sum.f90. If you don't have a Fortran compiler, you can install e.g. the gcc package (beside a compiler for C and C++ it also offers a compiler for Fortran called gfortran). You can find instructions on how to install gcc, for instance, on Ubuntu here and for OSX via brew here. After you installed gcc, you can compile the Fortran file using the following command:
gfortran -o fib_sum fib_sum.f90
`
An executable file named fib_sum should have been created. Let's call it.
./fib_sum
`
You should get the following output.
result: 49597426547377748
We obtain 49597426547377748. This is the correct result. If we chose the lower precision datatype real*8 instead of real*16, we should have received the incorrect result of 49597426547377776 again (you can try it out!). Now, let's also measure the time. On the command line we can do this as follows:
time ./fib_sum
`
We should receive the following (or similar) output:
result: 49597426547377748
./fib_sum 0.00s user 0.00s system 69% cpu 0.008 total
It took 0.008 seconds, which are 8 milliseconds. This is faster than our code using the decimal module. On the other side, it is slower than our Python code implementing our brute force and improved approach. However, it's not so clear how comparable the time measurement between the %timeit magic and the command time for the command line are.
Fortran should be quite fast actually. Nevertheless, we got the correct solution and could definitely speed up.