Problem 239: Twenty-two Foolish Primes

Detailed mathematical approach

Problem Summary

Discs numbered \(1\) to \(100\) are placed in a row in random order, all \(100!\) orders being equally likely, and disc \(k\) is in its natural position when it is the \(k\)-th disc of the row. Of the discs, \(25\) carry a prime number. We must find the probability that exactly \(22\) of these prime discs are away from their natural positions, so that exactly \(3\) of them are in place, while the \(75\) other discs may be in or out of place. The answer is to be rounded to \(12\) decimal places, in the form 0.abcdefghijkl; the statement gives no intermediate values. ^[1]

The whole solution in one sentence: choose which \(3\) of the \(25\) prime discs stay in place, count by inclusion–exclusion the orders of the other \(97\) discs that leave none of the remaining \(22\) prime discs in place, divide by \(100!\), and add up the resulting \(23\) alternating terms as ratios, so that no huge factorial is ever formed.

Mathematical Approach

We count by hand on four discs, prove the formula in general, test it on ten discs, and finally show that floating-point arithmetic settles all twelve decimals.

1. Fixed points and the twenty-five primes

Describe an order by the permutation \(\sigma\) in which \(\sigma(i)\) is the number of the disc at position \(i\). Disc \(k\) is in place exactly when \(\sigma(k)=k\), that is, when \(k\) is a fixed point of \(\sigma\). The primes up to \(100\) are 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89 and 97. We call the prime discs marked and let \(K\) be the number of marked discs in place. The event of the problem is \(K=3\), so

$$P(K=3)=\frac{A(100,25,3)}{100!},$$

where \(A(n,m,f)\) is the number of orders of \(n\) discs, \(m\) of them marked, with exactly \(f\) marked discs in place. Only the number of marked discs matters: renaming the discs and the positions by the same one-to-one map keeps fixed points fixed, so we may take the marked discs to be the first \(m\) ones, as the C++ check does.

2. Four discs by hand

Take the discs \(1\) to \(4\), whose primes are \(2\) and \(3\), and ask for exactly one prime disc in place. If disc \(2\) is in place, the discs \(1\), \(3\) and \(4\) fill the other three positions in \(3!=6\) ways: 1234, 1243, 3214, 3241, 4213 and 4231. In two of them, 1234 and 4231, disc \(3\) is in place too, so \(6-2=4\) rows have disc \(2\) in place and disc \(3\) away; in 1243 disc \(1\) is in place as well, which is allowed. Likewise the four rows 1432, 2134, 2431 and 4132 have disc \(3\) in place and disc \(2\) away. Hence \(A(4,2,1)=8\) and the probability is \(8/24=1/3\); of the other \(16\) rows, \(14\) have no prime disc in place and \(2\) have both.

These two moves, choosing the marked discs that stay and removing unwanted rows, are the two steps of the general solution. The example also shows why probabilities cannot simply be multiplied: each prime disc is in place with probability \(1/4\), but both are with probability \(2/24=1/12\), not \(1/16\), because once disc \(2\) is at home, disc \(3\) has only three positions left. Independence would give \(2\cdot\frac14\cdot\frac34=\frac38\) instead of \(\frac13\).

3. Choosing the marked discs that stay in place

Let \(D(N,b)\) be the number of orders of \(N\) discs in which \(b\) specified discs are all away from their positions, the other \(N-b\) discs being free. Then

$$A(n,m,f)=\binom{m}{f}\,D(n-f,\,m-f).$$

Proof. Each order determines the set \(F\) of marked discs that it leaves in place, so the orders with \(K=f\) split according to \(F\), which is one of \(\binom{m}{f}\) sets of size \(f\). Fix \(F\). In an order with this set, the discs of \(F\) occupy their own positions, so the other \(n-f\) discs fill the other \(n-f\) positions, which carry the same numbers, and no marked disc outside \(F\) is in place. Conversely, every arrangement of these \(n-f\) discs that keeps the \(m-f\) marked discs outside \(F\) away from home becomes, with the discs of \(F\) put back, an order whose set is exactly \(F\). This correspondence is one-to-one, so each \(F\) contributes \(D(n-f,m-f)\) orders. With every disc marked, this is the classical formula for the rencontres numbers, which count the permutations with exactly \(k\) fixed points: choose them, then derange the rest. ^[2] In Section 2, \(\binom21 D(3,1)=2\cdot4=8\).

4. Inclusion–exclusion

Before the general formula, return to four discs and require both prime discs to be away. Of the \(24\) orders, \(6\) have disc \(2\) in place and \(6\) have disc \(3\) in place; the \(2\) orders with both were subtracted twice, so they are added back once, which gives \(24-6-6+2=14\), the count found in Section 2. With more specified discs one subtracts the single restrictions, adds back the pairs, subtracts the triples, and so on, which gives

$$D(N,b)=\sum_{j=0}^{b}(-1)^j\binom{b}{j}(N-j)!.$$

Proof. If a set \(J\) of \(j\) specified discs must be in place, the other \(N-j\) discs can be arranged freely, in \((N-j)!\) ways, and there are \(\binom{b}{j}\) such sets. So the right-hand side adds, over all sets \(J\) of specified discs, \((-1)^{|J|}\) times the number of orders that leave all of \(J\) in place. An order that leaves exactly \(s\) specified discs in place is counted once for each subset \(J\) of these \(s\) discs, with the sign \((-1)^{|J|}\), so it contributes

$$\sum_{j=0}^{s}(-1)^j\binom{s}{j}=(1-1)^s$$

by the binomial theorem: \(1\) if \(s=0\) and \(0\) if \(s\ge1\). The sum therefore counts exactly the orders that leave no specified disc in place. This is the inclusion–exclusion principle for the elements lying in none of several sets, here the sets \(E_i\) of orders in which the \(i\)-th specified disc is in place, with its standard proof. ^[3] For instance \(D(8,2)=8!-2\cdot7!+6!\)\({}=40320-10080+720\)\({}=30960\). With \(b=N\), \(D(N,N)\) counts the derangements of \(N\) objects. ^[2]

5. The formula for Problem 239

With \(n=100\), \(m=25\) and \(f=3\), Sections 1, 3 and 4 give

$$\begin{aligned}P(K=3)&=\binom{25}{3}\,S,\\ S&=\sum_{j=0}^{22}(-1)^j\binom{22}{j}\frac{(97-j)!}{100!},\end{aligned}$$

where \(\binom{25}{3}=2300\) counts the choices of the three prime discs in place, and \(S\) is the probability that the prime discs in place are exactly three given ones, say \(2\), \(3\) and \(5\).

6. Worked example: ten discs

For the discs \(1\) to \(10\) the marked discs are \(2\), \(3\), \(5\) and \(7\), so \(n=10\) and \(m=4\). Exactly two prime discs away means \(f=2\), and

$$A(10,4,2)=\binom42 D(8,2)=6\cdot30960=185760,$$

a probability of \(185760/10!=43/840\approx0.051190\). Listing all \(3628800\) orders confirms this count and the whole distribution of \(K\):

• \(K=0\): \(2399760\) orders, probability \(1111/1680\approx0.661310\);
• \(K=1\): \(1025280\) orders, probability \(89/315\approx0.282540\);
• \(K=2\): \(185760\) orders, probability \(43/840\approx0.051190\);
• \(K=3\): \(17280=4\cdot(7!-6!)\) orders, probability \(1/210\approx0.004762\);
• \(K=4\): \(720=6!\) orders, probability \(1/5040\approx0.000198\).

The counts add up to \(10!\), and the mean of \(K\), each value weighted by its probability, is exactly \(2/5\). This agrees with a short argument: a given disc is in place in \((n-1)!\) of the \(n!\) orders, so with probability \(1/n\), and the expectation of a sum is the sum of the expectations, so \(\mathbb E[K]=m/n\). For \(m=n\) this is the classical fact that a random permutation has one fixed point on average. ^[2] In Problem 239, \(\mathbb E[K]=1/4\).

7. A triangle of differences

The numbers \(D(N,b)\) also obey a recurrence, which gives a second way to compute them:

$$D(N,0)=N!,\qquad D(N,b)=D(N,b-1)-D(N-1,b-1).$$

Indeed, the orders that keep the first \(b-1\) specified discs away split into those that keep the \(b\)-th one away too, \(D(N,b)\) of them, and those that leave it in place; deleting that disc and its position maps the latter one-to-one onto the orders of \(N-1\) discs that keep the first \(b-1\) specified discs away. For \(b=0,1,\dots,N\) the rows \(N=4\) and \(N=5\) are \(24, 18, 14, 11, 9\) and \(120, 96, 78, 64, 53, 44\), for example \(78=96-18\); the rows start with the factorials and end with the derangement numbers. The OEIS lists this triangle of successive differences of the factorials together with the formula of Section 4. ^[4] We used it to check the count of Section 5 with exact integers.

8. Terms as ratios

Write \(N=n-f\) and \(b=m-f\), and divide each term of the alternating sum by \(n!\):

$$\begin{aligned}T_j&=(-1)^j\binom{b}{j}\frac{(N-j)!}{n!},\\ P(K=f)&=\binom{m}{f}\sum_{j=0}^{b}T_j.\end{aligned}$$

The first term is \(T_0=(n-f)!/n!\), the reciprocal of \(n(n-1)\cdots(n-f+1)\), here \(1/970200\). Since \(\binom{b}{j+1}\) is \(\binom{b}{j}\) times \(\frac{b-j}{j+1}\) and \((N-j-1)!\) is \((N-j)!\) divided by \(N-j\),

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

The ratios avoid huge numbers: \(100!\) has \(158\) digits, far beyond \(64\)-bit or \(128\)-bit integers, and although a double can hold the factorials up to \(170!\), the ratios keep every term and every partial sum at most \(1\) in size for any \(n\), for one division and one multiplication per term. Here the ratio is \(-22/97\) for \(j=0\) and \(-7/64\) for \(j=1\), and its absolute value keeps falling, as the numerator \(22-j\) shrinks and the denominator \((j+1)(97-j)\) grows. So the terms alternate in sign, each less than a quarter of the previous one in size.

For such a sum the partial sums lie alternately above and below the total, and stopping after \(T_J\) changes the total by at most \(|T_{J+1}|\); this is the alternating series estimate. ^[5] In particular \(T_0+T_1\le S\le T_0\). Multiplied by \(2300\), the term \(j=9\) is already below \(3\cdot10^{-15}\), so the nine terms \(j=0,\dots,8\) settle the twelfth decimal, given the margin found in Section 9; the programs still add all twenty-three terms, which takes little work.

9. Rounding errors and the twelfth decimal

Python's float, Java's double and, on some machines, C++'s long double are IEEE binary64 numbers: unless it underflows, each operation returns the exact result times a factor \(1+\delta\) with \(|\delta|\le u=2^{-53}\approx1.11\cdot10^{-16}\), and no term here is below \(10^{-49}\) in size; the 80-bit long double of GCC on x86 Linux has \(u=2^{-64}\). ^{[6] [7]} The numerators and denominators of the ratios are integers of at most \(1672\), hence exact. \(T_0\) takes three divisions and each later term one division and one multiplication, so the computed \(T_j\) is the exact one times at most \(3+2j\le47\) such factors. The \(22\) additions bring at most \(22\) more factors into the contribution of each term, so the computed sum is \(\sum_j T_j(1+\theta_j)\), where each \(1+\theta_j\) is a product of at most \(69\) factors. Such a product differs from \(1\) by at most \(\gamma_{69}\), where \(\gamma_r=ru/(1-ru)\): from above, \((1+u)^r-1\) is at most \(ru+(ru)^2+\cdots=\gamma_r\), and from below, \(1-(1-u)^r\le ru\). The binomial coefficient is exact, since its evaluation meets only the integers \(23\), \(552\), \(276\), \(6900\) and \(2300\), and the final product adds one more factor. So, with \(\hat P\) the computed and \(P\) the exact probability,

$$\begin{aligned}\frac{|\hat P-P|}{P}&\le(1+u)\,\gamma_{69}\,\kappa+u,\\ \kappa&=\frac{\sum_j|T_j|}{|S|},\end{aligned}$$

where \(\kappa\) is the condition number of the sum, the factor by which cancellation can magnify rounding errors. ^[8] By Section 8, \(\sum_j|T_j|\le T_0/(1-22/97)\) and \(S\ge(75/97)\,T_0\), so \(\kappa\le(97/75)^2\approx1.673\), and the relative error is below \(1.3\cdot10^{-14}\), as is the absolute error, since \(P\lt1\). The rounding to twelve decimals could only fail if the exact value were that close to a point where its twelfth decimal changes. It is not: exact rational arithmetic shows that it is more than \(4\cdot10^{-13}\) away from every such point, and that the binary64 result is off by only \(2.4\cdot10^{-17}\) in relative terms, about a fifth of \(u\).

How the Code Works

The C++ program ^[9] begins in main with validate, described below, and exits with code \(1\) if a check fails. It then reads up to three optional arguments, the numbers of discs, of marked discs and of marked discs away from their positions, with the defaults \(100\), \(25\) and \(22\). Supplied values are clamped into range, but defaults are not, so a single argument below \(25\) makes it print zero. It sets fixed_marked to marked - misplaced_marked, which is \(3\) by default, and prints probability_exact_marked_fixed(n, marked, fixed_marked) with std::fixed and std::setprecision(12). With the arguments 10 4 2 it prints 0.051190476190, the value \(43/840\) of Section 6.

probability_exact_marked_fixed carries out Section 8 in long double. It returns \(0\) for impossible parameters, sets bad_marked to \(b\), divides \(1\) by \(n,n-1,\dots,n-f+1\) to obtain \(T_0\), and then, for \(j=0,\dots,b-1\), multiplies term by numer / denom, where numer is \(-(b-j)\) and denom is \((j+1)(n-f-j)\), adding each new term to alternating_sum. The sum is multiplied by combination_ld(marked, fixed_marked), which replaces \(k\) by \(\min(k,n-k)\) and, for \(i=1,\dots,k\), multiplies by \(n-k+i\) and divides by \(i\).

validate checks the tiny cases \((n,m,f)=(7,3,1)\), \((8,4,2)\), \((9,5,1)\) and \((10,4,2)\). count_exact_marked_fixed_inclusion evaluates the count of Sections 3 and 4 exactly in __int128, the 128-bit integers of GCC and Clang, which hold the factorials up to \(33!\), using combination_i128 and factorial_i128. brute_count_exact_marked_fixed marks the discs \(0,\dots,m-1\), which by Section 1 loses no generality, and visits all \(n!\) orders: for each disc at position \(0\) it lists the other discs in increasing order and lets std::next_permutation run through their \((n-1)!\) orders until it returns false. ^[10] The two counts must agree; they are \(1512\), \(3024\), \(120120\) and \(185760\). Next, the floating-point probability must match the exact count divided by \(n!\) within \(10^{-16}\), and the probabilities for \(f=0,\dots,7\) with \(n=20\) and \(m=7\) must add up to \(1\) within \(10^{-16}\), since every order has some number of marked discs in place. Finally, it enumerates the ten-disc case twice more and requires the two counts to agree.

The Python and Java programs contain only the floating-point part, in binary64, with the same loops in the same order: combination, probability_exact_marked_fixed (probabilityExactMarkedFixed in Java) and solve, whose parameters are default arguments in Python and fixed in Java. Python forms the denominator with exact integers, Java as the exact double product (j + 1.0) * (n - fixedMarked - j). Python prints with f"{answer:.12f}", Java with String.format(Locale.US, "%.12f", answer), whose explicit locale guarantees a decimal point. Neither runs any checks.

Complexity and Verification

The probability routine needs \(f\) divisions, four floating-point operations for each of the \(b\) further terms and \(\min(f,m-f)\) steps for the binomial coefficient: \(O(m)\) operations and \(O(1)\) memory, about a hundred operations here. In C++, validate dominates: its four cases visit \(4037040\) orders, and the two extra enumerations of the ten-disc case bring the total to \(11294640\). Each order costs a scan of the \(m\) marked positions and, on average, about three comparisons in std::next_permutation, ^[10] so a case takes \(O((m+1)\,n!)\) time and \(O(n)\) memory. On our machine each of the three programs finishes within a few hundredths of a second.

The C++ self-check passes on our machine, where long double has the binary64 format of double. There the tolerance \(10^{-16}\) is below the gap \(2^{-53}\approx1.1\cdot10^{-16}\) between \(1\) and the next smaller double, so the total-probability check passes only because the computed sum is exactly \(1\) for \(n=20\) and \(m=7\); with 80 bits it has a wide margin. For this article we also checked, outside the published programs:

• listing all orders of \(0\le n\le8\) discs, the empty row included, gives \(A(n,m,f)\) for all \(0\le f\le m\le n\), \(165\) cases in all;
• the recurrence of Section 7 holds for all \(N\lt60\) and all \(b\), and its triangle up to \(N=100\) gives the same integer \(\binom{25}{3}D(97,22)\) as the alternating sum;
• with exact fractions, the probabilities \(P(K=f)\) for \(n=100\), \(m=25\) and \(f=0,\dots,25\) add up to \(1\), with mean \(1/4\);
• the error figures of Section 9 come from exact rational arithmetic;
• the C++ program, compiled with Clang 22 and GCC 16 at several optimisation levels, Python 3.11 and Java 17 print the same twelve decimals.

Footnotes and References

The references below supply the background facts; the counts, the distributions, the error bound and the checks were computed for this article.

1. Project Euler 239 — Twenty-two Foolish Primes. Discs numbered 1 to 100 lie in a line in random order; the task asks for the probability that exactly 22 prime-numbered discs are away from their natural positions, non-prime discs being unrestricted, rounded to 12 decimal places as 0.abcdefghijkl.
2. Wikipedia — Rencontres numbers. The rencontres number Dₙ,ₖ counts the permutations of {1, …, n} with exactly k fixed points, the partial derangements; Dₙ,ₖ = C(n, k) Dₙ₋ₖ,₀ (choose the fixed points, derange the rest), Dₙ,₀ counts the derangements, and for n ≥ 1 the mean number of fixed points is 1.
3. Wikipedia — Inclusion–exclusion principle. In the formula for the union of finite sets A₁, …, Aₙ, an element of t ≥ 1 of them gets the total coefficient 1, by the binomial theorem 0 = (1 − 1)ᵗ; the complementary form counts the elements in none of the sets, where such an element gets 0. Derangements are an application.
4. OEIS A047920 — Triangular array formed from successive differences of factorial numbers. T(n, k) = Σⱼ (−1)ʲ C(k, j) (n − j)!, with t(n, k) = t(n, k − 1) − t(n − 1, k − 1) and rows 1; 1, 0; 2, 1, 1; 6, 4, 3, 2; …; T(n, k) counts the permutations of {1, …, n} that do not fix the points 1, …, k.
5. Wikipedia — Alternating series test. If the absolute values of the terms decrease monotonically to zero, the partial sums lie alternately above and below the limit, and Sₖ differs from it by at most the next term; a finite sum followed by zeros is a special case.
6. Wikipedia — Machine epsilon. For binary64 with rounding to nearest, the rounding machine epsilon 2⁻⁵³ ≈ 1.11 × 10⁻¹⁶ bounds the relative error of a rounding; for the extended format with a 64-bit significand, used for long double, it is 2⁻⁶⁴.
7. Wikipedia — Long double. The standard only requires long double to be at least as precise as double; with GCC on x86 Linux it is the 80-bit extended format, while Microsoft Visual C++ and 64-bit ARM under Windows and macOS give it the format of double.
8. Wikipedia — Pairwise summation. Σ|xᵢ| / |Σxᵢ| is the condition number of the summation problem, and the relative error bound of naive summation of n numbers grows as O(εn) times it, ε being the machine precision.
9. C++ — Euler239.cpp. The linked immutable C++ revision contains probability_exact_marked_fixed, the exact and brute-force counts and validate. It is the implementation record for this article.
10. cppreference — std::next_permutation. It turns a range into the lexicographically next permutation and returns true, or, after the last one, into the first and returns false; typical implementations average about 3 comparisons and 1.5 swaps per call.