Solve problem 351

python/e351.py

@@ -0,0 +1,108 @@
+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