Problem 317: Firecracker

Project Euler Problem 317 Solution

EulerSolve provides an optimized solution for Project Euler Problem 317, Firecracker, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A firecracker explodes at height $$h=100\ \text{m}$$ above the ground. Every tiny fragment leaves the explosion point with the same speed $$v=20\ \text{m/s}$$ but in an arbitrary direction. Gravity is constant: $$g=9.81\ \text{m/s}^2.$$ Ignoring air resistance, we must find the volume of the whole region visited by the fragments before they hit the ground. Mathematical Approach 1) Rotational symmetry reduces the problem to a 2D cross-section. Because the launch speed is the same in every direction, the reachable region is symmetric around the vertical axis through the explosion point. So it is enough to work in cylindrical coordinates: $$r=\sqrt{x^2+z^2},\qquad y=\text{height above ground}.$$ Once we know the boundary curve $y(r)$ in the $(r,y)$-plane, the full 3D region is obtained by revolving that curve around the $y$-axis. 2) One projectile trajectory. Let $\theta$ be the launch angle above the horizontal. Then the fragment satisfies $$r(t)=v\cos\theta\; t,$$ $$y(t)=h+v\sin\theta\; t-\frac12gt^2.$$ Eliminating $t$ using $t=r/(v\cos\theta)$ gives the familiar projectile parabola $$y(r,\theta)=h+r\tan\theta-\frac{gr^2}{2v^2\cos^2\theta}.$$ 3) Rewrite with $u=\tan\theta$. Since $$\frac1{\cos^2\theta}=1+\tan^2\theta=1+u^2,$$ we can write $$y(r,u)=h+ru-\frac{gr^2}{2v^2}(1+u^2).$$ For fixed $r$, this is a concave quadratic polynomial in $u$....

Detailed mathematical approach

Problem Summary

A firecracker explodes at a height of $h=100$ metres above level ground and breaks into a very large number of very small fragments. They leave the burst point in every direction, all with the same initial speed $v=20$ metres per second, fly without air resistance in a uniform gravitational field with $g=9.81$ metres per second squared, and stop on reaching the ground. We must find the volume, in cubic metres, of the region through which they move before reaching the ground, rounded to four decimal places. The statement gives no checkpoint value. ^[1]

The whole solution in one sentence: every fragment flies along a parabola, and all these flights together fill exactly the region between the ground and one downward parabola — their envelope — spun about the vertical axis through the burst.

Mathematical Approach

Nothing has to be simulated. We describe one flight, find which points a flight can reach, show that these fill a paraboloid, and integrate. Throughout $y$ is the height above the ground and $r\ge0$ the distance from the axis through the burst, at $r=0$, $y=h$.

1. One fragment

Let a fragment leave at the angle $\theta$ above the horizontal, $-90^\circ\le\theta\le90^\circ$. Horizontal and vertical motion are independent, so after $t$ seconds

$$\begin{aligned}r&=v\cos\theta\,t,\\ y&=h+v\sin\theta\,t-\tfrac12gt^2,\end{aligned}$$

the standard equations of projectile motion. ^[2] Take the small values $h=10$, $v=10$, $g=5$ and a horizontal shot, $\theta=0$: then $r=10t$ and $y=10-2.5t^2$, so after one second the fragment is $10$ metres out and $7.5$ metres up, and it lands at $t=2$, $20$ metres away. With the data of the problem and $\theta=45^\circ$ a fragment rises $10.1937$ metres above the burst, stays up $6.1814$ seconds and lands $87.4181$ metres away; the one thrown straight up stays up longest, $6.9929$ seconds. By energy conservation every fragment reaches the ground with the same speed $\sqrt{v^2+2gh}$, about $48.6004$ metres per second.

2. Which points can be reached

Fix a radius $r\gt0$ and ask which heights are reachable there. Eliminating $t$ with $t=r/(v\cos\theta)$ and using $1/\cos^2\theta=1+\tan^2\theta$ turns the flight into a parabola over $r$ in the slope $u=\tan\theta$:

$$y=h+ru-\frac{gr^2}{2v^2}\bigl(1+u^2\bigr).$$

With the small values above and $r=10$ this reads $y=17.5-2.5(u-2)^2$: the slopes $u=0,1,2,3$ give the heights $7.5$, $15$, $17.5$, $15$, so $u=2$ is the best slope there and $17.5$ metres the ceiling at that radius.

In general, since the fragments fly in every direction, $u$ runs over all real numbers: large positive $u$ is a steep upward shot, negative $u$ a downward one. For fixed $r$ the right-hand side is a quadratic in $u$ with negative leading coefficient $-gr^2/(2v^2)$, so it has a largest value, attained where its derivative $r-(gr^2/v^2)u$ vanishes, that is at $u=v^2/(gr)$, and it tends to $-\infty$ as $u\to\pm\infty$. Substituting that $u$ back and writing

$$a=h+\frac{v^2}{2g},\qquad b=\frac{g}{2v^2},$$

the largest reachable height at radius $r$ is simply

$$y_{\max}(r)=a-br^2.$$

Being continuous in $u$, attaining $y_{\max}(r)$ and falling to $-\infty$ on both sides, the quadratic takes every value $y\le y_{\max}(r)$ and none above it, so the flight curves pass through exactly those heights at this radius. Of them only $y\ge0$ is above ground, and Section 6 shows that a fragment is still in flight when it arrives there. Here $a\approx120.3874$ and $b=0.0122625$; at $r=40$, where $br^2=19.62$, the ceiling is about $100.7674$ metres, grazed by the slope $u=v^2/(gr)\approx1.0194$, an angle of about $45.5495^\circ$.

3. The same condition as a discriminant

Solving the quadratic for $u$ instead of maximising it gives the angles that hit a prescribed point $(r,y)$,

$$\begin{aligned}u&=\frac{v^2\pm\sqrt W}{gr},\\ W&=v^4-g\bigl(gr^2+2(y-h)v^2\bigr),\end{aligned}$$

the classical formula for the launch angle required to hit a target. ^[2] A real solution exists exactly when the discriminant $W$ is non-negative, and dividing $W\ge0$ by $2gv^2$ returns $y\le a-br^2$. The two signs are the familiar high and low trajectory: for $r=40$, $y=100$ they are $50.5934^\circ$ and $39.4066^\circ$.

4. The envelope and its focus

Let $D(r)$ be the envelope height $a-br^2$ minus the height of the flight of slope $u$. Completing the square in the formula of Section 2 gives the identity

$$D(r)=b\bigl(ru-v^2/g\bigr)^2,$$

so the curve lies above every flight path, and $D$ and its derivative vanish together exactly where $ru=v^2/g$: each upward path is tangent to the curve at the radius $r=v^2/(gu)$ of Section 2, while horizontal and downward paths stay strictly below it for $r\gt0$. Such a curve is the envelope of the family; the standard recipe for one — set the family equation and its derivative in the parameter to zero — is the maximisation of Section 2. ^[3] In ballistics it is the parabola of safety. The contact lies on the real flight only if it is above the ground: for $h=10$, $v=10$, $g=5$ and $u=\tfrac12$ it would sit at $r=40$, $20$ metres below ground, while that shot lands after $27.5959$ metres.

Read as $r^2=(a-y)/b$ and compared with the standard form of a parabola with vertical axis, the envelope has focal length $1/(4b)=v^2/(2g)$. Its vertex sits at height $a$, one focal length above the burst, so the focus is the burst point and the directrix is the horizontal line at height $h+v^2/g\approx140.7747$ metres. ^[4]

One corollary is needed below: the envelope meets the ground at

$$R=\sqrt{\frac ab}=\frac vg\sqrt{v^2+2gh},$$

about $99.0834$ metres. No fragment lands farther out, and the one whose path touches the envelope at ground level does land there; its angle $\arctan\bigl(v/\sqrt{v^2+2gh}\bigr)\approx22.3681^\circ$ is the optimal angle for a throw from a height, much flatter than $45^\circ$, as maximising the known landing distance from height $h$ ^[2] over $\theta$ confirms numerically.

5. Three dimensions

Nothing in this depended on the compass direction, because the initial speed is the same in every direction. Rotating about the axis, the set of points that some fragment visits at some time is described by the two horizontal coordinates $x$, $z$ and the height $y$ through

$$0\le y\le a-b\bigl(x^2+z^2\bigr).$$

This is a solid paraboloid dome of height $a\approx120.3874$ metres on a ground disc of radius $R\approx99.0834$ metres. The axis is covered too: the fragment thrown straight up rises to $a$ and falls back through everything below it.

6. Nothing is lost to the ground

One point still needs care, and it is the only subtle step: a fragment counts only until it reaches the ground. Could one pass underground and come back, so that part of the dome is reached only by flights already over? It cannot. Take a point off the axis, so $r\gt0$ and $0\le y\le y_{\max}(r)$. Section 2 gives it a finite slope $u$, hence an angle with $\cos\theta\gt0$ and an arrival time $t=r/(v\cos\theta)$; the height $y(s)=h+v\sin\theta\,s-\tfrac12gs^2$ then satisfies

$$\begin{aligned}y(s)&=\frac{(t-s)h+s\,y(t)}{t}\\&\quad+\frac g2\,s(t-s),\end{aligned}$$

because both sides are quadratics in $s$ with the same leading coefficient that agree at $s=0$ and at $s=t$. For $0\lt s\lt t$ the first term is non-negative, since $y(0)=h\ge0$ and $y(t)=y\ge0$, and the second is strictly positive: the fragment stays strictly above the ground and has certainly not landed before it arrives. Points on the axis are covered by the vertical shot of Section 5. The dome is therefore swept in full, and the region asked for is precisely that dome.

7. The volume of the dome

At height $y$ the dome cuts a disc of radius $r(y)$ with $r(y)^2=(a-y)/b$, so the disc method ^[5] gives

$$\begin{aligned}V&=\pi\int_0^a r(y)^2\,dy\\&=\frac\pi b\int_0^a(a-y)\,dy=\frac{\pi a^2}{2b}.\end{aligned}$$

The shell method, the form used in the solution's comments, is just as short:

$$\begin{aligned}V&=2\pi\int_0^R r\,(a-br^2)\,dr\\&=2\pi\Bigl[\frac{ar^2}2-\frac{br^4}4\Bigr]_0^R,\end{aligned}$$

and with $R^2=a/b$ this collapses to $\pi a^2/(2b)$ as well. In the data of the problem,

$$V=\frac{\pi v^2}{g}\Bigl(h+\frac{v^2}{2g}\Bigr)^2.$$

The whole problem is this one formula; the programs only evaluate it.

8. Four checks on the formula

Half a cylinder. The dome has height $a$ and base radius $R$, and $\pi a^2/(2b)=\pi R^2a/2$: it fills exactly half its circumscribing cylinder, the classical volume of a paraboloid of revolution. ^[6]
Pappus. In one vertical half-plane the region under the envelope has area $\int_0^R(a-br^2)\,dr=\tfrac23aR$ and first moment $\int_0^R r(a-br^2)\,dr=\tfrac14aR^2$ about the axis, which touches it without crossing, so its centroid lies $\tfrac38R$ out and Pappus's centroid theorem ^[5] gives $2\pi\cdot\tfrac38R\cdot\tfrac23aR=\tfrac12\pi R^2a$, the same value.
A burst on the ground. For $h=0$ the formula collapses to $\pi v^6/(4g^3)$, the check built into the C++ program; with $v=10$ and $g=5$ that is $2000\pi$.
A small worked example. For $h=10$, $v=10$, $g=5$ we get $a=20$, $b=0.025$, $R=\sqrt{800}\approx28.2843$ and $V=\pi\cdot400/0.05=8000\pi$, that is $25132.7412$.

Finally $v^2/g$ and $h+v^2/(2g)$ are lengths, so the formula is a length cubed, and it grows like $h^2$ for a high burst and like $v^6$ for a strong one.

9. Getting four decimals right

The value has $7$ digits before the decimal point, so four correct decimals mean $11$ correct significant digits. A double-precision number carries $15$ to $17$, with a unit roundoff of $2^{-53}\approx1.1\cdot10^{-16}$ per operation, ^[7] and the formula is a handful of multiplications and divisions of positive quantities, with no cancellation. Recomputing it in $60$-digit decimal arithmetic shows the double-precision value is off by about $3\cdot10^{-10}$, five orders of magnitude below the $10^{-4}$ the fourth decimal needs, and the rounding is not borderline: the exact value lies about $2.2\cdot10^{-5}$ from the nearest fifth-decimal half-way point, so no tie-breaking rule can change the printed digits. Finally $g=9.81$ must be used exactly as given; the standard value $9.80665$ would shift the volume by $4.6\cdot10^{-4}$ of its size.

How the Code Works

The C++ program ^[8] is built around the two constants of Section 2. A struct Options holds height, speed and gravity, preset to $100$, $20$ and $9.81$, and a flag run_checkpoints. parse_ld_after_prefix reads one option of the form --height=… with std::strtold and rejects a tail that is not consumed completely; parse_arguments loops over the arguments, also accepts --skip-checkpoints, and requires a non-negative height and a positive speed and gravity. Those are the right conditions for finite inputs, but an infinity slips through: --height=inf, or a literal that overflows, makes the program print an infinite volume.

solve(h, v, g) is Sections 2 and 7 verbatim: it forms a $=h+v^2/(2g)$, b $=g/(2v^2)$ and r_max2 $=a/b$, the squared ground radius of Section 4, and returns pi * a * a / (2 * b) with $\pi$ from std::acos(-1.0L); the shell integral stands in the comments above that line, and r_max2 is a by-product the formula never reads. Unless --skip-checkpoints was given, main next calls run_checkpoints, which compares solve(0, 10, 5) with $\pi v^6/(4g^3)$ and, if they differ by more than $10^{-12}$ in absolute value, makes main exit with status $2$. Only then is solve called with the requested parameters and the volume printed with std::fixed and std::setprecision(4).

The C++ arithmetic uses long double, whose width depends on the compiler and the platform: a $64$-bit significand with GCC and Clang on x86-64, the same as double on Apple silicon, quadruple precision on AArch64 Linux. The result is cast back to double for printing in every case, and Section 9 shows that double is already enough. Python keeps the parameters as the defaults of solve(h, v, g) and returns f"{ans:.4f}" computed with math.pi; Java hard-codes them in a parameterless solve() and formats with Math.PI, String.format and an explicit Locale.US, so the separator is a point and not a comma. Neither port has the command-line interface, its validation or the checkpoint.

Complexity and Verification

Evaluating the closed form takes a fixed number of arithmetic operations and a constant number of scalars, $O(1)$ time and $O(1)$ space; the only loop in any of the three programs is the C++ pass over its command-line arguments. All three print the same value immediately. The C++ checkpoint passes, the variants $(h,v,g)=(0,10,5)$ and $(10,10,5)$ print $6283.1853$ and $25132.7412$ as Section 8 predicts, and a negative height, a zero speed or a malformed option value makes it exit with status $1$.

The geometry was also tested independently of the formula, outside the published solutions. Flight paths for $200\,000$ launch angles, sampled every $5$ milliseconds, never rose above the envelope, the largest excess being $1.4\cdot10^{-14}$. For $200\,000$ random points of the bounding cylinder the discriminant test of Section 3 and the inequality $y\le a-br^2$ agreed in every case, and for $100\,000$ random points inside the dome both angles that reach the point were computed and the height along the way never became negative, the claim of Section 6. Both integrals of Section 7 by Simpson's rule, a sample of $2\,000\,000$ random points of the cylinder, and the Pappus computation all agree with the closed form, the sampling to one part in a thousand. A scan of the launch angle in steps of $0.001^\circ$ reproduces $R$ and picks $22.368^\circ$, the analytic optimum of Section 4 at that grid's resolution.

Footnotes and References

The references below supply the background facts used above. The envelope, the proof that the whole dome is swept before the fragments land, the volume formula and every number quoted here are derived or computed in this article.

Project Euler 317 — Firecracker. The official statement has a firecracker exploding 100 m above level ground into very small fragments moving in every direction with the same initial speed of 20 m/s, without air resistance and with g = 9.81 m/s², and asks for the volume in m³ of the region they move through before reaching the ground, to four decimal places.
Wikipedia — Projectile motion. The article derives x(t) = v₀t cos θ and y(t) = v₀t sin θ − gt²/2 and eliminates t to get the parabolic path y = x tan θ − gx²/(2v₀²cos²θ) of Section 2. For the angle required to hit a point (x, y) it states θ = arctan((v² ± √(v⁴ − g(gx² + 2yv²)))/(gx)), the two roots being the launch angles when they are not imaginary, the condition of Section 3; and for a launch from a height y₀ above the impact point it gives d = (v cos θ/g)(v sin θ + √((v sin θ)² + 2gy₀)), maximised over θ in Section 4. The title "Range of a projectile" redirects here.
Wikipedia — Envelope (mathematics). An envelope of a family F(t, x, y) = 0 is a curve tangent to every member; it follows from F = 0 with ∂F/∂t = 0, and when F is polynomial in t, equivalently from the vanishing of the discriminant in t. Example 5 throws projectiles with constant initial speed at different angles and obtains the concave parabola y = v²/(2g) − gx²/(2v²), the envelope of Sections 2 and 4.
Wikipedia — Parabola. A parabola is the set of points whose distance to a fixed point, the focus, equals the distance to a fixed line, the directrix; with the vertex at the origin, y = x²/(4f) has focal length f, focus (0, f) and directrix y = −f. Section 4 applies this to the envelope.
Wikipedia — Solid of revolution. Revolving a plane region about the y-axis gives π∫|f(y)² − g(y)²|dy by the disc method and 2π∫x|f(x) − g(x)|dx by the shell method; by Pappus's second centroid theorem the volume is the region's area times the circumference traced by its centroid, if the region does not cross the axis. Sections 7 and 8 use all three.
Wikipedia — Paraboloid. A paraboloidal dish of rim radius R and depth D holds πR²D/2, half the πR²D of the cylinder around it: the check of Section 8 with D = a.
Wikipedia — Double-precision floating-point format. Binary64 has a 53-bit significand, which gives 15 to 17 significant decimal digits; its round-to-nearest unit roundoff is 2⁻⁵³ ≈ 1.11·10⁻¹⁶, the precision budget of Section 9.
C++ — Euler317.cpp. The linked immutable C++ revision contains Options, parse_ld_after_prefix, parse_arguments, solve and run_checkpoints. It is the implementation record for this article; the sampling experiments above are extra checks, not tests in that file.

Problem 317 source code

C++

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

namespace {

struct Options {
    long double height = 100.0L;
    long double speed = 20.0L;
    long double gravity = 9.81L;
    bool run_checkpoints = 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;
    }
    char* end_ptr = nullptr;
    value = std::strtold(tail.c_str(), &end_ptr);
    return end_ptr != nullptr && *end_ptr == '\0';
}

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, "--height=", options.height) ||
            parse_ld_after_prefix(arg, "--speed=", options.speed) ||
            parse_ld_after_prefix(arg, "--gravity=", options.gravity)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.height >= 0.0L && options.speed > 0.0L && options.gravity > 0.0L;
}

long double solve(const long double h, const long double v, const long double g) {
    // Envelope in cylindrical coordinates:
    // y(r) = h + v^2/(2g) - g*r^2/(2v^2), for 0 <= r <= r_max where y=0.
    const long double a = h + (v * v) / (2.0L * g);
    const long double b = g / (2.0L * v * v);
    const long double r_max2 = a / b;

    // Volume of revolution: V = pi * integral_0^{r_max} r^2 * 2? No:
    // V = pi * integral_0^{y_max} r(y)^2 dy = pi * integral_0^{r_max} 2r y(r) dr.
    // Directly with y(r)=a-br^2:
    // V = pi * integral_0^{r_max} (a r^2 - b r^4)'? easiest via shells:
    // V = 2*pi * integral_0^{r_max} r*(a-br^2) dr = pi * a^2 / (2b).
    const long double pi = std::acos(-1.0L);
    return pi * a * a / (2.0L * b);
}

bool run_checkpoints() {
    // For h=0, the envelope is centered at the launch point and the volume becomes:
    // pi * (v^6) / (4 g^3).
    {
        const long double v = 10.0L;
        const long double g = 5.0L;
        const long double expected =
            std::acos(-1.0L) * v * v * v * v * v * v / (4.0L * g * g * g);
        const long double got = solve(0.0L, v, g);
        if (std::fabsl(got - expected) > 1e-12L) {
            std::cerr << "Checkpoint failed for h=0 closed form" << '\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 ans = solve(options.height, options.speed, options.gravity);
    std::cout << std::fixed << std::setprecision(4) << static_cast<double>(ans) << '\n';
    return 0;
}

Python

import math

def solve(h=100.0, v=20.0, g=9.81):
    a = h + (v * v) / (2.0 * g)
    b = g / (2.0 * v * v)
    ans = math.pi * a * a / (2.0 * b)
    return f"{ans:.4f}"

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

Java

import java.util.Locale;

public class Euler317 {
    public static String solve() {
        double h = 100.0;
        double v = 20.0;
        double g = 9.81;

        double a = h + (v * v) / (2.0 * g);
        double b = g / (2.0 * v * v);

        double ans = Math.PI * a * a / (2.0 * b);
        return String.format(Locale.US, "%.4f", ans);
    }

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