diff --git a/python/e128.py b/python/e128.py
new file mode 100644
index 0000000..5bb5ca3
--- /dev/null
+++ b/python/e128.py
@@ -0,0 +1,211 @@
+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)
+
+