Problem 363: Bézier Curves
View on Project EulerProject Euler Problem 363 Solution
EulerSolve provides an optimized solution for Project Euler Problem 363, Bézier Curves, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The curve is the cubic Bézier arc with control points $P_0=(1,0),\qquad P_1=(1,v),\qquad P_2=(v,1),\qquad P_3=(0,1).$ It starts at \((1,0)\), ends at \((0,1)\), and is meant to mimic a quarter of the unit circle. The parameter \(v\) is not arbitrary: the region bounded by the coordinate axes and the Bézier arc must have the same area as a quarter disk, namely \(\pi/4\). Once that area condition determines \(v\), we compute the Bézier arc length \(L\) and compare it with the exact quarter-circle length \(\pi/2\)....
Detailed mathematical approach
Problem Summary
The curve is the cubic Bézier arc with control points
$P_0=(1,0),\qquad P_1=(1,v),\qquad P_2=(v,1),\qquad P_3=(0,1).$
It starts at \((1,0)\), ends at \((0,1)\), and is meant to mimic a quarter of the unit circle. The parameter \(v\) is not arbitrary: the region bounded by the coordinate axes and the Bézier arc must have the same area as a quarter disk, namely \(\pi/4\). Once that area condition determines \(v\), we compute the Bézier arc length \(L\) and compare it with the exact quarter-circle length \(\pi/2\).
The quantity printed by the program is
$100\cdot \frac{L-\pi/2}{\pi/2}.$
Mathematical Approach
Step 1: Expand the Cubic Bézier Curve
A cubic Bézier curve has the Bernstein form
$B(t)=(1-t)^3P_0+3(1-t)^2tP_1+3(1-t)t^2P_2+t^3P_3,\qquad 0\le t\le 1.$
Substituting the four control points and separating coordinates gives
$x(t)=(1-t)^3+3(1-t)^2t+3v(1-t)t^2,$
$y(t)=3v(1-t)^2t+3(1-t)t^2+t^3.$
After expansion this becomes exactly the polynomial form used in the solution files:
$x(t)=1+3(v-1)t^2+(2-3v)t^3,$
$y(t)=3vt+(3-6v)t^2+(3v-2)t^3.$
Differentiating gives the velocity components
$x'(t)=6(v-1)t+3(2-3v)t^2,$
$y'(t)=3v+2(3-6v)t+3(3v-2)t^2.$
The implementations package these derivatives inside the speed function
$s(t)=\sqrt{(x'(t))^2+(y'(t))^2}.$
Step 2: Convert the Area Condition into an Equation for \(v\)
The enclosed region is bounded by the \(x\)-axis from \((0,0)\) to \((1,0)\), the Bézier arc from \((1,0)\) to \((0,1)\), and the \(y\)-axis back to the origin. Green's theorem allows us to write the area as the line integral
$A=\oint_C x\,dy.$
The two axis segments contribute nothing: on the \(x\)-axis we have \(dy=0\), and on the \(y\)-axis we have \(x=0\). Therefore only the Bézier arc remains:
$A(v)=\int_0^1 x(t)\,y'(t)\,dt.$
Multiplying the two polynomials gives
$\begin{aligned} x(t)y'(t)=&\,3v+(6-12v)t+(9v^2-6)t^2+(-45v^2+60v-18)t^3\\ &+(63v^2-87v+30)t^4+(-27v^2+36v-12)t^5. \end{aligned}$
Integrating term by term over \([0,1]\) simplifies to
$A(v)=\frac12+\frac{3v}{5}-\frac{3v^2}{20}.$
This is the exact closed form used by curve_area(v) in the C++ source.
Step 3: Solve the Quadratic Exactly
The problem requires the area to equal the area of a quarter unit disk:
$A(v)=\frac{\pi}{4}.$
Substituting the formula for \(A(v)\) gives
$\frac12+\frac{3v}{5}-\frac{3v^2}{20}=\frac{\pi}{4}.$
Multiplying by \(20\) and rearranging yields
$10+12v-3v^2=5\pi,$
$3v^2-12v+(5\pi-10)=0.$
Applying the quadratic formula gives
$v=\frac{12\pm\sqrt{144-12(5\pi-10)}}{6}=2\pm \frac{\sqrt{66-15\pi}}{3}.$
The plus sign gives a value greater than \(1\), so the interior control points would leave the unit square. The geometric branch used by every implementation is therefore
$\boxed{v=2-\frac{\sqrt{66-15\pi}}{3}}\approx 0.551778477804468.$
Step 4: Set Up the Arc-Length Integral
For a parametric plane curve, the arc length is
$L(v)=\int_0^1 \sqrt{(x'(t))^2+(y'(t))^2}\,dt=\int_0^1 s(t)\,dt.$
Here the integrand is a square root of a quartic polynomial in \(t\). The repository solutions do not try to derive an elementary antiderivative; instead they evaluate the integral numerically to high precision.
The final percentage is then
$\text{percent}=100\cdot \frac{L(v)-\pi/2}{\pi/2}.$
Step 5: Adaptive Simpson Integration
The numerical part is identical in all three languages. For any interval \([a,b]\), Simpson's rule approximates the integral of \(s(t)\) by
$S(a,b)=\frac{b-a}{6}\left(s(a)+4s\!\left(\frac{a+b}{2}\right)+s(b)\right).$
Adaptive Simpson then splits the interval at the midpoint \(m\), computes \(S(a,m)\) and \(S(m,b)\), and compares
$S(a,m)+S(m,b)$
with the old estimate \(S(a,b)\). If
$\left|S(a,m)+S(m,b)-S(a,b)\right|\le 15\varepsilon,$
the code accepts the corrected estimate
$S(a,m)+S(m,b)+\frac{S(a,m)+S(m,b)-S(a,b)}{15}.$
Otherwise it recurses on both halves with tolerance \(\varepsilon/2\). The C++ version uses long double, tolerance \(10^{-18}\), and maximum depth \(30\). The Python and Java versions use the same recursion pattern with tolerance \(10^{-15}\) and depth \(30\).
Step 6: Numerical Value and Sanity Checks
Using the exact value of \(v\) above, the numerical integration gives
$L\approx 1.570796911273925.$
Hence
$100\cdot \frac{L-\pi/2}{\pi/2}\approx 0.0000372091.$
This matches the formatted output of the local solution programs.
The C++ implementation also contains useful checkpoints. At \(v=0\), the area formula gives
$A(0)=\frac12,$
and the curve degenerates to the diagonal segment \(x+y=1\), whose length is
$\sqrt{2}.$
The code verifies both facts numerically before computing the final answer, and it also checks that the chosen \(v\) lies in the interval \((0,1)\) and reproduces the target area \(\pi/4\).
How the Code Works
The three implementation files share the same structure. solve_v() uses the closed-form quadratic solution, so the parameter search is not numerical. curve_speed(t, v) evaluates the norm of the derivative vector. curve_length(v) starts Simpson's rule on \([0,1]\), and adaptive_simpson(...) refines the partition until the error estimate is below the requested tolerance or the recursion depth reaches zero.
The C++ file adds curve_area(v) and run_checkpoints() explicitly, which makes the mathematical structure especially transparent: first validate the exact area formula, then integrate the speed, then form the relative error percentage.
Complexity Analysis
Deriving \(v\) from the quadratic equation is \(O(1)\) time and \(O(1)\) memory. The arc-length computation dominates the runtime. If the adaptive Simpson routine ends up using \(m\) accepted subintervals, then the time cost is \(O(m)\) function evaluations up to a constant factor, and the memory usage is \(O(d)\), where \(d\) is the recursion depth. In this repository, \(d\le 30\) by construction, so the memory footprint is effectively constant.
Footnotes and References
- Problem page: https://projecteuler.net/problem=363
- Bézier curve: Wikipedia - Bézier curve
- Green's theorem: Wikipedia - Green's theorem
- Simpson's rule: Wikipedia - Simpson's rule