Generics in Rust: little library for Bezier curves -- Part 2
Tags: rust, generic, beginners, numeric
Published: 2024-05-27 14:10:20

Some time ago, I decided to write a couple of posts about my experience with [generic programming](https://en.wikipedia.org/wiki/Generic_programming) in Rust. For me, when someone says generics, [C++ template metaprogramming](https://en.cppreference.com/w/cpp/language/template_metaprogramming) and [Alexander Stepanov](https://en.wikipedia.org/wiki/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](https://en.wikipedia.org/wiki/B%C3%A9zier_curve) (polynomials in the [Bernstein basis](https://en.wikipedia.org/wiki/Bernstein_polynomial)) 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](https://dev.to/iprosk/experimenting-with-generics-in-rust-little-library-for-bezier-curves-part-1-4093) 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](https://dev.to/iprosk/generic-constant-expressions-a-future-bright-side-of-nightly-rust-3bp7).
- [Experimenting with generics in Rust: little library for Bezier curves - part 1](https://dev.to/iprosk/experimenting-with-generics-in-rust-little-library-for-bezier-curves-part-1-4093).

And the [Github repo](https://github.com/sciprosk/bernstein) 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](https://dev.to/iprosk/generic-constant-expressions-a-future-bright-side-of-nightly-rust-3bp7), 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](https://hackmd.io/OZG_XiLFRs2Xmw5s39jRzA). 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](https://dev.to/iprosk/experimenting-with-generics-in-rust-little-library-for-bezier-curves-part-1-4093). 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](https://en.wikipedia.org/wiki/De_Casteljau%27s_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`](https://docs.rs/num/latest/num/trait.Num.html) 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](https://github.com/sciprosk/bernstein)):
```
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](https://hackmd.io/OZG_XiLFRs2Xmw5s39jRzA)

> 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.