# Bisection, linear recurrences and even Fibonacci numbers

So this code is yet another attempt at solving the second Project Euler problem to improve my handling of Python. The purpose of the code is to solve the problem below

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.

### Method

However, I wanted to do it as fast as possible meaning:

• Cache small values and use a fast lookup table (bisection lookup)
• For larger values iterate over the even fibonacci values directly using the linear recurrence relation E(n) = 4 E(n-1) + E(n - 2) with E(0)=0 and E(1)=2.

### Questions / Wanted feedback

In particular I wanted to know if my general definition of a recurrence relation could be improved. In particular values[:-1], values[-1] = values[1:], last feels quite unpythonic to me. Secondly I am wondering if my docstrings are clear enough. I tried to strictly follow the Google Docstring style.

As a side note # fmt: on and #fmt: off are strictly neccecary to make sure my formater does not format my lookup tables.

### Code

"""
This code solves Project Euler Problem 2:

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.

See https://projecteuler.net/index.php?section=problems&id=2 for details
"""

import bisect

EvenFibSum = int
Limit = int

PE_002_LIMIT = 4 * 10 ** 6

# fmt: off
EVEN_FIBS = [
0, 2, 8, 34, 144, 610, 2584, 10946, 46368, 196418, 832040, 3524578,
14930352, 63245986, 267914296, 1134903170, 4807526976, 20365011074,
86267571272, 365435296162, 1548008755920, 6557470319842, 27777890035288,
117669030460994, 498454011879264, 2111485077978050, 8944394323791464,
37889062373143906, 160500643816367088, 679891637638612258,
2880067194370816120, 12200160415121876738
]
# fmt: on

EVEN_FIBS_ = set(EVEN_FIBS)

# fmt: off
EVEN_FIBS_CUMSUM = [
0, 2, 10, 44, 188, 798, 3382, 14328, 60696, 257114, 1089154, 4613732,
19544084, 82790070, 350704366, 1485607536, 6293134512, 26658145586,
112925716858, 478361013020, 2026369768940, 8583840088782, 36361730124070,
154030760585064, 652484772464328, 2763969850442378, 11708364174233842,
49597426547377748, 210098070363744836, 889989708002357094,
3770056902373173214, 15970217317495049952
]
# fmt: on

def linear_reccurence(constants: list[int], initials: list[int]) -> int:
"""Returns a generator for a linear recucurence relation

A linear recurrence relation is of the form

A(n) = c0 * A(n-1) + c0 * A(n-2) + ... + ck * A(n - k);
A(0) = I0, A(1) = I1, ..., A(k) = Ik

and would correspond to constants = [ck ..., c1, c0] and initials =
[Ik, ..., I1, I0].  Note that the values here are stored in ascending order

Args:
constants: Defines the constants in the recucurence relation.
initials: The initial values for the recurrence relation.

Yields:
The next value in the recucurence relation.

Examples:
Returns the Fibonacci numbers F(n) = F(n-1) + F(n-2) with F(0)=0, F(1)=1.

>>> fibonacci = linear_reccurence([1, 1], [0, 1])
>>> print([next(fibonacci) for _ in range(10)])
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34]

Returns the Lucas numbers L(n) = L(n-1) + L(n-2) with L(0)=1, L(1)=3.

>>> lucas = linear_reccurence([1, 1], [1, 3])
>>> print([next(lucas) for _ in range(10)])
[1, 3, 4, 7, 11, 18, 29, 47, 76, 123]
"""

values = initials.copy()
for value in values:
yield value
while True:
last = sum(const * value for const, value in zip(constants, values))
values[:-1], values[-1] = values[1:], last
yield last

def PE_002(limit: Limit = PE_002_LIMIT) -> EvenFibSum:
"""Sums all even fibonacci numbers under some limit

Args:
limit: Sums all even fibonacci numbers less than this limit

Returns:
The sum of all even fibonacci numbers less than some limit

Examples:
>>> limits = [0, 2, 8, 10, 10**8]
>>> print([PE_002(lim) for lim in limits])
[0, 2, 10, 10, 82790070]

>>> print(PE_002(2**65-1))
15970217317495049952
"""

def _even_fib_sum_large(limit: Limit) -> EvenFibSum:
total = EVEN_FIBS_CUMSUM[-1] - EVEN_FIBS[-2] - EVEN_FIBS[-1]
# The linear recurrence relation for the even fibonacci numbers is
#   E(n) = 4 * E(n - 1) + 1 * E(n - 2); E(0) = 0, E(1) = 2
# See https://math.stackexchange.com/questions/94359/need-help-deriving-recurrence-relation-for-even-valued-fibonacci-numbers
even_fibonacci = linear_reccurence([1, 4], [EVEN_FIBS[-2], EVEN_FIBS[-1]])
while (even_fib := next(even_fibonacci)) < limit:
total += even_fib

def _even_fib_sum_small(limit: Limit) -> EvenFibSum:
# The code performs a lookup to find the largest index such that EVEN_FIBS[index] <= limit.
# The lookup is done in O(log n) using a basic bisection algorithm.
# The offset is added because bisection performs < and we need <=
offset = 0 if limit in EVEN_FIBS_ else 1
index = bisect.bisect_left(EVEN_FIBS, limit)
return EVEN_FIBS_CUMSUM[index - offset]

if limit > EVEN_FIBS[-1]:
return _even_fib_sum_large(limit)
return _even_fib_sum_small(limit)

if __name__ == "__main__":
import doctest
import argparse

doctest.testmod()

parser = argparse.ArgumentParser(
description="Solves Project Euler 2; Sums all even fibonacci numbers less than limit"
)
dest="limit",
nargs="?",
type=int,
default=PE_002_LIMIT,
help="Sums all even fibonacci numbers less than this number",
)
args = parser.parse_args()

print(PE_002(args.limit))


1. Type aliases

These two lines confused me at first:

EvenFibSum = int
Limit = int


After I read further down, it became clear that the purpose of these two lines was to create type aliases for clearer annotations. However, you could maybe consider adding a comment above these two lines to explain what's going on. Alternatively, consider using typing.Annotated or typing.NewType. (Personally, I'm a fan of typing.Annotated.)

For example, perhaps:

from typing import Annotated

EvenFibSum = Annotated[
int,
"A positive integer representing the sum "
"of all even Fibonacci numbers below some limit"
]

Limit = Annotated[
int,
"A positive integer representing the limit "
"below which all even Fibonacci numbers are to be summed"
]


2. [Edited]: Consider choosing more distinct names for EVEN_FIBS and EVEN_FIBS_, as discussed in the comments.

3. Consider moving your argparsing to a separate function.

Suggested refactoring of your if __name__ == '__main__' block:

def get_cmd_args():
from argparse import ArgumentParser

parser = ArgumentParser(
description="Solves Project Euler 2; Sums all even fibonacci numbers less than limit"
)

dest="limit",
nargs="?",
type=int,
default=PE_002_LIMIT,
help="Sums all even fibonacci numbers less than this number",
)

return parser.parse_args()

if __name__ == "__main__":
import doctest
doctest.testmod()
args = get_cmd_args()
print(PE_002(args.limit))


Your docstrings look great to me!

• Thanks for the great feedback! Could you perhaps add an example on how to implement the Annotations? I can not seem to figure out a good value for the metadata x in Annotated[T, x]. Also how would the code work if EVEN_FIBS is removed? Do note that I use EVEN_FIBS_ (underscore) as a lookup and EVEN_FIBS as a list. For instance I traverse this list in the bisection part. I thought variable_ was the standard way of defining the set / dict conterpart for a list? Lastly, would get_cmd_args be better as a "private method" (_get_cmd_args)? Aug 9 at 13:03
• My apologies, I hadn't noticed that EVEN_FIBS and EVEN_FIBS_ were different variables. I wasn't aware of that convention regarding trailing underscore names — personally, I think I would always go for more distinct names, but that's just my opinion. Aug 9 at 13:14
• I've added some suggestions regarding the Annotations, though it's obviously fairly opinionated! Can see the argument for get_cmd_args being a private method. Aug 9 at 13:29

Disclaimer 1: This is not a review, but an extended comment.

Disclaimer 2: ProjectEuler is not about programming. It is about math. As long as you understand what the problem is about, an actual computation is supposed to be very simple.

I wanted to do it as fast as possible

You can do it much faster.

Observation #1: a parity of Fibonacci numbers follows the pattern of

odd even odd
odd even odd
...


(more or less obvious, but still try to prove it) so the sum of even-valued terms is the sum of every third term.

Observation #2: since Fibonacci numbers (as any linear recurrence) have a nice closed-form representation, observe that the answer is a sum of two geometric progressions. With a little work you may express it in terms of a couple of larger Fibonacci numbers.

Observation \$3: computing an nth Fibonacci number (as any linear recurrence) does not need to be linear in n. A (matrix) exponentiation-by-squaring let you do it in a logarithmic time.

• Here is a code implementing most of your ideas pastebin.com/TFNQBHVx. The reason I decided to try to implement it using recurrence relations and bisections is more for mostly for learning. My only issue with your comment is that we need to figure out how many terms we need to sum to take advantage of the closed forms. This leads to rounding off errors for large values unfortunately. Secondly taking advantage of maths is very specific to this particular problem and hard to generalize. Aug 7 at 0:40
• This leads to rounding off errors for large values unfortunately - with a correct implementation it should not. Pay attention to the Observation #2.
– vnp
Aug 7 at 0:42
• @vpn Did you check out my pastebin? If you use Bennets formula you still need to solve F(3n) < limit for n to figure out how many terms to sum. In addition there are far faster algorithms for finding the nth fibonacci number than matrix multiplication. The pastebin uses the one from here stackoverflow.com/a/48171368/1048781. Aug 7 at 0:49
• Solving F(3n) < limit is more or less $O(\log \log {limit})$. The SO answer you linked is still $\log{n}$ as far as I can tell (n being properly defined of course).
– vnp
Aug 7 at 1:06
• Solving F(3n) < limit is what causes the rounding errors Aug 7 at 1:09