Problem 363: Bézier Curves

Project 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

Problem 363 source code

C++

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

namespace {

long double curve_area(const long double v) {
    // Integral of x(t) * y'(t) over t in [0,1].
    return 0.5L + (3.0L * v) / 5.0L - (3.0L * v * v) / 20.0L;
}

long double curve_speed(const long double t, const long double v) {
    const long double dx = 6.0L * (v - 1.0L) * t + 3.0L * (2.0L - 3.0L * v) * t * t;
    const long double dy = 3.0L * v + 2.0L * (3.0L - 6.0L * v) * t + 3.0L * (3.0L * v - 2.0L) * t * t;
    return std::sqrt(dx * dx + dy * dy);
}

long double adaptive_simpson(const long double v,
                             const long double a,
                             const long double b,
                             const long double eps,
                             const long double whole,
                             const long double fa,
                             const long double fb,
                             const long double fm,
                             const int depth) {
    const long double m = (a + b) * 0.5L;
    const long double lm = (a + m) * 0.5L;
    const long double rm = (m + b) * 0.5L;

    const long double flm = curve_speed(lm, v);
    const long double frm = curve_speed(rm, v);

    const long double left = (m - a) * (fa + 4.0L * flm + fm) / 6.0L;
    const long double right = (b - m) * (fm + 4.0L * frm + fb) / 6.0L;
    const long double combined = left + right;

    if (depth <= 0 || std::fabsl(combined - whole) <= 15.0L * eps) {
        return combined + (combined - whole) / 15.0L;
    }

    return adaptive_simpson(v, a, m, eps * 0.5L, left, fa, fm, flm, depth - 1) +
           adaptive_simpson(v, m, b, eps * 0.5L, right, fm, fb, frm, depth - 1);
}

long double curve_length(const long double v) {
    const long double a = 0.0L;
    const long double b = 1.0L;
    const long double fa = curve_speed(a, v);
    const long double fb = curve_speed(b, v);
    const long double m = 0.5L;
    const long double fm = curve_speed(m, v);
    const long double whole = (b - a) * (fa + 4.0L * fm + fb) / 6.0L;
    return adaptive_simpson(v, a, b, 1e-18L, whole, fa, fb, fm, 30);
}

long double solve_v() {
    // From area equation:
    //   1/2 + 3v/5 - 3v^2/20 = pi/4
    // => v = 2 +/- sqrt(66 - 15*pi) / 3
    const long double disc = 66.0L - 15.0L * std::acos(-1.0L);
    return 2.0L - std::sqrt(disc) / 3.0L;
}

bool run_checkpoints() {
    if (std::fabsl(curve_area(0.0L) - 0.5L) > 1e-18L) {
        std::cerr << "Checkpoint failed: area(v=0)\n";
        return false;
    }

    if (std::fabsl(curve_length(0.0L) - std::sqrt(2.0L)) > 1e-12L) {
        std::cerr << "Checkpoint failed: length(v=0)\n";
        return false;
    }

    const long double v = solve_v();
    if (!(v > 0.0L && v < 1.0L)) {
        std::cerr << "Checkpoint failed: v range\n";
        return false;
    }

    const long double target_area = std::acos(-1.0L) / 4.0L;
    if (std::fabsl(curve_area(v) - target_area) > 1e-18L) {
        std::cerr << "Checkpoint failed: area match\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const long double pi = std::acos(-1.0L);
    const long double v = solve_v();
    const long double L = curve_length(v);
    const long double percent = 100.0L * (L - pi / 2.0L) / (pi / 2.0L);

    std::cout << std::fixed << std::setprecision(10) << percent << '\n';
    return 0;
}

Python

import math

def curve_speed(t, v):
    dx = 6.0 * (v - 1.0) * t + 3.0 * (2.0 - 3.0 * v) * t * t
    dy = 3.0 * v + 2.0 * (3.0 - 6.0 * v) * t + 3.0 * (3.0 * v - 2.0) * t * t
    return math.sqrt(dx * dx + dy * dy)

def adaptive_simpson(v, a, b, eps, whole, fa, fb, fm, depth):
    m = (a + b) * 0.5
    lm = (a + m) * 0.5
    rm = (m + b) * 0.5

    flm = curve_speed(lm, v)
    frm = curve_speed(rm, v)

    left = (m - a) * (fa + 4.0 * flm + fm) / 6.0
    right = (b - m) * (fm + 4.0 * frm + fb) / 6.0
    combined = left + right

    if depth <= 0 or abs(combined - whole) <= 15.0 * eps:
        return combined + (combined - whole) / 15.0

    return adaptive_simpson(v, a, m, eps * 0.5, left, fa, fm, flm, depth - 1) + \
           adaptive_simpson(v, m, b, eps * 0.5, right, fm, fb, frm, depth - 1)

def curve_length(v):
    a = 0.0
    b = 1.0
    fa = curve_speed(a, v)
    fb = curve_speed(b, v)
    m = 0.5
    fm = curve_speed(m, v)
    whole = (b - a) * (fa + 4.0 * fm + fb) / 6.0
    return adaptive_simpson(v, a, b, 1e-15, whole, fa, fb, fm, 30)

def solve_v():
    disc = 66.0 - 15.0 * math.acos(-1.0)
    return 2.0 - math.sqrt(disc) / 3.0

def solve():
    pi = math.acos(-1.0)
    v = solve_v()
    L = curve_length(v)
    percent = 100.0 * (L - pi / 2.0) / (pi / 2.0)
    return f"{percent:.10f}"

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

Java

public class Euler363 {
    static double curveSpeed(double t, double v) {
        double dx = 6.0 * (v - 1.0) * t + 3.0 * (2.0 - 3.0 * v) * t * t;
        double dy = 3.0 * v + 2.0 * (3.0 - 6.0 * v) * t + 3.0 * (3.0 * v - 2.0) * t * t;
        return Math.sqrt(dx * dx + dy * dy);
    }

    static double adaptiveSimpson(double v, double a, double b, double eps, double whole, double fa, double fb,
            double fm, int depth) {
        double m = (a + b) * 0.5;
        double lm = (a + m) * 0.5;
        double rm = (m + b) * 0.5;

        double flm = curveSpeed(lm, v);
        double frm = curveSpeed(rm, v);

        double left = (m - a) * (fa + 4.0 * flm + fm) / 6.0;
        double right = (b - m) * (fm + 4.0 * frm + fb) / 6.0;
        double combined = left + right;

        if (depth <= 0 || Math.abs(combined - whole) <= 15.0 * eps) {
            return combined + (combined - whole) / 15.0;
        }

        return adaptiveSimpson(v, a, m, eps * 0.5, left, fa, fm, flm, depth - 1) +
                adaptiveSimpson(v, m, b, eps * 0.5, right, fm, fb, frm, depth - 1);
    }

    static double curveLength(double v) {
        double a = 0.0;
        double b = 1.0;
        double fa = curveSpeed(a, v);
        double fb = curveSpeed(b, v);
        double m = 0.5;
        double fm = curveSpeed(m, v);
        double whole = (b - a) * (fa + 4.0 * fm + fb) / 6.0;
        return adaptiveSimpson(v, a, b, 1e-15, whole, fa, fb, fm, 30);
    }

    static double solveV() {
        double disc = 66.0 - 15.0 * Math.acos(-1.0);
        return 2.0 - Math.sqrt(disc) / 3.0;
    }

    public static String solve() {
        double pi = Math.acos(-1.0);
        double v = solveV();
        double L = curveLength(v);
        double percent = 100.0 * (L - pi / 2.0) / (pi / 2.0);
        return String.format(java.util.Locale.US, "%.10f", percent);
    }

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