Project Euler

Bezier Curves

A cubic Bezier curve is defined by four control points P_0, P_1, P_2, P_3. Find the cubic Bezier curve that best approximates a quarter circle (from (0,1) to (1,0) along the unit circle) in the min...

Source sync Apr 19, 2026

Problem #0363

Level Level 13

Solved By 1,304

Languages C++, Python

Answer 0.0000372091

Length 403 words

optimizationgeometryalgebra

A cubic Bézier curve is defined by four points: and .

The curve is constructed as follows:

On the segments , , and the points and are drawn such that , with in .

On the segments and the points and are drawn such that
for the same value of .

On the segment the point is drawn such that for the same value of .

The Bézier curve defined by the points is the locus of as takes all possible positions on the segment .
(Please note that for all points the value of is the same.)

From the construction it is clear that the Bézier curve will be tangent to the segments in and in .

A cubic Bézier curve with and is used to approximate a quarter circle.
The value is chosen such that the area enclosed by the lines and the curve is equal to (the area of the quarter circle).

By how many percent does the length of the curve differ from the length of the quarter circle?
That is, if is the length of the curve, calculate
Give your answer rounded to 10 digits behind the decimal point.

Problem 363: Bezier Curves

Mathematical Foundation

Theorem 1 (Symmetric Optimal Bezier). By the symmetry of the quarter circle under reflection through the line $y = x$ , the optimal cubic Bezier approximation has control points of the form

P_0 = (0, 1), \quad P_1 = (k, 1), \quad P_2 = (1, k), \quad P_3 = (1, 0)

for a unique $k > 0$ .

Proof. The quarter circle from $(0,1)$ to $(1,0)$ is symmetric under the map $(x,y) \mapsto (y,x)$ . Any optimal approximant must share this symmetry (otherwise, reflecting it would give another approximant with the same error, and their average would have strictly smaller error by strict convexity of the max-norm on the error space). The symmetry constraint $P_1 = (k, 1)$ , $P_2 = (1, k)$ follows from requiring $B(t)$ and its reflection $\tilde{B}(t) = (B_y(1-t), B_x(1-t))$ to coincide. $\square$

Lemma 1 (Bezier Parametrization). With the above control points, the Bezier curve is

\mathbf{B}(t) = \bigl(x(t),\, y(t)\bigr)

where

x(t) = 3k t(1-t)^2 + 3t^2(1-t) + t^3, \quad y(t) = (1-t)^3 + 3(1-t)^2 t + 3k(1-t)t^2.

Proof. Direct substitution of $P_0, P_1, P_2, P_3$ into the Bernstein polynomial formula $\mathbf{B}(t) = \sum_{i=0}^{3} \binom{3}{i}(1-t)^{3-i}t^i P_i$ . $\square$

Definition (Radial Error). The radial error at parameter $t$ is

\varepsilon(t; k) = \sqrt{x(t)^2 + y(t)^2} - 1.

Theorem 2 (Chebyshev Equioscillation). The optimal $k$ is characterized by the Chebyshev equioscillation property: the radial error $\varepsilon(t; k^*)$ attains its maximum absolute value at (at least) 3 points in $[0, 1]$ with alternating signs.

Proof. This is a consequence of the Chebyshev equioscillation theorem applied to best uniform approximation. The approximation problem is: among all one-parameter families of curves (parametrized by $k$ ), find the one minimizing $\|\varepsilon(\cdot; k)\|_\infty$ . The equioscillation theorem (extended to nonlinear but unisolvent families) guarantees that at the optimum, the error equioscillates. The symmetry about $t = 1/2$ means the error is even about this point, so the extrema occur at $t = 0, 1/2, 1$ (boundary) and at symmetric interior points. $\square$

Lemma 2 (Numerical Solution). Solving the equioscillation conditions numerically (e.g., via Newton’s method on the system $\varepsilon'(t_1; k) = 0$ , $|\varepsilon(t_1; k)| = |\varepsilon(1/2; k)|$ ) yields

k^* \approx 0.5522847498.

The corresponding minimax radial error is $\varepsilon^* \approx 0.0000372091$ . $\square$

Editorial

Minimize the maximum radial deviation from the unit circle. By symmetry, control points are: P0 = (0,1), P1 = (k,1), P2 = (1,k), P3 = (1,0) We optimize k to minimize max |sqrt(x^2+y^2) - 1| over t in [0,1]. We binary search / Newton iteration on k. We then compute max and min of epsilon(t; k) on [0, 1]. Finally, by finding critical points of epsilon via derivative = 0.

Pseudocode

Binary search / Newton iteration on k
Compute max and min of epsilon(t; k) on [0, 1]
by finding critical points of epsilon via derivative = 0
k is too small or too large; adjust
else
equioscillation achieved
More precisely: use Newton's method
Let g(k) = max_t epsilon(t;k) + min_t epsilon(t;k)
Optimal k satisfies g(k) = 0 (symmetric equioscillation)

Complexity Analysis

Time: $O(I \cdot N)$ where $I$ is the number of Newton iterations (typically $< 20$ ) and $N$ is the number of function evaluations per iteration to find critical points (polynomial root-finding, $O(1)$ for a cubic).
Space: $O(1)$ .

Answer

\boxed{0.0000372091}

C++ project_euler/problem_363/solution.cpp

#include <bits/stdc++.h>
using namespace std;

/*
 * Problem 363: Bezier Curves
 *
 * Find optimal cubic Bezier approximation to a quarter circle.
 * Minimize the maximum radial deviation from the unit circle.
 *
 * By symmetry, control points are:
 *   P0 = (0,1), P1 = (k,1), P2 = (1,k), P3 = (1,0)
 *
 * We optimize k to minimize max |sqrt(x^2+y^2) - 1| over t in [0,1].
 *
 * Answer: 0.0000372091
 */

typedef long double ld;

// Evaluate Bezier curve at parameter t
pair<ld, ld> bezier(ld t, ld k) {
    ld s = 1.0L - t;
    ld x = 3.0L * s * s * t * k + 3.0L * s * t * t * 1.0L + t * t * t * 1.0L;
    ld y = s * s * s * 1.0L + 3.0L * s * s * t * 1.0L + 3.0L * s * t * t * k;
    return {x, y};
}

// Radial deviation at parameter t
ld radial_error(ld t, ld k) {
    auto [x, y] = bezier(t, k);
    return sqrtl(x * x + y * y) - 1.0L;
}

// Maximum absolute radial error over [0,1]
ld max_error(ld k) {
    ld mx = 0.0L;
    int N = 100000;
    for (int i = 0; i <= N; i++) {
        ld t = (ld)i / N;
        ld err = fabsl(radial_error(t, k));
        mx = max(mx, err);
    }
    return mx;
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);

    // Golden section search for optimal k
    ld lo = 0.5L, hi = 0.6L;
    ld gr = (sqrtl(5.0L) - 1.0L) / 2.0L;

    for (int iter = 0; iter < 200; iter++) {
        ld a = hi - gr * (hi - lo);
        ld b = lo + gr * (hi - lo);
        if (max_error(a) < max_error(b)) {
            hi = b;
        } else {
            lo = a;
        }
    }

    ld optimal_k = (lo + hi) / 2.0L;
    ld result = max_error(optimal_k);

    cout << fixed << setprecision(10) << result << endl;
    // Expected: 0.0000372091

    return 0;
}

Python project_euler/problem_363/solution.py

"""
Problem 363: Bezier Curves

Find optimal cubic Bezier approximation to a quarter circle.
Minimize the maximum radial deviation from the unit circle.

By symmetry, control points are:
  P0 = (0,1), P1 = (k,1), P2 = (1,k), P3 = (1,0)

We optimize k to minimize max |sqrt(x^2+y^2) - 1| over t in [0,1].

Answer: 0.0000372091
"""

import math
from scipy.optimize import minimize_scalar

def bezier(t, k):
    """Evaluate cubic Bezier at parameter t with symmetry parameter k."""
    s = 1.0 - t
    x = 3 * s * s * t * k + 3 * s * t * t * 1.0 + t ** 3
    y = s ** 3 + 3 * s * s * t * 1.0 + 3 * s * t * t * k
    return x, y

def radial_error(t, k):
    """Radial deviation from unit circle at parameter t."""
    x, y = bezier(t, k)
    return math.sqrt(x * x + y * y) - 1.0

def max_abs_error(k):
    """Maximum absolute radial error over t in [0, 1]."""
    N = 100000
    max_err = 0.0
    for i in range(N + 1):
        t = i / N
        err = abs(radial_error(t, k))
        if err > max_err:
            max_err = err
    return max_err

def solve():
    # Golden section search for optimal k
    result = minimize_scalar(max_abs_error, bounds=(0.5, 0.6), method='bounded',
                             options={'xatol': 1e-15})
    optimal_k = result.x
    min_max_dev = max_abs_error(optimal_k)
    return min_max_dev

if __name__ == "__main__":
    answer = solve()
    print(f"{answer:.10f}")
    # Expected: 0.0000372091