Detailed mathematical approach

Problem Summary

The silo has circular base radius \(R\). Grain is poured from a point whose vertical projection onto the base is \(P=(x,0)\), where \(x\in[0,R]\). Because the heap settles at angle of repose \(\alpha\), the empty region between the flat roof level and the conical surface has volume \(V(x)\). We must find every offset \(x\) for which \(V(x)\) is a perfect square, and then sum those offsets.

Mathematical Approach

1. Geometry of the Empty Volume

Let the base disk be

$D=\{(u,v)\in\mathbb{R}^2:u^2+v^2\le R^2\}.$

For a base point \(Q=(u,v)\), the horizontal distance to the pouring point is

$r=\sqrt{(u-x)^2+v^2}.$

The grain surface rises linearly with slope \(\tan(\alpha)\), so the local empty-space height above \(Q\) is \(r\tan(\alpha)\). Therefore the wasted volume is the surface integral

$V(x)=\tan(\alpha)\iint_D \sqrt{(u-x)^2+v^2}\,du\,dv.$

This is the quantity evaluated by the implementations.

2. Reduce the Double Integral to One Angular Integral

Use polar coordinates centered at \(P\):

$u=x+\rho\cos\varphi,\qquad v=\rho\sin\varphi.$

Along a fixed direction \(\varphi\), the ray leaves the disk when

$\left(x+\rho\cos\varphi\right)^2+\left(\rho\sin\varphi\right)^2=R^2.$

Solving this quadratic for the nonnegative root gives the radial intersection length

$s(\varphi)=-x\cos\varphi+\sqrt{R^2-x^2\sin^2\varphi}.$

The Jacobian contributes a factor \(\rho\), and the local height contributes another factor \(\rho\tan(\alpha)\). Hence

$V(x)=\tan(\alpha)\int_0^{2\pi}\int_0^{s(\varphi)} \rho^2\,d\rho\,d\varphi =\frac{\tan(\alpha)}{3}\int_0^{2\pi}s(\varphi)^3\,d\varphi.$

This is exactly the formula used in the C++, Python, and Java implementations.

3. Endpoint Checks and the Range of Square Targets

At the center, \(x=0\), every ray has the same length \(s(\varphi)=R\). So

$V(0)=\frac{\tan(\alpha)}{3}\int_0^{2\pi}R^3\,d\varphi=\frac{2\pi R^3}{3}\tan(\alpha).$

At the wall, \(x=R\), the ray length becomes

$s(\varphi)=-R\cos\varphi+R\lvert\cos\varphi\rvert,$

which is zero on the outward half-plane and equals \(-2R\cos\varphi\) on the inward half-plane. Therefore

$V(R)=\frac{\tan(\alpha)}{3}\int_{\pi/2}^{3\pi/2}\left(-2R\cos\varphi\right)^3\,d\varphi =\frac{32R^3}{9}\tan(\alpha).$

For the actual parameters \(R=6\) and \(\alpha=40^\circ\), this gives

$V(0)\approx 379.599730118848,\qquad V(R)\approx 644.428516744151.$

So the only possible square volumes are

$20^2,\,21^2,\,22^2,\,23^2,\,24^2,\,25^2.$

This is the key simplification: instead of searching over a continuum of offsets, we only need to solve six scalar equations \(V(x)=k^2\).

4. Worked Checkpoint Example

The sample checkpoint used by the implementations is \(R=3\) and \(\alpha=30^\circ\). Then

$V(0)=\frac{2\pi\cdot 3^3}{3}\tan(30^\circ)=\frac{18\pi}{\sqrt{3}}\approx 32.648388556,$

$V(R)=\frac{32\cdot 3^3}{9}\tan(30^\circ)=\frac{96}{\sqrt{3}}\approx 55.425625842.$

The only square targets in this interval are \(36\) and \(49\). Solving \(V(x)=36\) and \(V(x)=49\) numerically yields

$x\approx 1.114785284,\qquad x\approx 2.511167869,$

which matches the checkpoint values verified by the C++ implementation.

5. Numerical Integration and Root Search

No elementary antiderivative is used for the angular integral, so the implementation evaluates

$\int_0^{2\pi}s(\varphi)^3\,d\varphi$

with Simpson's rule using an even number \(N_q=4096\) of subintervals. If \(h=2\pi/N_q\) and \(f_i=f(ih)\), then

$\int_0^{2\pi}f(\varphi)\,d\varphi\approx \frac{h}{3}\left(f_0+f_{N_q}+4\sum_{j=1}^{N_q/2}f_{2j-1}+2\sum_{j=1}^{N_q/2-1}f_{2j}\right).$

After the endpoint values are known, the implementations enumerate every integer \(k\) with \(k^2\in[V(0),V(R)]\). For each target square they apply bisection on the interval \([0,R]\). The method relies on the smooth increase of \(V(x)\) from the center toward the wall, so each admissible square contributes one root. Using 120 bisection steps is far more than enough for nine correct decimal places.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they evaluate \(V(0)\) and \(V(R)\). Next they list all integer squares inside that volume range. For each target square they repeatedly bisect the offset interval and reevaluate the Simpson integral at the midpoint until the root is isolated. Finally they sum the resulting offsets and format the answer to nine decimal places. The C++ version also checks the sample values before solving the main instance.

Complexity Analysis

Let \(K\) be the number of square targets, \(B\) the number of bisection iterations, and \(N_q\) the Simpson subinterval count. One evaluation of \(V(x)\) costs \(O(N_q)\), so the full computation costs \(O(KBN_q)\) time and \(O(1)\) extra memory. For the actual problem, \(K=6\), \(B=120\), and \(N_q=4096\), so the runtime is dominated by a modest number of accurate integral evaluations.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=431
2. Angle of repose: Wikipedia — Angle of repose
3. Polar coordinates: Wikipedia — Polar coordinate system
4. Simpson's rule: Wikipedia — Simpson's rule
5. Bisection method: Wikipedia — Bisection method