We present new representations of Gauss--Legendre polynomials and their derivatives in the shifted power basis and in bases related to symmetric orthogonal Jacobi polynomials. Using these representations and certain recurrence relations, we propose $O(n^2+dn)$ methods for evaluating a Gauss--Legendre curve of degree $n$ in $\mathbb E^d$. We also propose algorithms for multipoint evaluation with computational complexity $O(Mdn+dn^2)$, where $M$ is the number of evaluation points. Numerical tests show that the approach using the shifted power basis is very fast, but it works well only for small $n$ (approximately up to 20), whereas the method employing Jacobi polynomials gives good numerical results even for large $n$ (around 100).
It is well-known that a $d$-dimensional polynomial Bézier curve of degree $n$ can be subdivided into two segments using the famous de Casteljau algorithm in $O(dn^2)$ time. We show that it is possible to do this in $O(dn\log{n})$ time using the fast Fourier transform and its inverse. The direct application of the method is numerically unstable for larger values of $n$. Therefore, we also introduce a modified version based on parameter scaling, which retains the same asymptotic computational complexity. This scaling approach can improve numerical accuracy in some cases (for $n\lesssim 50$), as demonstrated by our experiments. Moreover, the new method has a nice property: if a Bézier curve is extended by an additional control point, the subdivision can be updated in $O(d)$ time. The method can be adapted to subdivide rational Bézier curves and rectangular Bézier surfaces, as well as to compute the derivatives of Bézier curves. Although motivated by applications in CAGD, the proposed methods may also be useful in other areas involving polynomial reparameterization.