Problem 239: Twenty-two Foolish Primes

Project Euler Problem 239 Solution

EulerSolve provides an optimized solution for Project Euler Problem 239, Twenty-two Foolish Primes, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We consider a uniformly random permutation of the numbers $1$ through $100$. Among them, the marked objects are the $25$ primes in that range. In the language of this problem, a prime is "foolish" when it does not return to its own position under the permutation. The target event is therefore very specific: exactly $22$ of those $25$ prime-labelled items are misplaced, so exactly $3$ primes are fixed points. The implementations are written for general parameters $(n,\;m,\;f)$, but for Problem 239 the relevant values are $n=100$, $m=25$, and $f=3$. Mathematical Approach The whole problem is a fixed-point counting question on a distinguished subset. Nothing special is required from the non-prime labels; once the prime labels satisfy the condition, every remaining label may be arranged freely. That observation is what makes a clean inclusion-exclusion count possible. Prime-labelled fixed points are the only constrained objects Let $n$ be the total number of labels, let $m$ be the number of marked labels, and let $f$ be the number of marked labels that are allowed to stay fixed. For Problem 239, those marked labels are precisely the primes, so $$n=100,\qquad m=25,\qquad f=25-22=3.$$ We want permutations in which exactly $f$ marked labels are fixed points and the remaining $m-f$ marked labels are not....

Detailed mathematical approach

Problem Summary

We consider a uniformly random permutation of the numbers $1$ through $100$. Among them, the marked objects are the $25$ primes in that range. In the language of this problem, a prime is "foolish" when it does not return to its own position under the permutation.

The target event is therefore very specific: exactly $22$ of those $25$ prime-labelled items are misplaced, so exactly $3$ primes are fixed points. The implementations are written for general parameters $(n,\;m,\;f)$, but for Problem 239 the relevant values are $n=100$, $m=25$, and $f=3$.

Mathematical Approach

The whole problem is a fixed-point counting question on a distinguished subset. Nothing special is required from the non-prime labels; once the prime labels satisfy the condition, every remaining label may be arranged freely. That observation is what makes a clean inclusion-exclusion count possible.

Prime-labelled fixed points are the only constrained objects

Let $n$ be the total number of labels, let $m$ be the number of marked labels, and let $f$ be the number of marked labels that are allowed to stay fixed. For Problem 239, those marked labels are precisely the primes, so

$$n=100,\qquad m=25,\qquad f=25-22=3.$$

We want permutations in which exactly $f$ marked labels are fixed points and the remaining $m-f$ marked labels are not. The non-marked labels are unrestricted except for filling the remaining places.

Choose which primes stay fixed

The first choice is purely combinatorial: decide which $f$ marked labels are the ones that remain in their original positions. There are

$$\binom{m}{f}$$

ways to make that choice. After fixing those labels, both their positions and their values are no longer part of the count, so the problem shrinks to the remaining $n-f$ slots.

Inclusion-exclusion on the remaining marked labels

Now let

$$b=m-f$$

be the number of marked labels that must avoid their own positions. In the actual problem, $b=22$. If we ignore that restriction for a moment, there are $(n-f)!$ ways to arrange the remaining labels.

For inclusion-exclusion, choose a set of $j$ among those $b$ forbidden marked labels and pretend that all of them are fixed after all. Once those $j$ extra positions are forced, the remaining objects can be permuted arbitrarily, giving

$$ (n-f-j)! $$

arrangements. There are $\binom{b}{j}$ ways to choose the set, so the number of valid completions is

$$\sum_{j=0}^{b} (-1)^j \binom{b}{j}(n-f-j)!.$$

Multiplying by the earlier choice of which marked labels stay fixed gives the exact count

$$A(n,m,f)=\binom{m}{f}\sum_{j=0}^{m-f} (-1)^j \binom{m-f}{j}(n-f-j)!.$$

The closed formula for Problem 239

Since all $n!$ permutations are equally likely, the desired probability is just $A(n,m,f)/n!$. For this problem that becomes

$$\boxed{P=\binom{25}{3}\sum_{j=0}^{22} (-1)^j \binom{22}{j}\frac{(97-j)!}{100!}.}$$

This is the exact probability that exactly three primes are fixed and the other twenty-two are foolish. The implementations evaluate this normalized sum directly rather than forming the huge integer count first.

The multiplicative recurrence used by the implementations

The solution files do not repeatedly compute factorials and binomial coefficients from scratch inside the alternating sum. Instead they use the normalized term

$$T_j = (-1)^j \binom{b}{j}\frac{(n-f-j)!}{n!},$$

so that

$$P=\binom{m}{f}\sum_{j=0}^{b} T_j.$$

The first term is

$$T_0=\frac{(n-f)!}{n!}=\prod_{i=0}^{f-1}\frac{1}{n-i},$$

and adjacent terms satisfy the simple ratio

$$T_{j+1}=T_j\cdot\frac{-(b-j)}{(j+1)(n-f-j)}.$$

That recurrence is exactly what the C++, Python, and Java implementations update inside their main loop. It keeps the computation compact and avoids manufacturing enormous intermediate factorial values only to divide them away again a moment later.

Worked example

A smaller example shows the same logic more transparently. Suppose $n=8$, $m=3$, and we want exactly $f=1$ marked label fixed. Then $b=2$. First choose which marked label is fixed: $\binom{3}{1}=3$ choices.

After that, $7$ labels remain. The other two marked labels must avoid their own positions. Inclusion-exclusion gives

$$3\left(7!-\binom{2}{1}6!+\binom{2}{2}5!\right)=3(5040-1440+120)=11160.$$

Since $8!=40320$, the probability is

$$\frac{11160}{40320}=\frac{31}{112}.$$

Problem 239 is exactly the same argument, just with $100$ total labels, $25$ marked primes, $3$ fixed primes, and $22$ forbidden prime fixed points.

How the Code Works

The implementations start from the general parameters $n$, $m$, and the number of misplaced marked labels. From that they derive $f=m-r$, the number of marked fixed points that must occur. For Problem 239 this means converting "22 foolish primes" into "3 fixed prime-labelled items."

Next, they compute $\binom{m}{f}$ multiplicatively, one factor at a time. The alternating sum is then evaluated through the recurrence for $T_j$: the code builds $T_0$ by dividing by $n$, $n-1$, and so on for exactly $f$ factors, and each later term is obtained from the previous one by multiplying with the rational ratio above. The running sum therefore mirrors the inclusion-exclusion formula term by term.

Finally, the sum is multiplied by $\binom{m}{f}$ and formatted as a decimal with twelve digits after the point. The C++ implementation also performs small-instance checks: it compares the inclusion-exclusion count against exact enumeration on tiny cases and verifies on a sample input that summing the probabilities over all possible numbers of fixed marked labels gives total mass $1$.

Complexity Analysis

The binomial coefficient $\binom{m}{f}$ is computed in $O(\min(f,m-f))$ arithmetic steps, and the alternating sum has exactly $b+1=m-f+1$ terms. So the overall running time is $O(m)$, with constant extra space.

For the actual Project Euler instance this is tiny: $f=3$ and $b=22$, so the probability comes from only $23$ inclusion-exclusion terms. The mathematical work is in deriving the formula; the numerical evaluation itself is very light.

Footnotes and References

Problem page: Project Euler 239
Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
Fixed point of a permutation: Wikipedia - Fixed point (mathematics)
Rencontres numbers and partial derangements: Wikipedia - Rencontres numbers
Prime numbers up to 100: Wikipedia - Prime number

Problem 239 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <thread>
#include <vector>

namespace {

long double combination_ld(int n, int k) {
    if (k < 0 || k > n) return 0.0L;
    k = std::min(k, n - k);
    long double value = 1.0L;
    for (int i = 1; i <= k; ++i) {
        value *= static_cast<long double>(n - k + i);
        value /= static_cast<long double>(i);
    }
    return value;
}

long double probability_exact_marked_fixed(int n, int marked, int fixed_marked) {
    if (n < 0 || marked < 0 || marked > n) return 0.0L;
    if (fixed_marked < 0 || fixed_marked > marked) return 0.0L;

    const int bad_marked = marked - fixed_marked;

    long double term = 1.0L;
    for (int i = 0; i < fixed_marked; ++i) {
        term /= static_cast<long double>(n - i);
    }

    long double alternating_sum = term;
    for (int j = 0; j < bad_marked; ++j) {
        const long double numer = -static_cast<long double>(bad_marked - j);
        const long double denom = static_cast<long double>(j + 1) *
                                  static_cast<long double>(n - fixed_marked - j);
        term *= numer / denom;
        alternating_sum += term;
    }

    return combination_ld(marked, fixed_marked) * alternating_sum;
}

__int128 combination_i128(int n, int k) {
    if (k < 0 || k > n) return 0;
    k = std::min(k, n - k);
    __int128 value = 1;
    for (int i = 1; i <= k; ++i) {
        value = value * static_cast<__int128>(n - k + i) / static_cast<__int128>(i);
    }
    return value;
}

__int128 factorial_i128(int n) {
    __int128 value = 1;
    for (int i = 2; i <= n; ++i) {
        value *= static_cast<__int128>(i);
    }
    return value;
}

__int128 count_exact_marked_fixed_inclusion(int n, int marked, int fixed_marked) {
    if (n < 0 || marked < 0 || marked > n) return 0;
    if (fixed_marked < 0 || fixed_marked > marked) return 0;

    const int bad_marked = marked - fixed_marked;
    __int128 alternating_sum = 0;
    for (int j = 0; j <= bad_marked; ++j) {
        __int128 term = combination_i128(bad_marked, j) * factorial_i128(n - fixed_marked - j);
        if ((j & 1) == 0) {
            alternating_sum += term;
        } else {
            alternating_sum -= term;
        }
    }
    return combination_i128(marked, fixed_marked) * alternating_sum;
}

std::uint64_t brute_count_exact_marked_fixed(int n, int marked, int fixed_marked, int threads) {
    if (threads < 1) threads = 1;
    if (n <= 0) return (marked == 0 && fixed_marked == 0) ? 1ULL : 0ULL;

    std::atomic<int> next_first{0};
    std::vector<std::uint64_t> partial(static_cast<std::size_t>(threads), 0ULL);
    std::vector<std::thread> pool;
    pool.reserve(static_cast<std::size_t>(threads));

    for (int t = 0; t < threads; ++t) {
        pool.emplace_back([&, t]() {
            std::uint64_t local = 0;
            while (true) {
                const int first = next_first.fetch_add(1, std::memory_order_relaxed);
                if (first >= n) break;

                std::vector<int> rest;
                rest.reserve(static_cast<std::size_t>(n - 1));
                for (int v = 0; v < n; ++v) {
                    if (v != first) rest.push_back(v);
                }

                do {
                    int fixed = 0;
                    if (marked > 0 && first == 0) ++fixed;
                    for (int pos = 1; pos < marked; ++pos) {
                        if (rest[static_cast<std::size_t>(pos - 1)] == pos) ++fixed;
                    }
                    if (fixed == fixed_marked) ++local;
                } while (std::next_permutation(rest.begin(), rest.end()));
            }
            partial[static_cast<std::size_t>(t)] = local;
        });
    }

    for (std::thread& th : pool) {
        th.join();
    }

    std::uint64_t total = 0;
    for (std::uint64_t value : partial) total += value;
    return total;
}

bool validate() {
    struct TinyCase {
        int n;
        int marked;
        int fixed_marked;
    };

    const std::vector<TinyCase> tiny_cases = {
        {7, 3, 1},
        {8, 4, 2},
        {9, 5, 1},
        {10, 4, 2},
    };

    for (const TinyCase& c : tiny_cases) {
        const __int128 exact_count = count_exact_marked_fixed_inclusion(c.n, c.marked, c.fixed_marked);
        const std::uint64_t brute_count = brute_count_exact_marked_fixed(c.n, c.marked, c.fixed_marked, 1);
        if (exact_count != static_cast<__int128>(brute_count)) {
            std::cerr << "Validation failed (inclusion vs brute) at n=" << c.n
                      << ", marked=" << c.marked << ", fixed=" << c.fixed_marked
                      << ": inclusion=" << static_cast<long long>(exact_count)
                      << ", brute=" << brute_count << "\n";
            return false;
        }

        const long double prob_formula = probability_exact_marked_fixed(c.n, c.marked, c.fixed_marked);
        const long double prob_exact = static_cast<long double>(brute_count) /
                                       static_cast<long double>(factorial_i128(c.n));
        if (std::fabsl(prob_formula - prob_exact) > 1e-16L) {
            std::cerr << "Validation failed (probability mismatch) at n=" << c.n
                      << ", marked=" << c.marked << ", fixed=" << c.fixed_marked
                      << ": formula=" << std::setprecision(22) << prob_formula
                      << ", exact=" << prob_exact << "\n";
            return false;
        }
    }

    const int n_mass = 20;
    const int marked_mass = 7;
    long double mass = 0.0L;
    for (int fixed = 0; fixed <= marked_mass; ++fixed) {
        mass += probability_exact_marked_fixed(n_mass, marked_mass, fixed);
    }
    if (std::fabsl(mass - 1.0L) > 1e-16L) {
        std::cerr << "Validation failed (total mass): got=" << std::setprecision(22) << mass << "\n";
        return false;
    }

    unsigned hw = std::thread::hardware_concurrency();
    if (hw == 0) hw = 2;
    const int threads = static_cast<int>(std::min<unsigned>(hw, 8));
    if (threads > 1) {
        const TinyCase thread_case{10, 4, 2};
        const std::uint64_t single =
            brute_count_exact_marked_fixed(thread_case.n, thread_case.marked, thread_case.fixed_marked, 1);
        const std::uint64_t multi =
            brute_count_exact_marked_fixed(thread_case.n, thread_case.marked, thread_case.fixed_marked, threads);
        if (single != multi) {
            std::cerr << "Validation failed (thread consistency): single=" << single
                      << ", multi=" << multi << "\n";
            return false;
        }
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    if (!validate()) return 1;

    int n = 100;
    int marked = 25;
    int misplaced_marked = 22;

    if (argc > 1) n = std::max(1, std::atoi(argv[1]));
    if (argc > 2) marked = std::clamp(std::atoi(argv[2]), 0, n);
    if (argc > 3) misplaced_marked = std::clamp(std::atoi(argv[3]), 0, marked);

    const int fixed_marked = marked - misplaced_marked;
    const long double answer = probability_exact_marked_fixed(n, marked, fixed_marked);
    std::cout << std::fixed << std::setprecision(12) << answer << '\n';
    return 0;
}

Python

def combination(n, k):
    if k < 0 or k > n:
        return 0.0
    k = min(k, n - k)
    value = 1.0
    for i in range(1, k + 1):
        value *= (n - k + i)
        value /= i
    return value

def probability_exact_marked_fixed(n, marked, fixed_marked):
    if n < 0 or marked < 0 or marked > n:
        return 0.0
    if fixed_marked < 0 or fixed_marked > marked:
        return 0.0
        
    bad_marked = marked - fixed_marked
    
    term = 1.0
    for i in range(fixed_marked):
        term /= (n - i)
        
    alternating_sum = term
    for j in range(bad_marked):
        numer = -(bad_marked - j)
        denom = (j + 1) * (n - fixed_marked - j)
        term *= numer / denom
        alternating_sum += term
        
    return combination(marked, fixed_marked) * alternating_sum

def solve(n=100, marked=25, misplaced_marked=22):
    fixed_marked = marked - misplaced_marked
    answer = probability_exact_marked_fixed(n, marked, fixed_marked)
    return f"{answer:.12f}"

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

Java

import java.util.Locale;

public class Euler239 {
    static double combination(int n, int k) {
        if (k < 0 || k > n)
            return 0.0;
        k = Math.min(k, n - k);
        double value = 1.0;
        for (int i = 1; i <= k; ++i) {
            value *= (n - k + i);
            value /= i;
        }
        return value;
    }

    static double probabilityExactMarkedFixed(int n, int marked, int fixedMarked) {
        if (n < 0 || marked < 0 || marked > n)
            return 0.0;
        if (fixedMarked < 0 || fixedMarked > marked)
            return 0.0;

        int badMarked = marked - fixedMarked;

        double term = 1.0;
        for (int i = 0; i < fixedMarked; ++i) {
            term /= (n - i);
        }

        double alternatingSum = term;
        for (int j = 0; j < badMarked; ++j) {
            double numer = -(badMarked - j);
            double denom = (j + 1.0) * (n - fixedMarked - j);
            term *= numer / denom;
            alternatingSum += term;
        }

        return combination(marked, fixedMarked) * alternatingSum;
    }

    public static String solve() {
        int n = 100;
        int marked = 25;
        int misplacedMarked = 22;

        int fixedMarked = marked - misplacedMarked;
        double answer = probabilityExactMarkedFixed(n, marked, fixedMarked);

        return String.format(Locale.US, "%.12f", answer);
    }

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