Performance: 58% faster Project Euler 070

project_euler/problem_070/sol1.py

@@ -26,22 +26,78 @@
 
 References:
 Finding totients
-https://en.wikipedia.org/wiki/Euler's_totient_function#Euler's_product_formula
+https://en.wikipedia.org/wiki/Euler%27s_totient_function#Euler%27s_product_formula
 """
 from __future__ import annotations
 
+from math import isqrt
+
 import numpy as np
 
 
-def get_totients(max_one: int) -> list[int]:
+def calculate_prime_numbers(max_number: int) -> list[int]:
+    """
+    Returns prime numbers below max_number.
+    See: https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
+
+    >>> calculate_prime_numbers(10)
+    [2, 3, 5, 7]
+    >>> calculate_prime_numbers(2)
+    []
+    """
+    if max_number <= 2:
+        return []
+
+    # List containing a bool value for every odd number below max_number/2
+    is_prime = [True] * (max_number // 2)
+
+    for i in range(3, isqrt(max_number - 1) + 1, 2):
+        if is_prime[i // 2]:
+            # Mark all multiple of i as not prime using list slicing
+            is_prime[i**2 // 2 :: i] = [False] * (
+                # Same as: (max_number - (i**2)) // (2 * i) + 1
+                # but faster than len(is_prime[i**2 // 2 :: i])
+                len(range(i**2 // 2, max_number // 2, i))
+            )
+
+    return [2] + [2 * i + 1 for i in range(1, max_number // 2) if is_prime[i]]
+
+
+def np_calculate_prime_numbers(max_number: int) -> list[int]:
+    """
+    Returns prime numbers below max_number.
+    See: https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
+
+    >>> np_calculate_prime_numbers(10)
+    [2, 3, 5, 7]
+    >>> np_calculate_prime_numbers(2)
+    []
+    """
+    if max_number <= 2:
+        return []
+
+    # List containing a bool value for every odd number below max_number/2
+    is_prime = np.ones(max_number // 2, dtype=bool)
+
+    for i in range(3, isqrt(max_number - 1) + 1, 2):
+        if is_prime[i // 2]:
+            # Mark all multiple of i as not prime using list slicing
+            is_prime[i**2 // 2 :: i] = False
+
+    primes = np.where(is_prime)[0] * 2 + 1
+    primes[0] = 2
+    return primes.tolist()
+
+
+def slow_get_totients(max_one: int) -> list[int]:
     """
     Calculates a list of totients from 0 to max_one exclusive, using the
     definition of Euler's product formula.
 
-    >>> get_totients(5)
+    >>> slow_get_totients(5)
     [0, 1, 1, 2, 2]
 
-    >>> get_totients(10)
+    >>> slow_get_totients(10)
     [0, 1, 1, 2, 2, 4, 2, 6, 4, 6]
     """
     totients = np.arange(max_one)
@@ -54,6 +110,66 @@ def get_totients(max_one: int) -> list[int]:
     return totients.tolist()
 
 
+def slicing_get_totients(max_one: int) -> list[int]:
+    """
+    Calculates a list of totients from 0 to max_one exclusive, using the
+    definition of Euler's product formula.
+
+    >>> slicing_get_totients(5)
+    [0, 1, 1, 2, 2]
+
+    >>> slicing_get_totients(10)
+    [0, 1, 1, 2, 2, 4, 2, 6, 4, 6]
+    """
+    totients = np.arange(max_one)
+
+    for i in range(2, max_one):
+        if totients[i] == i:
+            totients[i::i] -= totients[i::i] // i
+
+    return totients.tolist()
+
+
+def get_totients(limit) -> list[int]:
+    """
+    Calculates a list of totients from 0 to max_one exclusive, using the
+    definition of Euler's product formula.
+
+    >>> get_totients(5)
+    [0, 1, 1, 2, 2]
+
+    >>> get_totients(10)
+    [0, 1, 1, 2, 2, 4, 2, 6, 4, 6]
+    """
+    totients = np.arange(limit)
+    primes = calculate_prime_numbers(limit)
+
+    for i in primes:
+        totients[i::i] -= totients[i::i] // i
+
+    return totients.tolist()
+
+
+def np_get_totients(limit) -> list[int]:
+    """
+    Calculates a list of totients from 0 to max_one exclusive, using the
+    definition of Euler's product formula.
+
+    >>> np_get_totients(5)
+    [0, 1, 1, 2, 2]
+
+    >>> np_get_totients(10)
+    [0, 1, 1, 2, 2, 4, 2, 6, 4, 6]
+    """
+    totients = np.arange(limit)
+    primes = np_calculate_prime_numbers(limit)
+
+    for i in primes:
+        totients[i::i] -= totients[i::i] // i
+
+    return totients.tolist()
+
+
 def has_same_digits(num1: int, num2: int) -> bool:
     """
     Return True if num1 and num2 have the same frequency of every digit, False
@@ -71,6 +187,51 @@ def has_same_digits(num1: int, num2: int) -> bool:
     return sorted(str(num1)) == sorted(str(num2))
 
 
+def slow_solution(max_n: int = 10000000) -> int:
+    """
+    Finds the value of n from 1 to max such that n/φ(n) produces a minimum.
+
+    >>> slow_solution(100)
+    21
+
+    >>> slow_solution(10000)
+    4435
+    """
+    totients = slow_get_totients(max_n + 1)
+
+    return common_solution(totients, max_n)
+
+
+def slicing_solution(max_n: int = 10000000) -> int:
+    """
+    Finds the value of n from 1 to max such that n/φ(n) produces a minimum.
+
+    >>> slicing_solution(100)
+    21
+
+    >>> slicing_solution(10000)
+    4435
+    """
+    totients = slicing_get_totients(max_n + 1)
+
+    return common_solution(totients, max_n)
+
+
+def py_solution(max_n: int = 10000000) -> int:
+    """
+    Finds the value of n from 1 to max such that n/φ(n) produces a minimum.
+
+    >>> py_solution(100)
+    21
+
+    >>> py_solution(10000)
+    4435
+    """
+    totients = get_totients(max_n + 1)
+
+    return common_solution(totients, max_n)
+
+
 def solution(max_n: int = 10000000) -> int:
     """
     Finds the value of n from 1 to max such that n/φ(n) produces a minimum.
@@ -81,10 +242,23 @@ def solution(max_n: int = 10000000) -> int:
     >>> solution(10000)
     4435
     """
+    totients = np_get_totients(max_n + 1)
+
+    return common_solution(totients, max_n)
 
+
+def common_solution(totients: list[int], max_n: int = 10000000) -> int:
+    """
+    Finds the value of n from 1 to max such that n/φ(n) produces a minimum.
+
+    >>> common_solution(get_totients(101), 100)
+    21
+
+    >>> common_solution(get_totients(10001), 10000)
+    4435
+    """
     min_numerator = 1  # i
     min_denominator = 0  # φ(i)
-    totients = get_totients(max_n + 1)
 
     for i in range(2, max_n + 1):
         t = totients[i]
@@ -96,5 +270,29 @@ def solution(max_n: int = 10000000) -> int:
     return min_numerator
 
 
+def benchmark() -> None:
+    """
+    Benchmark
+    """
+    # Running performance benchmarks...
+    # Solution    : 49.389978999999585
+    # Py Solution : 56.19136740000067
+    # Slicing Sol : 70.83823779999875
+    # Slow Sol    : 118.29514729999937
+
+    from timeit import timeit
+
+    print("Running performance benchmarks...")
+
+    print(f"Solution    : {timeit('solution()', globals=globals(), number=10)}")
+    print(f"Py Solution : {timeit('py_solution()', globals=globals(), number=10)}")
+    print(f"Slicing Sol : {timeit('slicing_solution()', globals=globals(), number=10)}")
+    print(f"Slow Sol    : {timeit('slow_solution()', globals=globals(), number=10)}")
+
+
 if __name__ == "__main__":
-    print(f"{solution() = }")
+    print(f"Solution    : {solution()}")
+    print(f"Py Solution : {py_solution()}")
+    print(f"Slicing Sol : {slicing_solution()}")
+    print(f"Slow Sol    : {slow_solution()}")
+    benchmark()