Generics in Rust: little library for Bezier curves -- Part 2

Some time ago, I decided to write a couple of posts about my experience with generic programming in Rust. For me, when someone says generics, C++ template metaprogramming and Alexander Stepanov immediately pop in my mind. Rust is different. So it was interesting to see what's going on out there, which motivated my original interest.

Learning by doing is the best for quick hands-on, so I decided to write a little generic numeric library for manipulating Bezier curves (polynomials in the Bernstein basis) that is (i) uses static dispatch (and no heap allocation calls), and (ii) can be used with different types for specifying Bezier control polygon: reals, rationals, complex, and, in general, something that implements standard vector-space operations.

What I learned is that writing generic libraries in Rust is not a piece of cake. There are two major facts that contribute into this: (i) Rust's explicit safety-oriented type system, and (ii) the absence of decent support of generic constant expressions in stable Rust. So what I am writing about here is based on Rust Nightly.

So before we take of, there is a couple of previous posts

• Generic constant expressions: a future bright side of nightly Rust.
• Experimenting with generics in Rust: little library for Bezier curves - part 1.

And the Github repo that contains examples from these posts, some tests, and maybe even some demos, if I will manage to add them in nearby future.

First steps

First, make sure we use the unstable rust build. I will just set it up for all repositories with rustup default nightly, but it is possible to apply it to a folder with rustup override set nightly. After that:

PS > rustc --version
rustc 1.80.0-nightly (1ba35e9bb 2024-05-25)

and we also make sure that we include the following lines in our crate

#![allow(incomplete_features)]
#![feature(generic_const_exprs)]

As I briefly outlined here, Rust currently does not support generic constant expression in its type system. It is not that something is wrong with it in general, it is just not implemented due to technical difficulties. This is considered to be an unstable feature.

To describe a generic Bezier curve c(u) = (x(u), y(u), z(u), ...), I basically wrap a primitive array type [T; N] into a struct

pub struct Bernstein<T, U, const N: usize> {
    coef: [T; N],
    segm: (U, U),
}

that contains a generic type parameter T for Bezier control polygon, type parameter U for curve parametrization, and a generic constexpression parameter N for the number of basis polynomials (or just the size of the Bezier control polygon). For example, a cubic Bezier curve should have four points in its control polygon, i. e. N = 4. The curve is always parameterized on 0 <= u <= 1, so right now segm is defaulted to (0, 1).

Some more details can be found in my previous post. In this post, I would like to discuss how implement generic methods on this type that would leverage generic constant expressions on the size of the control polygon N.

Implementing eval-method

The first method to do is, of course, to implement eval() method to find a point on the curve at some value of the parameter u, so that we can write something like this

let p0 = Complex::new(0.0, 0.0);
let p1 = Complex::new(2.5, 1.0);
let p2 = Complex::new(-0.5, 1.0);
let p3 = Complex::new(2.0, 0.0);

// Define cubic Bezier curve in the complex plane.
let c: Bernstein<Complex<f32>, f32, 4> = Bernstein::new([p0, p1, p2, p3]);

let p = c.eval(0.5); // point on a curve of type Complex<f32>

This part is easy, and can be done even in stable Rust. I just use the De Casteljau's algorithm

impl<T, U, const N: usize> Bernstein<T, U, N> where
    T: Copy + Add<T, Output = T> + Sub<T, Output = T> + Mul<U, Output = T>,
    U: Copy + Num,
{
    pub fn eval(&self, u: U) -> T {
        // -- snippet --
        // De Casteljau's algorithm
    }
}

Trait bounds on types are quite transparent. I require Copy trait to make my life easier when manipulating mathematical expressions, type U should be a number, which is required by num::Num trait. This is especially useful because Num trait requires generic One and Zero traits, which provide methods such as U::zero() and U::one(). Type T is required to implement vector space operations of addition, subtraction, and right-hand-side multiplication by a variable of type U with the result of being of type T (Mul<U, Output = T>).

Implementing diff and integ methods

The next step is to implement generic diff() and integ() methods to find the parametric derivative of the Bezier curve dc(u)/du, and the integral with respect to parameter u. That's where generic constant expression come into play.

The problem is that our methods should take an array of control points [T; N] of size N as an input, and return an array of size N - 1 for diff() or N + 1 for integ() as output nicely wrapped into our custom Bernstein type so that the signatures of the functions should be like these:

fn diff(&self) -> Bernstein<T, U, {N - 1}> {}

// c: T -- is the initial point to fix the constant of integration.
fn integ(&self, c: T) -> Bernstein<T, U, {N + 1}> {}

And stable Rust does not allow us to do that. Using generic_const_exprs feature, it becomes possible, as we shall see shortly.

Another difficulty is related to Rust's explicit type system. In these methods, the size of the array N becomes a part of mathematical expressions in the diff() and integ() algorithms. Rust requires the size of the array to be of a machine-dependent pointer size usize (which totally make sense but it is not generic). Converting from usize to other type is not considered to be a safe operation so I have to rely on third party traits for that purpose, such as num::FromPrimitive trait that is implemented for usize in the num crate. Otherwise, multiplying and expression of type, let's say, f64, by N is not defined.

Having this in mind, let's discuss diff() method (implementing integ() is similar and may be found in the repo):

impl<T, U, const N: usize> Bernstein<T, U, N> where
    T: Copy + Add<T, Output = T> + Sub<T, Output = T> + Mul<U, Output = T>,
    U: Copy + Num + FromPrimitive,
{
    pub fn diff(&self) -> Bernstein<T, U, {N - 1}> where
        [(); N - 1]:
    {
        let coef: [T; N - 1] = array::from_fn(
            |i| -> T {
                (self.coef[i + 1] - self.coef[i]) *
                (U::from_usize(N - 1).unwrap() / (self.segm.1 - self.segm.0))
            }
        );

        Bernstein {
            segm: self.segm,
            coef: coef,
        }
    }
}

Here, there is a couple of new details. First, I require type U to be bounded by FromPrimitive trait that allows to convert from usize to U in a generic environment by calling U::from_usize(N - 1).unwrap(). Second is that there is new bound [(); N - 1]: which is required by generic_const_exprs feature

We currently use where [(); expr]: as a way to add additional const wf bounds. Once we have started experimenting with this it is probably worth it to add a more intuitive way to add const wf bounds.

The bounds on T type is basically the same as in the eval() method.

Now, we can find a derivative type from a basic type, by using for example the following

// Define cubic Bezier curve in the complex plane.
let c: Bernstein<Complex<f32>, f32, 4> = Bernstein::new([p0, p1, p2, p3]);

// Get the derivative, or hodograph curve at u = 0.2.
let d = c.diff().eval(0.2); // `d` is of type Complex<f32>

Generic product of two polynomials

The next example is a product of two polynomials of order N - 1 and M - 1 (the size of arrays of coefficients is N and M respectively). This is a little bit more involved since it has to take the array [T; N] as an input, multiply it by [T; M] and the type of output should be [T; M + N - 1]. For example, multiplying a polynomial of the third order (N = 4) by a second-order polynomial (N = 3) should give a quintic polynomial (N = 6).

The implementation may looks like this:

impl<T, U, const N: usize, const M: usize>
Mul<Bernstein<T, U, {M}>> for Bernstein<T, U, {N}> where
    T: Copy + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Mul<U, Output = T>,
    U: Num + FromPrimitive,
    [(); N]:,
    [(); M]:,
    [(); N + M - 1]:
{
    type Output = Bernstein<T, U, {N + M - 1}>;

    fn mul(self, rhs: Bernstein<T, U, {M}>) -> Self::Output {
        let mut coef = [self.coef[0] - self.coef[0]; N + M - 1];

        // -- snippet --
        // actual algorithm

        Bernstein {
            coef: coef,
            segm: self.segm,
        }
    }
}

The required operation is specified by Mul<Bernstein<T, U, {M}>> in the impl, and the resulting type should be type Output = Bernstein<T, U, {N + M - 1}>.

Another subtle moment is that I have to initialize the array in the body of the function mul() because Rust does not allow to use uninitialized variables (remember old plain C89? -- who cared about initializing all the variables). One way to do it is to put a trait bound on T to implement T::zero() that can be used as an initial value. In this case, I chose a workaround instead (may change it later) which is to require subtraction Sub<Output = T> and use self.coef[0] - self.coef[0] as a kind of generic zero.

Note that generic-const-exprs require additional trait bounds to be imposed for each of the array types we use [(); N]:, [(); M]:, [(); N + M - 1]:.

Now, it possible to write

let p: Bernstein<f64, f64, 3> = Bernstein::new([0.0, 2.0, 1.0]);
let q: Bernstein<f64, f64, 4> = Bernstein::new([1.0, 2.0, 0.0, 0.0]);

// Quintic polynomial with real coefficient 
let c = p * q;

Summary

Generic constant expressions in Rust give the flexibility of implementing generic types which size is known at compile time. So far, using an unstable Rust nightly 1.80, I didn't notice any issues with using generic-const-exprs feature.