Project Euler Problem 734

Project Euler Problem 734

Solution

Answer: 557988060

Let $P_n$ be the set of primes $\le n$, and define

$$f(x)= \begin{cases} 1 & x\in P_n\ 0 & \text{otherwise} \end{cases}$$

We want the number of $k$-tuples of primes whose bitwise-OR is also prime.

The key observation is that bitwise-OR convolution becomes simple under the subset zeta transform.

For a mask $m$, define

$$F(m)=\sum_{s\subseteq m} f(s)$$

which counts primes whose set bits are contained in $m$.

Then:

• $F(m)^k$ counts ordered $k$-tuples of primes whose OR is a subset of $m$.
• Applying Möbius inversion on subsets recovers the number of tuples whose OR is exactly $m$.

Thus:

1. Build the prime-indicator array.
2. Apply the OR-zeta transform.
3. Raise each transformed value to the $k$-th power.
4. Apply inverse OR-zeta transform.
5. Sum over prime masks.

The following Python code computes the answer exactly.

from math import isqrt

MOD = 1_000_000_007
n = 10**6
k = 999983

# ------------------------------------------------------------
# Sieve of Eratosthenes
# ------------------------------------------------------------
sieve = bytearray(b'\x01') * (n + 1)
sieve[0] = sieve[1] = 0

for p in range(2, isqrt(n) + 1):
    if sieve[p]:
        sieve[p * p:n + 1:p] = b'\x00' * (((n - p * p) // p) + 1)

primes = [i for i, v in enumerate(sieve) if v]

# ------------------------------------------------------------
# OR convolution via subset zeta transform
# ------------------------------------------------------------
BITS = 20                  # since 2^20 > 10^6
SIZE = 1 << BITS

# f[x] = 1 if x is prime else 0
F = [0] * SIZE
for p in primes:
    F[p] = 1

# ------------------------------------------------------------
# OR-zeta transform
# After this:
# F[m] = number of primes s with s ⊆ m (bitwise)
# ------------------------------------------------------------
for b in range(BITS):
    bit = 1 << b
    for mask in range(SIZE):
        if mask & bit:
            F[mask] += F[mask ^ bit]
            F[mask] %= MOD

# ------------------------------------------------------------
# k-th OR convolution power
# ------------------------------------------------------------
for i in range(SIZE):
    F[i] = pow(F[i], k, MOD)

# ------------------------------------------------------------
# Möbius inversion for OR transform
# ------------------------------------------------------------
for b in range(BITS):
    bit = 1 << b
    for mask in range(SIZE):
        if mask & bit:
            F[mask] -= F[mask ^ bit]
            F[mask] %= MOD

# ------------------------------------------------------------
# Sum counts whose OR is prime
# ------------------------------------------------------------
answer = sum(F[p] for p in primes) % MOD

print(answer)

Verification against the examples:

• $T(5,2)=5$
• $T(100,3)=3355$
• $T(1000,10)\equiv 2071632 \pmod{10^9+7}$

all match the problem statement.

Running the computation for $T(10^6,999983)$ gives:

Answer: 557988060