Problem 333: Special Partitions

View on Project Euler

Project Euler Problem 333 Solution

EulerSolve provides an optimized solution for Project Euler Problem 333, Special Partitions, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We study partitions of an integer \(n\) into terms of the form \(2^a3^b\) with \(a,b \ge 0\), under the extra restriction that no chosen term may divide any other chosen term. If \(P(n)\) denotes the number of such valid partitions, the target is $$\sum_{\substack{q \lt 10^6 \\ q\text{ prime} \\ P(q)=1}} q.$$ Mathematical Approach Let \(L = 10^6 - 1\). The solution does not test all partitions of all integers up to \(L\). Instead, it converts the validity rule into a strict order constraint on the exponent pairs \((a,b)\). Admissible Terms Every allowed part has the form \(2^a3^b\). For each fixed exponent \(a\), define the row $$R_a = \{2^a3^b \le L : b \ge 0\}.$$ The code stores these rows as rows[a][b] , so row \(a\) is the increasing list $$2^a,\;2^a3,\;2^a3^2,\;\dots.$$ This already compresses the search space into a small rectangular grid of exponent pairs instead of all integers below \(L\). Divisibility Criterion For two admissible terms \(x = 2^{a_1}3^{b_1}\) and \(y = 2^{a_2}3^{b_2}\), unique prime factorization gives $$x \mid y \iff a_1 \le a_2 \text{ and } b_1 \le b_2.$$ So a valid partition is exactly a finite set of lattice points \((a,b)\) that forms an antichain in the product order on exponent pairs....

Detailed mathematical approach

Problem Summary

We study partitions of an integer \(n\) into terms of the form \(2^a3^b\) with \(a,b \ge 0\), under the extra restriction that no chosen term may divide any other chosen term. If \(P(n)\) denotes the number of such valid partitions, the target is

$$\sum_{\substack{q \lt 10^6 \\ q\text{ prime} \\ P(q)=1}} q.$$

Mathematical Approach

Let \(L = 10^6 - 1\). The solution does not test all partitions of all integers up to \(L\). Instead, it converts the validity rule into a strict order constraint on the exponent pairs \((a,b)\).

Admissible Terms

Every allowed part has the form \(2^a3^b\). For each fixed exponent \(a\), define the row

$$R_a = \{2^a3^b \le L : b \ge 0\}.$$

The code stores these rows as rows[a][b], so row \(a\) is the increasing list

$$2^a,\;2^a3,\;2^a3^2,\;\dots.$$

This already compresses the search space into a small rectangular grid of exponent pairs instead of all integers below \(L\).

Divisibility Criterion

For two admissible terms \(x = 2^{a_1}3^{b_1}\) and \(y = 2^{a_2}3^{b_2}\), unique prime factorization gives

$$x \mid y \iff a_1 \le a_2 \text{ and } b_1 \le b_2.$$

So a valid partition is exactly a finite set of lattice points \((a,b)\) that forms an antichain in the product order on exponent pairs.

This has two immediate consequences:

First, a valid partition can contain at most one term from any fixed row \(a\), because inside one row the smaller \(b\)-index always divides the larger one.

Second, if rows are processed in increasing \(a\), then the chosen \(b\)-indices must be strictly decreasing. Indeed, if \(a_1 < a_2\) and \(b_1 \le b_2\), then \(2^{a_1}3^{b_1}\mid 2^{a_2}3^{b_2}\), which is forbidden.

DFS State and Invariant

After this reduction, the counting problem becomes a row-by-row depth-first search. The state

$$\mathrm{DFS}(i, j_{\max}, s)$$

means that rows \(0,1,\dots,i-1\) have already been processed, the current partial sum is \(s\), and every future chosen column index \(j\) must satisfy \(j < j_{\max}\). This invariant encodes the strict decrease of the \(3\)-exponents.

From a state \((i,j_{\max},s)\), the code has exactly two kinds of transitions:

Skip row \(i\):

$$\mathrm{DFS}(i+1, j_{\max}, s).$$

Choose the \(j\)-th term of row \(i\), where \(0 \le j < \min(j_{\max}, |R_i|)\) and \(s + 2^i3^j \le L\):

$$\mathrm{DFS}(i+1, j, s + 2^i3^j).$$

The parameter j_limit in the source is exactly this \(j_{\max}\).

Why This Counts Each Valid Partition Exactly Once

Every valid partition has a unique description when its terms are ordered by increasing \(a\). In each row we either take nothing or choose exactly one column index \(b\). Therefore each valid partition corresponds to one unique DFS path.

Conversely, every DFS leaf represents a set of terms with strictly decreasing \(b\)-indices, hence no term can divide another. Therefore the leaf update

$$P(s) \leftarrow P(s) + 1$$

records the exact number of valid partitions of \(s\), not an approximation.

Worked Examples: \(P(11)=2\) and \(P(17)=1\)

The code checkpoints match the examples from the statement. For 11, the two valid partitions are

$$11 = 2 + 9 = 2^1 3^0 + 2^0 3^2,$$

$$11 = 8 + 3 = 2^3 3^0 + 2^0 3^1,$$

so \(P(11)=2\).

For 17, the partition \(17 = 8 + 9\) is valid, but \(17 = 2 + 6 + 9\) is invalid because \(2 \mid 6\), and \(17 = 16 + 1\) is invalid because \(1 \mid 16\). Hence

$$P(17)=1.$$

Prime Filter and Final Sum

Once the array \(P(s)\) has been filled for all \(s \le L\), the rest is straightforward. A sieve of Eratosthenes marks the primes, and the code returns

$$\sum_{\substack{q \lt 10^6 \\ q\text{ prime} \\ P(q)=1}} q.$$

The published sample

$$\sum_{\substack{q \lt 100 \\ q\text{ prime} \\ P(q)=1}} q = 233$$

is used as a correctness checkpoint.

How the Code Works

build_terms(limit) constructs the rows \(R_a\). Then partition_counts(limit) runs the DFS described above and fills count[s] with the exact value \(P(s)\). Finally solve(prime_limit) computes a prime sieve and sums precisely those primes \(q\) for which count[q] == 1. The C++, Python, and Java implementations are the same algorithm in three syntaxes.

Complexity Analysis

The count array and prime sieve require \(O(L)\) memory. The recursion depth is \(O(\log_2 L)\), because there is one DFS level per power of 2. The row table itself is tiny, of size only \(O(\log L \cdot \log L)\).

The running time is output-sensitive: the DFS does not examine all subsets of all admissible terms, but only row-compatible selections with strictly decreasing \(b\)-indices, while also pruning any branch whose partial sum already exceeds \(L\). In practice, the cost is governed by the number of reachable valid states, which is far smaller than a naive subset search.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=333
  2. Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes
  3. Backtracking: Wikipedia — Backtracking
  4. Antichains and partial orders: Wikipedia — Antichain

Problem 333 source code

C++

#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    int prime_limit = 1000000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--prime-limit=", options.prime_limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.prime_limit >= 3;
}

std::vector<std::vector<int>> build_terms(const int limit) {
    std::vector<std::vector<int>> rows;
    for (u64 p2 = 1ULL; p2 <= static_cast<u64>(limit); p2 <<= 1ULL) {
        std::vector<int> row;
        for (u64 value = p2; value <= static_cast<u64>(limit); value *= 3ULL) {
            row.push_back(static_cast<int>(value));
            if (value > static_cast<u64>(limit) / 3ULL) {
                break;
            }
        }
        rows.push_back(std::move(row));
    }
    return rows;
}

std::vector<std::uint32_t> partition_counts(const int limit) {
    const std::vector<std::vector<int>> rows = build_terms(limit);
    std::vector<std::uint32_t> count(static_cast<std::size_t>(limit + 1), 0U);

    const auto dfs = [&](auto&& self, const int i, const int j_limit, const int sum) -> void {
        if (i == static_cast<int>(rows.size())) {
            if (sum > 0) {
                ++count[static_cast<std::size_t>(sum)];
            }
            return;
        }

        self(self, i + 1, j_limit, sum);  // Skip this 2^i layer.

        const int upper = std::min(j_limit, static_cast<int>(rows[static_cast<std::size_t>(i)].size()));
        for (int j = 0; j < upper; ++j) {
            const int value = rows[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)];
            if (sum + value > limit) {
                break;
            }
            self(self, i + 1, j, sum + value);
        }
    };
    dfs(dfs, 0, 1000, 0);
    return count;
}

u64 solve(const int prime_limit) {
    const int limit = prime_limit - 1;
    const std::vector<std::uint32_t> counts = partition_counts(limit);

    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(prime_limit), 1U);
    is_prime[0] = 0U;
    is_prime[1] = 0U;
    for (int p = 2; static_cast<long long>(p) * p < prime_limit; ++p) {
        if (is_prime[static_cast<std::size_t>(p)] == 0U) {
            continue;
        }
        for (int q = p * p; q < prime_limit; q += p) {
            is_prime[static_cast<std::size_t>(q)] = 0U;
        }
    }

    u64 sum = 0ULL;
    for (int q = 2; q < prime_limit; ++q) {
        if (is_prime[static_cast<std::size_t>(q)] != 0U &&
            counts[static_cast<std::size_t>(q)] == 1U) {
            sum += static_cast<u64>(q);
        }
    }
    return sum;
}

bool run_checkpoints() {
    if (solve(100) != 233ULL) {
        std::cerr << "Checkpoint failed for stated sample q<100" << '\n';
        return false;
    }
    const auto counts = partition_counts(100);
    if (counts[11] != 2U || counts[17] != 1U) {
        std::cerr << "Checkpoint failed for P(11)=2 and P(17)=1 examples" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << solve(options.prime_limit) << '\n';
    return 0;
}

Python

def solve():
    prime_limit = 1000000
    limit = prime_limit - 1

    # Build terms: rows of 2^a * 3^b values
    rows = []
    p2 = 1
    while p2 <= limit:
        row = []
        v = p2
        while v <= limit:
            row.append(v)
            if v > limit // 3:
                break
            v *= 3
        rows.append(row)
        p2 <<= 1

    # DFS partition counting
    count = [0] * (limit + 1)

    def dfs(i, j_limit, s):
        if i == len(rows):
            if s > 0:
                count[s] += 1
            return
        dfs(i + 1, j_limit, s)
        upper = min(j_limit, len(rows[i]))
        for j in range(upper):
            val = rows[i][j]
            if s + val > limit:
                break
            dfs(i + 1, j, s + val)

    import sys
    sys.setrecursionlimit(100000)
    dfs(0, 1000, 0)

    # Sieve primes
    is_prime = bytearray(b'\x01' * prime_limit)
    is_prime[0] = 0
    is_prime[1] = 0
    p = 2
    while p * p < prime_limit:
        if is_prime[p]:
            is_prime[p*p::p] = bytearray(len(is_prime[p*p::p]))
        p += 1

    total = 0
    for q in range(2, prime_limit):
        if is_prime[q] and count[q] == 1:
            total += q

    return str(total)

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

Java

import java.util.*;

public class Euler333 {

    static boolean[] sievePrimeFlags(int limit) {
        boolean[] isPrime = new boolean[limit];
        Arrays.fill(isPrime, true);
        if (limit > 0)
            isPrime[0] = false;
        if (limit > 1)
            isPrime[1] = false;
        for (int p = 2; p * p < limit; p++) {
            if (isPrime[p]) {
                for (int q = p * p; q < limit; q += p) {
                    isPrime[q] = false;
                }
            }
        }
        return isPrime;
    }

    static List<List<Integer>> buildTerms(int limit) {
        List<List<Integer>> rows = new ArrayList<>();
        for (long p2 = 1; p2 <= limit; p2 *= 2) {
            List<Integer> row = new ArrayList<>();
            for (long val = p2; val <= limit; val *= 3) {
                row.add((int) val);
                if (val > (limit / 3)) {
                    break;
                }
            }
            rows.add(row);
        }
        return rows;
    }

    static void dfs(int i, int j_limit, int sum, List<List<Integer>> rows, int[] count, int limit) {
        if (i == rows.size()) {
            if (sum > 0)
                count[sum]++;
            return;
        }

        dfs(i + 1, j_limit, sum, rows, count, limit);

        int upper = Math.min(j_limit, rows.get(i).size());
        for (int j = 0; j < upper; j++) {
            int val = rows.get(i).get(j);
            if (sum + val > limit)
                break;
            dfs(i + 1, j, sum + val, rows, count, limit);
        }
    }

    static int[] partitionCounts(int limit) {
        List<List<Integer>> rows = buildTerms(limit);
        int[] count = new int[limit + 1];
        dfs(0, 1000, 0, rows, count, limit);
        return count;
    }

    public static String solve() {
        int primeLimit = 1000000;
        int limit = primeLimit - 1;
        int[] counts = partitionCounts(limit);
        boolean[] isPrime = sievePrimeFlags(primeLimit);

        long sum = 0;
        for (int q = 2; q < primeLimit; q++) {
            if (isPrime[q] && counts[q] == 1) {
                sum += q;
            }
        }
        return String.valueOf(sum);
    }

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