Problem 394: Eating Pie

Detailed mathematical approach

Problem Summary

Jeff eats a circular pie, beginning with one cut along a radius. Then, while at least a given fraction \(F\) of the pie is left, he repeats a round: two slices from the centre to two points of the remaining border, every point equally likely, divide the remains into three pieces, and going anticlockwise from the initial cut he eats the first two. As soon as less than \(F\) remains he stops repeating and eats the rest. For \(x\ge1\) let \(E(x)\) be the expected number of times the round is repeated when \(F=1/x\). The statement records \(E(1)=1\), \(E(2)\approx1.2676536759\) and \(E(7.5)\approx2.1215732071\), and asks for \(E(40)\) rounded to ten decimal places. ^[1]

The whole solution in one sentence: one round leaves a random fraction of the pie with density \(2(1-w)\), so averaging over one round turns the expected count into an equation whose solution is \(\frac79+\frac23\ln x+\frac{2}{9x^3}\).

Mathematical Approach

1. What is left is always a sector

Measure angles anticlockwise from the initial cut, so the pie occupies \(0\) to \(2\pi\). In the first round the new cuts sit at angles \(a\), \(b\); write \(u=\min(a,b)\), \(v=\max(a,b)\). The three pieces are the sectors \(0\) to \(u\), \(u\) to \(v\), \(v\) to \(2\pi\), also the order met walking anticlockwise from the initial cut. Jeff eats the first two, so the survivor is \(v\) to \(2\pi\).

That survivor still has the initial cut as one of its straight edges, so the remains are always a sector \(c\) to \(2\pi\). The new cuts land between \(c\) and \(2\pi\), and walking anticlockwise from the initial cut one crosses the eaten gap and then meets \(c\) to \(u\), \(u\) to \(v\), \(v\) to \(2\pi\); the survivor is again the last, touching the initial cut from the clockwise side. The rule is the same every round, so the state is one number: the fraction \(p\) of the whole pie still there.

2. The shrinking factor of one round

An example first: cuts one quarter and three quarters along the remaining arc give pieces of relative size \(\frac14\), \(\frac12\), \(\frac14\), and the round keeps the last, a quarter of what was there. In general let the sector span an angle \(\theta\) and write the cuts as the fractions \(U_1\), \(U_2\) along its arc, independent and uniform on \([0,1]\). By Section 1 the survivor spans \(\theta\bigl(1-\max(U_1,U_2)\bigr)\), so the surviving fraction of the previous remains is

$$W=1-\max(U_1,U_2)=\min(1-U_1,1-U_2).$$

Each \(1-U_i\) is again uniform on \([0,1]\), so \(W\) is the smaller of two independent uniform numbers, and \(W\) exceeds \(w\) exactly when both \(1-U_i\) do:

$$P(W\le w)=1-(1-w)^2=2w-w^2,\quad 0\le w\le1 .$$

Differentiating, the density of \(W\) on \([0,1]\) is \(2(1-w)\): the Beta\((1,2)\) law, the case \(k=1\), \(n=2\) of the \(k\)-th smallest of \(n\) uniform numbers. ^[2] Its mean \(\int_0^1 2w(1-w)\,dw=\frac13\) says a round eats two thirds of what was there on average.

The law of \(W\) does not depend on \(\theta\): a sector behaves like a whole pie, only scaled, and that invariance keeps the state one-dimensional. All three pieces share this law: the first exceeds \(w\) when both cuts lie beyond \(w\) and the last when both lie before \(1-w\), each with probability \((1-w)^2\), while the middle exceeds \(w\) when \(|U_1-U_2|\gt w\), two triangles of area \(\frac12(1-w)^2\) in the unit square of cut pairs. So keeping any position fixed in advance — not, say, always the largest piece — gives the same \(E(x)\).

3. One number is enough: the recursion

Jeff performs a round exactly when \(p\ge F\), and a round replaces \(p\) by \(pW\) with a fresh independent \(W\). With \(F=1/x\) the test \(p\ge1/x\) is \(px\ge1\), which invites the rescaled state \(y=px\): the process starts at \(y=x\) and each round, performed whenever \(y\ge1\), multiplies \(y\) by \(W\). Let \(E(y)\) be the expected number of rounds still to come from state \(y\). Then \(E(y)=0\) for \(y\lt1\), and for \(y\ge1\) one round happens at once and leaves the state \(yW\), so averaging over the density of \(W\) gives

$$E(x)=1+2\int_0^1(1-w)E(xw)\,dw .$$

The quantity asked for is this same \(E\) at \(y=x\), since the process starts at \(p=1\). At \(x=1\) the factor \(E(w)\) vanishes for every \(w\lt1\), the integral is \(0\) and \(E(1)=1\): the statement's first checkpoint, for free.

4. Counting by tails: a worked example at \(x=2\)

Let \(N\) be the number of rounds, \(W_1,W_2,\dots\) the successive factors and \(Z_k=W_1\cdots W_k\). After \(k\) rounds the state is \(xZ_k\), so round \(k+1\) happens exactly when \(xZ_k\ge1\). A trial of \(n\) rounds makes exactly the \(n\) statements \(N\ge1,\dots,N\ge n\) true, so averaging their indicators gives \(\mathbb{E}[N]=\sum_{j\ge1}P(N\ge j)\), that is

$$E(x)=\sum_{k\ge0}P\!\left(Z_k\ge\frac1x\right),$$

the empty product \(Z_0=1\) contributing the leading term \(1\).

The series also settles the regularity the next sections need. The factors are independent with mean \(\frac13\), so \(\mathbb{E}[Z_k]=3^{-k}\), and on the event \(Z_k\ge1/x\) we have \(xZ_k\ge1\), whence

$$P\!\left(Z_k\ge\frac1x\right)\le x\,\mathbb{E}[Z_k]=x\,3^{-k}.$$

On \([1,B]\) the terms are dominated by \(B\,3^{-k}\), so the series converges uniformly and \(E\) is finite; each term with \(k\ge1\) is continuous because \(Z_k\) has a density, so \(E\) is continuous on \([1,\infty)\).

Now take \(x=2\). The term \(k=1\) is \(P(W\ge\frac12)=(1-\frac12)^2=\frac14\). For \(k=2\), condition on \(W_1=w\) and use \(P(W\ge u)=(1-u)^2\):

$$P(Z_2\ge c)=\int_{c}^{1}2(1-w)\Bigl(1-\frac{c}{w}\Bigr)^{2}dw .$$

With \(c=\frac12\) the integrand expands to \(4-\frac{5}{2w}+\frac{1}{2w^2}-2w\), whose antiderivative is \(4w-\frac52\ln w-\frac{1}{2w}-w^2\); evaluating it between \(\frac12\) and \(1\) leaves \(\frac74-\frac52\ln2\approx0.0171320486\). The same conditioning gives \(0.000513055562\) and \(0.0000084826627\) for \(k=3\) and \(k=4\). The running totals are \(1\), \(1.25\), \(1.26713204860\), \(1.26764510416\), \(1.26765358682\), and seven terms already reach \(1.26765367593\). This agrees with the statement's \(E(2)\approx1.2676536759\), a strong check on the model and the indexing of rounds. The series is exact but needs more terms the larger \(x\) is; we want a closed form of fixed cost.

5. From the recursion to a Volterra equation

Substitute \(t=xw\) in the integral of Section 3, so that \(w=t/x\) and \(dw=dt/x\). Since \(E(t)=0\) below \(t=1\), the lower limit may be raised to \(1\):

$$E(x)=1+\frac{2}{x^2}\int_1^x(x-t)E(t)\,dt .$$

The unknown now appears as \(E(t)\) under a fixed lower limit. This is a linear Volterra integral equation of the second kind, \(\varphi(x)=f(x)+\int_1^xK(x,t)\varphi(t)\,dt\), with \(f\equiv1\) and the continuous kernel \(K(x,t)=2(x-t)/x^{2}\). ^[3]

It has at most one continuous solution. Suppose \(E_1\) and \(E_2\) both solve it and put \(D=E_1-E_2\); subtracting kills the \(1\), so \(D(x)=\frac{2}{x^2}\int_1^x(x-t)D(t)\,dt\). Fix \(B\gt1\) and let \(c\) be the largest \(|D|\) on \([1,B]\). There

$$\begin{aligned}|D(x)|&\le\frac{2c}{x^2}\int_1^x(x-t)\,dt\\&=c\Bigl(1-\frac1x\Bigr)^{2}.\end{aligned}$$

The last factor increases with \(x\), so \(c\le c\bigl(1-\frac1B\bigr)^{2}\), which forces \(c=0\); as \(B\) was arbitrary, any continuous solution is the answer. The equation also gives the smoothness Section 6 needs: \(E\) is continuous by Section 4, so \(I(x)=\int_1^x(x-t)E(t)\,dt\) is twice differentiable, and hence so is \(E=1+2I/x^{2}\).

6. Two differentiations give an Euler equation

With \(I\) as above, \(x^2\bigl(E(x)-1\bigr)=2I(x)\). Both the upper limit and the integrand depend on \(x\), so the Leibniz rule applies, ^[4] and the boundary term \((x-x)E(x)\) it produces is \(0\):

$$I'(x)=\int_1^xE(t)\,dt,\qquad I''(x)=E(x).$$

Differentiating \(x^2(E-1)=2I\) once and then once more gives

$$2x\bigl(E(x)-1\bigr)+x^2E'(x)=2\int_1^xE(t)\,dt,$$

$$2\bigl(E(x)-1\bigr)+4xE'(x)+x^2E''(x)=2E(x).$$

In the second line the two copies of \(E(x)\) cancel, and what survives is a differential equation with no integral left in it:

$$x^2E''(x)+4xE'(x)=2,\qquad x\gt1 .$$

The same identities at \(x=1\) supply the starting data: the undifferentiated one gives \(E(1)=1\), since \(I(1)=0\), and the once-differentiated one becomes \(2\bigl(E(1)-1\bigr)+E'(1)=0\), so \(E'(1)=0\). Derivatives at \(x=1\) are one-sided, from the right: the extension of Section 3 jumps there.

7. Solving the equation

Dropping the right-hand side leaves \(x^2E''+4xE'=0\), a Cauchy–Euler or equidimensional equation. Trying \(E=x^{m}\) turns it into the indicial equation \(m(m-1)+4m=m(m+3)=0\), ^[5] so the homogeneous solutions are the constants, from \(m=0\), and \(x^{-3}\), from \(m=-3\). For the rest, set \(Y=E'\), which satisfies the first-order equation \(x^2Y'+4xY=2\). Multiplying by \(x^2\) makes the left side an exact derivative:

$$\bigl(x^4Y\bigr)'=2x^2,\qquad x^4Y=\tfrac23x^3+C_1 .$$

So \(E'(x)=\frac{2}{3x}+C_1x^{-4}\), and integrating once more

$$E(x)=\frac23\ln x-\frac{C_1}{3x^3}+C_2 .$$

Now \(E'(1)=0\) forces \(\frac23+C_1=0\), that is \(C_1=-\frac23\); then \(E(1)=1\) forces \(\frac29+C_2=1\), that is \(C_2=\frac79\). The closed form is

$$E(x)=\frac79+\frac23\ln x+\frac{2}{9x^3},\qquad x\ge1 .$$

8. Why this really solves the original equation

Differentiating can admit functions the original equation did not, so the candidate must go back in. Two elementary antiderivatives, both vanishing at \(x=1\), are

$$\int_1^xE(t)\,dt=\frac{6x^3\ln x+x^3-1}{9x^2},$$

$$\int_1^xtE(t)\,dt=\frac{3x^3\ln x+2x^3-2}{9x}.$$

Taking \(x\) times the first minus the second,

$$\int_1^x(x-t)E(t)\,dt=\frac{3x^{3}\ln x-x^{3}+1}{9x},$$

and multiplying by \(2/x^{2}\) yields \(\frac23\ln x-\frac29+\frac{2}{9x^{3}}\), which is exactly \(E(x)-1\). So

$$\frac{2}{x^2}\int_1^x(x-t)E(t)\,dt=E(x)-1,$$

the equation of Section 5 holds, and by the uniqueness proved there nothing else does.

9. A second derivation: the renewal transform

The closed form hides where \(\frac23\), \(\frac79\) and \(\frac29\) come from; taking logarithms exposes all three at once. Put \(X_i=-\ln W_i\ge0\) and \(S_k=X_1+\cdots+X_k\), with \(S_0=0\). The condition \(Z_k\ge1/x\) of Section 4 becomes \(S_k\le\ln x\), so with \(T=\ln x\) and

$$M(T)=\sum_{k\ge0}P(S_k\le T)$$

we have \(E(x)=M(\ln x)\). The \(X_i\) are independent, identically distributed waiting times and \(M\) their renewal function, counting the renewal at time \(0\) too. ^[6] Splitting the first waiting time off each term with \(k\ge1\) gives the renewal equation

$$M(T)=1+\int_0^T f(\tau)M(T-\tau)\,d\tau,$$

where \(f\) is the density of \(X\); substituting \(w=e^{-t}\) in Section 2 turns it into

$$f(t)=2\bigl(e^{-t}-e^{-2t}\bigr),\qquad t\ge0 .$$

A Laplace transform weights a function exponentially and integrates; for \(\lambda\gt0\),

$$\widetilde M(\lambda)=\int_0^{\infty}e^{-\lambda T}M(T)\,dT .$$

One elementary integral does all the work. With \(\varphi(s)=\mathbb{E}\bigl[W^{s}\bigr]\), for \(s\gt-1\),

$$\varphi(s)=\int_0^1 2(1-w)w^{s}\,dw=\frac{2}{(s+1)(s+2)} .$$

Since \(W^{s}=e^{-sX}\), this \(\varphi\) is also the Laplace transform of \(X\). The renewal integral is a convolution, so transforming both sides with \(\mathcal{L}\{1\}=1/\lambda\) ^[7] gives \(\widetilde M=\frac1\lambda+\varphi\widetilde M\), hence

$$1-\varphi(\lambda)=\frac{\lambda(\lambda+3)}{(\lambda+1)(\lambda+2)},$$

$$\begin{aligned}\widetilde M(\lambda)&=\frac{1}{\lambda\bigl(1-\varphi(\lambda)\bigr)}\\&=\frac{(\lambda+1)(\lambda+2)}{\lambda^{2}(\lambda+3)}.\end{aligned}$$

Partial fractions split the result into three standard pieces,

$$\frac{2/3}{\lambda^{2}}+\frac{7/9}{\lambda}+\frac{2/9}{\lambda+3},$$

and the table entries \(\mathcal{L}\{t\}=1/\lambda^{2}\), \(\mathcal{L}\{1\}=1/\lambda\) and \(\mathcal{L}\{e^{-3t}\}=1/(\lambda+3)\) ^[7] invert them term by term:

$$M(T)=\frac23T+\frac79+\frac29e^{-3T} .$$

Setting \(T=\ln x\), so that \(e^{-3T}=x^{-3}\), reproduces Section 7 exactly, all three constants at once. The exponents match too: the denominator \(\lambda^{2}(\lambda+3)\) vanishes at \(0\) and \(-3\), the roots of \(m(m+3)=0\), and the double zero at \(0\) produces the term linear in \(T\), logarithmic in \(x\). This rational function extends \(\varphi\) beyond the range \(s\gt-1\) where its integral converges, and \(-3\) lies outside that range, so the match concerns the extension, not the expectation.

10. Reading the formula, and why the small term matters

For this law the behaviour of \(\varphi\) at \(\lambda=0\) explains the first two terms. The moments

$$\mu=\mathbb{E}[X]=\tfrac32,\qquad m_2=\mathbb{E}\bigl[X^{2}\bigr]=\tfrac72,$$

which follow from \(\varphi'(0)=-\mu\) and \(\varphi''(0)=m_2\), give \(1-\varphi(\lambda)=\mu\lambda-\frac12m_2\lambda^{2}+\cdots\) and so

$$\widetilde M(\lambda)=\frac{1}{\mu\lambda^{2}}+\frac{m_2}{2\mu^{2}\lambda}+\cdots .$$

Comparing with the exact partial fractions, the two leading coefficients are

$$\frac1\mu=\frac23,\qquad \frac{m_2}{2\mu^{2}}=\frac79 ,$$

the slope and the constant of the closed form. In words, each round divides the pie by \(e^{3/2}\approx4.4817\) in the geometric mean, so rounds accumulate at rate \(\frac23\) per unit of \(\ln x\). Such an expansion says nothing about the remainder, but Section 9 inverted the whole rational transform, so here it is exactly \(\frac{2}{9}x^{-3}\).

That term must not be thrown away. At \(x=40\) it is about \(3.47\cdot10^{-6}\), far above the half unit \(5\cdot10^{-11}\) of the tenth decimal, so using only \(\frac23\ln x+\frac79\) would corrupt the sixth. The value there also lies more than \(4\cdot10^{-11}\) from a halfway point, so the three languages, whose round-to-nearest rules break ties differently, cannot print different digits.

How the Code Works

All three programs evaluate the closed form of Section 7; the mathematics above is the whole algorithm. The C++ file ^[8] is the fullest. An Options struct carries the argument x, default \(40\), and the flag run_checkpoints. parse_ld_after_prefix tests the prefix and converts the tail with std::stold, rejecting an empty tail or a throwing conversion. parse_arguments walks argv, accepts --skip-checkpoints and --x=, reports anything else on standard error, and requires options.x to be at least \(1\), the statement's condition.

That is a lower bound, not a full check: std::stold is never asked how much of the tail it consumed, so --x=2junk reads as \(2\) and --x=1e as \(1\), while --x=inf passes the comparison and the program prints inf. Values below \(1\) and --x=nan are rejected silently, no comparison with a NaN holding.

expected_steps is the one line that matters: in long double it adds 7.0L / 9.0L, (2.0L / 3.0L) * std::log(x) and 2.0L / (9.0L * x * x * x) — the three terms of Section 7, every literal a floating-point constant so that nothing is divided in integers — and a comment records the recursion of Section 3. nearly_equal compares with a relative tolerance against the larger magnitude, and run_checkpoints tests \(E(1)=1\) to \(10^{-14}\) and \(E(2)\), \(E(7.5)\) against the statement's ten-digit constants to \(10^{-10}\); that looser tolerance is necessary, not sloppy, since the true values differ from those constants by \(2.9\cdot10^{-11}\) and \(1.0\cdot10^{-11}\), the statement's own rounding. main runs the checkpoints unless switched off, exits with status \(2\) on failure, and prints with std::fixed and std::setprecision(10).

Python keeps only expected_steps and solve, which formats expected_steps(40.0) with "{:.10f}"; Java has expectedSteps and solve, calling String.format("%.10f", …) then .replace(',', '.'). Three differences between the files are real.

• Precision. C++ evaluates in long double and casts to double before printing; Python and Java work in binary64 throughout. What the cast discards is platform-dependent: here long double has the same \(53\)-bit significand as double, so nothing, and where it is wider the extra bits sit far below the last requested digit. All three reproduce the statement's checkpoints to ten decimals.
• Locale. String.format uses the virtual machine's default formatting locale when none is given, and localises the decimal separator and the digits. ^[9] A German or French default gives a comma, which .replace(',', '.') repairs — but not every locale: Arabic (Egypt) yields non-ASCII digits and another separator, both left untouched. Python's "{:.10f}" is locale-independent, and a C++ stream keeps the classic locale unless the program changes it.
• Scope. Only the C++ version takes an argument and only it packages the checkpoints; the other two hard-code \(x=40\) and check nothing.

Complexity and Verification

Time and memory are both \(O(1)\): one logarithm, a few multiplications and divisions, a fixed set of variables — Sections 5 to 9 traded a series and an integral equation for a formula.

Simulation is hopeless at this precision. For renewal counts the variance grows like \(\sigma^{2}T/\mu^{3}\) with \(T=\ln x\), ^[6] and here \(\sigma^{2}=m_2-\mu^{2}=\frac54\) and \(\mu^{3}=\frac{27}{8}\), so the rate is \(\frac{10}{27}\ln x\). Solving the second-moment renewal equation symbolically gives variance \(0.87616\) at \(x=10\), where \(2\cdot10^{6}\) trials gave \(0.875\). A sample mean would need about \(3.5\cdot10^{20}\) trials to push its standard error below half a unit in the tenth decimal.

Only the C++ version packages checks — the three checkpoints, run unless --skip-checkpoints is given. It passes them under GCC and Clang, and its --x= option reproduces \(E(1)\), \(E(2)\) and \(E(7.5)\) digit for digit; calling Python's expected_steps and Java's expectedSteps at those arguments gives the same digits, though neither file tests itself. All three print one identical line. The checks below were made for this article.

• The geometric model. Simulating the literal cuts — angles, three pieces, eat the first two anticlockwise — for \(300\,000\) trials at each of \(x=1,2,7.5,10,25\) matched the closed form within two standard errors every time; the largest gap was \(1.63\cdot10^{-3}\) against a standard error of \(0.88\cdot10^{-3}\), at \(x=2\).
• The series of Section 4. Its first seven terms, in exact arithmetic, sum to within \(3.4\cdot10^{-12}\) of \(\frac79+\frac23\ln2+\frac{1}{36}\).
• The substitution of Section 8, symbolically, in the Volterra form and in the untransformed recursion of Section 3: both residuals simplify to exactly \(0\).
• The integral equation of Section 5, marched with the trapezoidal rule ^[3] over \([1,40]\); the kernel vanishes at \(t=x\), so the scheme is explicit. Its error against the closed form fell from \(4.8\cdot10^{-8}\) to \(1.9\cdot10^{-10}\) as the uniform grid went from \(2\cdot10^{4}\) to \(3.2\cdot10^{5}\) steps — the factor \(16\) per fourfold refinement a second-order rule should give.
• The renewal equation of Section 9, solved the same way in \(T\) over \([0,\ln40]\), converging to \(\frac23T+\frac79+\frac29e^{-3T}\) at the same rate.
• The exponents. The denominator \(\lambda^{2}(\lambda+3)\) and the indicial polynomial \(m(m+3)\) give the same two numbers.

Footnotes and References

The references below supply the standard background: order statistics, Volterra equations, differentiation under the integral sign, Cauchy–Euler equations, renewal theory and Laplace transforms. The recursion, closed form, finiteness and uniqueness arguments and transform calculation are derived above.

1. Project Euler 394 — Eating Pie. The official statement describes the circular pie, the radial first cut and the repeated round that splits the remains into three pieces at two equally likely border points, the first two counterclockwise from the initial cut being eaten; it defines E(x) for F = 1/x, records E(1) = 1, E(2) ≈ 1.2676536759 and E(7.5) ≈ 2.1215732071, and asks for E(40) to ten decimals.
2. Wikipedia — Order statistic. The article gives the k-th order statistic of n independent standard uniform variables as Beta(k, n+1−k), with density n!/((k−1)!(n−k)!) · u^(k−1) (1−u)^(n−k); the case k = 1, n = 2 is the density 2(1−w) of Section 2.
3. Wikipedia — Volterra integral equation. The article gives the normal form x(t) = f(t) + ∫ K(t,s) x(s) ds of a linear second-kind equation, the shape reached in Section 5, and its trapezoidal discretisation, used in the verification.
4. Wikipedia — Leibniz integral rule. The article gives the derivative of an integral with variable limits as f(x, b(x)) b′(x) − f(x, a(x)) a′(x) plus the integral of the partial derivative of f, which Section 6 applies twice.
5. Wikipedia — Cauchy–Euler equation. The article treats the equidimensional equation and gives, for the second-order form with coefficients 1, a, b, the indicial equation m² + (a−1)m + b = 0 from the trial x^m; with a = 4, b = 0 that is m(m+3) = 0.
6. Wikipedia — Renewal theory. The article defines the renewal function m(t) as the expected number of renewals after time zero and gives m(t) = F(t) + ∫ f(s) m(t−s) ds. Section 9 counts the renewal at time zero too, so M = 1 + m and the equation becomes M(t) = 1 + ∫ f(s) M(t−s) ds. The elementary renewal theorem there explains the reciprocal-mean slope, and the central limit theorem for renewal counts supplies the variance rate t σ²/μ³.
7. Wikipedia — Laplace transform. The properties table gives the convolution rule, and the table of selected transforms lists the unit step as 1/s, the ramp t as 1/s² and the decay e^(−αt) as 1/(s+α) — the three entries Section 9 inverts.
8. C++ — Euler394.cpp. The linked immutable revision contains expected_steps, nearly_equal, run_checkpoints, parse_arguments and the two command-line options discussed above. It is the implementation record for this article.
9. Java — java.util.Formatter. The documentation states that with no locale given the default formatting locale is used, and that number localisation replaces each digit relative to the locale's zero digit and substitutes a locale-specific decimal separator.