Problem 242: Odd Triplets

Project Euler Problem 242 Solution

EulerSolve provides an optimized solution for Project Euler Problem 242, Odd Triplets, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each pair $(n,k)$, let $f(n,k)$ be the number of $k$-element subsets of $\{1,2,\dots,n\}$ whose element sum is odd. An odd triplet is a triplet $(n,k,f(n,k))$ in which all three entries are odd. The quantity to evaluate is $$A(N)=\#\{(n,k):1\le k\le n\le N,\ n,\ k,\ f(n,k)\text{ are all odd}\},$$ and the implementations compute $A(10^{12})$. The search space is far too large for direct enumeration, so the key is to understand only the parity of $f(n,k)$, not its exact size. Mathematical Approach The derivation starts from a straightforward combinatorial count, turns that count into a parity statement over modulo 2 arithmetic, and then collapses the whole problem to a short prefix sum involving binary popcounts. Separating odd and even elements Let $$o=\left\lceil\frac n2\right\rceil,\qquad e=\left\lfloor\frac n2\right\rfloor.$$ Any $k$-subset of $\{1,\dots,n\}$ has odd total exactly when it contains an odd number of odd elements. Therefore $$f(n,k)=\sum_{\substack{j=0\\ j\text{ odd}}}^{k}\binom{o}{j}\binom{e}{k-j}.$$ Because an odd triplet must have odd first and second coordinates, we only need $n=2m+1$ and $k=2r+1$....

Detailed mathematical approach

Problem Summary

For each pair $(n,k)$, let $f(n,k)$ be the number of $k$-element subsets of $\{1,2,\dots,n\}$ whose element sum is odd. An odd triplet is a triplet $(n,k,f(n,k))$ in which all three entries are odd.

The quantity to evaluate is

$$A(N)=\#\{(n,k):1\le k\le n\le N,\ n,\ k,\ f(n,k)\text{ are all odd}\},$$

and the implementations compute $A(10^{12})$. The search space is far too large for direct enumeration, so the key is to understand only the parity of $f(n,k)$, not its exact size.

Mathematical Approach

The derivation starts from a straightforward combinatorial count, turns that count into a parity statement over modulo 2 arithmetic, and then collapses the whole problem to a short prefix sum involving binary popcounts.

Separating odd and even elements

Let

$$o=\left\lceil\frac n2\right\rceil,\qquad e=\left\lfloor\frac n2\right\rfloor.$$

Any $k$-subset of $\{1,\dots,n\}$ has odd total exactly when it contains an odd number of odd elements. Therefore

$$f(n,k)=\sum_{\substack{j=0\\ j\text{ odd}}}^{k}\binom{o}{j}\binom{e}{k-j}.$$

Because an odd triplet must have odd first and second coordinates, we only need $n=2m+1$ and $k=2r+1$. Then $o=m+1$ and $e=m$, so

$$f(2m+1,2r+1)=\sum_{\substack{j=0\\ j\text{ odd}}}^{2r+1}\binom{m+1}{j}\binom{m}{2r+1-j}.$$

A parity-generating polynomial modulo 2

To study only parity, encode subset size by $x$ and the parity of the subset sum by $y$. Modulo 2, the relevant polynomial is

$$P_n(x,y)=\prod_{i=1}^{n}\left(1+x\,y^{i\bmod 2}\right)=(1+x)^e(1+xy)^o,$$

with the rule $y^2=1$. The coefficient of $x^k y$ is exactly $f(n,k)\bmod 2$, because each chosen odd number contributes one factor of $y$, and only the parity of the number of chosen odd elements matters.

For $n=2m+1$, this becomes

$$P_{2m+1}(x,y)=(1+x)^m(1+xy)^{m+1}.$$

Now split into the parity of $m$.

Why only $n\equiv 1 \pmod 4$ can contribute

If $m=2t+1$ is odd, then modulo 2 we can square factors:

$$ (1+x)^{2t+1}=(1+x)(1+x^2)^t,\qquad (1+xy)^{2t+2}=(1+x^2)^{t+1}. $$

Hence

$$P_{4t+3}(x,y)=(1+x)(1+x^2)^{2t+1},$$

which contains no $y$-term at all. Therefore every $f(4t+3,2r+1)$ is even, so these values of $n$ never produce odd triplets.

If instead $m=2t$ is even, then

$$ (1+x)^{2t}=(1+x^2)^t,\qquad (1+xy)^{2t+1}=(1+x^2)^t(1+xy), $$

and therefore

$$P_{4t+1}(x,y)=(1+x^2)^{2t}(1+xy).$$

The coefficient of $x^{2r+1}y$ is now the coefficient of $x^{2r}$ in $(1+x^2)^{2t}$, namely

$$f(4t+1,2r+1)\equiv \binom{2t}{r}\pmod 2.$$

Equivalently, with $m=(n-1)/2$,

$$f(2m+1,2r+1)\text{ is odd}\iff m\text{ is even and }\binom{m}{r}\text{ is odd}.$$

Lucas' theorem and the surviving values of $k$

Modulo 2, Lucas' theorem says

$$\binom{u}{v}\equiv 1\pmod 2\iff (v\ \&\ \sim u)=0.$$

In words: every 1-bit of $v$ must sit under a 1-bit of $u$. So if $m=2t$, the admissible values of $r$ are exactly the binary submasks of $m$.

The number of such $r$ is the number of odd entries in row $m$ of Pascal's triangle:

$$\#\{r:\binom{m}{r}\text{ odd}\}=2^{\operatorname{popcount}(m)}.$$

Since $m=2t$ only appends a trailing zero in binary, $\operatorname{popcount}(m)=\operatorname{popcount}(t)$. Therefore each $t$ contributes

$$2^{\operatorname{popcount}(t)}$$

odd triplets.

Collapsing the global count

The contributing values of $n$ are exactly $n=4t+1$. For a bound $N$, the largest possible $t$ is

$$t_{\max}=\left\lfloor\frac{N-1}{4}\right\rfloor.$$

So the entire Project Euler question reduces to

$$A(N)=\sum_{t=0}^{t_{\max}}2^{\operatorname{popcount}(t)}.$$

This is already a complete characterization of all odd triplets: every contributing $n$ has the form $4t+1$, and for that $n$, the valid odd $k$ are precisely those with $k=2r+1$ where $r$ is a binary submask of $2t$.

Worked example: counting odd triplets up to $N=10$

Here

$$t_{\max}=\left\lfloor\frac{10-1}{4}\right\rfloor=2.$$

So we only inspect $t=0,1,2$:

$$ 2^{\operatorname{popcount}(0)}=1,\qquad 2^{\operatorname{popcount}(1)}=2,\qquad 2^{\operatorname{popcount}(2)}=2. $$

The total is $1+2+2=5$. Concretely, the contributing values are:

$t=0 \Rightarrow n=1$, giving $k=1$.

$t=1 \Rightarrow n=5$, where row 2 of Pascal's triangle has odd entries at $r=0,2$, giving $k=1,5$.

$t=2 \Rightarrow n=9$, where row 4 has odd entries at $r=0,4$, giving $k=1,9$.

That matches the small checkpoint used by the reference implementation: there are exactly 5 odd triplets with $n\le 10$.

Evaluating the prefix sum in $O(\log N)$

Define

$$S(T)=\sum_{t=0}^{T}2^{\operatorname{popcount}(t)}.$$

The implementation scans the bits of $T$ from most significant to least significant. Suppose the current bit is $b$, this bit of $T$ is 1, and there are already $s$ ones in the higher prefix. If we set bit $b$ to 0 and allow the lower $b$ bits to vary freely, every lower pattern $x$ contributes

$$2^{s+\operatorname{popcount}(x)}.$$

Summing over all $x\in[0,2^b)$ gives

$$2^s\sum_{x=0}^{2^b-1}2^{\operatorname{popcount}(x)}=2^s\cdot 3^b,$$

because each lower bit independently contributes either a factor $1$ or a factor $2$, so the total product is $(1+2)^b=3^b$. After processing every set bit, one final term $2^s$ accounts for $T$ itself.

This transforms the whole problem into a constant-memory bit walk over at most 64 positions.

How the Code Works

Reducing the bound

The C++, Python, and Java implementations first replace the original upper bound $N$ by

$$T=\left\lfloor\frac{N-1}{4}\right\rfloor,$$

because only $n=4t+1$ can contribute. From that point on, the task is just to compute $S(T)=\sum_{t=0}^{T}2^{\operatorname{popcount}(t)}$.

Prefix accumulation with powers of three

All three implementations precompute $3^b$ for bit positions $0$ through $63$. They then scan the bits of $T$ from high to low, keep track of how many 1-bits have already appeared in the prefix, and whenever a set bit is encountered they add the contribution $2^s3^b$ for the branch in which that bit is turned off and the remaining suffix is arbitrary.

After the scan finishes, the remaining prefix itself contributes $2^s$. That exactly reproduces the mathematical prefix-sum formula above, and it is why the final solver is extremely short in every language.

Validation strategy in the reference implementation

The C++ implementation also contains explicit correctness checks on small inputs. It verifies a direct sample value $f(5,3)=4$, checks that the total number of odd triplets up to $n=10$ is 5, compares the closed parity rule against direct combinatorial evaluation for all odd $n\le 1001$, and tests the prefix-sum routine against a brute-force accumulation on a large range of small $t$.

The Python and Java implementations keep only the compact production calculation, but they implement exactly the same mathematics.

Complexity Analysis

The actual solver performs one scan over the bit representation of $T=\lfloor(N-1)/4\rfloor$. For a 64-bit bound this means at most 64 iterations, so the time complexity is $O(\log N)$ and the memory usage is $O(1)$.

The optional brute-force and parity-validation routines in the C++ version are only small safety checks; they are not part of the main Project Euler computation.

Footnotes and References

Problem page: https://projecteuler.net/problem=242
Lucas' theorem: Wikipedia - Lucas' theorem
Binomial coefficient: Wikipedia - Binomial coefficient
Pascal's triangle modulo 2 and the Sierpinski pattern: Wikipedia - Sierpinski triangle
Generating functions: Wikipedia - Generating function

Problem 242 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 kDefaultLimit = 1'000'000'000'000ULL;

struct Options {
    u64 limit = kDefaultLimit;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }

    std::string result;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        result.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }

    std::reverse(result.begin(), result.end());
    return result;
}

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const 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 (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--n=", parsed_u64)) {
            options.limit = parsed_u64;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload) {
    if (!allow_multithreading || workload < 8ULL) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}

bool binom_is_odd(const u64 n, const u64 k) {
    return k <= n && ((k & ~n) == 0ULL);
}

u128 binom_exact(const u64 n, u64 k) {
    if (k > n) {
        return 0;
    }
    if (k > n - k) {
        k = n - k;
    }

    u128 value = 1;
    for (u64 i = 1; i <= k; ++i) {
        value = (value * static_cast<u128>(n - k + i)) / static_cast<u128>(i);
    }
    return value;
}

u128 f_exact(const u64 n, const u64 k) {
    const u64 odd_count = (n + 1ULL) / 2ULL;
    const u64 even_count = n / 2ULL;

    u128 total = 0;
    for (u64 odd_taken = 1ULL; odd_taken <= k && odd_taken <= odd_count; odd_taken += 2ULL) {
        const u64 even_taken = k - odd_taken;
        if (even_taken > even_count) {
            continue;
        }
        total += binom_exact(odd_count, odd_taken) * binom_exact(even_count, even_taken);
    }
    return total;
}

bool f_is_odd_direct(const u64 n, const u64 k) {
    if (k > n) {
        return false;
    }

    const u64 odd_count = (n + 1ULL) / 2ULL;
    const u64 even_count = n / 2ULL;

    bool parity = false;
    for (u64 odd_taken = 1ULL; odd_taken <= k && odd_taken <= odd_count; odd_taken += 2ULL) {
        const u64 even_taken = k - odd_taken;
        if (even_taken > even_count) {
            continue;
        }

        if (binom_is_odd(odd_count, odd_taken) && binom_is_odd(even_count, even_taken)) {
            parity = !parity;
        }
    }
    return parity;
}

bool f_is_odd_formula(const u64 n, const u64 k) {
    if ((n & 1ULL) == 0ULL || (k & 1ULL) == 0ULL || k > n) {
        return false;
    }

    const u64 m = (n - 1ULL) / 2ULL;
    if ((m & 1ULL) != 0ULL) {
        return false;
    }

    const u64 r = (k - 1ULL) / 2ULL;
    return binom_is_odd(m, r);
}

u64 count_odd_triplets_bruteforce(const u64 limit) {
    u64 total = 0;
    for (u64 n = 1ULL; n <= limit; n += 2ULL) {
        for (u64 k = 1ULL; k <= n; k += 2ULL) {
            if (f_is_odd_direct(n, k)) {
                ++total;
            }
        }
    }
    return total;
}

const std::vector<u128>& powers_of_three() {
    static const std::vector<u128> table = [] {
        std::vector<u128> values(65, 0);
        values[0] = 1;
        for (std::size_t i = 1; i < values.size(); ++i) {
            values[i] = values[i - 1] * static_cast<u128>(3);
        }
        return values;
    }();
    return table;
}

u128 weighted_popcount_prefix_sum(const u64 n) {
    const auto& pow3 = powers_of_three();

    u128 result = 0;
    unsigned ones_so_far = 0U;

    for (int bit = 63; bit >= 0; --bit) {
        if (((n >> bit) & 1ULL) == 0ULL) {
            continue;
        }

        result += (static_cast<u128>(1) << ones_so_far) * pow3[static_cast<std::size_t>(bit)];
        ++ones_so_far;
    }

    result += (static_cast<u128>(1) << ones_so_far);
    return result;
}

u128 solve(const u64 limit) {
    if (limit == 0ULL) {
        return 0;
    }

    // n = 2m+1 and only even m contribute, so m = 2t.
    const u64 m_max = (limit - 1ULL) / 2ULL;
    const u64 t_max = m_max / 2ULL;
    return weighted_popcount_prefix_sum(t_max);
}

bool run_checkpoints(const Options& options) {
    bool ok = true;

    if (f_exact(5, 3) != static_cast<u128>(4)) {
        std::cerr << "Checkpoint failed: f(5,3) should be 4.\n";
        ok = false;
    }

    const u64 sample_count = count_odd_triplets_bruteforce(10ULL);
    if (sample_count != 5ULL) {
        std::cerr << "Checkpoint failed: odd-triplets up to n=10 should be 5.\n";
        ok = false;
    }

    const u64 parity_check_max_n = 1001ULL;
    const std::size_t workload = static_cast<std::size_t>((parity_check_max_n + 1ULL) / 2ULL);
    const unsigned thread_count =
        choose_thread_count(options.allow_multithreading, options.requested_threads, workload);

    std::atomic<bool> parity_ok(true);
    std::atomic<u64> next_n(1ULL);

    auto worker = [&]() {
        while (parity_ok.load(std::memory_order_relaxed)) {
            const u64 n = next_n.fetch_add(2ULL, std::memory_order_relaxed);
            if (n > parity_check_max_n) {
                break;
            }

            for (u64 k = 1ULL; k <= n; k += 2ULL) {
                if (f_is_odd_direct(n, k) != f_is_odd_formula(n, k)) {
                    parity_ok.store(false, std::memory_order_relaxed);
                    return;
                }
            }
        }
    };

    std::vector<std::thread> pool;
    pool.reserve(thread_count);
    for (unsigned i = 0U; i < thread_count; ++i) {
        pool.emplace_back(worker);
    }
    for (auto& th : pool) {
        th.join();
    }

    if (!parity_ok.load()) {
        std::cerr << "Checkpoint failed: parity identity mismatch.\n";
        ok = false;
    }

    u128 running = 0;
    constexpr u64 kPrefixCheckMax = 200000ULL;
    for (u64 x = 0ULL; x <= kPrefixCheckMax; ++x) {
        running += (static_cast<u128>(1) << __builtin_popcountll(x));
        if (weighted_popcount_prefix_sum(x) != running) {
            std::cerr << "Checkpoint failed: weighted prefix mismatch at x=" << x << ".\n";
            ok = false;
            break;
        }
    }

    if (ok) {
        std::cout << "Checkpoints passed (threads=" << thread_count << ").\n";
    }

    return ok;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints(options)) {
        return 1;
    }

    const u128 answer = solve(options.limit);
    std::cout << "Answer: " << to_string_u128(answer) << '\n';
    return 0;
}

Python

def solve():
    limit = 1000000000000
    
    # Precompute powers of 3
    pow3 = [1] * 65
    for i in range(1, 65):
        pow3[i] = pow3[i-1] * 3
        
    def weighted_popcount_prefix_sum(n):
        result = 0
        ones_so_far = 0
        
        for bit in range(63, -1, -1):
            if ((n >> bit) & 1) == 0:
                continue
                
            result += (1 << ones_so_far) * pow3[bit]
            ones_so_far += 1
            
        result += (1 << ones_so_far)
        return result
        
    m_max = (limit - 1) // 2
    t_max = m_max // 2
    
    ans = weighted_popcount_prefix_sum(t_max)
    return str(ans)

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

Java

public class Euler242 {
    static long[] pow3 = new long[65];
    static {
        pow3[0] = 1;
        for (int i = 1; i < 65; i++) {
            pow3[i] = pow3[i - 1] * 3;
        }
    }

    static long weightedPopcountPrefixSum(long n) {
        long result = 0;
        int onesSoFar = 0;

        for (int bit = 63; bit >= 0; --bit) {
            if (((n >> bit) & 1) == 0) {
                continue;
            }

            result += (1L << onesSoFar) * pow3[bit];
            onesSoFar++;
        }

        result += (1L << onesSoFar);
        return result;
    }

    public static String solve() {
        long limit = 1_000_000_000_000L;
        long mMax = (limit - 1) / 2;
        long tMax = mMax / 2;

        long ans = weightedPopcountPrefixSum(tMax);
        return String.valueOf(ans);
    }

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