Problem 363: Bézier Curves

Detailed mathematical approach

Problem Summary

A cubic Bézier curve is built from four points \(P_0,P_1,P_2,P_3\) by repeated interpolation. For a fixed ratio \(t\) in \([0,1]\), mark the three points dividing \(P_0P_1\), \(P_1P_2\), \(P_2P_3\) in that ratio, then the two dividing the segments between those, then the point \(B\) dividing the last one. As \(t\) runs through \([0,1]\), \(B\) traces the curve, which leaves \(P_0\) along \(P_0P_1\) and reaches \(P_3\) along \(P_2P_3\). ^[1]

Here \(P_0=(1,0)\), \(P_1=(1,v)\), \(P_2=(v,1)\), \(P_3=(0,1)\), so the curve joins \((1,0)\) to \((0,1)\) and imitates a quarter of the unit circle. The positive number \(v\) is fixed by one demand: the region bounded by the segment \(OP_0\), the curve and the segment \(OP_3\) must have the quarter disc's area \(\pi/4\). For that \(v\), with \(L\) the length of the curve, we must give

$$100\cdot\frac{L-\pi/2}{\pi/2}$$

rounded to ten digits after the decimal point.

The whole solution in one sentence: the construction makes \(x\) and \(y\) cubic polynomials in \(t\), so the area condition becomes a quadratic equation in \(v\) that we solve in closed form, and the length is then one smooth integral that adaptive Simpson quadrature evaluates to far more digits than the ten we need.

Mathematical Approach

1. From the construction to a polynomial

The point dividing a segment \(XY\) in the ratio \(t\) is \(X+t(Y-X)\). Writing the three levels in that form, \(Q_i=P_i+t(P_{i+1}-P_i)\), then \(R_i=Q_i+t(Q_{i+1}-Q_i)\), then \(B=R_0+t(R_1-R_0)\), and expanding gives

$$\begin{aligned}B(t)={}&(1-t)^3P_0+3(1-t)^2tP_1\\&+3(1-t)t^2P_2+t^3P_3.\end{aligned}$$

Every step is an affine combination, so the expansion may be done one coordinate at a time. This is de Casteljau's algorithm, and the coefficients are the cubic Bernstein basis polynomials, non-negative on \([0,1]\) and summing to \(1\); the curve is a weighted average of the control points and never leaves their convex hull. ^[2]

A case small enough to check by hand is \(v=1\) at \(t=\tfrac12\), where every mark is a midpoint. From \((1,0)\), \((1,1)\), \((1,1)\), \((0,1)\) we get \(Q=(1,\tfrac12)\), \((1,1)\), \((\tfrac12,1)\), then \(R=(1,\tfrac34)\), \((\tfrac34,1)\), and finally \(B=(\tfrac78,\tfrac78)\) — what the coefficients \(\tfrac18,\tfrac38,\tfrac38,\tfrac18\) give.

2. The curve of this problem

Substituting the control points and collecting powers of \(t\) gives exactly the polynomials all three programs use:

$$x(t)=1+3(v-1)t^2+(2-3v)t^3,$$

$$y(t)=3vt+(3-6v)t^2+(3v-2)t^3,$$

$$x'(t)=3t\bigl((2-3v)t-2(1-v)\bigr),$$

$$y'(t)=3(t-1)\bigl((3v-2)t-v\bigr).$$

Reflection in the line \(y=x\) swaps \(P_0\) with \(P_3\) and \(P_1\) with \(P_2\), reversing the list of control points, which replaces \(t\) by \(1-t\). The curve is therefore symmetric in that line: \(y(t)=x(1-t)\), as the polynomials confirm, and the speed satisfies \(s(t)=s(1-t)\). At the ends \(x'(0)=0\) with \(y'(0)=3v\), and \(y'(1)=0\) with \(x'(1)=-3v\): the tangent at \(P_0\) is vertical and the one at \(P_3\) horizontal, the tangency the statement mentions.

For \(0\lt v\lt1\) the brackets settle two signs. Each is linear in \(t\) and negative at both ends of \([0,1]\) — the first takes \(-2(1-v)\) and \(-v\), the second \(-v\) and \(2v-2\) — hence negative throughout. So \(x'(t)\lt0\) for \(0\lt t\le1\), and \(y'(t)\gt0\) for \(0\le t\lt1\) because \(3(t-1)\) is negative there too. The arc is a strictly monotone graph inside the unit square and never touches itself.

3. The enclosed area as a line integral

Green's theorem turns a double integral over a region into a line integral around its positively oriented, piecewise smooth, simple boundary; choosing the integrand so that the double integral becomes the area gives \(A=\oint x\,dy\). ^[3] Our boundary, counter-clockwise, runs from \(O\) to \(P_0\) along the \(x\)-axis, along the arc to \(P_3\), and back to \(O\) along the \(y\)-axis; by Section 2 it is simple, so the theorem applies. On the first piece \(y\) is constant, so \(dy=0\), and on the third \(x=0\): only the arc is left,

$$A(v)=\int_0^1x(t)\,y'(t)\,dt.$$

4. The area is a quadratic in \(v\)

Multiplying the two polynomials of Section 2 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 \(t^k\) over \([0,1]\) divides its coefficient by \(k+1\), and collecting powers of \(v\) leaves a very short result:

$$A(v)=\frac12+\frac{3v}{5}-\frac{3v^2}{20}.$$

Two endpoint values check it, the first by geometry, the second by an elementary integral. For \(v=0\) the inner control points collapse onto the outer ones, \(x(t)+y(t)=1\) identically, and the region is the triangle with legs \(1\): area \(\tfrac12=A(0)\). For \(v=1\) we have \(P_1=P_2=(1,1)\), so \(x(t)=1-t^3\), \(y'(t)=3(1-t)^2\) and

$$\begin{aligned}\int_0^1(1-t^3)\,3(1-t)^2\,dt&=3\cdot\frac{19}{60}\\&=\frac{19}{20},\end{aligned}$$

which is \(\tfrac12+\tfrac35-\tfrac3{20}\) as promised. Since \(A'(v)=\tfrac35-\tfrac{3v}{10}\) is positive for \(v\lt2\), the area grows strictly from \(\tfrac12\) to \(\tfrac{19}{20}\) as \(v\) goes from \(0\) to \(1\), and \(\pi/4\approx0.7853981634\) lies strictly between: exactly one admissible \(v\) exists in \((0,1)\).

5. Solving for \(v\)

Multiplying \(A(v)=\pi/4\) by \(20\) gives \(10+12v-3v^2=5\pi\), that is

$$\begin{aligned}3v^2-12v+5\pi-10&=0,\\v&=2\pm\frac{\sqrt{66-15\pi}}{3}.\end{aligned}$$

Here \(66-15\pi\approx18.8761101962\) and \(\sqrt{66-15\pi}/3\) is about \(1.4482215222\), so the roots are near \(0.5517784778\) and \(3.4482215222\). The closed form is exact: with \(d=2-v\) the area is \(\tfrac{11}{10}-\tfrac{3}{20}d^2\), and \(d^2=(66-15\pi)/9\) makes \(\tfrac{3}{20}d^2\) equal \(\tfrac{11}{10}-\pi/4\).

The larger root describes nothing of the kind asked for. Subtracting the polynomials of Section 2 and factoring,

$$x(t)-y(t)=(1-2t)\bigl(1-(3v-2)t(1-t)\bigr),$$

so besides \(t=\tfrac12\) the curve meets the line \(y=x\) wherever \(t(1-t)=1/(3v-2)\). As \(t(1-t)\) never exceeds \(\tfrac14\), that happens exactly when \(v\ge2\) — and the larger root is above \(2\). Its two solutions are swapped by \(t\mapsto1-t\), and by the mirror symmetry of Section 2 the two curve points then coincide: the path crosses itself. It also leaves the unit square, since \(y'\) now vanishes at \(t_2=v/(3v-2)\approx0.4132247\) and \(y\) climbs to about \(1.8429355\) before coming back. Around such a path the line integral is a signed area, not the area of the single region the statement describes. The curve of the problem is therefore

$$v=2-\frac{\sqrt{66-15\pi}}{3}=0.5517784778\ldots$$

6. The length is an elliptic integral

The length of a parametrised arc is the integral of its speed:

$$L=\int_0^1\sqrt{x'(t)^2+y'(t)^2}\,dt.$$

Under the root is a quartic. Centred at \(u=t-\tfrac12\), where the speed symmetry kills the odd powers, it is \(q=au^4+bu^2+c\) with

$$\begin{aligned}a&=18(3v-2)^2,\\b&=-9(v^2-8v+4),\\c&=\tfrac98(v-2)^2.\end{aligned}$$

Here \(a\) and \(c\) are positive multiples of squares, and \(b\) is positive exactly for \(v\gt4-2\sqrt3\approx0.5358984\), which our \(v\) satisfies. So \(q\) is positive for every real \(t\) and grows with \(u^2\): the speed is least at the midpoint, \(\sqrt c=\sqrt{132-30\pi}/4\), about \(1.5360708885\), and greatest at the ends, \(3v\approx1.6553354334\). Such a positive, nearly constant integrand is what makes the quadrature cheap.

An integral of a rational function of \(t\) and the square root of a quartic is of elliptic type, and these have in general no elementary expression. ^[4] Two of the degenerate cases that reference lists are ruled out here, as \(q\) has no repeated root and the integrand is an odd power of \(\sqrt q\), so no elementary antiderivative is to be expected and the programs integrate numerically.

7. Why the length must exceed \(\pi/2\)

The sign of the answer is not a numerical accident. Reflect the arc in the \(x\)-axis, in the \(y\)-axis, and in both. The four copies join end to end — the arc runs from \((1,0)\) to \((0,1)\), its mirror in the \(y\)-axis from \((-1,0)\) to \((0,1)\) — and form a closed curve \(\Gamma\) of length \(4L\) through \((\pm1,0)\) and \((0,\pm1)\). By Section 2 each copy is a monotone graph in its own quadrant, so \(\Gamma\) is simple and its interior meets the first quadrant in our region alone: it encloses the area \(4\cdot\pi/4=\pi\).

The isoperimetric inequality states that a closed plane curve of length \(P\) enclosing an area \(A\) satisfies \(4\pi A\le P^2\), with equality only for a circle. ^[5] For \(\Gamma\) this reads \(4\pi\cdot\pi\le(4L)^2\), so \(\pi^2\le4L^2\) and \(L\ge\pi/2\). Equality would force \(\Gamma\) to be a circle, and a circle through those four points is the unit circle, making the arc the quarter circle. It is not: a unit circle has curvature constantly \(1\), whereas here

$$\kappa=\frac{x'y''-y'x''}{(x'^2+y'^2)^{3/2}}$$

equals \(2(1-v)/(3v^2)\approx0.9814584841\) at both ends. Equality is therefore impossible: \(L\gt\pi/2\) strictly, and the percentage asked for is positive.

The two conditions are genuinely different. Measured radially, the arc runs inside the circle in the middle, by about \(0.000268491\) at \(t=\tfrac12\), and outside it nearer the ends, by about \(0.000172081\); the area condition balances those deviations exactly, and nothing makes the length balance too.

8. Simpson's rule and the adaptive scheme

Simpson's rule approximates an integral from the ends and the midpoint \(m\) of an interval,

$$S(a,b)=\frac{b-a}{6}\bigl(f(a)+4f(m)+f(b)\bigr),$$

and is exact for polynomials of degree at most \(3\). ^[6] Let \(n\) count whole panels of that kind, each of width \(h=(b-a)/n\) with its own ends and midpoint. Adding up the one-panel error \(-h^5f^{(4)}(\xi)/2880\) gives the composite value \(S_n\) the error

$$I-S_n=-\frac{(b-a)h^4}{2880}f^{(4)}(\xi),$$

so halving the panel width divides the leading error by about \(16\). For a check by hand take \(f(t)=t^4\) on \([0,1]\), whose integral is \(\tfrac15\): one panel gives \(\tfrac5{24}\) and two give \(\tfrac{77}{384}\), errors of \(-\tfrac1{120}\) and \(-\tfrac1{1920}\) whose ratio is exactly \(16\). On our speed the ratio of consecutive errors falls to \(16.0046\) by \(n=64\).

Writing \(S_1=S(a,b)\) and \(S_2=S(a,m)+S(m,b)\), the finer estimate is about \(16\) times better, so its error is roughly one fifteenth of what the refinement changed: \(I-S_2\approx(S_2-S_1)/15\). Lyness's criterion therefore accepts an interval when \(|S_2-S_1|\le15\varepsilon\) and returns \(S_2+(S_2-S_1)/15\), a combination of the five samples which is the five-point Newton–Cotes rule and so exact up to degree \(5\); otherwise the interval is halved and each half treated with the tolerance \(\varepsilon/2\). ^[7] The test compares two estimates rather than the true error, so it is a heuristic: near the rounding noise their computed difference may fall to zero and pass, or keep failing however small the interval becomes. The depth limit then guarantees termination, but the value returned carries no promise of meeting the tolerance.

9. How much precision ten decimals need

The quantity asked for is \(p=(200/\pi)L-100\), so \(dp/dL=200/\pi\approx63.662\): a budget of \(0.5\cdot10^{-10}\) in \(p\) allows about \(7.85\cdot10^{-13}\) in \(L\), roughly twelve decimal places of it, since \(L-\pi/2\) cancels about six leading digits. One unit in the last place of \(L\) in 53-bit arithmetic is \(2.2\cdot10^{-16}\), some \(3500\) times smaller, so double precision suffices if the quadrature error is controlled. The parameter needs care too: \(dL/dv\approx0.4340266455\), so \(v\) must be right to \(1.8\cdot10^{-12}\), and the closed form in double precision is off by \(6.9\cdot10^{-17}\), under one unit in its last place. Such a budget is not a proof of correct rounding, for which the error interval of \(p\) would have to lie inside one rounding interval.

How the Code Works

All three files share four numerical helpers — solve_v, curve_speed, adaptive_simpson and curve_length, which Java spells solveV, curveSpeed, adaptiveSimpson and curveLength. Python and Java wrap the final arithmetic and formatting in solve; the C++ does that in main and adds curve_area and run_checkpoints. ^[8] solve_v returns the closed form of Section 5, using acos(-1) for \(\pi\), so no iterative root search is needed anywhere. curve_speed evaluates the derivative polynomials of Section 2 and returns the square root of \(dx^2+dy^2\), and curve_length evaluates it at \(0\), \(\tfrac12\) and \(1\), forms the one-panel Simpson value whole, and calls adaptive_simpson on \([0,1]\) with the tolerance eps and depth \(30\).

adaptive_simpson implements Section 8. It computes the midpoint m and the quarter points lm, rm, evaluates the speed only there — fa, fm, fb arrive as arguments — and forms left, right, combined. If depth has run out or Lyness's test passes it returns combined + (combined - whole) / 15; otherwise it recurses on the halves, handing each the Simpson value just computed as its new whole and eps * 0.5 as its tolerance. With the \(10^{-15}\) of the Python and Java files it accepts \(944\) subintervals after \(3777\) speed evaluations, never passes level \(10\), and ends on an almost uniform grid.

run_checkpoints tests four things before the main computation: \(A(0)=\tfrac12\) to \(10^{-18}\); curve_length(0) equal to \(\sqrt2\) to \(10^{-12}\); \(0\lt v\lt1\); and \(A(v)\) equal to \(\pi/4\) to \(10^{-18}\). The second has an exact target: at \(v=0\) the speed is \(6\sqrt2\,t(1-t)\), a quadratic, which Simpson integrates exactly, so the first panel is accepted after five evaluations and the result is \(\sqrt2\) to one unit in the last place. The fourth puts the computed root back into the same area formula, so it checks consistency rather than the derivation, and its \(10^{-18}\) is demanding: one unit in the last place near \(\pi/4\) is \(1.11\cdot10^{-16}\) with a 53-bit significand but only \(5.42\cdot10^{-20}\) with a 64-bit one. In binary64 the residual came out exactly \(0\); with a 64-bit significand it was one unit in the last place, nonzero but inside the tolerance. main accepts --skip-checkpoints and exits with \(1\) on any other argument, \(2\) on a failed checkpoint.

One difference is deliberate: the C++ carries long double with the tolerance \(10^{-18}\), Python and Java double precision with \(10^{-15}\). An emulation with a 64-bit significand accepts \(5823\) subintervals and reaches level \(13\) without touching the depth limit, so on an 80-bit extended type the tighter tolerance costs little; where long double has only a 53-bit significand it is out of reach, and one measured run accepted \(209530\) subintervals, \(60938\) of them stopped by the depth limit and the rest passing on a gap that had rounded to zero. Such counts depend on the platform and the compiler settings; the ten printed decimals did not.

Complexity and Verification

Finding \(v\) is \(O(1)\). The quadrature cost is set by the tolerance, not by any input size. A recursion with \(M\) accepted leaves makes \(2M-1\) calls, each evaluating the speed at two new points on top of the three samples of the first call, so it spends exactly \(4M+1\) evaluations — \(3777\) for \(M=944\), as measured — each one square root and a few arithmetic operations. While truncation error dominates rounding and the depth cap stays idle, a leaf of width \(h=2^{-d}\) is accepted once its Simpson gap, of order \(h^5\), drops below the \(15\varepsilon/2^{d}\) allotted to it, which needs \(2^{4d}\) to exceed a constant times \(1/\varepsilon\); so \(M\) grows like \(\varepsilon^{-1/4}\). The measured \(188\) leaves at \(\varepsilon=10^{-12}\) against \(944\) at \(10^{-15}\) fit that scale; at far smaller tolerances rounding and the cap take over. No array is kept, so memory is the recursion stack: a depth argument of \(30\) permits \(31\) nested calls, a fixed bound.

Every claim above was checked by computation. Symbolically: the construction of Section 1 equals the Bernstein cubic for arbitrary control points; \(x,y,x',y'\) and their factorisations match all three sources; \(y(t)=x(1-t)\); the product of Section 4, the factorisation of \(x-y\) and the centred quartic of Section 6 are identities; and \(A(v)-\pi/4=0\) exactly. By hand the construction gives \((\tfrac78,\tfrac78)\), \(A(0)=\tfrac12\) and \(A(1)=\tfrac{19}{20}\). Numerically, integrating \(x\,y'\) at the chosen \(v\) returns \(\pi/4\) to 50 digits, and the shoelace area of a 400,000-point polygon along the real boundary differs from \(\pi/4\) by \(2.0\cdot10^{-12}\), the expected size for a polygon.

The published routine was then measured against independent quadrature. At its tolerance \(10^{-15}\) it differs from a 50-digit tanh–sinh evaluation by \(2.8\cdot10^{-16}\) in \(L\) and \(1.1\cdot10^{-14}\) in the percentage, about \(4500\) times smaller than the budget of Section 9, and tanh–sinh and Gauss–Legendre agree to all 60 digits asked of them. The tolerances \(10^{-12}\), \(10^{-15}\) and \(10^{-18}\), in double precision and with a 64-bit significand, all give the same ten decimals, and the three programs print byte-identical output. The sign is not left to arithmetic: Section 7 proves \(L\gt\pi/2\), so a negative or zero result would mean a bug rather than rounding.

Footnotes and References

The references supply the background facts used above; the area polynomial, the arc's monotonicity and symmetry, the rejected root, the reflection argument and the measured counts are derived or computed here.

1. Project Euler 363 — Bézier Curves. The official statement defines the curve by the repeated-ratio construction through Q₀, Q₁, Q₂, R₀, R₁ and B, notes its tangency to P₀P₁ and P₂P₃, and sets P₀ = (1, 0), P₁ = (1, v), P₂ = (v, 1), P₃ = (0, 1) with v > 0 making the area enclosed by OP₀, OP₃ and the curve π/4; it asks for 100(L − π/2)/(π/2) to ten decimals.
2. Wikipedia — Bézier curve. The article writes the cubic with the Bernstein weights (1 − t)³, 3(1 − t)²t, 3(1 − t)t², t³, names de Casteljau's algorithm behind the construction, and states the convex-hull property and the tangency to the first and last section of the control polygon.
3. Wikipedia — Green's theorem. For a positively oriented, piecewise smooth, simple closed curve bounding a region, the line integral of L dx + M dy around it equals the double integral of ∂M/∂x − ∂L/∂y; with that integrand 1 the area becomes ∮ x dy.
4. Wikipedia — Elliptic integral. An elliptic integral is one of R(t, √P(t)) with R rational and P of degree 3 or 4 without repeated roots, and in general such integrals have no elementary expression; its exceptions are repeated roots, an even dependence on the root, and pseudo-elliptic cases. This is background for Section 6, not a proof about that quartic.
5. Wikipedia — Isoperimetric inequality. A closed plane curve of length L enclosing an area A satisfies 4πA ≤ L², with equality only for a circle; the proofs quoted there cover arbitrary rectifiable curves, all Section 7 needs.
6. Wikipedia — Simpson's rule. The rule approximates the integral over [a, b] by (b − a)/6 times f(a) + 4f((a + b)/2) + f(b), is exact for polynomials of degree three or less, and has the single-interval error −(b − a)⁵f⁗(ξ)/2880. Its composite error −(b − a)h⁴f⁗(ξ)/180 counts h as the node-to-node spacing, half a panel width; Section 8 takes h as the panel width, turning 180 into 2880.
7. Wikipedia — Adaptive Simpson's method. The method estimates the error of a refinement by S(a, m) + S(m, b) − S(a, b), accepts it when that difference is below 15ε (Lyness's criterion), adds one fifteenth of it as a Richardson correction that makes the result exact up to degree five, and otherwise halves the tolerance for each half; a depth counter is its usual guard against an unreachable tolerance.
8. C++ — Euler363.cpp. The linked immutable revision contains curve_area, curve_speed, adaptive_simpson, curve_length, solve_v and run_checkpoints. It is the implementation record for this article.