Problem 449: Chocolate Covered Candy

View on Project Euler

Project Euler Problem 449 Solution

EulerSolve provides an optimized solution for Project Euler Problem 449, Chocolate Covered Candy, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The candy core is the ellipsoid of revolution $$\frac{x^2+y^2}{a^2}+\frac{z^2}{b^2}=1,$$ equivalently \(b^2x^2+b^2y^2+a^2z^2=a^2b^2\). A chocolate layer of thickness \(1\) mm is added along the outward surface normal , so the outer boundary is the parallel surface at distance \(1\). The required quantity \(C(a,b)\) is the volume of that outer body minus the volume of the original ellipsoid. The important geometric point is that a normal offset is not obtained by simply replacing \((a,b)\) with \((a+1,b+1)\). We must derive the offset surface explicitly. Mathematical Approach Step 1: Parametrize the meridian ellipse Because the body is rotationally symmetric around the \(z\)-axis, it is enough to work in the meridian plane \((\rho,z)\), where \(\rho=\sqrt{x^2+y^2}\). The inner ellipse is $$\frac{\rho^2}{a^2}+\frac{z^2}{b^2}=1.$$ A convenient parametrization of its upper half is $$\rho(\theta)=a\sin\theta,\qquad z(\theta)=b\cos\theta,\qquad 0\le\theta\le\frac{\pi}{2}.$$ Rotating this arc around the \(z\)-axis recovers the whole ellipsoid....

Detailed mathematical approach

Problem Summary

The candy core is the ellipsoid of revolution

$$\frac{x^2+y^2}{a^2}+\frac{z^2}{b^2}=1,$$

equivalently \(b^2x^2+b^2y^2+a^2z^2=a^2b^2\). A chocolate layer of thickness \(1\) mm is added along the outward surface normal, so the outer boundary is the parallel surface at distance \(1\). The required quantity \(C(a,b)\) is the volume of that outer body minus the volume of the original ellipsoid.

The important geometric point is that a normal offset is not obtained by simply replacing \((a,b)\) with \((a+1,b+1)\). We must derive the offset surface explicitly.

Mathematical Approach

Step 1: Parametrize the meridian ellipse

Because the body is rotationally symmetric around the \(z\)-axis, it is enough to work in the meridian plane \((\rho,z)\), where \(\rho=\sqrt{x^2+y^2}\). The inner ellipse is

$$\frac{\rho^2}{a^2}+\frac{z^2}{b^2}=1.$$

A convenient parametrization of its upper half is

$$\rho(\theta)=a\sin\theta,\qquad z(\theta)=b\cos\theta,\qquad 0\le\theta\le\frac{\pi}{2}.$$

Rotating this arc around the \(z\)-axis recovers the whole ellipsoid.

Step 2: Compute the outward unit normal

Write the ellipse implicitly as

$$F(\rho,z)=\frac{\rho^2}{a^2}+\frac{z^2}{b^2}-1=0.$$

An outward normal is proportional to \(\nabla F\), so on the parametrized curve we get the normal direction

$$\left(\frac{\rho}{a^2},\frac{z}{b^2}\right)=\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$

Define

$$E(\theta)=\sqrt{\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2}}.$$

Then the outward unit normal in the meridian plane is

$$\mathbf{n}(\theta)=\frac{1}{E(\theta)}\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$

Step 3: Offset the profile by distance \(1\)

Moving one unit along the outward normal sends the meridian point \((\rho(\theta),z(\theta))\) to the outer profile \((R(\theta),Z(\theta))\):

$$R(\theta)=a\sin\theta+\frac{\sin\theta}{aE(\theta)}=\sin\theta\left(a+\frac{1}{aE(\theta)}\right),$$

$$Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}=\cos\theta\left(b+\frac{1}{bE(\theta)}\right).$$

This is exactly the geometry evaluated by the implementations.

Step 4: Differentiate the outer height

For the volume integral we need \(dZ\). First differentiate \(E(\theta)\):

$$E'(\theta)=\frac{\sin\theta\cos\theta}{E(\theta)}\left(\frac{1}{a^2}-\frac{1}{b^2}\right).$$

Now differentiate \(Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}\). After simplification,

$$-Z'(\theta)=\sin\theta\left(b+\frac{1}{bE(\theta)}+\frac{\cos^2\theta}{b}\frac{\frac{1}{a^2}-\frac{1}{b^2}}{E(\theta)^3}\right).$$

The integrand used numerically is therefore \(R(\theta)^2(-Z'(\theta))\).

Step 5: Volume by disks of revolution

The outer surface is symmetric with respect to the plane \(z=0\), so we integrate the upper half and double. Using the disk method for a solid of revolution,

$$V_{\text{outer}}=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta.$$

The inner ellipsoid volume is the standard formula

$$V_{\text{inner}}=\frac{4\pi a^2 b}{3}.$$

Hence the chocolate volume is

$$\boxed{C(a,b)=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta-\frac{4\pi a^2 b}{3}.}$$

Step 6: Numerical checks

When \(a=b=1\), the core is the unit sphere. Then \(E(\theta)=1\), so

$$R(\theta)=2\sin\theta,\qquad Z(\theta)=2\cos\theta,$$

which means the outer body is a sphere of radius \(2\). Therefore

$$C(1,1)=\frac{4\pi}{3}(2^3-1^3)=\frac{28\pi}{3}.$$

For the second checkpoint, numerical evaluation gives

$$C(2,1)\approx 60.35475635.$$

Applying the same integral to the target parameters yields

$$C(3,1)\approx 103.37870096,$$

which is the value rounded to eight decimal places by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. They evaluate the smooth meridian integrand derived above, apply Simpson's rule on \([0,\pi/2]\), and then refine the interval recursively with adaptive Simpson integration until the local error estimate is below the requested tolerance or the recursion depth limit is reached.

Once the integral is stable, the implementation multiplies by \(2\pi\), subtracts the ellipsoid volume \(\frac{4\pi a^2b}{3}\), and formats the result to eight digits after the decimal point. The two published checkpoints are used as numerical sanity checks for the geometry and the quadrature.

Complexity Analysis

Adaptive Simpson integration is data-dependent. If the interval is split into \(m\) accepted subintervals, the total work is \(O(m)\) function evaluations up to a constant factor, and the extra memory is \(O(d)\), where \(d\) is the recursion depth. For this problem the integrand is smooth on \([0,\pi/2]\), so convergence is fast and memory usage is negligible.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=449
  2. Ellipsoid: Wikipedia - Ellipsoid
  3. Parallel curve and offset geometry: Wikipedia - Parallel curve
  4. Simpson's rule: Wikipedia - Simpson's rule
  5. Adaptive Simpson's method: Wikipedia - Adaptive Simpson's method

Problem 449 source code

C++

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

namespace {

struct Options {
    long double a = 3.0L;
    long double b = 1.0L;
    bool run_checkpoints = true;
};

bool parse_ld_after_prefix(const std::string& arg, const std::string& prefix, long double& out) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        out = 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, "--a=", options.a)) {
            continue;
        }
        if (parse_ld_after_prefix(arg, "--b=", options.b)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.a > 0.0L && options.b > 0.0L;
}

long double integrand(long double theta, long double a, long double b) {
    const long double s = std::sin(theta);
    const long double c = std::cos(theta);
    const long double aa = a * a;
    const long double bb = b * b;
    const long double e = std::sqrt((s * s) / aa + (c * c) / bb);

    const long double r = s * (a + 1.0L / (a * e));
    const long double minus_z_prime =
        s * (b + 1.0L / (b * e) + (c * c / b) * (1.0L / aa - 1.0L / bb) / (e * e * e));
    return r * r * minus_z_prime;
}

long double simpson(long double l, long double r, long double a, long double b) {
    const long double m = (l + r) * 0.5L;
    return (r - l) *
           (integrand(l, a, b) + 4.0L * integrand(m, a, b) + integrand(r, a, b)) / 6.0L;
}

long double adaptive_simpson(long double l,
                             long double r,
                             long double eps,
                             long double whole,
                             int depth,
                             long double a,
                             long double b) {
    const long double m = (l + r) * 0.5L;
    const long double left = simpson(l, m, a, b);
    const long double right = simpson(m, r, a, b);
    const long double delta = left + right - whole;

    if (depth <= 0 || std::fabsl(delta) <= 15.0L * eps) {
        return left + right + delta / 15.0L;
    }
    return adaptive_simpson(l, m, eps * 0.5L, left, depth - 1, a, b) +
           adaptive_simpson(m, r, eps * 0.5L, right, depth - 1, a, b);
}

long double compute_chocolate_volume(long double a, long double b) {
    constexpr long double pi = 3.141592653589793238462643383279502884L;
    const long double half_pi = pi * 0.5L;

    const long double whole = simpson(0.0L, half_pi, a, b);
    const long double integral = adaptive_simpson(0.0L, half_pi, 1e-14L, whole, 24, a, b);

    const long double outer = 2.0L * pi * integral;
    const long double inner = 4.0L * pi * a * a * b / 3.0L;
    return outer - inner;
}

bool run_checkpoints() {
    constexpr long double pi = 3.141592653589793238462643383279502884L;
    const long double c11 = compute_chocolate_volume(1.0L, 1.0L);
    const long double expected11 = 28.0L * pi / 3.0L;
    if (std::fabsl(c11 - expected11) > 1e-12L) {
        std::cerr << "Checkpoint failed: (a,b)=(1,1)\n";
        return false;
    }

    const long double c21 = compute_chocolate_volume(2.0L, 1.0L);
    if (std::fabsl(c21 - 60.35475635L) > 5e-9L) {
        std::cerr << "Checkpoint failed: (a,b)=(2,1)\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 answer = compute_chocolate_volume(options.a, options.b);
    std::cout << std::fixed << std::setprecision(8) << answer << '\n';
    return 0;
}

Python

import math

def integrand(theta, a, b):
    s = math.sin(theta)
    c = math.cos(theta)
    aa = a * a
    bb = b * b
    e = math.sqrt((s * s) / aa + (c * c) / bb)
    
    r = s * (a + 1.0 / (a * e))
    minus_z_prime = s * (b + 1.0 / (b * e) + (c * c / b) * (1.0 / aa - 1.0 / bb) / (e * e * e))
    return r * r * minus_z_prime

def simpson(l, r, a, b):
    m = (l + r) * 0.5
    return (r - l) * (integrand(l, a, b) + 4.0 * integrand(m, a, b) + integrand(r, a, b)) / 6.0

def adaptive_simpson(l, r, eps, whole, depth, a, b):
    m = (l + r) * 0.5
    left = simpson(l, m, a, b)
    right = simpson(m, r, a, b)
    delta = left + right - whole
    
    if depth <= 0 or abs(delta) <= 15.0 * eps:
        return left + right + delta / 15.0
    return adaptive_simpson(l, m, eps * 0.5, left, depth - 1, a, b) + \
           adaptive_simpson(m, r, eps * 0.5, right, depth - 1, a, b)

def compute_chocolate_volume(a, b):
    pi = 3.141592653589793238462643383279502884
    half_pi = pi * 0.5
    whole = simpson(0.0, half_pi, a, b)
    integral = adaptive_simpson(0.0, half_pi, 1e-14, whole, 24, a, b)
    
    outer = 2.0 * pi * integral
    inner = 4.0 * pi * a * a * b / 3.0
    return outer - inner

def solve():
    ans = compute_chocolate_volume(3.0, 1.0)
    return f"{ans:.8f}"

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

Java

public class Euler449 {
    private static double integrand(double theta, double a, double b) {
        double s = Math.sin(theta);
        double c = Math.cos(theta);
        double aa = a * a;
        double bb = b * b;
        double e = Math.sqrt((s * s) / aa + (c * c) / bb);

        double r = s * (a + 1.0 / (a * e));
        double minusZPrime = s * (b + 1.0 / (b * e) + (c * c / b) * (1.0 / aa - 1.0 / bb) / (e * e * e));
        return r * r * minusZPrime;
    }

    private static double simpson(double l, double r, double a, double b) {
        double m = (l + r) * 0.5;
        return (r - l) * (integrand(l, a, b) + 4.0 * integrand(m, a, b) + integrand(r, a, b)) / 6.0;
    }

    private static double adaptiveSimpson(double l, double r, double eps, double whole, int depth, double a, double b) {
        double m = (l + r) * 0.5;
        double left = simpson(l, m, a, b);
        double right = simpson(m, r, a, b);
        double delta = left + right - whole;

        if (depth <= 0 || Math.abs(delta) <= 15.0 * eps) {
            return left + right + delta / 15.0;
        }
        return adaptiveSimpson(l, m, eps * 0.5, left, depth - 1, a, b) +
                adaptiveSimpson(m, r, eps * 0.5, right, depth - 1, a, b);
    }

    private static double computeChocolateVolume(double a, double b) {
        double pi = 3.141592653589793238462643383279502884;
        double halfPi = pi * 0.5;

        double whole = simpson(0.0, halfPi, a, b);
        double integral = adaptiveSimpson(0.0, halfPi, 1e-14, whole, 24, a, b);

        double outer = 2.0 * pi * integral;
        double inner = 4.0 * pi * a * a * b / 3.0;
        return outer - inner;
    }

    public static String solve() {
        double ans = computeChocolateVolume(3.0, 1.0);
        return String.format(java.util.Locale.US, "%.8f", ans);
    }

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