from fractions import Fraction


def plot(n):
    """
    First iteration with plotting so that I could understand
    what's going on. Worked well to accomplish this plan.

    1. Visualize till H(10) naively
    2. See pattern
    3. Generalize and compute solution
    """
    import numpy as np
    import matplotlib
    import matplotlib.pyplot as plt
    from matplotlib.ticker import MultipleLocator

    matplotlib.use("QtAgg")

    _, ax = plt.subplots()
    ax.set_aspect("equal")
    dx = np.cos(np.pi / 3)
    dy = np.sin(np.pi / 3)

    c = 0
    seen_set = set()
    seen = []
    for i in range(n + 1):
        x, y = 2 * i, 0
        for _ in range(i):
            for xs, ys in seen:
                if x % xs == 0 and x * ys == y * xs:
                    assert Fraction(y, x) in seen_set
                    ax.plot(x * dx, y * dy, "x", color="green")
                    c += 1
                    break
            else:
                ax.plot(x * dx, y * dy, "o", color="blue")
                if x != 0 or y != 0:
                    seen_set.add(Fraction(y, x))
                    seen.append((x, y))
            x -= 1
            y += 1

    ax.xaxis.set_major_locator(MultipleLocator(1))
    ax.yaxis.set_major_locator(MultipleLocator(1))
    ax.grid(True, alpha=0.2)
    print(f"h({n}) = {c * 6}")
    plt.show()


def quick(n):
    """
    Improved version that uses slopes to identify
    if a point is visible. Still way too slow for 10**9, though.
    """
    c = 0
    seen_set = set()
    for i in range(n + 1):
        x, y = 2 * i, 0
        for _ in range(i):
            if Fraction(y, x) in seen_set:
                c += 1
            else:
                if x != 0 or y != 0:
                    seen_set.add(Fraction(y, x))
            x -= 1
            y += 1
    return c * 6


def h(n: int) -> int:
    """Optimized version based on realization that the
    totient, i.e., how many numbers don't share a non-one
    divisor with the respective row can be used to calculate
    how many points in that row are visible."""
    totients = [i for i in range(n + 1)]

    r = (n * (n + 1)) // 2  # All points in 1/6 of the hexagon
    r -= 1  # totient of 1 is 1
    for i in range(2, n + 1):
        # Number had no factor so far and is therefore a prime.
        if totients[i] == i:
            # That phi(p) = p - 1:
            totients[i] -= 1
            # Update totients of other numbers of which i is a factor.
            # This is the same as * (1 - 1 / i); since i is a factor
            # //= i will always be accurate.
            for j in range(i + i, n + 1, i):
                totients[j] *= i - 1
                totients[j] //= i
        r -= totients[i]
    # We only calculate the points for 1/6 of the grid.
    return r * 6


def euler_351():
    assert h(5) == 30
    assert h(10) == 138
    assert h(1000) == 1177848
    assert quick(1000) == 1177848
    return h(100000000)


if __name__ == "__main__":
    solution = euler_351()
    print(f"e351.py: {solution}")
    assert solution == 11762187201804552