Problem 431: Square Space Silo

View on Project Euler

Project Euler Problem 431 Solution

EulerSolve provides an optimized solution for Project Euler Problem 431, Square Space Silo, with C++, Python, Java, and a step-by-step mathematical explanation.

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....

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

Problem 431 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = long long;

constexpr long double PI = 3.141592653589793238462643383279502884L;

struct Options {
    long double radius = 6.0L;
    long double alpha_deg = 40.0L;
    int simpson_n = 4096;
    bool run_checkpoints = 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_ld_after_prefix(const std::string& arg, const std::string& prefix, long double& 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::stold(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_ld_after_prefix(arg, "--radius=", options.radius) ||
            parse_ld_after_prefix(arg, "--alpha=", options.alpha_deg) ||
            parse_int_after_prefix(arg, "--simpson-n=", options.simpson_n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.radius > 0.0L && options.alpha_deg > 0.0L && options.alpha_deg < 90.0L &&
           options.simpson_n >= 128 && (options.simpson_n % 2 == 0);
}

long double integral_distance(long double x, long double radius, int simpson_n) {
    // I(x) = ∬_disk dist((u,v),(x,0)) dudv
    //      = (1/3) ∫_0^{2π} s(phi)^3 dphi,
    // where s(phi) is radial intersection in polar coordinates centered at (x,0).
    const long double a = 0.0L;
    const long double b = 2.0L * PI;
    const long double h = (b - a) / static_cast<long double>(simpson_n);

    long double acc = 0.0L;
    for (int i = 0; i <= simpson_n; ++i) {
        const long double phi = a + h * static_cast<long double>(i);
        const long double sinp = std::sin(phi);
        const long double cosp = std::cos(phi);
        const long double rad = std::sqrt(std::max(0.0L, radius * radius - x * x * sinp * sinp));
        const long double s = -x * cosp + rad;
        const long double f = s * s * s;

        int coeff = 2;
        if (i == 0 || i == simpson_n) {
            coeff = 1;
        } else if (i & 1) {
            coeff = 4;
        }
        acc += static_cast<long double>(coeff) * f;
    }

    return (h / 3.0L) * acc / 3.0L;
}

long double wasted_volume(long double x, long double radius, long double alpha_deg, int simpson_n) {
    const long double tan_alpha = std::tan(alpha_deg * PI / 180.0L);
    return tan_alpha * integral_distance(x, radius, simpson_n);
}

std::vector<long double> solve_x_values(long double radius, long double alpha_deg, int simpson_n) {
    const long double v0 = wasted_volume(0.0L, radius, alpha_deg, simpson_n);
    const long double vr = wasted_volume(radius, radius, alpha_deg, simpson_n);

    const int k_lo = static_cast<int>(std::ceill(std::sqrt(v0) - 1e-14L));
    const int k_hi = static_cast<int>(std::floorl(std::sqrt(vr) + 1e-14L));

    std::vector<long double> xs;
    for (int k = k_lo; k <= k_hi; ++k) {
        const long double target = static_cast<long double>(k) * static_cast<long double>(k);
        if (target < v0 - 1e-12L || target > vr + 1e-12L) {
            continue;
        }

        long double lo = 0.0L;
        long double hi = radius;
        for (int it = 0; it < 120; ++it) {
            const long double mid = (lo + hi) * 0.5L;
            const long double vm = wasted_volume(mid, radius, alpha_deg, simpson_n);
            if (vm < target) {
                lo = mid;
            } else {
                hi = mid;
            }
        }
        xs.push_back((lo + hi) * 0.5L);
    }

    return xs;
}

bool close_to(long double a, long double b, long double tol) {
    return std::fabsl(a - b) <= tol;
}

bool run_checkpoints(const int simpson_n) {
    const long double v0 = wasted_volume(0.0L, 3.0L, 30.0L, simpson_n);
    if (!close_to(v0, 32.648388556L, 5e-10L)) {
        std::cerr << "Checkpoint failed: V(0) for R=3, alpha=30\n";
        return false;
    }

    const auto xs = solve_x_values(3.0L, 30.0L, simpson_n);
    if (xs.size() != 2U) {
        std::cerr << "Checkpoint failed: expected two square solutions for sample case\n";
        return false;
    }
    if (!close_to(xs[0], 1.114785284L, 5e-10L)) {
        std::cerr << "Checkpoint failed: sample x for V=36\n";
        return false;
    }
    if (!close_to(xs[1], 2.511167869L, 5e-10L)) {
        std::cerr << "Checkpoint failed: sample x for V=49\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(options.simpson_n)) {
        return 2;
    }

    const auto xs = solve_x_values(options.radius, options.alpha_deg, options.simpson_n);
    long double sum = 0.0L;
    for (long double x : xs) {
        sum += x;
    }

    std::cout << std::fixed << std::setprecision(9) << static_cast<double>(sum) << '\n';
    return 0;
}

Python

import sys
import math

PI = math.pi

def integral_distance(x, radius, simpson_n):
    a = 0.0
    b = 2.0 * PI
    h = (b - a) / float(simpson_n)
    
    acc = 0.0
    for i in range(simpson_n + 1):
        phi = a + h * i
        sinp = math.sin(phi)
        cosp = math.cos(phi)
        
        val = radius * radius - x * x * sinp * sinp
        rad = math.sqrt(max(0.0, val))
        
        s = -x * cosp + rad
        f = s * s * s
        
        coeff = 2
        if i == 0 or i == simpson_n:
            coeff = 1
        elif i % 2 == 1:
            coeff = 4
            
        acc += coeff * f
        
    return (h / 3.0) * acc / 3.0

def wasted_volume(x, radius, alpha_deg, simpson_n):
    tan_alpha = math.tan(alpha_deg * PI / 180.0)
    return tan_alpha * integral_distance(x, radius, simpson_n)

def solve_x_values(radius, alpha_deg, simpson_n):
    v0 = wasted_volume(0.0, radius, alpha_deg, simpson_n)
    vr = wasted_volume(radius, radius, alpha_deg, simpson_n)
    
    k_lo = int(math.ceil(math.sqrt(v0) - 1e-14))
    k_hi = int(math.floor(math.sqrt(vr) + 1e-14))
    
    xs = []
    for k in range(k_lo, k_hi + 1):
        target = float(k) * float(k)
        if target < v0 - 1e-12 or target > vr + 1e-12:
            continue
            
        lo = 0.0
        hi = radius
        for _ in range(120):
            mid = (lo + hi) * 0.5
            vm = wasted_volume(mid, radius, alpha_deg, simpson_n)
            if vm < target:
                lo = mid
            else:
                hi = mid
        xs.append((lo + hi) * 0.5)
        
    return xs

def solve():
    radius = 6.0
    alpha_deg = 40.0
    simpson_n = 4096
    
    xs = solve_x_values(radius, alpha_deg, simpson_n)
    ans = sum(xs)
    return f"{ans:.9f}"

if __name__ == '__main__':
    print(solve())

Java

public class Euler431 {
    static final double PI = Math.PI;

    static double integralDistance(double x, double radius, int simpsonN) {
        double a = 0.0;
        double b = 2.0 * PI;
        double h = (b - a) / simpsonN;

        double acc = 0.0;
        for (int i = 0; i <= simpsonN; i++) {
            double phi = a + h * i;
            double sinp = Math.sin(phi);
            double cosp = Math.cos(phi);

            double val = radius * radius - x * x * sinp * sinp;
            double rad = Math.sqrt(Math.max(0.0, val));

            double s = -x * cosp + rad;
            double f = s * s * s;

            int coeff = 2;
            if (i == 0 || i == simpsonN) {
                coeff = 1;
            } else if (i % 2 == 1) {
                coeff = 4;
            }

            acc += coeff * f;
        }

        return (h / 3.0) * acc / 3.0;
    }

    static double wastedVolume(double x, double radius, double alphaDeg, int simpsonN) {
        double tanAlpha = Math.tan(alphaDeg * PI / 180.0);
        return tanAlpha * integralDistance(x, radius, simpsonN);
    }

    static java.util.List<Double> solveXValues(double radius, double alphaDeg, int simpsonN) {
        double v0 = wastedVolume(0.0, radius, alphaDeg, simpsonN);
        double vr = wastedVolume(radius, radius, alphaDeg, simpsonN);

        int kLo = (int) Math.ceil(Math.sqrt(v0) - 1e-14);
        int kHi = (int) Math.floor(Math.sqrt(vr) + 1e-14);

        java.util.List<Double> xs = new java.util.ArrayList<>();
        for (int k = kLo; k <= kHi; k++) {
            double target = (double) k * (double) k;
            if (target < v0 - 1e-12 || target > vr + 1e-12)
                continue;

            double lo = 0.0;
            double hi = radius;
            for (int it = 0; it < 120; it++) {
                double mid = (lo + hi) * 0.5;
                double vm = wastedVolume(mid, radius, alphaDeg, simpsonN);
                if (vm < target)
                    lo = mid;
                else
                    hi = mid;
            }
            xs.add((lo + hi) * 0.5);
        }

        return xs;
    }

    public static String solve() {
        double radius = 6.0;
        double alphaDeg = 40.0;
        int simpsonN = 4096;

        java.util.List<Double> xs = solveXValues(radius, alphaDeg, simpsonN);
        double sum = 0.0;
        for (double x : xs) {
            sum += x;
        }

        return String.format(java.util.Locale.US, "%.9f", sum);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}