perf(project_euler): optimize problem_073/sol1 with Mobius inversion

Replaces the O(n^2 log n) gcd-per-pair loop with Mobius inversion,
reducing complexity to O(n log log n). Benchmarks: 2719ms -> 50ms (54x).

Fixes #8594

Author: Ankur-Kataria

project_euler/problem_073/sol1.py

@@ -14,15 +14,28 @@
 
 How many fractions lie between 1/3 and 1/2 in the sorted set
 of reduced proper fractions for d ≤ 12,000?
-"""
 
-from math import gcd
+References:
+    - https://en.wikipedia.org/wiki/M%C3%B6bius_function
+    - https://en.wikipedia.org/wiki/Farey_sequence
+"""
 
 
 def solution(max_d: int = 12_000) -> int:
     """
-    Returns number of fractions lie between 1/3 and 1/2 in the sorted set
-    of reduced proper fractions for d ≤ max_d
+    Returns the number of fractions lying between 1/3 and 1/2 in the sorted
+    set of reduced proper fractions for d ≤ max_d.
+
+    Uses Möbius inversion to avoid a per-fraction gcd check.  Let
+    f(j) = #{m ∈ ℤ : j/3 < m < j/2} = max(0, (j-1)//2 - j//3).
+
+    By Möbius inversion on the coprimality condition:
+        count = Σ_{e=1}^{max_d} μ(e) · Σ_{j=1}^{⌊max_d/e⌋} f(j)
+
+    Building prefix sums of f reduces the inner sum to an O(1) lookup,
+    giving overall complexity O(max_d log log max_d) for the Möbius sieve
+    plus O(max_d) for the final summation — compared with O(max_d² log max_d)
+    for the naive gcd-per-pair approach.
 
     >>> solution(4)
     0
@@ -32,19 +45,41 @@ def solution(max_d: int = 12_000) -> int:
 
     >>> solution(8)
     3
+
+    >>> solution(12_000)
+    7295372
     """
+    # Compute the Möbius function for 1..max_d via a linear sieve.
+    # μ(n) = 1 if n=1; (-1)^k if n is a product of k distinct primes; 0 otherwise.
+    mobius = [0] * (max_d + 1)
+    mobius[1] = 1
+    primes: list[int] = []
+    is_composite = bytearray(max_d + 1)
+
+    for i in range(2, max_d + 1):
+        if not is_composite[i]:
+            primes.append(i)
+            mobius[i] = -1
+        for prime in primes:
+            product = i * prime
+            if product > max_d:
+                break
+            is_composite[product] = 1
+            if i % prime == 0:
+                mobius[product] = 0  # prime² divides product
+                break
+            mobius[product] = -mobius[i]  # new prime factor
+
+    # prefix[j] = Σ_{k=1}^{j} f(k),  where f(k) = #{m : k/3 < m < k/2}
+    prefix = [0] * (max_d + 1)
+    for j in range(1, max_d + 1):
+        prefix[j] = prefix[j - 1] + max(0, (j - 1) // 2 - j // 3)
 
-    fractions_number = 0
-    for d in range(max_d + 1):
-        n_start = d // 3 + 1
-        n_step = 1
-        if d % 2 == 0:
-            n_start += 1 - n_start % 2
-            n_step = 2
-        for n in range(n_start, (d + 1) // 2, n_step):
-            if gcd(n, d) == 1:
-                fractions_number += 1
-    return fractions_number
+    return sum(
+        mobius[e] * prefix[max_d // e]
+        for e in range(1, max_d + 1)
+        if mobius[e]
+    )
 
 
 if __name__ == "__main__":