Recursive const generics and specialization

I am writing a square matrix type utilizing const generics.

#[derive(Clone, Debug, Eq, PartialEq, Hash)]
pub struct SquareMatrix<T, const DIM: usize>
where
    [(); DIM.pow(2)]:,
{
    data: [T; DIM.pow(2)],
}

I want to compute the determinant of the matrix. For 1x1, 2x2, 3x3 I use a hardcoded implementation.

pub trait Determinant<T> {
    fn det(&self) -> T;
}

impl<T> Determinant<T> for SquareMatrix<T, 1>
where
    T: Mul<Output = T> + Add<Output = T> + Sub<Output = T> + Neg<Output = T> + Copy + Default,
{
    fn det(&self) -> T {
        self[[0, 0]]
    }
}

impl<T> SquareMatrix<T, 2>
where
    T: Mul<Output = T> + Add<Output = T> + Sub<Output = T> + Neg<Output = T> + Copy + Default,
{
    fn det(&self) -> T {
        self[[0, 0]] * self[[1, 1]] - self[[1, 0]] * self[[0, 1]]
    }
}

// ...

For bigger matrices, I want to (using specialization) use a recursive algorithm that splits up the matrix up into smaller matrices.

impl<T, const DIM: usize> Determinant<T> for SquareMatrix<T, DIM>
where
    [(); DIM.pow(2)]:,
    T: Mul<Output = T> + Add<Output = T> + Sub<Output = T> + Neg<Output = T> + Copy + Default,
    SquareMatrix<T, { DIM - 1 }>: Determinant<T>,
{
    default fn det(&self) -> T {
        let mut val = None;

        for a in 0..DIM {
            let sub_matrix: SquareMatrix<T, { DIM - 1 }> = { 
                // Actually create the submatrix here, but that doesn't seem relevant to my error
                todo!()
            };

            let det = sub_matrix.det();
            if let Some(val) = val.as_mut() {
                *val = *val + det;
            } else {
                val = Some(det);
            }
        }

        val.unwrap()
    }
}
SquareMatrix<T, { DIM - 1 }>: Determinant<T>,
                  ^^^^^^^ attempt to compute `0_usize - 1_usize`, which would overflow

I have tried implement the Determinant trait for SquareMatrix<T, 0> as well, but that didn't seem to help.

Any attempts to remedy have only resulted in different compiler errors (such as Determinant not being implemented for SquareMatrix<T, { DIM - 1 }>.

I am aware that the generic_const_exprs and min_specialization are highly experimental and unfinished. Is this possible to fix?

Sourcecode on the playground

Edit:
Simplified version of the sourcecode without generic type parameters

1 Like

Your link is broken, the gist hash wasn't fully copied.


Edit: Thanks for fixing it :slightly_smiling_face:

1 Like

I don't think the compiler is smart enough to do this kind of thing yet, but here are a few notes regardless:

  • You can use [[T; DIM]; DIM] instead of [T; DIM.pow(2)] to get rid of the pow(2) bounds.
  • The recursive determinant algorithm you have is exponential, so you should probably avoid it anyway. You can use LU decomposition (using Gaussian elimination) instead.
1 Like

It's a shame the compiler isn't there yet, though this does seem to be a very niche usecase.

Using a faster algorithm is deffintly a better idea.

(I don't have a practical useage for this program, I only use it as an excuse to play around with fancy compiler features :wink: )

Thanks!

2 Likes