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