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)