euler/python/e128.py

from lib_prime import primes, is_prime_rabin_miller
from functools import lru_cache


def first_approach():
    """ This was my original naiv approach that worked, but was way to slow
    because we computed the whole hex tile grid. """
    ps = set(primes(1000000))
    hextiles = {(0, 0): 1}
    add_ring(1, hextiles)
    add_ring(2, hextiles)

    # Visualization of how the coord system for the hex tile grid works.
    #                 
    #     2 1 0 1 2  
    #  4      8          
    #  3    9  19        
    #  2 10   2  18      
    #  1    3   7        
    #  0 11   1  17      
    #  1    4   6        
    #  2 12   5  16      
    #  3   13  15        
    #  4     14          
    #                 
    #                 

    pds = [1, 2, 8]
    for n in range(0, 100):
        ring_coord = ring_start(n)
        print("  ", ring_coord, hextiles[ring_coord])
        if n % 10 == 0:
            print(n)
        add_ring(n + 1, hextiles)
        for coord in ring_coords(n):
            pdv = pd(coord, hextiles, ps)
            if pdv == 3:
                v = hextiles[coord]
                assert v == first_number_ring(n) or v == first_number_ring(n + 1) - 1
                pds.append(hextiles[coord])
    print(len(pds))
    print(sorted(pds))
    target = 10
    if len(pds) >= target:
        return sorted(pds)[target - 1]
    return None


def ring_start(n):
    return (-n * 2, 0)


@lru_cache
def first_number_ring(n):
    if n == 0:
        return 1
    elif n == 1:
        return 2
    else:
        return first_number_ring(n - 1) + (n - 1) * 6


def get_nbvs_ring_start(n):
    #    8        1
    #  9   19   2   6   
    #    2        0  
    #  3   7    3   5
    #    l        4 
    v0 = first_number_ring(n)
    v1 = first_number_ring(n + 1)
    v2 = v1 + 1
    v3 = v0 + 1
    v4 = first_number_ring(n - 1)
    v5 = v1 - 1
    v6 = first_number_ring(n + 2) - 1
    nbvs = [v1, v2, v3, v4, v5, v6]
    # print(v0, nbvs)
    return nbvs


def pd_ring_start(n):
    v0 = first_number_ring(n)
    nbvs = get_nbvs_ring_start(n)
    r = 0
    for v in nbvs:
        if is_prime_rabin_miller(abs(v0 - v)):
            r += 1
    return r


def get_nbvs_ring_prev(n):
    """ Get neighbors and values for tile before top tile. """
    v0 = first_number_ring(n + 1) - 1
    v1 = first_number_ring(n + 2) - 1
    v2 = first_number_ring(n)
    v3 = first_number_ring(n - 1)
    v4 = v2 - 1
    v5 = first_number_ring(n + 1) - 2
    v6 = v1 - 1
    nbvs = [v1, v2, v3, v4, v5, v6]
    # print(v0, nbvs)
    return nbvs


def pd_ring_prev(n):
    v0 = first_number_ring(n + 1) - 1
    nbvs = get_nbvs_ring_prev(n)
    r = 0
    for v in nbvs:
        if is_prime_rabin_miller(abs(v0 - v)):
            r += 1
    return r


def add_ring(n, numbers):
    """ Adds ring n coords and values to numbers. """
    if n == 0:
        return
    first_coord = ring_start(n)
    current_number = first_number_ring(n)
    numbers[first_coord] = current_number
    current_coord = tuple(first_coord)

    for ro, co in [(1, -1), (2, 0), (1, 1), (-1, 1), (-2, 0), (-1, -1)]:
        for _ in range(n):
            current_coord = (current_coord[0] + ro, current_coord[1] + co)
            current_number += 1
            numbers[current_coord] = current_number

    # Reset first coord which is overriden.
    numbers[first_coord] = first_number_ring(n)


def get_neighbor_values(coord, numbers):
    neighbors = []

    for ro, co in [(-2, 0), (-1, 1), (1, 1), (2, 0), (1, -1), (-1, -1)]:
        nc = (coord[0] + ro, coord[1] + co)
        neighbors.append(numbers[nc])
    return neighbors


def ring_coords(n):
    """ Returns coords for ring n. """
    if n == 0:
        yield (0, 0)
    current_coord = ring_start(n)
    for ro, co in [(1, -1), (2, 0), (1, 1), (-1, 1), (-2, 0), (-1, -1)]:
        for _ in range(n):
            current_coord = (current_coord[0] + ro, current_coord[1] + co)
            yield current_coord


def pd(coord, numbers, primes):
    prime_delta = 0
    v = numbers[coord] 
    for nbv in get_neighbor_values(coord, numbers):
        if abs(v - nbv) in primes:
            prime_delta += 1
    return prime_delta


def test_get_nbvs_functions():
    """ Make sure that the function to compute the neighbor values works
    correctly by comparing the values to the original grid based approach. """
    hextiles = {(0, 0): 1}
    add_ring(1, hextiles)
    add_ring(2, hextiles)

    for n in range(2, 30):
        add_ring(n + 1, hextiles)
        nbvs1 = get_nbvs_ring_start(n)
        ring_coord = ring_start(n)
        nbvs2 = get_neighbor_values(ring_coord, hextiles)
        assert sorted(nbvs1) == sorted(nbvs2)

        nbvs1 = get_nbvs_ring_prev(n)
        ring_coord = (ring_start(n)[0] + 1, ring_start(n)[1] + 1)
        nbvs2 = get_neighbor_values(ring_coord, hextiles)
        assert sorted(nbvs1) == sorted(nbvs2)


def euler_128():
    # Conjecture: PD3s only occur at ring starts or at ring start minus one.
    # Under this assumption, we can significantly reduce the search space.
    # The only challenge is to find a way to directly compute the relevant
    # coords and neighbor values.

    test_get_nbvs_functions()

    target = 2000
    pds = [1, 2]
    for n in range(2, 80000):
        if pd_ring_start(n) == 3:
            v0 = first_number_ring(n)
            pds.append(v0)

        if pd_ring_prev(n) == 3:
            v0 = first_number_ring(n + 1) - 1
            pds.append(v0)

    assert len(pds) > target
    return sorted(pds)[target - 1]


if __name__ == "__main__":
    solution = euler_128()
    print("e128.py: " + str(solution))
    assert(solution == 14516824220)