Solve e816 because I have to get easy ones for my Euler goal.

python/e816.py

@@ -0,0 +1,39 @@
+from math import sqrt
+
+
+def distance(p1, p2) -> float:
+    return (p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2
+
+
+def points(k: int):
+    s = 290797
+    mod = 50515093
+    # s = 5
+    # mod = 23
+    points = []
+    for _ in range(k):
+        sn = (s ** 2) % mod
+        p = (s, sn)
+        points.append(p)
+        s = (sn ** 2) % mod
+    return points
+
+
+def euler_816():
+    ps = sorted(points(2000000), key=lambda p: p[0])
+    d_min = float("inf")
+    sweep_dist = 5  # heuristic
+    for i in range(len(ps)):
+        for j in range(i, min(i + sweep_dist, len(ps))):
+            if i == j:
+                continue
+            d  = distance(ps[i], ps[j])
+            if d < d_min:
+                d_min = d
+    return round(sqrt(d_min), 9)
+
+
+if __name__ == "__main__":
+    solution = euler_816()
+    print("e816.py: " + str(solution))
+    assert(solution == 20.880613018)