diff --git a/maths/prime_numbers.py b/maths/prime_numbers.py
index 38cc6670385d..dbf489c35af6 100644
--- a/maths/prime_numbers.py
+++ b/maths/prime_numbers.py
@@ -97,9 +97,9 @@ def benchmark():
     from timeit import timeit
 
     setup = "from __main__ import slow_primes, primes, fast_primes"
-    print(timeit("slow_primes(1_000_000_000_000)", setup=setup, number=1_000_000))
-    print(timeit("primes(1_000_000_000_000)", setup=setup, number=1_000_000))
-    print(timeit("fast_primes(1_000_000_000_000)", setup=setup, number=1_000_000))
+    print(timeit("list(slow_primes(1_000))", setup=setup, number=1_000))
+    print(timeit("list(primes(1_000))", setup=setup, number=1_000))
+    print(timeit("list(fast_primes(1_000))", setup=setup, number=1_000))
 
 
 if __name__ == "__main__":
diff --git a/maths/prime_sieve_eratosthenes.py b/maths/prime_sieve_eratosthenes.py
index 32eef9165bba..5183ba430152 100644
--- a/maths/prime_sieve_eratosthenes.py
+++ b/maths/prime_sieve_eratosthenes.py
@@ -11,6 +11,10 @@
 https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
 """
 
+from math import isqrt
+
+import numpy as np
+
 
 def prime_sieve_eratosthenes(num: int) -> list[int]:
     """
@@ -45,10 +49,53 @@ def prime_sieve_eratosthenes(num: int) -> list[int]:
     return [prime for prime in range(2, num + 1) if primes[prime]]
 
 
+def np_prime_sieve_eratosthenes(max_number: int) -> list[int]:
+    """
+    Returns prime numbers below max_number.
+    See: https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
+
+    >>> np_prime_sieve_eratosthenes(10)
+    [2, 3, 5, 7]
+    >>> np_prime_sieve_eratosthenes(2)
+    [2]
+    >>> np_prime_sieve_eratosthenes(1)
+    []
+    """
+    if max_number < 2:
+        return []
+
+    # List containing a bool value for every odd number below max_number/2
+    is_prime = np.ones((max_number + 1) // 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 benchmark():
+    """
+    Benchmarks
+    """
+    from timeit import timeit
+
+    print("Running performance benchmarks...")
+
+    functions = ["prime_sieve_eratosthenes", "np_prime_sieve_eratosthenes"]
+    for func in functions:
+        print(f"{func} : {timeit(f'{func}(10_000)', globals=globals(), number=10_000)}")
+
+
 if __name__ == "__main__":
     import doctest
 
     doctest.testmod()
 
     user_num = int(input("Enter a positive integer: ").strip())
-    print(prime_sieve_eratosthenes(user_num))
+    print(np_prime_sieve_eratosthenes(user_num))
+
+    benchmark()
diff --git a/project_euler/problem_187/sol1.py b/project_euler/problem_187/sol1.py
index 8944776fef50..2bfb5b8f189e 100644
--- a/project_euler/problem_187/sol1.py
+++ b/project_euler/problem_187/sol1.py
@@ -13,6 +13,8 @@
 
 from math import isqrt
 
+from maths.prime_sieve_eratosthenes import np_prime_sieve_eratosthenes
+
 
 def slow_calculate_prime_numbers(max_number: int) -> list[int]:
     """
@@ -38,15 +40,15 @@ def slow_calculate_prime_numbers(max_number: int) -> list[int]:
     return [i for i in range(2, max_number) if is_prime[i]]
 
 
-def calculate_prime_numbers(max_number: int) -> list[int]:
+def py_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)
+    >>> py_calculate_prime_numbers(10)
     [2, 3, 5, 7]
 
-    >>> calculate_prime_numbers(2)
+    >>> py_calculate_prime_numbers(2)
     []
     """
 
@@ -100,7 +102,7 @@ def while_solution(max_number: int = 10**8) -> int:
     10
     """
 
-    prime_numbers = calculate_prime_numbers(max_number // 2)
+    prime_numbers = py_calculate_prime_numbers(max_number // 2)
 
     semiprimes_count = 0
     left = 0
@@ -114,6 +116,31 @@ def while_solution(max_number: int = 10**8) -> int:
     return semiprimes_count
 
 
+def for_solution(max_number: int = 10**8) -> int:
+    """
+    Returns the number of composite integers below max_number have precisely two,
+    not necessarily distinct, prime factors.
+
+    >>> for_solution(30)
+    10
+    """
+
+    prime_numbers = py_calculate_prime_numbers(max_number // 2)
+
+    semiprimes_count = 0
+    right = len(prime_numbers) - 1
+    for left in range(len(prime_numbers)):
+        if left > right:
+            break
+        for r in range(right, left - 2, -1):
+            if prime_numbers[left] * prime_numbers[r] < max_number:
+                break
+        right = r
+        semiprimes_count += right - left + 1
+
+    return semiprimes_count
+
+
 def solution(max_number: int = 10**8) -> int:
     """
     Returns the number of composite integers below max_number have precisely two,
@@ -123,7 +150,7 @@ def solution(max_number: int = 10**8) -> int:
     10
     """
 
-    prime_numbers = calculate_prime_numbers(max_number // 2)
+    prime_numbers = np_prime_sieve_eratosthenes((max_number - 1) // 2)
 
     semiprimes_count = 0
     right = len(prime_numbers) - 1
@@ -146,7 +173,8 @@ def benchmark() -> None:
     # Running performance benchmarks...
     # slow_solution : 108.50874730000032
     # while_sol     : 28.09581200000048
-    # solution      : 25.063097400000515
+    # for_sol       : 25.063097400000515
+    # solution      : 5.219610300000568
 
     from timeit import timeit
 
@@ -154,6 +182,7 @@ def benchmark() -> None:
 
     print(f"slow_solution : {timeit('slow_solution()', globals=globals(), number=10)}")
     print(f"while_sol     : {timeit('while_solution()', globals=globals(), number=10)}")
+    print(f"for_sol       : {timeit('for_solution()', globals=globals(), number=10)}")
     print(f"solution      : {timeit('solution()', globals=globals(), number=10)}")