fn system(mut my_dir: Local<Dir3>) {
    my_dir *= some_quaternion;
    method_that_assumes_the_dir_is_well_normalised(my_dir);
}

/// Performs an approximate and fast normalisation, assuming the object is already close to normalised.
/// Useful for mitigating the accumulation of numerical error.
fn approx_renormalise(&mut self) {
  // We numerically approximate the inverse square root by a taylor series around 1
  // As we expect the error (x) to be small
  // inverse_sqrt(length_squared) = (1 + x)^(-1/2) = 1 - 1/2 x + O(x^2)
  // inverse_sqrt(length_squared) ~~ 1 - 1/2 (length_squared - 1) = 1/2 (3 - length_squared)

  // Iterative calls to this method quickly converge to a normalised value,
  // so long as the denormalisation is not large O(1/10).
  // One iteration can be described as:
  // l_sq <- l_sq * (1 - 1/2 (l_sq - 1))^2;
  // Rewriting in terms of the error x:
  // 1 + x <- (1 + x) * (1 - 1/2 x)^2
  // 1 + x <- (1 + x) * (1 - x + 1/4 x^2)
  // 1 + x <- 1 - x + 1/4 x^2 + x - x^2 + 1/4 x^3
  // x <- -1/4 x^2 (3 - x)
  // if we have a reasonable error, say in a range of (-1/3, 1/3), then:
  // |-1/4 x^2 (3 - x)| <= (3/4 + 1/4 * |x|) * x^2 <= (3/4 + 1/4 * 1/3) * x^2 < x^2 < 1/3 x
  // Therefore the sequence of iterates converges to 0 error, at least geometrically, but practically it is a second order method.

  let length_squared = self.length_square();
  self *= 0.5 * (3 - length_squared);
}