Project Euler 431

题目

Square Space Silo

Fred the farmer arranges to have a new storage silo installed on his farm and having an obsession for all things square he is absolutely devastated when he discovers that it is circular. Quentin, the representative from the company that installed the silo, explains that they only manufacture cylindrical silos, but he points out that it is resting on a square base. Fred is not amused and insists that it is removed from his property.

Quick thinking Quentin explains that when granular materials are delivered from above a conical slope is formed and the natural angle made with the horizontal is called the angle of repose. For example if the angle of repose, \(\alpha = 30\) degrees, and grain is delivered at the centre of the silo then a perfect cone will form towards the top of the cylinder. In the case of this silo, which has a diameter of 6m, the amount of space wasted would be approximately \(32.648388556 \text{m}^3\). However, if grain is delivered at a point on the top which has a horizontal distance of \(x\) metres from the centre then a cone with a strangely curved and sloping base is formed. He shows Fred a picture.

We shall let the amount of space wasted in cubic metres be given by \(V(x)\). If \(x = 1.114785284\), which happens to have three squared decimal places, then the amount of space wasted, \(V(1.114785284) \approx 36\). Given the range of possible solutions to this problem there is exactly one other option: \(V(2.511167869) \approx 49\). It would be like knowing that the square is king of the silo, sitting in splendid glory on top of your grain.

Fred’s eyes light up with delight at this elegant resolution, but on closer inspection of Quentin’s drawings and calculations his happiness turns to despondency once more. Fred points out to Quentin that it’s the radius of the silo that is \(6\) metres, not the diameter, and the angle of repose for his grain is \(40\) degrees. However, if Quentin can find a set of solutions for this particular silo then he will be more than happy to keep it.

If Quick thinking Quentin is to satisfy frustratingly fussy Fred the farmer’s appetite for all things square then determine the values of \(x\) for all possible square space wastage options and calculate \(\sum x\) correct to \(9\) decimal places.

解决方案

把粮仓顶面当作平面区域\(D=\{(u,v)\mid u^2+v^2\le r^2\}.\)取坐标使圆心 \(O=(0,0)\)，投料点 \(P=(x,0)\)。

休止角定义为颗粒表面与水平面的夹角为 \(\alpha\)，因此在顶面任一点 \(Q=(u,v)\) 处，空腔高度（顶面到粮面）为\(z(Q)=\tan\alpha\cdot |PQ| =\tan\alpha\cdot \sqrt{(u-x)^2+v^2}.\)故浪费体积就是

\[ V(x)=\iint_D z(Q)dA=\tan\alpha\iint_D |PQ|dA. \]

以 \(P\) 为极点，那么有\(Q=P+(\rho\cos\theta,\rho\sin\theta), dA=\rho d\rho d\theta.\)对每个 \(\theta\)，射线与圆 \(u^2+v^2=r^2\) 的交点满足\((x+\rho\cos\theta)^2+(\rho\sin\theta)^2=r^2,\)即\(\rho^2+2x\cos\theta\rho+x^2-r^2=0.\)取正根，那么有：\(\rho_0(\theta)=-x\cos\theta+\sqrt{r^2-x^2\sin^2\theta}.\)于是

\[ \begin{aligned} V(x) &=\tan\alpha\int_0^{2\pi}\int_0^{\rho_0(\theta)} \rho\cdot \rho d\rho d\theta\\ &=\frac{\tan\alpha}{3}\int_0^{2\pi}\rho_0(\theta)^3d\theta\\ &=\frac{2\tan\alpha}{3}\int_0^\pi \rho_0(\theta)^3d\theta. \end{aligned} \]

可以通过计算证明，\(V(x)\)在\(x\in[0,r]\)是单调递增的。为计算所有\(V(x)\)为完全平方数下的\(x\)，先算两端的浪费体积，得到\(V(0),V(r)\)的值，所以可能的完全平方数范围是\(w^2=\lceil\sqrt{V(0)}\rceil^2,\dots,\lfloor\sqrt{V(r)}\rfloor^2.\)即只需解这些方程：\(V(x)=w^2.\)为此，我们对每个目标平方 \(w^2\) 用二分法解 \(V(x)=w^2\)。

代码

import math

from scipy.integrate import quad
from scipy.optimize import brentq

R = 6.0
A = math.radians(40.0)
T = math.tan(A)

def volume_quad(x: float) -> float:
    def f(theta: float) -> float:
        s = math.sin(theta)
        c = math.cos(theta)
        d = R * R - x * x * s * s
        if d < 0.0:
            d = 0.0
        rho = -x * c + math.sqrt(d)
        return rho * rho * rho
    v, err, *_ = quad(f, 0.0, math.pi, full_output=1, epsabs=1e-13, epsrel=1e-13, limit=400)
    return (2.0 * T / 3.0) * v

def solve():
    vf = volume_quad
    v0 = vf(0.0)
    vr = vf(R)
    n0 = math.ceil(v0 ** 0.5)
    n1 = math.floor(vr ** 0.5)
    roots = []
    lo = 0.0
    for n in range(n0, n1 + 1):
        target = float(n * n)
        x = brentq(lambda z: vf(z) - target, lo, R, xtol=1e-13, rtol=1e-13, maxiter=200)
        roots.append(x)
        lo = x
    return sum(roots)

if __name__ == "__main__":
    total = solve()
    print(f"{total:.9f}")