The code is now generic over the dimension. There is a slight performance
loss that I need to investigate; perhaps some things are not getting
specialised? Maybe it is better to be more explicit about staging and
splice in the functions with fixed dimensions.