Problem 869: Prime Guessing

View on Project Euler

Project Euler Problem 869 Solution

EulerSolve provides an optimized solution for Project Euler Problem 869, Prime Guessing, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let \(P_N\) be the set of primes not exceeding \(N\), with \(N=10^8\). One prime is chosen uniformly at random. We then try to guess its binary expansion one bit at a time, starting from the least significant bit. After each guess, the true next bit is effectively revealed because the candidate set is restricted to primes sharing the observed suffix. A correct guess scores one point. The goal is to maximize the expected total score, so the problem is really about making the best local decision at every partially revealed binary suffix. Mathematical Approach The natural state space is a binary trie built from reversed binary representations. Each node corresponds to a suffix of already revealed low bits, so every decision can be expressed in terms of how many remaining primes continue with bit \(0\) and how many continue with bit \(1\). Step 1: Encode Each Prime by Low Bits First Write a prime \(p\) in binary as $$p=\sum_{d=0}^{\ell(p)-1} b_d(p)\,2^d,\qquad b_d(p)\in\{0,1\},\qquad b_{\ell(p)-1}(p)=1.$$ The game reveals bits in the order \(b_0(p),b_1(p),b_2(p),\dots\), so a state is determined by an already known suffix \(s\) of low bits. A trie over these reversed bit strings stores exactly the primes compatible with each suffix....

Detailed mathematical approach

Problem Summary

Let \(P_N\) be the set of primes not exceeding \(N\), with \(N=10^8\). One prime is chosen uniformly at random. We then try to guess its binary expansion one bit at a time, starting from the least significant bit.

After each guess, the true next bit is effectively revealed because the candidate set is restricted to primes sharing the observed suffix. A correct guess scores one point. The goal is to maximize the expected total score, so the problem is really about making the best local decision at every partially revealed binary suffix.

Mathematical Approach

The natural state space is a binary trie built from reversed binary representations. Each node corresponds to a suffix of already revealed low bits, so every decision can be expressed in terms of how many remaining primes continue with bit \(0\) and how many continue with bit \(1\).

Step 1: Encode Each Prime by Low Bits First

Write a prime \(p\) in binary as

$$p=\sum_{d=0}^{\ell(p)-1} b_d(p)\,2^d,\qquad b_d(p)\in\{0,1\},\qquad b_{\ell(p)-1}(p)=1.$$

The game reveals bits in the order \(b_0(p),b_1(p),b_2(p),\dots\), so a state is determined by an already known suffix \(s\) of low bits. A trie over these reversed bit strings stores exactly the primes compatible with each suffix.

If \(S(s)\) denotes the set of primes whose reversed binary expansion begins with \(s\), then reaching node \(s\) means that the hidden prime is known to lie in \(S(s)\).

Step 2: Separate Continuing Primes from Finished Primes

At a trie node \(s\), let \(c(s)=|S(s)|\). Some primes end exactly at that node because their entire binary expansion has already been consumed. Let \(t(s)\) be that terminal count, and define

$$q(s)=c(s)-t(s).$$

Only those \(q(s)\) primes require another guess. If \(q(s)=0\), then the game is already over for every prime in that state, so the remaining expected score is \(0\).

Now look at the two child states \(s^{(0)}\) and \(s^{(1)}\), obtained by appending the next revealed bit. Let \(c_0(s)\) and \(c_1(s)\) be the numbers of compatible primes in those two children. Because every continuing prime must go to exactly one child, we have

$$c_0(s)+c_1(s)=q(s).$$

Step 3: The Optimal Guess Is the Local Majority Bit

Suppose we are at state \(s\) and must guess the next bit. Choosing \(0\) or \(1\) changes only the score of the current round. It does not change which child state is reached, because that depends on the hidden prime, not on our guess.

Therefore the optimal policy at node \(s\) is simply to guess the more common next bit among the continuing primes. The best immediate success probability is

$$\frac{\max(c_0(s),c_1(s))}{q(s)}.$$

This greedy choice is globally optimal because the future subproblem is the same regardless of which bit we guessed; only the current point is affected by the choice.

Step 4: Dynamic Programming Recurrence on Trie Nodes

Let \(E(s)\) be the maximum expected additional score once suffix \(s\) is known. For each child, some primes may stop immediately after that next bit. Let \(q_0(s)\) and \(q_1(s)\) be the numbers of primes that still continue beyond \(s^{(0)}\) and \(s^{(1)}\), respectively.

If \(q(s)=0\), then

$$E(s)=0.$$

If \(q(s)>0\), then

$$E(s)=\frac{\max(c_0(s),c_1(s))}{q(s)}+\frac{q_0(s)}{q(s)}E\!\left(s^{(0)}\right)+\frac{q_1(s)}{q(s)}E\!\left(s^{(1)}\right).$$

The first term is the expected point earned on the current guess. The remaining terms are weighted by the probabilities that the game continues into the two child states.

The required value for the problem is simply

$$E(\varnothing),$$

where \(\varnothing\) denotes the root, meaning that no bits have been revealed yet.

Worked Example: \(N=10\)

The primes are \(2,3,5,7\). Their reversed binary strings are

$$2\to 01,\qquad 3\to 11,\qquad 5\to 101,\qquad 7\to 111.$$

At the root, the next-bit counts are \(c_0=1\) and \(c_1=3\), so the optimal first guess is \(1\), worth \(3/4\) in expectation.

At the node reached after observing suffix \(1\), the next-bit counts are \(c_0=1\) and \(c_1=2\), so the optimal local score there is \(2/3\). From both of its relevant children, the continuation contributes exactly one more expected point.

Thus

$$E(1)=\frac{2}{3}+\frac{1}{3}\cdot 1+\frac{1}{3}\cdot 1=\frac{4}{3},$$

and the root value is

$$E(\varnothing)=\frac{3}{4}+\frac{1}{4}\cdot 1+\frac{3}{4}\cdot \frac{4}{3}=2.$$

This matches the small checkpoint used by the implementation. A second checkpoint in the implementations is \(N=30\), for which the same recurrence gives \(2.9\).

How the Code Works

The C++, Python, and Java implementations all follow the same two-phase plan. First they generate every prime up to \(10^8\) using an odd-only sieve, inserting \(2\) separately and then scanning only odd candidates.

Each prime is inserted into a binary trie from least significant bit to most significant bit. Every node records how many primes pass through it, whether a prime ends exactly there, and where the \(0\)-child and \(1\)-child are located.

Once the trie has been built, a depth-first traversal evaluates the recurrence bottom-up. At each node, the implementation reads the child counts, chooses the heavier child for the immediate guess, subtracts terminal primes when computing continuation probabilities, and combines those pieces into the expected score for that node.

Finally the root value is printed with fixed decimal precision. The three language versions differ only in storage details; mathematically they all compute the same trie dynamic program.

Complexity Analysis

Let \(N=10^8\). The sieve phase runs in \(O(N\log\log N)\) time. Inserting the primes into the trie costs

$$O\!\left(\sum_{p\in P_N} \ell(p)\right)=O\!\bigl(\pi(N)\log N\bigr),$$

because each prime contributes one trie step per binary digit. The final depth-first pass is linear in the number of trie nodes, so it is also \(O\!\bigl(\pi(N)\log N\bigr)\).

Memory usage is dominated by the odd-only sieve array together with the trie itself. In asymptotic form, the method uses \(O(N)\) sieve storage plus \(O\!\bigl(\pi(N)\log N\bigr)\) trie storage.

Footnotes and References

  1. Project Euler problem page: https://projecteuler.net/problem=869
  2. Trie data structure: Wikipedia - Trie
  3. Binary numeral system: Wikipedia - Binary number
  4. Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes
  5. Expected value and dynamic programming ideas: Wikipedia - Expected value and Wikipedia - Dynamic programming

Problem 869 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>

struct Node {
    int child[2];
    std::uint32_t cnt_term;
    Node() : child{-1, -1}, cnt_term(0) {}
};

class PrimeGuessing {
public:
    explicit PrimeGuessing(int limit) : limit_(limit) {
        nodes_.reserve(32'000'000);
        nodes_.push_back(Node());
        build();
    }

    long double expectation() {
        return dfs(0);
    }

private:
    static constexpr std::uint32_t TERM_BIT = 0x80000000u;
    static constexpr std::uint32_t CNT_MASK = 0x7fffffffu;

    int limit_;
    std::vector<Node> nodes_;

    static inline std::uint32_t count_of(const Node& n) {
        return n.cnt_term & CNT_MASK;
    }

    static inline std::uint32_t term_of(const Node& n) {
        return n.cnt_term >> 31;
    }

    inline void inc_count(int idx) {
        Node& n = nodes_[idx];
        std::uint32_t c = (n.cnt_term & CNT_MASK) + 1u;
        n.cnt_term = (n.cnt_term & TERM_BIT) | c;
    }

    inline void set_term(int idx) {
        nodes_[idx].cnt_term |= TERM_BIT;
    }

    void insert_prime(std::uint32_t p) {
        int u = 0;
        inc_count(u);

        int len = 32 - __builtin_clz(p);
        for (int d = 0; d < len; ++d) {
            int b = (p >> d) & 1u;
            int v = nodes_[u].child[b];
            if (v == -1) {
                v = static_cast<int>(nodes_.size());
                nodes_[u].child[b] = v;
                nodes_.push_back(Node());
            }
            u = v;
            inc_count(u);
        }
        set_term(u);
    }

    void build() {
        if (limit_ >= 2) insert_prime(2);
        if (limit_ < 3) return;

        int m = (limit_ >> 1) + 1;
        std::vector<std::uint8_t> is_prime(m, 1);
        is_prime[0] = 0;

        int r = static_cast<int>(std::sqrt(static_cast<long double>(limit_)));
        for (int p = 3; p <= r; p += 2) {
            if (!is_prime[p >> 1]) continue;
            int step = p << 1;
            for (int x = p * p; x <= limit_; x += step) {
                is_prime[x >> 1] = 0;
            }
        }

        for (int p = 3; p <= limit_; p += 2) {
            if (is_prime[p >> 1]) insert_prime(static_cast<std::uint32_t>(p));
        }
    }

    long double dfs(int u) {
        const Node& cur = nodes_[u];
        std::uint32_t total = count_of(cur);
        std::uint32_t term = term_of(cur);
        std::uint32_t cont = total - term;
        if (cont == 0) return 0.0L;

        std::uint32_t c0 = 0, c1 = 0;
        std::uint32_t n0 = 0, n1 = 0;
        long double v0 = 0.0L, v1 = 0.0L;

        int ch0 = cur.child[0];
        if (ch0 != -1) {
            c0 = count_of(nodes_[ch0]);
            n0 = c0 - term_of(nodes_[ch0]);
            if (n0 > 0) v0 = dfs(ch0);
        }

        int ch1 = cur.child[1];
        if (ch1 != -1) {
            c1 = count_of(nodes_[ch1]);
            n1 = c1 - term_of(nodes_[ch1]);
            if (n1 > 0) v1 = dfs(ch1);
        }

        long double denom = static_cast<long double>(cont);
        long double score = static_cast<long double>(std::max(c0, c1)) / denom;
        score += static_cast<long double>(n0) / denom * v0;
        score += static_cast<long double>(n1) / denom * v1;
        return score;
    }
};

int main() {
    {
        PrimeGuessing g(10);
        assert(std::fabsl(g.expectation() - 2.0L) < 1e-15L);
    }
    {
        PrimeGuessing g(30);
        assert(std::fabsl(g.expectation() - 2.9L) < 1e-15L);
    }

    PrimeGuessing g(100'000'000);
    std::cout << std::fixed << std::setprecision(8) << g.expectation() << '\n';
    return 0;
}

Python

import math

def solve():
    limit = 100000000

    class Node:
        __slots__ = ['child', 'cnt', 'term']
        def __init__(self):
            self.child = [-1, -1]; self.cnt = 0; self.term = 0

    nodes = [Node()]

    def insert(p):
        u = 0; nodes[u].cnt += 1
        l = p.bit_length()
        for d in range(l):
            b = (p >> d) & 1
            v = nodes[u].child[b]
            if v == -1:
                v = len(nodes); nodes[u].child[b] = v; nodes.append(Node())
            u = v; nodes[u].cnt += 1
        nodes[u].term = 1

    # Sieve and insert primes
    if limit >= 2: insert(2)
    m = (limit >> 1) + 1
    is_prime = bytearray(b'\x01' * m); is_prime[0] = 0
    r = int(math.isqrt(limit))
    for p in range(3, r+1, 2):
        if not is_prime[p >> 1]: continue
        for x in range(p*p, limit+1, 2*p): is_prime[x >> 1] = 0
    for p in range(3, limit+1, 2):
        if is_prime[p >> 1]: insert(p)

    def dfs(u):
        cur = nodes[u]
        cont = cur.cnt - cur.term
        if cont == 0: return 0.0
        c0 = c1 = n0 = n1 = 0; v0 = v1 = 0.0
        ch0 = cur.child[0]
        if ch0 != -1:
            c0 = nodes[ch0].cnt; n0 = c0 - nodes[ch0].term
            if n0 > 0: v0 = dfs(ch0)
        ch1 = cur.child[1]
        if ch1 != -1:
            c1 = nodes[ch1].cnt; n1 = c1 - nodes[ch1].term
            if n1 > 0: v1 = dfs(ch1)
        d = cont
        return max(c0, c1) / d + n0 / d * v0 + n1 / d * v1

    return f'{dfs(0):.8f}'

if __name__ == '__main__':
    print(solve())

Java

import java.util.Arrays;

public class Euler869 {

    static final int TERM_BIT = 0x80000000;
    static final int CNT_MASK = 0x7fffffff;

    static class PrimeGuessing {
        int limit;
        int[] child0;
        int[] child1;
        int[] cntTerm;
        int nodeCount;

        PrimeGuessing(int limit) {
            this.limit = limit;
            int capacity = 32000000;
            child0 = new int[capacity];
            child1 = new int[capacity];
            cntTerm = new int[capacity];
            Arrays.fill(child0, -1);
            Arrays.fill(child1, -1);
            nodeCount = 1; // root node
            build();
        }

        int addNode() {
            if (nodeCount == child0.length) {
                int newCap = child0.length * 2;
                child0 = Arrays.copyOf(child0, newCap);
                child1 = Arrays.copyOf(child1, newCap);
                cntTerm = Arrays.copyOf(cntTerm, newCap);
                for (int i = nodeCount; i < newCap; i++) {
                    child0[i] = -1;
                    child1[i] = -1;
                }
            }
            return nodeCount++;
        }

        void insertPrime(int p) {
            int u = 0;
            cntTerm[u] = (cntTerm[u] & TERM_BIT) | ((cntTerm[u] & CNT_MASK) + 1);

            int len = 32 - Integer.numberOfLeadingZeros(p);
            for (int d = 0; d < len; ++d) {
                int b = (p >> d) & 1;
                int v;
                if (b == 0) {
                    v = child0[u];
                    if (v == -1) {
                        v = addNode();
                        child0[u] = v;
                    }
                } else {
                    v = child1[u];
                    if (v == -1) {
                        v = addNode();
                        child1[u] = v;
                    }
                }

                u = v;
                cntTerm[u] = (cntTerm[u] & TERM_BIT) | ((cntTerm[u] & CNT_MASK) + 1);
            }
            cntTerm[u] |= TERM_BIT;
        }

        void build() {
            if (limit >= 2)
                insertPrime(2);
            if (limit < 3)
                return;

            int m = (limit >> 1) + 1;
            byte[] isPrime = new byte[m];
            Arrays.fill(isPrime, (byte) 1);
            isPrime[0] = 0;

            int r = (int) Math.sqrt(limit);
            for (int p = 3; p <= r; p += 2) {
                if (isPrime[p >> 1] == 0)
                    continue;
                int step = p << 1;
                for (int x = p * p; x <= limit && x > 0; x += step) {
                    isPrime[x >> 1] = 0;
                }
            }

            for (int p = 3; p <= limit; p += 2) {
                if (isPrime[p >> 1] == 1) {
                    insertPrime(p);
                }
            }
        }

        double dfs(int u) {
            int ct = cntTerm[u];
            int total = ct & CNT_MASK;
            int term = (ct >>> 31);
            int cont = total - term;
            if (cont == 0)
                return 0.0;

            int c0 = 0, n0 = 0;
            double v0 = 0.0;

            int ch0 = child0[u];
            if (ch0 != -1) {
                int ct0 = cntTerm[ch0];
                c0 = ct0 & CNT_MASK;
                n0 = c0 - (ct0 >>> 31);
                if (n0 > 0)
                    v0 = dfs(ch0);
            }

            int c1 = 0, n1 = 0;
            double v1 = 0.0;

            int ch1 = child1[u];
            if (ch1 != -1) {
                int ct1 = cntTerm[ch1];
                c1 = ct1 & CNT_MASK;
                n1 = c1 - (ct1 >>> 31);
                if (n1 > 0)
                    v1 = dfs(ch1);
            }

            double denom = (double) cont;
            double score = Math.max(c0, c1) / denom;
            score += (n0 / denom) * v0;
            score += (n1 / denom) * v1;
            return score;
        }

        double expectation() {
            return dfs(0);
        }
    }

    public static String solve() {
        PrimeGuessing pg = new PrimeGuessing(100000000);
        return String.format(java.util.Locale.US, "%.8f", pg.expectation());
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}