diff --git a/project_euler/problem_810/__init__.py b/project_euler/problem_810/__init__.py
new file mode 100644
index 000000000000..e69de29bb2d1
diff --git a/project_euler/problem_810/sol1.py b/project_euler/problem_810/sol1.py
new file mode 100644
index 000000000000..c47dae8a5e42
--- /dev/null
+++ b/project_euler/problem_810/sol1.py
@@ -0,0 +1,151 @@
+"""
+Project Euler Problem 810: https://projecteuler.net/problem=810
+
+We use x ⊕ y for the bitwise XOR of x and y.
+Define the XOR-product of x and y, denoted by x ⊗ y,
+similar to a long multiplication in base 2,
+except the intermediate results are XORed instead of usual integer addition.
+
+For example, 7 ⊗ 3 = 9, or in base 2, 111_2 ⊗ 11_2 = 1001_2:
+        111
+    ⊗   11
+    -------
+        111
+       111
+    -------
+       1001
+An XOR-Prime is an integer n greater than 1 that is not an
+XOR-product of two integers greater than 1.
+The above example shows that 9 is not an XOR-prime.
+Similarly, 5 = 3 ⊗ 3 is not an XOR-prime.
+The first few XOR-primes are 2, 3, 7, 11, 13, ... and the 10th XOR-prime is 41.
+
+Find the 5,000,000 th XOR-prime.
+
+References:
+http://en.wikipedia.org/wiki/M%C3%B6bius_function
+
+"""
+
+import math
+
+
+def get_divisors(num: int) -> set[int]:
+    """
+    Return all positive divisors of num.
+    >>> get_divisors(12)
+    {1, 2, 3, 4, 6, 12}
+    """
+    divisors = {1}
+    for i in range(2, int(math.sqrt(num)) + 1):
+        if num % i == 0:
+            divisors.add(i)
+            divisors.add(num // i)
+    divisors.add(num)
+    return divisors
+
+
+def xor_multiply(op_a: int, op_b: int) -> int:
+    """
+    Perform XOR-based multiplication (polynomial multiplication mod 2).
+    >>> xor_multiply(3, 5)
+    15
+    """
+    result = 0
+    while op_b:
+        if op_b & 1:
+            result ^= op_a
+        op_a <<= 1
+        op_b >>= 1
+    return result
+
+
+def mobius_table(lim: int, k: int = 2) -> list[int]:
+    """
+    Compute a modified Mobius function table up to `lim`.
+    >>> mobius_table(10)[:6]
+    [0, 1, -1, -1, 0, -1]
+    """
+    mob = [0] + [1] * lim
+    is_prime = [True] * (lim + 1)
+    is_prime[0] = is_prime[1] = False
+
+    for p in range(2, lim + 1):
+        if is_prime[p]:
+            mob[p] *= -1
+            for mul in range(2 * p, lim + 1, p):
+                is_prime[mul] = False
+                mob[mul] *= -1
+
+            p_pow = p**k
+            if p_pow <= lim:
+                for mul in range(p_pow, lim + 1, p_pow):
+                    mob[mul] = 0
+    return mob
+
+
+def cnt_irred_pol(num: int) -> int:
+    """
+    Return the number of monic irreducible polynomials of degree num over GF(2).
+    >>> cnt_irred_pol(3)
+    2
+    """
+    mob = mobius_table(num)
+    total = sum(mob[d] * (2 ** (num // d)) for d in get_divisors(num))
+    return total // num
+
+
+def xor_prime_func(tgt_idx: int) -> int:
+    """
+    Find the N-th XOR-prime (irreducible polynomial) index approximation.
+    >>> xor_prime_func(10)
+    41
+    """
+    total, degree = 0, 1
+
+    while True:
+        cnt = cnt_irred_pol(degree)
+        if total + cnt > tgt_idx:
+            break
+        total += cnt
+        degree += 1
+
+    lim = 1 << (degree + 1)
+    is_prime = [True] * lim
+    is_prime[0] = is_prime[1] = False
+
+    for even in range(4, lim, 2):
+        is_prime[even] = False
+
+    cnt = 1
+    for num in range(3, lim, 2):
+        if not is_prime[num]:
+            continue
+
+        cnt += 1
+        if cnt == tgt_idx:
+            return num
+
+        mul = num
+        while True:
+            prod = xor_multiply(mul, num)
+            if prod >= lim:
+                break
+            is_prime[prod] = False
+            mul += 2
+
+    raise ValueError("Could not compute the XOR-prime.")
+
+
+def solution(nth: int = 5000000) -> int:
+    """
+    Compute the Nth XOR prime and print timing.
+    >>> solution(10)
+    41
+    """
+    result = xor_prime_func(nth)
+    return result
+
+
+if __name__ == "__main__":
+    print(f"{solution() = }")