Problem 304

Comments

The key to this problem is in the modulo. Ends up that the Fibonacci sequence is cyclical for every modulo. It’s just a matter of finding what the period – referred to as the Pisano period – of this cycle is. Without this, the values of the Fibonacci numbers would be prohibitively high.

I was able to find the period by looking at the OEIS page for the Pisano periods. It included several methods on how to find this number. I ended up using a small python algorithm shared there. I decided not to include this in my final code since it is much more code to find the number than it is to just figure it out once and save it to a variable.

My program is much slower than I would like it to be. Most of the time is taken up in just finding and saving all the Fibonacci numbers. When not saved to a list, this process seems to be much faster. I was unable to find something that worked any faster. However, there are a couple of things that I believe will speed this up. First, only find the Fibonacci numbers that are required instead of the full 900,788,112 of them to complete the cycle. Second, I am reading and learning of a way to extend the sieve of Eratosthenes past the first \(\sqrt{n}\) numbers. I believe I can find the first group of primes as I have here, then use those to sieve for another set of them \(>10^{14}\).

Code

modulo = 1234567891011

def solve(n):

    top = int(2 * 10**7)
    sieve = [False, False] + [True] * (top - 1)
    primes = []
    for p, isprime in enumerate(sieve):
        if not isprime:
            continue
        primes.append(p)
        num = p * p
        while num <= top:
            sieve[num] = False
            num += p

    period = 900788112
    a = i = 0
    b = 1
    fibs = []
    while i < period:
        fibs.append(a)
        a, b, i = b, (b + a) % modulo, i + 1

    def next_prime(n):
        k = n - (n % 2 == 0)
        while True:
            k += 2
            stop = int(k ** 0.5)
            for p in primes:
                if p > stop:
                    return k
                if k % p == 0:
                    break
            else:
                raise AssertionError('not enough primes')

    total, p = 0, 10**14
    for k in range(1, n + 1):
        p = next_prime(p)
        total += fibs[p % period]
    return total % modulo

if __name__ == '__main__':
    import sys
    n = eval(sys.argv[1])
    print(solve(n))

Tests

import pytest
from problem import solve

_test_solve = (
        (1, 428562224098),
        (2, 1081741301637),
        (3, 526706959321),
        (4, 262434673500),
)

@pytest.mark.parametrize('n,expect', _test_solve)
def test_solve(n, expect):
    assert expect == solve(n)