Fast subdivision of B\'{e}zier curves
Pith reviewed 2026-05-18 16:31 UTC · model grok-4.3
The pith
Bézier curve subdivision can be performed in O(d n log n) time using the fast Fourier transform and its inverse.
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
Core claim
The central claim is that subdivision of a polynomial Bézier curve reduces to a small number of forward and inverse FFTs on the control-point vectors, delivering O(d n log n) arithmetic operations instead of the O(d n²) cost of repeated linear interpolations in de Casteljau’s algorithm. The same FFT-based reduction also permits an O(d) update when the curve is extended by one additional control point and extends without essential change to rational curves, rectangular surfaces, and derivative evaluation.
What carries the argument
An FFT-based reduction that maps the control polygon to frequency-domain coefficients, performs the subdivision by pointwise scaling, and recovers the two new polygons via inverse transforms.
If this is right
- Subdivision of an extended curve can be refreshed in linear time rather than recomputed from scratch.
- The same technique yields an asymptotically faster algorithm for subdividing rational Bézier curves.
- Rectangular Bézier surface patches can be subdivided in each parameter direction with the same complexity gain.
- First and higher derivatives of a Bézier curve become available at reduced cost by reusing the same frequency-domain representation.
Where Pith is reading between the lines
- The O(d) update property could support incremental editing of high-degree curves in interactive design tools.
- Because FFT routines are highly optimized and parallelizable on current hardware, the method may translate directly into faster GPU-based curve rendering pipelines.
- The numerical-stability fix suggests that similar FFT reductions might be made robust for other polynomial-basis conversions in computational geometry.
Load-bearing premise
A slight change to the direct FFT procedure can preserve both the linearithmic complexity and sufficient numerical accuracy for degrees encountered in practice.
What would settle it
Compare the maximum coordinate-wise difference between the two output polygons produced by the modified FFT method and by de Casteljau on a random degree-64 curve in double precision; the claim fails if this difference grows faster than machine epsilon times a small constant.
read the original abstract
It is well-known that a $d$-dimensional polynomial B\'{e}zier curve of degree $n$ can be subdivided into two segments using the famous de Casteljau algorithm in $O(dn^2)$ time. Can this problem be solved more efficiently? In this paper, we show that it is possible to do this in $O(dn\log{n})$ time using the fast Fourier transform and its inverse. Experiments show that the direct application of the new method performs well only for small values of $n$, as the algorithm is numerically unstable. However, a slightly modified version -- which still has $O(dn\log{n})$ computational complexity -- offers good numerical quality, which is confirmed by numerical experiments conducted in \textsf{Python}. Moreover, the new method has a nice property: if a B\'{e}zier curve is extended by an additional control point, the subdivision can be updated in $O(d)$ time. A similar idea can be applied to speed up the subdivision of rational B\'{e}zier curves and rectangular B\'{e}zier surfaces, as well as to compute the derivatives of B\'{e}zier curves more efficiently.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The manuscript claims that subdividing a d-dimensional degree-n Bézier curve can be performed in O(d n log n) time via the FFT and its inverse, improving on the O(d n^2) de Casteljau algorithm. The direct FFT method is numerically unstable for practical n, but a slightly modified version retains the O(d n log n) bound while delivering acceptable accuracy, as shown by Python experiments. The approach also supports O(d)-time updates when extending the curve by one control point and extends analogously to rational curves, rectangular surfaces, and derivative evaluation.
Significance. If the central complexity claim for the stable variant holds with a rigorous operation count and the numerical results generalize, the result would be a meaningful advance for geometric algorithms in computer graphics, replacing a quadratic primitive with a near-linear one that also admits fast incremental updates. The empirical validation and the extensions to related problems strengthen the contribution, provided the stability modification is shown not to hide super-logarithmic costs.
major comments (2)
- Abstract: The statement that the 'slightly modified version -- which still has O(dn log n) computational complexity' is load-bearing for the main result, yet the manuscript provides no explicit recurrence, operation count, or pseudocode for the modification. If the fix involves iterative refinement, adaptive precision, or data-dependent reordering, the worst-case cost may exceed O(n log n); a concrete accounting in the complexity section is required to secure the claim.
- Section on the FFT-based procedure (likely §3): The derivation that subdivision reduces to FFT-based coefficient conversion must be checked for hidden factors. Standard FFT polynomial multiplication is O(n log n), but converting between Bernstein and monomial bases or enforcing stability may introduce additional linear algebra steps whose cost is not shown to remain O(n log n) in the worst case.
minor comments (2)
- Abstract: The phrase 'O(dn log n)' should be written with consistent spacing and a clear base for the logarithm (FFT literature conventionally uses base 2).
- Experiments: Include wall-clock timings and maximum absolute error tables for both the direct and modified methods across a range of n up to at least 1024, with direct comparison to de Casteljau.
Simulated Author's Rebuttal
We thank the referee for the careful and constructive report. The comments highlight the need for explicit complexity accounting for the stable variant and verification of all steps in the FFT derivation. We will revise the manuscript to include these details while preserving the core claims.
read point-by-point responses
-
Referee: Abstract: The statement that the 'slightly modified version -- which still has O(dn log n) computational complexity' is load-bearing for the main result, yet the manuscript provides no explicit recurrence, operation count, or pseudocode for the modification. If the fix involves iterative refinement, adaptive precision, or data-dependent reordering, the worst-case cost may exceed O(n log n); a concrete accounting in the complexity section is required to secure the claim.
Authors: We agree that an explicit recurrence and operation count for the modified algorithm are required to rigorously support the O(d n log n) bound. The modification consists of a constant number of additional FFT and inverse-FFT passes (plus O(n) post-processing) that do not change the asymptotic cost. In the revised manuscript we will add a dedicated complexity subsection containing the recurrence, a full operation count, and pseudocode for the stable procedure. revision: yes
-
Referee: Section on the FFT-based procedure (likely §3): The derivation that subdivision reduces to FFT-based coefficient conversion must be checked for hidden factors. Standard FFT polynomial multiplication is O(n log n), but converting between Bernstein and monomial bases or enforcing stability may introduce additional linear algebra steps whose cost is not shown to remain O(n log n) in the worst case.
Authors: The Bernstein-to-monomial conversion is performed directly by the FFT (exploiting the explicit binomial structure) without general matrix multiplication or other linear-algebra primitives. Stability is restored by a fixed number of additional FFT passes whose cost is absorbed into the O(n log n) term. We will expand the relevant section with an itemized complexity breakdown of every algorithmic step to make these bounds explicit and to confirm the absence of hidden super-logarithmic factors. revision: yes
Circularity Check
No circularity; derivation uses independent standard FFT properties
full rationale
The paper derives an O(dn log n) subdivision algorithm for degree-n Bézier curves by invoking the fast Fourier transform and its inverse, which are external, well-established algorithmic primitives whose complexity and correctness are independent of the Bézier problem and do not reduce to any fitted parameter, self-definition, or prior self-citation within the paper. The abstract's reference to a 'slightly modified version' that retains the same complexity for numerical stability is presented as an assertion supported by experiments, without any quoted reduction showing that the claimed bound is forced by construction from the inputs or from a self-referential definition. No load-bearing step equates the output to the input by renaming, fitting, or imported uniqueness; the derivation remains self-contained against external benchmarks.
Axiom & Free-Parameter Ledger
axioms (1)
- standard math The discrete Fourier transform and its inverse correctly evaluate and interpolate polynomials of degree at most n in O(n log n) arithmetic operations.
Reference graph
Works this paper leans on
-
[1]
Chudy, New algorithms for Bernstein polynomials, their dual bases, and B-spline func- tions, Ph.D
F. Chudy, New algorithms for Bernstein polynomials, their dual bases, and B-spline func- tions, Ph.D. thesis, University of Wrocław (2022)
work page 2022
- [2]
-
[3]
T. H. Cormen, C. E. Leiserson, R. L. Rivest, C. S., Introduction to Algorithms, 2nd ed., The MIT Press, 2001
work page 2001
-
[4]
Farin, Curves and surfaces for Computer-Aided Geometric Design
G. Farin, Curves and surfaces for Computer-Aided Geometric Design. A practical guide, 5th ed., Academic Press, Boston, 2002
work page 2002
-
[5]
R. T. Farouki, The Bernstein polynomial basis: A centennial retrospective, Computer Aided Geometric Design 29 (2012) 379–419
work page 2012
-
[6]
Fuda, Numerical stability of barycentric interpolation, Ph.D
C. Fuda, Numerical stability of barycentric interpolation, Ph.D. thesis, Universit` a della Svizzera italiana (2024)
work page 2024
-
[7]
C. Fuda, A. Ramanantoanina, K. Hormann, A comprehensive comparison of algorithms for evaluating rational B´ ezier curves, Dolomites Research Notes on Approximation 17 (2023) 56–79
work page 2023
- [8]
discussion (0)
Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.