109 lines
3.0 KiB
Python
109 lines
3.0 KiB
Python
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
|