Problem 398: Cutting Rope
View on Project EulerProject Euler Problem 398 Solution
EulerSolve provides an optimized solution for Project Euler Problem 398, Cutting Rope, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Choose uniformly from all ordered compositions of \(n\) into \(m\) positive integers: $$\mathcal{C}_{n,m}=\left\{(x_1,\dots,x_m)\in \mathbb{Z}_{>0}^m : x_1+\cdots+x_m=n\right\}.$$ After sorting one sampled composition into nondecreasing order, write $$X_{(1)}\le X_{(2)}\le \cdots \le X_{(m)}.$$ The task solved by the code is to compute the expected value of the second shortest piece, namely \(\mathbb{E}[X_{(2)}]\), for the Project Euler parameters \(n=10^7\) and \(m=100\). Mathematical Approach 1) Total number of equally likely cuts By the standard stars-and-bars argument, the number of ordered compositions of \(n\) into \(m\) positive parts is $$\left|\mathcal{C}_{n,m}\right|=\binom{n-1}{m-1}.$$ So every probability in the program is a composition count divided by \(\binom{n-1}{m-1}\). 2) Tail-sum identity for an integer-valued random variable Because \(X_{(2)}\) is a positive integer, we may use the tail-sum formula $$\mathbb{E}[X_{(2)}]=\sum_{k\ge 1}\Pr(X_{(2)}\ge k).$$ This reduces the whole problem to counting, for each \(k\), how many compositions have second smallest part at least \(k\). 3) Interpret the event \(X_{(2)}\ge k\) The condition \(X_{(2)}\ge k\) means that at most one part is smaller than \(k\). Equivalently, the composition belongs to one of two disjoint cases: Case A: all \(m\) parts are at least \(k\)....
Detailed mathematical approach
Problem Summary
A rope of length \(n\) carries \(n-1\) marked points, one at every integer distance from the two ends. We choose \(m-1\) of them uniformly at random, all \(\binom{n-1}{m-1}\) choices being equally likely, and cut there, leaving \(m\) segments whose lengths are positive integers adding up to \(n\). Sort those lengths into nondecreasing order and take the second entry: \(E(n,m)\) is its expected value. Sorting with repetitions is exactly the tie rule of the statement, which declares the second-shortest length equal to the shortest when segments tie. The statement gives the checkpoints \(E(3,2)=2\) and \(E(8,3)=\tfrac{16}{7}\), and asks for \(E(10^7,100)\), rounded to five decimal places. [1]
The whole solution in one sentence: for every length \(k\) we count the cuts that leave at most one segment shorter than \(k\) β a two-term expression in binomial coefficients β divide by the number of cuts to get a probability, and add those probabilities over all \(k\), carrying the two terms as running ratios so that no large number is ever built.
Mathematical Approach
1. Cuts are compositions
Write the cut positions as \(0\lt c_1\lt\cdots\lt c_{m-1}\lt n\) and read off the lengths \(x_1=c_1\), \(x_2=c_2-c_1\), β¦, \(x_m=n-c_{m-1}\). Every \(x_i\) is a positive integer with \(x_1+\cdots+x_m=n\), and the partial sums of such a list give the cut positions back. A choice of cuts is therefore the same thing as a composition of \(n\) into \(m\) parts, an ordered tuple of \(m\) positive integers summing to \(n\), and there are \(\binom{n-1}{m-1}\) of those, one per \((m-1)\)-element subset of the marked points. [2] Because the cut points are chosen uniformly, all these compositions are equally likely, and every probability below is a count divided by \(\binom{n-1}{m-1}\).
Take \(n=8\), \(m=3\). There are \(\binom{7}{2}=21\) compositions, and they fall into five groups according to which multiset of lengths they produce:
- \(\{1,1,6\}\): \(3\) orderings, sorted list \(1,1,6\), second entry \(1\);
- \(\{1,2,5\}\): \(6\) orderings, second entry \(2\);
- \(\{1,3,4\}\): \(6\) orderings, second entry \(3\);
- \(\{2,2,4\}\): \(3\) orderings, second entry \(2\);
- \(\{2,3,3\}\): \(3\) orderings, second entry \(3\).
The second entry is \(1\) for \(3\) compositions, \(2\) for \(6+3=9\) and \(3\) for \(6+3=9\), which gives
$$\begin{aligned}E(8,3)&=\frac{3\cdot1+9\cdot2+9\cdot3}{21}\\&=\frac{48}{21}=\frac{16}{7}.\end{aligned}$$
The first group shows the tie rule at work: in \(1,1,6\) the second-shortest length counts as \(1\). For \(m=2\) the rule makes the "second-shortest" of two segments the longer one, hence \(E(3,2)=2\): the only cuts give \(1,2\) and \(2,1\).
2. The second-shortest length as an order statistic
Let \(X_1,\dots,X_m\) be the segment lengths in the order they lie on the rope, and let
$$X_{(1)}\le X_{(2)}\le\cdots\le X_{(m)}$$
be the same numbers sorted. The \(k\)-th entry of a sorted sample is its \(k\)-th order statistic; repeated values simply appear repeatedly. [3] We want \(E(n,m)=\mathbb{E}\bigl[X_{(2)}\bigr]\), with no special case for ties.
3. Expectation as a sum of tail probabilities
\(X_{(2)}\) takes only positive integer values, and for such a random variable
$$\mathbb{E}[X]=\sum_{k\ge1}\Pr(X\ge k).$$
The proof is a rearrangement of a double sum: [4]
$$\begin{aligned}\sum_{k\ge1}\Pr(X\ge k)&=\sum_{k\ge1}\sum_{x\ge k}\Pr(X=x)\\&=\sum_{x\ge1}x\Pr(X=x),\end{aligned}$$
because the value \(x\) is counted once for each of \(k=1,\dots,x\). The sum is finite here, as Section 7 shows \(X_{(2)}\le\lfloor(n-1)/(m-1)\rfloor\). The whole task is now to count, for each \(k\), the compositions whose second-smallest part is at least \(k\).
4. What "second-smallest at least \(k\)" means
If \(X_{(2)}\ge k\) then \(X_{(2)},\dots,X_{(m)}\) are all at least \(k\), so at most one part β the smallest β is below \(k\). Conversely, if at most one part is below \(k\), the sorted list has at least \(m-1\) entries \(\ge k\), and these occupy positions \(2,\dots,m\), so \(X_{(2)}\ge k\). The event therefore splits into two disjoint cases:
- Case A: every part is at least \(k\);
- Case B: exactly one part lies in \(\{1,\dots,k-1\}\) and the other \(m-1\) parts are at least \(k\).
For \(k=1\), case B is empty and case A is everything, so \(\Pr(X_{(2)}\ge1)=1\), as it must be.
5. Counting case A
Demanding \(x_i\ge k\) for all \(i\) is a change of variables. Put \(y_i=x_i-(k-1)\); then \(y_i\ge1\) and
$$y_1+\cdots+y_m=n-m(k-1).$$
Every composition with all parts \(\ge k\) gives exactly one composition of \(n-m(k-1)\) into \(m\) positive parts and vice versa, so by Section 1 the count is
$$A_k=\binom{n-m(k-1)-1}{m-1}.$$
Throughout we use the convention \(\binom{a}{r}=0\) when \(a\lt r\), including negative \(a\). This matches the counting problem: all parts can be \(\ge k\) only if \(n\ge mk\), which is the same as \(n-m(k-1)-1\ge m-1\), so the binomial vanishes exactly when no such composition exists.
6. Counting case B, and the hockey-stick identity
In case B exactly one position holds the short part, so that position can be chosen in \(m\) ways with no double counting. Say it holds the value \(t\) with \(1\le t\le k-1\). The other \(m-1\) parts are all \(\ge k\) and sum to \(n-t\), so the substitution of Section 5 applied to \(m-1\) parts counts them as \(\binom{n-t-(m-1)(k-1)-1}{m-2}\). With \(M=n-(m-1)(k-1)-1\) that top is \(M-t\), and letting \(t\) run from \(1\) to \(k-1\) sweeps the tops \(M-1,\dots,M-k+1\):
$$B_k=m\sum_{i=M-k+1}^{M-1}\binom{i}{m-2}.$$
A consecutive run of binomial coefficients with a fixed lower index collapses, by the hockey-stick identity [5]
$$\sum_{i=r}^{U}\binom{i}{r}=\binom{U+1}{r+1},$$
whose proof is to classify the \((r+1)\)-element subsets of \(\{1,\dots,U+1\}\) by their largest element: if that element is \(i+1\), the other \(r\) members come from \(\{1,\dots,i\}\) in \(\binom{i}{r}\) ways, and letting \(i\) run over \(r,\dots,U\) counts every subset once. Applying the identity to both ends of our run with \(r=m-2\) gives
$$\begin{aligned}B_k&=m\binom{M}{m-1}\\&\quad-m\binom{M-k+1}{m-1},\end{aligned}$$
and \(M-k+1=n-m(k-1)-1\), the very top that appeared in case A, so case B is expressed through the same two binomial coefficients.
7. The tail probability, and where the sum stops
Abbreviate the two counts
$$\begin{aligned}T_k&=\binom{n-(m-1)(k-1)-1}{m-1},\\ U_k&=\binom{n-m(k-1)-1}{m-1}.\end{aligned}$$
Then \(A_k=U_k\) and \(B_k=m\,(T_k-U_k)\), and adding the two disjoint cases gives the count \(m\,T_k-(m-1)\,U_k\), hence
$$\Pr\bigl(X_{(2)}\ge k\bigr)=\frac{m\,T_k-(m-1)\,U_k}{\binom{n-1}{m-1}}.$$
Both tops start at \(n-1\) for \(k=1\) and fall linearly: the top of \(T_k\) by \(m-1\) at every step, the top of \(U_k\) by \(m\). By the rule of Section 5, \(T_k\ne0\) exactly while \(k\le(n-1)/(m-1)\) and \(U_k\ne0\) exactly while \(k\le n/m\); since \(m\le n\), the first bound is the larger, so
$$\begin{aligned}E(n,m)&=\sum_{k=1}^{K}\Pr\bigl(X_{(2)}\ge k\bigr),\\ K&=\Bigl\lfloor\frac{n-1}{m-1}\Bigr\rfloor,\end{aligned}$$
so also \(X_{(2)}\le K\) always β directly obvious, since \(m-1\) segments have length at least \(X_{(2)}\) and one more at least \(1\). For \(n=10^7\), \(m=100\) that is \(K=101010\) levels, with \(U_k\ne0\) only for the first \(\lfloor n/m\rfloor=100000\).
8. The checkpoint \(E(8,3)\) from the formula
With \(n=8\), \(m=3\) the denominator is \(\binom{7}{2}=21\).
- \(k=1\): \(T_1=U_1=\binom72=21\), count \(3\cdot21-2\cdot21=21\), probability \(1\).
- \(k=2\): \(T_2=\binom52=10\), \(U_2=\binom42=6\), count \(3\cdot10-2\cdot6=18\), probability \(\tfrac{18}{21}=\tfrac67\).
- \(k=3\): \(T_3=\binom32=3\), \(U_3=\binom12=0\), count \(9\), probability \(\tfrac{9}{21}=\tfrac37\).
- \(k=4\): \(T_4=\binom12=0\), so the sum stops, matching \(K=\lfloor7/2\rfloor=3\).
Hence \(E(8,3)=1+\tfrac67+\tfrac37=\tfrac{16}{7}\), the value found by enumeration in Section 1. The cases match that enumeration too: at \(k=2\), \(A_2=\binom42=6\) counts the orderings of \(\{2,2,4\}\) and \(\{2,3,3\}\) while \(B_2=3(10-6)=12\) counts those of \(\{1,2,5\}\) and \(\{1,3,4\}\); at \(k=3\), case A is empty, since all parts \(\ge3\) would need \(n\ge9\).
9. Binomial coefficients too large to build
The formula is exact, but evaluating it literally calls for arbitrary-precision integers: \(\binom{10^7-1}{99}\) has \(538\) decimal digits, and \(101010\) numbers of that size would have to be added. The published programs instead divide by the denominator as they go and track only the normalised quantities
$$R_k=\frac{T_k}{\binom{n-1}{m-1}},\qquad S_k=\frac{U_k}{\binom{n-1}{m-1}},$$
so that \(R_1=S_1=1\) and \(\Pr(X_{(2)}\ge k)=m\,R_k-(m-1)\,S_k\). Two binomial coefficients with the same lower index have a ratio that is a short product: for integers \(d\ge1\), \(r\ge1\) and \(a-d\ge r\), writing \(\binom{a}{r}=a(a-1)\cdots(a-r+1)/r!\) and cancelling \(r!\) gives
$$\frac{\binom{a-d}{r}}{\binom{a}{r}}=\prod_{j=0}^{r-1}\frac{a-d-j}{a-j}.$$
The hypothesis \(a-d\ge r\) keeps every numerator positive β for \(a=2\), \(r=2\), \(d=1\) the last factor would be \(0\) β and \(d\ge1\) makes each numerator smaller than its denominator, so every factor lies strictly between \(0\) and \(1\) and the running ratios only shrink. Once the next top falls below \(r\), the next ratio is \(0\) and no product is needed. Stepping from \(k\) to \(k+1\) uses \(d=m-1\) for \(R\) and \(d=m\) for \(S\), costing \(r=m-1\) divisions and \(r\) multiplications in the product plus one into the running ratio; no number above \(n\) is ever formed. For \(n=8\), \(m=3\) the first update multiplies \(R\) by \(\tfrac57\cdot\tfrac46=\tfrac{10}{21}\) and \(S\) by \(\tfrac47\cdot\tfrac36=\tfrac{6}{21}\), the normalised counts behind \(\Pr(X_{(2)}\ge2)=\tfrac67\) in Section 8.
10. A size estimate as a sanity check
It is worth knowing roughly how large the answer should be before trusting a program that prints it. Keep \(m\) fixed and let \(n\) grow. For \(a\) small beside \(M=n-1\) we have \(\binom{M-a}{r}\big/\binom{M}{r}\approx(1-a/M)^r\), so with \(u=(k-1)/n\)
$$\begin{aligned}&\Pr\bigl(X_{(2)}\ge k\bigr)\\&\approx m\bigl(1-(m-1)u\bigr)^{m-1}\\&\quad-(m-1)(1-mu)^{m-1},\end{aligned}$$
where each power must be read as \(0\) once its base turns negative, since by Section 5 the matching count is then \(0\). Summing over \(k\) is a Riemann sum of spacing \(1/n\) in \(u\), and each power integrates exactly over its own range, \(u\le1/(m-1)\) and \(u\le1/m\):
$$\begin{aligned}E(n,m)&\approx n\,c(m),\\ c(m)&=\frac{1}{m-1}-\frac{m-1}{m^2}\\&=\frac{2m-1}{m^2(m-1)}.\end{aligned}$$
Equivalently \(c(m)=\frac1m\bigl(\frac1m+\frac1{m-1}\bigr)\), the familiar shape of a second-smallest spacing. To six decimals the exact values are \(E(1000,10)=21.505398\), \(E(10^4,25)=33.125817\) and \(E(10^6,50)=808.643061\), each under half a unit above \(n\,c(m)\). None of this enters the computation; it only says the answer is of order \(n/m^2\), enough to catch a gross error.
How the Code Works
The C++ program [6] is organised around the two running ratios. An Options struct holds the defaults \(n=10^7\) and \(m=100\) together with a flag for the checkpoints, and parse_arguments recognises --n=, --m= and --skip-checkpoints, then requires \(n\ge1\), \(m\ge2\), \(m\le n\). Its helpers parse_ll_after_prefix and parse_int_after_prefix use std::stoll and std::stoi but never check that the whole tail was consumed, so --n=8junk is silently read as \(n=8\).
ratio_binom_drop(a, drop, choose_r) is the product of Section 9: writing \(r\) for choose_r, it multiplies the \(r\) factors \((a-\mathtt{drop}-j)/(a-j)\) for \(j=0,\dots,r-1\) and returns \(\binom{a-\mathtt{drop}}{r}/\binom{a}{r}\).
expected_second_shortest(n, m) runs the sum of Section 7. It sets \(r=m-1\), starts the tops a_r, a_s at \(n-1\) and the ratios ratio_r, ratio_s at \(1\), which is the level \(k=1\). Each pass first adds m * ratio_r - (m - 1) * ratio_s, the current tail probability, to answer, then advances each ratio separately. The guard a_r - (m - 1) >= r is exactly the condition of Section 5 for the next \(T_k\) to be nonzero, so when it fails the code sets ratio_r to \(0\) instead of multiplying; ratio_s is advanced the same way with the drop \(m\). A pass that advances neither sets has_next to false and breaks β correctly, since both ratios are then \(0\). These guards also supply the hypothesis \(a-d\ge r\) of Section 9.
run_checkpoints verifies the two values of the statement before any long computation: \(E(3,2)\) against \(2\) within relative tolerance \(10^{-14}\) and \(E(8,3)\) against \(16/7\) within \(10^{-13}\), both through nearly_equal, which compares the difference against the tolerance times the larger magnitude. main parses the options, runs the checkpoints unless skipped, and prints with std::fixed and std::setprecision(5).
The Python and Java files mirror the same two routines β Java names them ratioBinomDrop and expectedSecondShortest β with four differences worth knowing. Integer types: C++ keeps \(n\) and the two tops in long long, Java in long, Python in arbitrary-precision integers; at the target every integer stays below \(n\), so overflow is no concern, though the C++ options allow larger \(n\). Floating-point width: the two ratios, the per-level tail and the accumulator are long double in C++ and binary64 in Python and Java. [7] [8] long double need only be at least as wide as double, [9] and on the Apple-silicon machine used here it is binary64 too (LDBL_MANT_DIG\({}=53\)); where it is wider, the extra precision is carried through the loop and then dropped by the explicit static_cast<double> in the print statement. Output: both ports hard-code \(n\) and \(m\) in solve() and format with "{:.5f}" respectively String.format("%.5f"); the Java call ends with .replace(',', '.') because String.format without a locale uses the default one, which sets separator and digits alike, [10] and under an ar-EG default we measured a checkpoint printed in Arabic-Indic digits that the comma replacement does not repair. Checks: neither port runs the checkpoints.
Complexity and Verification
Put \(K=\lfloor(n-1)/(m-1)\rfloor\), \(L=\lfloor n/m\rfloor\) and \(H=K+L-2\). The loop makes \(K\) passes and advances its ratios \(K-1\) and \(L-1\) times, so the helper is called \(H\) times; each call does \(m-1\) divisions and \(m-1\) multiplications, its result is multiplied into the running ratio once, and every pass multiplies the two ratios by \(m\) and by \(m-1\). Counting the source's floating-point operations:
$$\begin{aligned}N_{\mathrm{div}}&=(m-1)H,\\ N_{\mathrm{mul}}&=mH+2K.\end{aligned}$$
Since \(K\approx n/(m-1)\) and \(L\approx n/m\), the division count is close to \(2n-n/m\): a larger \(m\) makes every step dearer but shortens the loop by the same factor. The work is \(O(n)\) for every admissible pair and \(\Theta(n)\) for fixed \(m\) as \(n\) grows, but not \(\Theta(n)\) uniformly β at \(m=n\) we get \(K=L=1\), \(H=0\) and a single pass. Storage is a constant number of scalars, with no array. For \(n=10^7\), \(m=100\): \(K=101010\), \(L=100000\), \(H=201008\), hence \(19\,899\,792\) divisions and \(20\,302\,820\) multiplications. Measured here, the optimised C++ build takes about \(0.02\) s and OpenJDK 17 about \(0.04\) s, while CPython 3.11 needs about \(1.0\) s for the same \(4\cdot10^7\) operations, most of it interpreter work rather than arithmetic.
The accuracy of the floating-point run needs its own argument, because five decimal places of a number of this size ask for about nine correct significant digits, and the errors of the two products build up over the whole loop. One update of a ratio performs \(2m-1\) roundings, so after \(k-1\) updates its relative error is at most about \((k-1)(2m-1)\varepsilon\) with \(\varepsilon=2^{-53}\approx1.11\cdot10^{-16}\), [7] [8] which bounds the absolute error of the sum by
$$\begin{aligned}&\varepsilon(2m-1)\sum_{k=1}^{K}(k-1)\,w_k,\\&w_k=m\,R_k+(m-1)S_k,\end{aligned}$$
which comes to about \(4.4\cdot10^{-6}\) at the target, plus \(2.3\cdot10^{-8}\) from the \(K\) additions β the same order as the \(5\cdot10^{-6}\) half-width of a rounded output, so the bound by itself certifies nothing. Underflow, which a relative-error bound does not cover, is harmless: the exact ratios sink below \(10^{-300}\) and finally below \(10^{-524}\), and those levels together contribute less than \(10^{-298}\).
What settles the question is an exact computation. Write the two tops of Section 7 as \(a_k=n-(m-1)(k-1)-1\) and \(b_k=n-m(k-1)-1\), so \(T_k=\binom{a_k}{r}\) and \(U_k=\binom{b_k}{r}\) with \(r=m-1\), and put \(F(a)=a(a-1)\cdots(a-r+1)\) for \(a\ge r\) and \(F(a)=0\) otherwise, giving \(\binom{a}{r}=F(a)/r!\). The factorials cancel:
$$\begin{aligned}N&=m\sum_{k=1}^{K}F(a_k)\\&\quad-(m-1)\sum_{k=1}^{K}F(b_k),\\ E(n,m)&=\frac{N}{F(n-1)}.\end{aligned}$$
At \(n=8\), \(m=3\) the two sums are \(42+20+6=68\) and \(42+12+0=54\) while \(F(7)=42\), giving \((3\cdot68-2\cdot54)/42=\tfrac{96}{42}=\tfrac{16}{7}\). Evaluated in exact integer arithmetic at \(n=10^7\), \(m=100\) and compared with the binary64 run, this ratio agrees to \(12\) significant digits: the absolute difference is \(4.9\cdot10^{-10}\), a relative error of \(2.4\cdot10^{-13}\), so the roundings plainly do not line up the way the worst case allows. The exact value also lies far from a rounding boundary, so the fifth decimal is unambiguous.
Beyond the C++ program's two checkpoints, these checks were run for this article. A brute force over all \(\binom{n-1}{m-1}\) cut choices, sorting each and averaging the second entry, agrees with the tail formula in all \(153\) cases with \(2\le n\le18\), \(2\le m\le n\). The case counts of Sections 5 and 6 were compared with direct enumeration for every triple with \(2\le n\le16\), \(2\le m\le n\), \(1\le k\le n+2\): all \(1600\) comparisons of \(A_k\), \(B_k\) and \(A_k+B_k\) passed, including the levels past the support, where \(\binom{a}{r}=0\) does the work. The published loop re-run with exact rationals instead of floating point reproduced the tail-sum formula for all \(3081\) pairs with \(2\le n\le79\), \(2\le m\le n\), which exercises the guards and the termination test on every shape of input in that range, and the same instrumented replay confirmed the operation counts above. Finally the three programs were built and run at the target and printed identical output, which the exact integer evaluation confirms; \(E(20,4)=\tfrac{53}{17}\) and \(E(30,5)=\tfrac{23512}{7917}\) were checked by brute force too. None of these extra checks is packaged in the published solutions.
Footnotes and References
The references below supply the background facts used above. The two case counts, their collapse to \(m\,T_k-(m-1)\,U_k\), the ratio recurrence, the error bound, the operation counts and the size estimate are derived explicitly in this article.
- Project Euler 398 β Cutting Rope. The official statement places n β 1 points at unit distance inside a rope of length n, cuts at m β 1 randomly chosen points to create m segments, defines E(n, m) as the expected length of the second-shortest segment with ties resolved so that it equals the shortest length, gives E(3, 2) = 2 and E(8, 3) = 16/7, and asks for E(10β·, 100) rounded to five decimal places.
- Wikipedia β Stars and bars (combinatorics). Its first theorem states that the number of k-tuples of positive integers summing to n equals the number of (k β 1)-element subsets of an (n β 1)-element set, that is C(n β 1, k β 1), and notes that these tuples are the compositions of n. Sections 1, 5 and 6 use this count.
- Wikipedia β Order statistic. The k-th order statistic of a sample is its k-th smallest value, written Xβββ and obtained by sorting into nondecreasing order. Section 2 uses this definition, whose handling of repeated values is the tie rule of the statement.
- Wikipedia β Expected value. The article expresses the expectation of a nonnegative integer-valued random variable through its tail probabilities as the sum of Pr(X > k) over k β₯ 0, the same series as the sum of Pr(X β₯ k) over k β₯ 1 used in Section 3.
- Wikipedia β Hockey-stick identity. It states that the sum of C(i, r) for i from r to U equals C(U + 1, r + 1), proved by classifying the (r + 1)-element subsets of a (U + 1)-element set by their largest element. Section 6 applies it twice, with the convention that a binomial coefficient whose lower index exceeds the upper one is zero.
- C++ β Euler398.cpp. The linked immutable revision contains
ratio_binom_drop,expected_second_shortestand the checkpoints for E(3, 2) and E(8, 3), and is the implementation record for this article. The exact-arithmetic comparison and the brute-force checks above are not packaged in it. - Wikipedia β Double-precision floating-point format. The format has a 53-bit significand, which "gives from 15 to 17 significant decimal digits precision (2β»β΅Β³ β 1.11 Γ 10β»ΒΉβΆ)" β the unit roundoff used in the error bound.
- Python β Floating-Point Arithmetic: Issues and Limitations. It states that almost all platforms map Python floats to IEEE-754 binary64 values with 53 bits of precision, and that errors are "on the order of no more than 1 part in 2**53 per operation" while "every float operation can suffer a new rounding error".
- cppreference β Fundamental types. It describes
long doubleas an extended-precision type that "does not necessarily map to types mandated by IEEE-754", listing implementations using the 80-bit x87 format, the 128-bit format, double-double or plain binary64 β which is why its width has to be measured, not assumed. - Java β java.util.Formatter. The
'f'conversion applies the localization algorithm, in which digits are replaced relative to the locale's zero digit and "a locale-specific decimal separator is substituted", and "if there is no explicit locale specified β¦ then the default locale is used".
Problem 398 source code
C++
#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>
namespace {
struct Options {
long long n = 10000000;
int m = 100;
bool run_checkpoints = true;
};
bool parse_ll_after_prefix(const std::string& arg, const std::string& prefix, long long& value) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
try {
value = std::stoll(tail);
} catch (...) {
return false;
}
return 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;
}
try {
value = std::stoi(tail);
} catch (...) {
return false;
}
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_ll_after_prefix(arg, "--n=", options.n) ||
parse_int_after_prefix(arg, "--m=", options.m)) {
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return options.n >= 1 && options.m >= 2 && options.m <= options.n;
}
long double ratio_binom_drop(long long a, int drop, int choose_r) {
// Computes C(a-drop, r) / C(a, r) by telescoping product.
long double ratio = 1.0L;
for (int j = 0; j < choose_r; ++j) {
ratio *= static_cast<long double>(a - drop - j) / static_cast<long double>(a - j);
}
return ratio;
}
long double expected_second_shortest(long long n, int m) {
// P(X_(2) >= k) = [ m*C(n-(m-1)(k-1)-1,m-1) - (m-1)*C(n-m(k-1)-1,m-1) ] / C(n-1,m-1).
// Then E[X_(2)] = sum_{k>=1} P(X_(2) >= k).
const int r = m - 1;
long long a_r = n - 1; // numerator top for first binomial ratio.
long long a_s = n - 1; // numerator top for second binomial ratio.
long double ratio_r = 1.0L;
long double ratio_s = 1.0L;
long double answer = 0.0L;
while (true) {
const long double tail = static_cast<long double>(m) * ratio_r -
static_cast<long double>(m - 1) * ratio_s;
answer += tail;
bool has_next = false;
if (a_r - (m - 1) >= r) {
ratio_r *= ratio_binom_drop(a_r, m - 1, r);
a_r -= (m - 1);
has_next = true;
} else {
ratio_r = 0.0L;
}
if (a_s - m >= r) {
ratio_s *= ratio_binom_drop(a_s, m, r);
a_s -= m;
has_next = true;
} else {
ratio_s = 0.0L;
}
if (!has_next) {
break;
}
}
return answer;
}
bool nearly_equal(long double a, long double b, long double rel_tol = 1e-12L) {
const long double scale = std::max(std::fabsl(a), std::fabsl(b));
if (scale == 0.0L) {
return true;
}
return std::fabsl(a - b) <= rel_tol * scale;
}
bool run_checkpoints() {
if (!nearly_equal(expected_second_shortest(3, 2), 2.0L, 1e-14L)) {
std::cerr << "Checkpoint failed: E(3,2)\n";
return false;
}
const long double e83 = expected_second_shortest(8, 3);
if (!nearly_equal(e83, 16.0L / 7.0L, 1e-13L)) {
std::cerr << "Checkpoint failed: E(8,3)\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;
}
const long double answer = expected_second_shortest(options.n, options.m);
std::cout << std::fixed << std::setprecision(5) << static_cast<double>(answer) << '\n';
return 0;
}
Python
def ratio_binom_drop(a, drop, choose_r):
ratio = 1.0
for j in range(choose_r):
ratio *= float(a - drop - j) / float(a - j)
return ratio
def expected_second_shortest(n, m):
r = m - 1
a_r = n - 1
a_s = n - 1
ratio_r = 1.0
ratio_s = 1.0
answer = 0.0
while True:
tail = float(m) * ratio_r - float(m - 1) * ratio_s
answer += tail
has_next = False
if a_r - (m - 1) >= r:
ratio_r *= ratio_binom_drop(a_r, m - 1, r)
a_r -= (m - 1)
has_next = True
else:
ratio_r = 0.0
if a_s - m >= r:
ratio_s *= ratio_binom_drop(a_s, m, r)
a_s -= m
has_next = True
else:
ratio_s = 0.0
if not has_next:
break
return answer
def solve():
n = 10000000
m = 100
ans = expected_second_shortest(n, m)
return "{:.5f}".format(ans)
if __name__ == '__main__':
print(solve())
Java
public class Euler398 {
private static double ratioBinomDrop(long a, int drop, int chooseR) {
double ratio = 1.0;
for (int j = 0; j < chooseR; ++j) {
ratio *= (double) (a - drop - j) / (double) (a - j);
}
return ratio;
}
private static double expectedSecondShortest(long n, int m) {
int r = m - 1;
long aR = n - 1;
long aS = n - 1;
double ratioR = 1.0;
double ratioS = 1.0;
double answer = 0.0;
while (true) {
double tail = (double) (m) * ratioR - (double) (m - 1) * ratioS;
answer += tail;
boolean hasNext = false;
if (aR - (m - 1) >= r) {
ratioR *= ratioBinomDrop(aR, m - 1, r);
aR -= (m - 1);
hasNext = true;
} else {
ratioR = 0.0;
}
if (aS - m >= r) {
ratioS *= ratioBinomDrop(aS, m, r);
aS -= m;
hasNext = true;
} else {
ratioS = 0.0;
}
if (!hasNext) {
break;
}
}
return answer;
}
public static String solve() {
long n = 10000000L;
int m = 100;
return String.format("%.5f", expectedSecondShortest(n, m)).replace(',', '.');
}
public static void main(String[] args) {
System.out.println(solve());
}
}