Minor changes suggested by clippy

src/camera.rs

@@ -16,12 +16,12 @@ struct ImageSampler {
 
     film_width: f64,
     film_height: f64,
-    camera_location: Vec3,
+    camera_location: Vec3<f64>,
     film_distance: f64,
 }
 
 impl ImageSampler {
-    pub fn new(width: usize, height: usize, camera_location: Vec3) -> ImageSampler {
+    pub fn new(width: usize, height: usize, camera_location: Vec3<f64>) -> ImageSampler {
         let (film_width, film_height) = {
             let width = width as f64;
             let height = height as f64;
@@ -173,12 +173,12 @@ mod tests {
                 }) => location,
                 None => panic!(),
             };
-            let expected_x: f64 =
-                ImageSampler::scale(200, 800, target.film_width) - target.film_width * 0.5;
-            assert!((point_on_film_plane.x() - expected_x).abs() < 0.5 / 800.0);
-            let expected_y =
-                -ImageSampler::scale(100, 600, target.film_height) + target.film_height * 0.5;
-            assert!((point_on_film_plane.y() - expected_y).abs() < 0.5 / 600.0);
+            let expected_x = target.film_width * (200.0/800.0 - 0.5);
+            assert!(point_on_film_plane.x() - expected_x < target.film_width / 800.0);
+            assert!(point_on_film_plane.x() - expected_x >= 0.0);
+            let expected_y = -target.film_height * (100.0/600.0 - 0.5);
+            assert!(expected_y - point_on_film_plane.y() < target.film_height / 600.0);
+            assert!(expected_y - point_on_film_plane.y() >= 0.0);
         }
     }
 }

src/colour/colour_rgb.rs

@@ -4,7 +4,7 @@ use std::ops::{Add, Mul};
 
 #[derive(Copy, Clone, Debug, Default)]
 pub struct ColourRgbF {
-    pub values: Vec3,
+    pub values: Vec3<f64>,
 }
 
 impl ColourRgbF {
@@ -34,7 +34,7 @@ impl ColourRgbF {
         }
     }
 
-    pub fn from_vec3(v: &Vec3) -> ColourRgbF {
+    pub fn from_vec3(v: &Vec3<f64>) -> ColourRgbF {
         ColourRgbF { values: *v }
     }
 
@@ -50,7 +50,7 @@ impl ColourRgbF {
         self.values.z()
     }
 
-    pub fn as_vec3(&self) -> &Vec3 {
+    pub fn as_vec3(&self) -> &Vec3<f64> {
         &self.values
     }
 }
@@ -114,7 +114,7 @@ mod tests {
         use quickcheck_macros::quickcheck;
         impl Arbitrary for ColourRgbF {
             fn arbitrary<G: Gen>(g: &mut G) -> ColourRgbF {
-                let values = <Vec3 as Arbitrary>::arbitrary(g);
+                let values = <Vec3<f64> as Arbitrary>::arbitrary(g);
                 ColourRgbF { values }
             }
         }

src/colour/colour_xyz.rs

@@ -5,7 +5,7 @@ use super::{ColourRgbF, Photon};
 /// A CIE XYZ Colour Value
 #[derive(Clone, Debug, Default, PartialEq)]
 pub struct ColourXyz {
-    pub values: Vec3,
+    pub values: Vec3<f64>,
 }
 
 impl ColourXyz {

src/colour/spectrum.rs

@@ -82,6 +82,7 @@ impl Spectrum {
         Spectrum {
             shortest_wavelength: rgb_reference_spectrum::SHORTEST_WAVELENGTH,
             longest_wavelength: rgb_reference_spectrum::LONGEST_WAVELENGTH,
+            #[allow(clippy::collapsible_else_if)]
             samples: if colour.red() <= colour.green() && colour.red() <= colour.blue() {
                 if colour.green() <= colour.blue() {
                     izip![
@@ -178,6 +179,7 @@ impl Spectrum {
 mod rgb_reference_spectrum {
     pub const SHORTEST_WAVELENGTH: f64 = 380.0;
     pub const LONGEST_WAVELENGTH: f64 = 720.0;
+    #[allow(clippy::excessive_precision)]
     pub mod reflection {
         pub const WHITE: [f64; 32] = [
             1.0618958571272863e+00,

src/image.rs

@@ -109,21 +109,21 @@ pub trait NormalizedAsByte {
 
 impl NormalizedAsByte for f32 {
     fn normalized_to_byte(self) -> u8 {
-        (self * (std::u8::MAX as f32)) as u8
+        (self * (u8::MAX as f32)) as u8
     }
 
     fn byte_to_normalized(byte: u8) -> f32 {
-        (byte as f32) / (std::u8::MAX as f32)
+        (byte as f32) / (u8::MAX as f32)
     }
 }
 
 impl NormalizedAsByte for f64 {
     fn normalized_to_byte(self) -> u8 {
-        (self * (std::u8::MAX as f64)) as u8
+        (self * (u8::MAX as f64)) as u8
     }
 
     fn byte_to_normalized(byte: u8) -> f64 {
-        (byte as f64) / (std::u8::MAX as f64)
+        (byte as f64) / (u8::MAX as f64)
     }
 }
 

src/integrators/simple_random_integrator.rs

@@ -54,7 +54,7 @@ impl Integrator for SimpleRandomIntegrator {
     }
 }
 
-pub fn test_lighting_environment(w_o: &Vec3, wavelength: f64) -> f64 {
+pub fn test_lighting_environment(w_o: &Vec3<f64>, wavelength: f64) -> f64 {
     //let sun_direction = Vec3::new(1.0, 1.0, -1.0).normalize();
     //if w_o.dot(&sun_direction) >= 0.99 {
     //    300.0

src/integrators/whitted_integrator.rs

@@ -8,7 +8,7 @@ use crate::util::algebra_utils::try_change_of_basis_matrix;
 use super::Integrator;
 
 pub struct DirectionalLight {
-    pub direction: Vec3,
+    pub direction: Vec3<f64>,
     pub spectrum: Spectrum,
 }
 

src/main.rs

@@ -79,7 +79,7 @@ fn update_texture(image: &ImageRgbU8, texture: &mut Texture) {
         .update(
             Rect::new(0, 0, image.get_width() as u32, image.get_height() as u32),
             image.get_pixel_data(),
-            (image.get_width() * ImageRgbU8::num_channels()) as usize,
+            image.get_width() * ImageRgbU8::num_channels(),
         )
         .expect("Couldn't update texture.");
 }

src/materials/lambertian_material.rs

@@ -25,15 +25,15 @@ impl LambertianMaterial {
 }
 
 impl Material for LambertianMaterial {
-    fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3, &Vec3, &Photon) -> Photon + 'a> {
-        Box::new(move |_w_o: &Vec3, _w_i: &Vec3, photon_in: &Photon| {
+    fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3<f64>, &Vec3<f64>, &Photon) -> Photon + 'a> {
+        Box::new(move |_w_o: &Vec3<f64>, _w_i: &Vec3<f64>, photon_in: &Photon| {
             let mut result = self.colour.scale_photon(photon_in);
             result.intensity *= self.diffuse_strength;
             result
         })
     }
 
-    fn sample(&self, _w_i: &Vec3, _photon: &Photon) -> MaterialSampleResult {
+    fn sample(&self, _w_i: &Vec3<f64>, _photon: &Photon) -> MaterialSampleResult {
         let mut rng = thread_rng();
         let mut w_o = Vec3::new(
             2.0 * rng.sample::<f64, _>(Open01) - 1.0,

src/materials/mod.rs

@@ -18,14 +18,16 @@ pub mod smooth_transparent_dialectric;
 pub use smooth_transparent_dialectric::SmoothTransparentDialectric;
 
 pub struct MaterialSampleResult {
-    pub direction: Vec3,
+    pub direction: Vec3<f64>,
     pub pdf: f64,
 }
 
-pub trait Material: Debug + Sync + Send {
-    fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3, &Vec3, &Photon) -> Photon + 'a>;
+type BsdfFunc<'a> = Box<dyn Fn(&Vec3<f64>, &Vec3<f64>, &Photon) -> Photon + 'a>;
 
-    fn sample(&self, _w_i: &Vec3, _photon: &Photon) -> MaterialSampleResult {
+pub trait Material: Debug + Sync + Send {
+    fn bsdf<'a>(&'a self) -> BsdfFunc<'a>;
+
+    fn sample(&self, _w_i: &Vec3<f64>, _photon: &Photon) -> MaterialSampleResult {
         let distribution = CosineWeightedHemisphere::new();
         let direction = distribution.value();
         let pdf = distribution.pdf(direction);

src/materials/phong_material.rs

@@ -14,8 +14,8 @@ pub struct PhongMaterial {
 }
 
 impl Material for PhongMaterial {
-    fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3, &Vec3, &Photon) -> Photon + 'a> {
-        Box::new(move |w_o: &Vec3, w_i: &Vec3, photon_in: &Photon| {
+    fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3<f64>, &Vec3<f64>, &Photon) -> Photon + 'a> {
+        Box::new(move |w_o: &Vec3<f64>, w_i: &Vec3<f64>, photon_in: &Photon| {
             if w_i.z() < 0.0 || w_o.z() < 0.0 {
                 Photon {
                     wavelength: photon_in.wavelength,

src/materials/reflective_material.rs

@@ -13,8 +13,8 @@ pub struct ReflectiveMaterial {
 }
 
 impl Material for ReflectiveMaterial {
-    fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3, &Vec3, &Photon) -> Photon + 'a> {
-        Box::new(move |w_o: &Vec3, w_i: &Vec3, photon_in: &Photon| {
+    fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3<f64>, &Vec3<f64>, &Photon) -> Photon + 'a> {
+        Box::new(move |w_o: &Vec3<f64>, w_i: &Vec3<f64>, photon_in: &Photon| {
             if w_i.z() <= 0.0 || w_o.z() <= 0.0 {
                 Photon {
                     wavelength: photon_in.wavelength,
@@ -39,7 +39,7 @@ impl Material for ReflectiveMaterial {
         })
     }
 
-    fn sample(&self, w_o: &Vec3, _photon: &Photon) -> MaterialSampleResult {
+    fn sample(&self, w_o: &Vec3<f64>, _photon: &Photon) -> MaterialSampleResult {
         MaterialSampleResult {
             direction: Vec3::new(-w_o.x(), -w_o.y(), w_o.z()),
             pdf: 1.0,

src/materials/smooth_transparent_dialectric.rs

@@ -6,13 +6,13 @@ use rand::random;
 
 #[derive(Debug)]
 struct FresnelResult {
-    reflection_direction: Vec3,
+    reflection_direction: Vec3<f64>,
     reflection_strength: f64,
-    transmission_direction: Vec3,
+    transmission_direction: Vec3<f64>,
     transmission_strength: f64,
 }
 
-fn fresnel(w_i: &Vec3, eta1: f64, eta2: f64) -> FresnelResult {
+fn fresnel(w_i: &Vec3<f64>, eta1: f64, eta2: f64) -> FresnelResult {
     let normal = if w_i.z() > 0.0 {
         Vec3::unit_z()
     } else {
@@ -70,8 +70,8 @@ impl SmoothTransparentDialectric {
 }
 
 impl Material for SmoothTransparentDialectric {
-    fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3, &Vec3, &Photon) -> Photon + 'a> {
-        Box::new(move |w_o: &Vec3, w_i: &Vec3, photon_in: &Photon| {
+    fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3<f64>, &Vec3<f64>, &Photon) -> Photon + 'a> {
+        Box::new(move |w_o: &Vec3<f64>, w_i: &Vec3<f64>, photon_in: &Photon| {
             let (eta1, eta2) = if w_i.z() >= 0.0 {
                 (1.0, self.eta.intensity_at_wavelength(photon_in.wavelength))
             } else {
@@ -88,7 +88,7 @@ impl Material for SmoothTransparentDialectric {
         })
     }
 
-    fn sample(&self, w_i: &Vec3, photon: &Photon) -> MaterialSampleResult {
+    fn sample(&self, w_i: &Vec3<f64>, photon: &Photon) -> MaterialSampleResult {
         let (eta1, eta2) = if w_i.z() >= 0.0 {
             (1.0, self.eta.intensity_at_wavelength(photon.wavelength))
         } else {

src/math/mat2.rs

@@ -1,16 +1,18 @@
+use super::Float;
+
 #[derive(PartialEq, Debug, Copy, Clone)]
-pub struct Mat2 {
-    pub elements: [[f64; 2]; 2],
+pub struct Mat2<T: Float> {
+    pub elements: [[T; 2]; 2],
 }
 
-impl Mat2 {
-    pub fn new(m00: f64, m01: f64, m10: f64, m11: f64) -> Mat2 {
+impl<T: Float> Mat2<T> {
+    pub fn new(m00: T, m01: T, m10: T, m11: T) -> Mat2<T> {
         Mat2 {
             elements: [[m00, m01], [m10, m11]],
         }
     }
 
-    pub fn determinant(&self) -> f64 {
+    pub fn determinant(&self) -> T {
         self.elements[0][0] * self.elements[1][1] - self.elements[0][1] * self.elements[1][0]
     }
 }

src/math/mat3.rs

@@ -1,38 +1,32 @@
-use super::Mat2;
-use super::Vec3;
+use super::{Float, Mat2, Vec3};
 
 use std::ops::{Mul, MulAssign};
 
 #[derive(PartialEq, Debug, Copy, Clone)]
-pub struct Mat3 {
-    elements: [[f64; 3]; 3],
+pub struct Mat3<T: Float> {
+    elements: [[T; 3]; 3],
 }
 
-impl Mat3 {
-    pub fn new(
-        m00: f64,
-        m01: f64,
-        m02: f64,
-        m10: f64,
-        m11: f64,
-        m12: f64,
-        m20: f64,
-        m21: f64,
-        m22: f64,
-    ) -> Mat3 {
+impl<T: Float> Mat3<T> {
+    #[allow(clippy::too_many_arguments)]
+    pub fn new(m00: T, m01: T, m02: T, m10: T, m11: T, m12: T, m20: T, m21: T, m22: T) -> Mat3<T> {
         Mat3 {
             elements: [[m00, m01, m02], [m10, m11, m12], [m20, m21, m22]],
         }
     }
 
-    pub fn identity() -> Mat3 {
+    pub fn identity() -> Mat3<T> {
         Mat3 {
-            elements: [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
+            elements: [
+                [T::one(), T::zero(), T::zero()],
+                [T::zero(), T::one(), T::zero()],
+                [T::zero(), T::zero(), T::one()],
+            ],
         }
     }
 
-    pub fn from_rows(r0: &Vec3, r1: &Vec3, r2: &Vec3) -> Mat3 {
-        let mut elements = [[0.0; 3]; 3];
+    pub fn from_rows(r0: &Vec3<T>, r1: &Vec3<T>, r2: &Vec3<T>) -> Mat3<T> {
+        let mut elements = [[T::zero(); 3]; 3];
         for (row, v) in elements.iter_mut().zip([r0, r1, r2].iter()) {
             for (it, val) in row.iter_mut().zip(v.coords.iter()) {
                 *it = *val;
@@ -41,26 +35,26 @@ impl Mat3 {
         Mat3 { elements }
     }
 
-    pub fn get_element(&self, row: usize, column: usize) -> f64 {
+    pub fn get_element(&self, row: usize, column: usize) -> T {
         self.elements[row][column]
     }
 
-    pub fn get_row(&self, row: usize) -> Vec3 {
+    pub fn get_row(&self, row: usize) -> Vec3<T> {
         Vec3 {
             coords: self.elements[row],
         }
     }
 
-    pub fn get_column(&self, column: usize) -> Vec3 {
-        let mut coords = [0.0; 3];
+    pub fn get_column(&self, column: usize) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
         for (coord, row) in coords.iter_mut().zip(self.elements.iter()) {
             *coord = row[column];
         }
         Vec3 { coords }
     }
 
-    pub fn transpose(&self) -> Mat3 {
-        let mut elements = [[0.0; 3]; 3];
+    pub fn transpose(&self) -> Mat3<T> {
+        let mut elements = [[T::zero(); 3]; 3];
         for i in 0..3 {
             for j in 0..3 {
                 elements[i][j] = self.elements[j][i];
@@ -69,8 +63,8 @@ impl Mat3 {
         Mat3 { elements }
     }
 
-    pub fn first_minor(&self, row: usize, column: usize) -> f64 {
-        let mut elements = [[0.0; 2]; 2];
+    pub fn first_minor(&self, row: usize, column: usize) -> T {
+        let mut elements = [[T::zero(); 2]; 2];
         let mut i_dst = 0;
         let mut j_dst = 0;
         for i_src in 0..3 {
@@ -89,12 +83,12 @@ impl Mat3 {
         minor_matrix.determinant()
     }
 
-    pub fn cofactor(&self, row: usize, column: usize) -> f64 {
-        ((-1i64).pow((row + column) as u32) as f64) * self.first_minor(row, column)
+    pub fn cofactor(&self, row: usize, column: usize) -> T {
+        T::from((-1i32).pow((row + column) as u32)) * self.first_minor(row, column)
     }
 
-    pub fn cofactor_matrix(&self) -> Mat3 {
-        let mut elements = [[0.0; 3]; 3];
+    pub fn cofactor_matrix(&self) -> Mat3<T> {
+        let mut elements = [[T::zero(); 3]; 3];
         for i in 0..3 {
             for j in 0..3 {
                 elements[i][j] = self.cofactor(i, j);
@@ -103,14 +97,14 @@ impl Mat3 {
         Mat3 { elements }
     }
 
-    pub fn determinant(&self) -> f64 {
+    pub fn determinant(&self) -> T {
         self.elements[0][0] * self.first_minor(0, 0) - self.elements[0][1] * self.first_minor(0, 1)
             + self.elements[0][2] * self.first_minor(0, 2)
     }
 
-    pub fn try_inverse(&self) -> Option<Mat3> {
+    pub fn try_inverse(&self) -> Option<Mat3<T>> {
         let determinant = self.determinant();
-        if determinant == 0.0 {
+        if determinant == T::zero() {
             None
         } else {
             Some(self.cofactor_matrix().transpose() * determinant)
@@ -118,11 +112,11 @@ impl Mat3 {
     }
 }
 
-impl Mul<Mat3> for Mat3 {
+impl<T: Float> Mul<Mat3<T>> for Mat3<T> {
     type Output = Self;
 
     fn mul(self, rhs: Self) -> Self {
-        let mut elements = [[0.0; 3]; 3];
+        let mut elements = [[T::zero(); 3]; 3];
         for row in 0..3 {
             for column in 0..3 {
                 elements[row][column] = self.get_row(row).dot(&rhs.get_column(column));
@@ -132,10 +126,10 @@ impl Mul<Mat3> for Mat3 {
     }
 }
 
-impl MulAssign<Mat3> for Mat3 {
+impl<T: Float> MulAssign<Mat3<T>> for Mat3<T> {
     fn mul_assign(&mut self, rhs: Self) {
         for row in 0..3 {
-            let mut new_row = [0.0; 3];
+            let mut new_row = [T::zero(); 3];
             for column in 0..3 {
                 new_row[column] = self.get_row(row).dot(&rhs.get_column(column));
             }
@@ -144,11 +138,11 @@ impl MulAssign<Mat3> for Mat3 {
     }
 }
 
-impl Mul<Vec3> for Mat3 {
-    type Output = Vec3;
+impl<T: Float> Mul<Vec3<T>> for Mat3<T> {
+    type Output = Vec3<T>;
 
-    fn mul(self, rhs: Vec3) -> Vec3 {
-        let mut coords = [0.0; 3];
+    fn mul(self, rhs: Vec3<T>) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
         for (coord, row) in coords.iter_mut().zip(self.elements.iter()) {
             *coord = Vec3 { coords: *row }.dot(&rhs);
         }
@@ -156,23 +150,23 @@ impl Mul<Vec3> for Mat3 {
     }
 }
 
-impl Mul<&Vec3> for Mat3 {
-    type Output = Vec3;
+impl<T: Float> Mul<&Vec3<T>> for Mat3<T> {
+    type Output = Vec3<T>;
 
-    fn mul(self, rhs: &Vec3) -> Vec3 {
-        let mut coords = [0.0; 3];
+    fn mul(self, rhs: &Vec3<T>) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
         for (coord, row) in coords.iter_mut().zip(self.elements.iter()) {
-            *coord = Vec3 { coords: *row }.dot(&rhs);
+            *coord = Vec3 { coords: *row }.dot(rhs);
         }
         Vec3 { coords }
     }
 }
 
-impl Mul<f64> for Mat3 {
-    type Output = Mat3;
+impl<T: Float> Mul<T> for Mat3<T> {
+    type Output = Mat3<T>;
 
-    fn mul(self, rhs: f64) -> Mat3 {
-        let mut elements = [[0.0; 3]; 3];
+    fn mul(self, rhs: T) -> Mat3<T> {
+        let mut elements = [[T::zero(); 3]; 3];
         for i in 0..3 {
             for j in 0..3 {
                 elements[i][j] = self.elements[i][j] * rhs;
@@ -283,7 +277,7 @@ mod tests {
 
     #[test]
     fn inverse_of_identity_is_identity() {
-        assert!(Mat3::identity().try_inverse() == Some(Mat3::identity()));
+        assert!(Mat3::<f64>::identity().try_inverse() == Some(Mat3::identity()));
     }
 
     #[test]

src/math/mat4.rs

@@ -1,31 +1,32 @@
-use super::Vec4;
+use super::{Float,Vec4};
 
 use std::ops::{Mul, MulAssign};
 
 #[derive(PartialEq, Debug)]
-pub struct Mat4 {
-    elements: [[f64; 4]; 4],
+pub struct Mat4<T:Float> {
+    elements: [[T; 4]; 4],
 }
 
-impl Mat4 {
+impl<T:Float> Mat4<T> {
+    #[allow(clippy::too_many_arguments)]
     pub fn new(
-        m00: f64,
-        m01: f64,
-        m02: f64,
-        m03: f64,
-        m10: f64,
-        m11: f64,
-        m12: f64,
-        m13: f64,
-        m20: f64,
-        m21: f64,
-        m22: f64,
-        m23: f64,
-        m30: f64,
-        m31: f64,
-        m32: f64,
-        m33: f64,
-    ) -> Mat4 {
+        m00: T,
+        m01: T,
+        m02: T,
+        m03: T,
+        m10: T,
+        m11: T,
+        m12: T,
+        m13: T,
+        m20: T,
+        m21: T,
+        m22: T,
+        m23: T,
+        m30: T,
+        m31: T,
+        m32: T,
+        m33: T,
+    ) -> Mat4<T> {
         Mat4 {
             elements: [
                 [m00, m01, m02, m03],
@@ -36,8 +37,8 @@ impl Mat4 {
         }
     }
 
-    pub fn from_rows(r0: &Vec4, r1: &Vec4, r2: &Vec4, r3: &Vec4) -> Mat4 {
-        let mut elements = [[0.0; 4]; 4];
+    pub fn from_rows(r0: &Vec4<T>, r1: &Vec4<T>, r2: &Vec4<T>, r3: &Vec4<T>) -> Mat4<T> {
+        let mut elements = [[T::zero(); 4]; 4];
         for (row, v) in elements.iter_mut().zip([r0, r1, r2, r3].iter()) {
             for (it, val) in row.iter_mut().zip(v.coords.iter()) {
                 *it = *val;
@@ -46,18 +47,18 @@ impl Mat4 {
         Mat4 { elements }
     }
 
-    pub fn get_element(&self, row: usize, column: usize) -> f64 {
+    pub fn get_element(&self, row: usize, column: usize) -> T {
         self.elements[row][column]
     }
 
-    pub fn get_row(&self, row: usize) -> Vec4 {
+    pub fn get_row(&self, row: usize) -> Vec4<T> {
         Vec4 {
             coords: self.elements[row],
         }
     }
 
-    pub fn get_column(&self, column: usize) -> Vec4 {
-        let mut coords = [0.0; 4];
+    pub fn get_column(&self, column: usize) -> Vec4<T> {
+        let mut coords = [T::zero(); 4];
         for (coord, row) in coords.iter_mut().zip(self.elements.iter()) {
             *coord = row[column];
         }
@@ -65,11 +66,11 @@ impl Mat4 {
     }
 }
 
-impl Mul<Mat4> for Mat4 {
+impl<T:Float> Mul<Mat4<T>> for Mat4<T> {
     type Output = Self;
 
     fn mul(self, rhs: Self) -> Self {
-        let mut elements = [[0.0; 4]; 4];
+        let mut elements = [[T::zero(); 4]; 4];
         for row in 0..4 {
             for column in 0..4 {
                 elements[row][column] = self.get_row(row).dot(&rhs.get_column(column));
@@ -79,10 +80,10 @@ impl Mul<Mat4> for Mat4 {
     }
 }
 
-impl MulAssign<Mat4> for Mat4 {
+impl<T:Float> MulAssign<Mat4<T>> for Mat4<T> {
     fn mul_assign(&mut self, rhs: Self) {
         for row in 0..4 {
-            let mut new_row = [0.0; 4];
+            let mut new_row = [T::zero(); 4];
             for column in 0..4 {
                 new_row[column] = self.get_row(row).dot(&rhs.get_column(column));
             }
@@ -91,11 +92,11 @@ impl MulAssign<Mat4> for Mat4 {
     }
 }
 
-impl Mul<Vec4> for Mat4 {
-    type Output = Vec4;
+impl<T:Float> Mul<Vec4<T>> for Mat4<T> {
+    type Output = Vec4<T>;
 
-    fn mul(self, rhs: Vec4) -> Vec4 {
-        let mut coords = [0.0; 4];
+    fn mul(self, rhs: Vec4<T>) -> Vec4<T> {
+        let mut coords = [T::zero(); 4];
         for (coord, row) in coords.iter_mut().zip(self.elements.iter()) {
             *coord = Vec4 { coords: *row }.dot(&rhs);
         }

src/math/mod.rs

@@ -1,3 +1,6 @@
+mod number;
+pub use number::*;
+
 mod vec2;
 pub use vec2::*;
 

src/math/number.rs

@@ -0,0 +1,61 @@
+use std::{
+    cmp::PartialOrd,
+    iter::Sum,
+    ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign},
+};
+
+pub trait HasZero {
+    fn zero() -> Self;
+}
+
+pub trait HasOne {
+    fn one() -> Self;
+}
+
+pub trait Float:
+    Copy
+    + HasZero
+    + HasOne
+    + Add<Output = Self>
+    + AddAssign
+    + Div<Output = Self>
+    + Mul<Output = Self>
+    + MulAssign
+    + Sub<Output = Self>
+    + SubAssign
+    + Neg<Output = Self>
+    + Sum
+    + PartialOrd
+    + From<f32>
+    + From<f64>
+    + From<u8>
+    + From<i8>
+    + From<u16>
+    + From<i16>
+    + From<u32>
+    + From<i32>
+{
+    fn abs(self) -> Self;
+    fn sqrt(self) -> Self;
+}
+
+impl HasZero for f64 {
+    fn zero() -> Self {
+        0.0
+    }
+}
+
+impl HasOne for f64 {
+    fn one() -> Self {
+        1.0
+    }
+}
+
+impl Float for f64 {
+    fn abs(self) -> Self {
+        self.abs()
+    }
+    fn sqrt(self) -> Self {
+        self.sqrt()
+    }
+}

src/math/vec2.rs

@@ -1,92 +1,103 @@
+use super::Float;
+
 use itertools::izip;
 
 use std::ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign};
 
 #[derive(PartialEq, Debug, Copy, Clone)]
-pub struct Vec2 {
-    coords: [f64; 2],
+pub struct Vec2<T: Float> {
+    coords: [T; 2],
 }
 
-impl Vec2 {
-    pub fn new(x: f64, y: f64) -> Self {
+impl<T: Float> Vec2<T> {
+    pub fn new(x: T, y: T) -> Self {
         Vec2 { coords: [x, y] }
     }
 
-    pub fn x(&self) -> f64 {
+    pub fn x(&self) -> T {
         self.coords[0]
     }
 
-    pub fn y(&self) -> f64 {
+    pub fn y(&self) -> T {
         self.coords[1]
     }
 
-    pub fn dot(&self, rhs: &Vec2) -> f64 {
+    pub fn dot(&self, rhs: &Vec2<T>) -> T {
         self.coords
             .iter()
-            .zip(rhs.coords.iter())
+            .copied()
+            .zip(rhs.coords.iter().copied())
             .map(|(a_elem, b_elem)| a_elem * b_elem)
             .sum()
     }
 
-    pub fn perp(&self, rhs: &Vec2) -> f64 {
+    pub fn perp(&self, rhs: &Vec2<T>) -> T {
         self.x() * rhs.y() - self.y() * rhs.x()
     }
 }
 
-impl Add for Vec2 {
+impl<T: Float> Add for Vec2<T> {
     type Output = Self;
 
     fn add(self, rhs: Self) -> Self {
-        let mut coords = [0.0; 2];
-        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+        let mut coords = [T::zero(); 2];
+        for (r, a, b) in izip!(
+            coords.iter_mut(),
+            self.coords.iter().copied(),
+            rhs.coords.iter().copied()
+        ) {
             *r = a + b;
         }
         Vec2 { coords }
     }
 }
 
-impl AddAssign for Vec2 {
+impl<T: Float> AddAssign for Vec2<T> {
     fn add_assign(&mut self, rhs: Self) {
-        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter()) {
+        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
             *a += b;
         }
     }
 }
 
-impl Sub for Vec2 {
+impl<T: Float> Sub for Vec2<T> {
     type Output = Self;
 
     fn sub(self, rhs: Self) -> Self {
-        let mut coords = [0.0; 2];
-        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+        let mut coords = [T::zero(); 2];
+        for (r, a, b) in izip!(
+            coords.iter_mut(),
+            self.coords.iter().copied(),
+            rhs.coords.iter().copied()
+        ) {
             *r = a - b;
         }
         Vec2 { coords }
     }
 }
 
-impl SubAssign for Vec2 {
+impl<T: Float> SubAssign for Vec2<T> {
     fn sub_assign(&mut self, rhs: Self) {
-        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter()) {
+        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
             *a -= b;
         }
     }
 }
 
-impl Mul<f64> for Vec2 {
+impl<T: Float> Mul<T> for Vec2<T> {
     type Output = Self;
 
-    fn mul(self, rhs: f64) -> Vec2 {
-        let mut coords = [0.0; 2];
-        for (r, a) in coords.iter_mut().zip(self.coords.iter()) {
+    fn mul(self, rhs: T) -> Vec2<T> {
+        let mut coords = [T::zero(); 2];
+        for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
             *r = a * rhs;
         }
         Vec2 { coords }
     }
 }
 
-impl MulAssign<f64> for Vec2 {
-    fn mul_assign(&mut self, rhs: f64) {
+impl<T: Float> MulAssign<T> for Vec2<T> {
+    fn mul_assign(&mut self, rhs: T) {
         for a in self.coords.iter_mut() {
             *a *= rhs;
         }
@@ -98,8 +109,8 @@ mod tests {
     use super::*;
     use quickcheck::{Arbitrary, Gen};
 
-    impl Arbitrary for Vec2 {
-        fn arbitrary<G: Gen>(g: &mut G) -> Vec2 {
+    impl Arbitrary for Vec2<f64> {
+        fn arbitrary<G: Gen>(g: &mut G) -> Vec2<f64> {
             Vec2::new(f64::arbitrary(g), f64::arbitrary(g))
         }
     }

src/math/vec3.rs

@@ -1,87 +1,75 @@
-use super::Mat3;
+use super::{Float, Mat3};
 
 use itertools::izip;
 
 use std::ops::{Add, AddAssign, Div, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign};
 
 #[derive(Copy, Clone, PartialEq, Debug, Default)]
-pub struct Vec3 {
-    pub coords: [f64; 3],
+pub struct Vec3<T: Float> {
+    pub coords: [T; 3],
 }
 
-impl Vec3 {
-    pub fn new(x: f64, y: f64, z: f64) -> Self {
+impl<T: Float> Vec3<T> {
+    pub fn new(x: T, y: T, z: T) -> Self {
         Vec3 { coords: [x, y, z] }
     }
 
-    pub fn from_slice(v: &[f64]) -> Self {
-        let mut coords = [0.0; 3];
+    pub fn from_slice(v: &[T]) -> Self {
+        let mut coords = [T::zero(); 3];
         coords.clone_from_slice(v);
         Vec3 { coords }
     }
 
-    /*pub fn from_iterator<I>(values: I) -> Vec3
-    where
-        I: Iterator<Item = f64>,
-    {
+    pub fn zeros() -> Vec3<T> {
         Vec3 {
-            coords: [
-                values.next().unwrap(),
-                values.next().unwrap(),
-                values.next().unwrap(),
-            ],
-        }
-    }*/
-
-    pub fn zeros() -> Vec3 {
-        Vec3 {
-            coords: [0.0, 0.0, 0.0],
+            coords: [T::zero(), T::zero(), T::zero()],
         }
     }
 
-    pub fn unit_x() -> Vec3 {
+    pub fn unit_x() -> Vec3<T> {
         Vec3 {
-            coords: [1.0, 0.0, 0.0],
+            coords: [T::one(), T::zero(), T::zero()],
         }
     }
 
-    pub fn unit_y() -> Vec3 {
+    pub fn unit_y() -> Vec3<T> {
         Vec3 {
-            coords: [0.0, 1.0, 0.0],
+            coords: [T::zero(), T::one(), T::zero()],
         }
     }
 
-    pub fn unit_z() -> Vec3 {
+    pub fn unit_z() -> Vec3<T> {
         Vec3 {
-            coords: [0.0, 0.0, 1.0],
+            coords: [T::zero(), T::zero(), T::one()],
         }
     }
 
-    pub fn x(&self) -> f64 {
+    pub fn x(&self) -> T {
         self.coords[0]
     }
 
-    pub fn y(&self) -> f64 {
+    pub fn y(&self) -> T {
         self.coords[1]
     }
 
-    pub fn z(&self) -> f64 {
+    pub fn z(&self) -> T {
         self.coords[2]
     }
 
-    pub fn as_slice(&self) -> &[f64] {
+    pub fn as_slice(&self) -> &[T] {
         &self.coords
     }
 
-    pub fn dot(&self, rhs: &Vec3) -> f64 {
+    pub fn dot(&self, rhs: &Vec3<T>) -> T {
         self.coords
             .iter()
-            .zip(rhs.coords.iter())
+            .copied()
+            .zip(rhs.coords.iter().copied())
             .map(|(a_elem, b_elem)| a_elem * b_elem)
             .sum()
     }
 
-    pub fn cross(&self, rhs: &Vec3) -> Vec3 {
+    pub fn cross(&self, rhs: &Vec3<T>) -> Vec3<T> {
         let x = self.y() * rhs.z() - self.z() * rhs.y();
         let y = self.z() * rhs.x() - self.x() * rhs.z();
         let z = self.x() * rhs.y() - self.y() * rhs.x();
@@ -92,18 +80,18 @@ impl Vec3 {
         Vec3::new(self.x().abs(), self.y().abs(), self.z().abs())
     }
 
-    pub fn norm_squared(&self) -> f64 {
+    pub fn norm_squared(&self) -> T {
         self.dot(self)
     }
 
-    pub fn norm(&self) -> f64 {
+    pub fn norm(&self) -> T {
         self.norm_squared().sqrt()
     }
 
     pub fn normalize(&self) -> Self {
-        let mut coords = [0.0; 3];
-        let inverse_norm = 1.0 / self.norm();
-        for (r, a) in coords.iter_mut().zip(self.coords.iter()) {
+        let mut coords = [T::zero(); 3];
+        let inverse_norm = T::one() / self.norm();
+        for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
             *r = a * inverse_norm;
         }
         Vec3 { coords }
@@ -127,150 +115,172 @@ impl Vec3 {
     }
 
     pub fn component_mul(&self, rhs: &Self) -> Self {
-        let mut coords = [0.0; 3];
-        for (elem, lhs_elem, rhs_elem) in
-            izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter())
-        {
+        let mut coords = [T::zero(); 3];
+        for (elem, lhs_elem, rhs_elem) in izip!(
+            coords.iter_mut(),
+            self.coords.iter().copied(),
+            rhs.coords.iter().copied()
+        ) {
             *elem = lhs_elem * rhs_elem;
         }
         Vec3 { coords }
     }
 }
 
-impl Index<usize> for Vec3 {
-    type Output = f64;
+impl<T: Float> Index<usize> for Vec3<T> {
+    type Output = T;
 
-    fn index(&self, i: usize) -> &f64 {
+    fn index(&self, i: usize) -> &T {
         &self.coords[i]
     }
 }
 
-impl IndexMut<usize> for Vec3 {
-    fn index_mut(&mut self, i: usize) -> &mut f64 {
+impl<T: Float> IndexMut<usize> for Vec3<T> {
+    fn index_mut(&mut self, i: usize) -> &mut T {
         &mut self.coords[i]
     }
 }
 
-impl Add<Vec3> for &Vec3 {
-    type Output = Vec3;
+impl<T: Float> Add<Vec3<T>> for &Vec3<T> {
+    type Output = Vec3<T>;
 
-    fn add(self, rhs: Vec3) -> Vec3 {
-        let mut coords = [0.0; 3];
-        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+    fn add(self, rhs: Vec3<T>) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
+        for (r, a, b) in izip!(
+            coords.iter_mut(),
+            self.coords.iter().copied(),
+            rhs.coords.iter().copied()
+        ) {
             *r = a + b;
         }
         Vec3 { coords }
     }
 }
 
-impl Add<&Vec3> for &Vec3 {
-    type Output = Vec3;
+impl<T: Float> Add<&Vec3<T>> for &Vec3<T> {
+    type Output = Vec3<T>;
 
-    fn add(self, rhs: &Vec3) -> Vec3 {
-        let mut coords = [0.0; 3];
-        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+    fn add(self, rhs: &Vec3<T>) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
+        for (r, a, b) in izip!(
+            coords.iter_mut(),
+            self.coords.iter().copied(),
+            rhs.coords.iter().copied()
+        ) {
             *r = a + b;
         }
         Vec3 { coords }
     }
 }
 
-impl Add for Vec3 {
+impl<T: Float> Add for Vec3<T> {
     type Output = Self;
 
     fn add(self, rhs: Self) -> Self {
-        let mut coords = [0.0; 3];
-        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+        let mut coords = [T::zero(); 3];
+        for (r, a, b) in izip!(
+            coords.iter_mut(),
+            self.coords.iter().copied(),
+            rhs.coords.iter().copied()
+        ) {
             *r = a + b;
         }
         Vec3 { coords }
     }
 }
 
-impl AddAssign for Vec3 {
+impl<T: Float> AddAssign for Vec3<T> {
     fn add_assign(&mut self, rhs: Self) {
-        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter()) {
+        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
             *a += b;
         }
     }
 }
 
-impl Neg for Vec3 {
-    type Output = Vec3;
-    fn neg(self) -> Vec3 {
+impl<T: Float> Neg for Vec3<T> {
+    type Output = Vec3<T>;
+    fn neg(self) -> Vec3<T> {
         Vec3::new(-self.x(), -self.y(), -self.z())
     }
 }
 
-impl Sub for &Vec3 {
-    type Output = Vec3;
+impl<T: Float> Sub for &Vec3<T> {
+    type Output = Vec3<T>;
 
-    fn sub(self, rhs: Self) -> Vec3 {
-        let mut coords = [0.0; 3];
-        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+    fn sub(self, rhs: Self) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
+        for (r, a, b) in izip!(
+            coords.iter_mut(),
+            self.coords.iter().copied(),
+            rhs.coords.iter().copied()
+        ) {
             *r = a - b;
         }
         Vec3 { coords }
     }
 }
 
-impl Sub for Vec3 {
-    type Output = Vec3;
+impl<T: Float> Sub for Vec3<T> {
+    type Output = Vec3<T>;
 
-    fn sub(self, rhs: Self) -> Vec3 {
-        let mut coords = [0.0; 3];
-        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+    fn sub(self, rhs: Self) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
+        for (r, a, b) in izip!(
+            coords.iter_mut(),
+            self.coords.iter().copied(),
+            rhs.coords.iter().copied()
+        ) {
             *r = a - b;
         }
         Vec3 { coords }
     }
 }
 
-impl SubAssign for Vec3 {
+impl<T: Float> SubAssign for Vec3<T> {
     fn sub_assign(&mut self, rhs: Self) {
-        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter()) {
+        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
             *a -= b;
         }
     }
 }
 
-impl Mul<f64> for &Vec3 {
-    type Output = Vec3;
+impl<T: Float> Mul<T> for &Vec3<T> {
+    type Output = Vec3<T>;
 
-    fn mul(self, rhs: f64) -> Vec3 {
-        let mut coords = [0.0; 3];
-        for (r, a) in coords.iter_mut().zip(self.coords.iter()) {
+    fn mul(self, rhs: T) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
+        for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
             *r = a * rhs;
         }
         Vec3 { coords }
     }
 }
 
-impl Mul<f64> for Vec3 {
-    type Output = Vec3;
+impl<T: Float> Mul<T> for Vec3<T> {
+    type Output = Vec3<T>;
 
-    fn mul(self, rhs: f64) -> Vec3 {
-        let mut coords = [0.0; 3];
-        for (r, a) in coords.iter_mut().zip(self.coords.iter()) {
+    fn mul(self, rhs: T) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
+        for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
             *r = a * rhs;
         }
         Vec3 { coords }
     }
 }
 
-impl MulAssign<f64> for Vec3 {
-    fn mul_assign(&mut self, rhs: f64) {
+impl<T: Float> MulAssign<T> for Vec3<T> {
+    fn mul_assign(&mut self, rhs: T) {
         for a in self.coords.iter_mut() {
             *a *= rhs;
         }
     }
 }
 
-impl Mul<Mat3> for &Vec3 {
-    type Output = Vec3;
+impl<T:Float> Mul<Mat3<T>> for &Vec3<T> {
+    type Output = Vec3<T>;
 
-    fn mul(self, rhs: Mat3) -> Vec3 {
-        let mut coords = [0.0; 3];
+    fn mul(self, rhs: Mat3<T>) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
         for i in 0..3 {
             coords[i] = self.dot(&rhs.get_column(i));
         }
@@ -278,11 +288,11 @@ impl Mul<Mat3> for &Vec3 {
     }
 }
 
-impl Mul<Mat3> for Vec3 {
+impl<T:Float> Mul<Mat3<T>> for Vec3<T> {
     type Output = Self;
 
-    fn mul(self, rhs: Mat3) -> Self {
-        let mut coords = [0.0; 3];
+    fn mul(self, rhs: Mat3<T>) -> Self {
+        let mut coords = [T::zero(); 3];
         for i in 0..3 {
             coords[i] = self.dot(&rhs.get_column(i));
         }
@@ -290,9 +300,9 @@ impl Mul<Mat3> for Vec3 {
     }
 }
 
-impl MulAssign<Mat3> for Vec3 {
-    fn mul_assign(&mut self, rhs: Mat3) {
-        let mut coords = [0.0; 3];
+impl<T:Float> MulAssign<Mat3<T>> for Vec3<T> {
+    fn mul_assign(&mut self, rhs: Mat3<T>) {
+        let mut coords = [T::zero(); 3];
         for i in 0..3 {
             coords[i] = self.dot(&rhs.get_column(i));
         }
@@ -300,40 +310,40 @@ impl MulAssign<Mat3> for Vec3 {
     }
 }
 
-impl Mul<Vec3> for f64 {
-    type Output = Vec3;
+impl Mul<Vec3<f64>> for f64 {
+    type Output = Vec3<f64>;
 
-    fn mul(self, rhs: Vec3) -> Vec3 {
+    fn mul(self, rhs: Vec3<f64>) -> Vec3<f64> {
         rhs * self
     }
 }
 
-impl Mul<&Vec3> for f64 {
-    type Output = Vec3;
+impl Mul<&Vec3<f64>> for f64 {
+    type Output = Vec3<f64>;
 
-    fn mul(self, rhs: &Vec3) -> Vec3 {
+    fn mul(self, rhs: &Vec3<f64>) -> Vec3<f64> {
         rhs * self
     }
 }
 
-impl Div<f64> for &Vec3 {
-    type Output = Vec3;
+impl<T: Float> Div<T> for &Vec3<T> {
+    type Output = Vec3<T>;
 
-    fn div(self, rhs: f64) -> Vec3 {
-        let mut coords = [0.0; 3];
-        for (r, a) in coords.iter_mut().zip(self.coords.iter()) {
+    fn div(self, rhs: T) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
+        for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
             *r = a / rhs;
         }
         Vec3 { coords }
     }
 }
 
-impl Div<f64> for Vec3 {
-    type Output = Vec3;
+impl<T: Float> Div<T> for Vec3<T> {
+    type Output = Vec3<T>;
 
-    fn div(self, rhs: f64) -> Vec3 {
-        let mut coords = [0.0; 3];
-        for (r, a) in coords.iter_mut().zip(self.coords.iter()) {
+    fn div(self, rhs: T) -> Vec3<T> {
+        let mut coords = [T::zero(); 3];
+        for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
             *r = a / rhs;
         }
         Vec3 { coords }
@@ -345,9 +355,9 @@ mod tests {
     use super::*;
     use quickcheck::{Arbitrary, Gen};
 
-    impl Arbitrary for Vec3 {
-        fn arbitrary<G: Gen>(g: &mut G) -> Vec3 {
-            Vec3::new(f64::arbitrary(g), f64::arbitrary(g), f64::arbitrary(g))
+    impl<T: Arbitrary + Float> Arbitrary for Vec3<T> {
+        fn arbitrary<G: Gen>(g: &mut G) -> Vec3<T> {
+            Vec3::new(T::arbitrary(g), T::arbitrary(g), T::arbitrary(g))
         }
     }
 

src/math/vec4.rs

@@ -1,98 +1,109 @@
+use super::Float;
+
 use itertools::izip;
 
 use std::ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign};
 
 #[derive(PartialEq, Debug)]
-pub struct Vec4 {
-    pub coords: [f64; 4],
+pub struct Vec4<T: Float> {
+    pub coords: [T; 4],
 }
 
-impl Vec4 {
-    pub fn new(x: f64, y: f64, z: f64, w: f64) -> Self {
+impl<T: Float> Vec4<T> {
+    pub fn new(x: T, y: T, z: T, w: T) -> Self {
         Vec4 {
             coords: [x, y, z, w],
         }
     }
 
-    pub fn x(&self) -> f64 {
+    pub fn x(&self) -> T {
         self.coords[0]
     }
 
-    pub fn y(&self) -> f64 {
+    pub fn y(&self) -> T {
         self.coords[1]
     }
 
-    pub fn z(&self) -> f64 {
+    pub fn z(&self) -> T {
         self.coords[2]
     }
 
-    pub fn w(&self) -> f64 {
+    pub fn w(&self) -> T {
         self.coords[3]
     }
 
-    pub fn dot(&self, rhs: &Vec4) -> f64 {
+    pub fn dot(&self, rhs: &Vec4<T>) -> T {
         self.coords
             .iter()
-            .zip(rhs.coords.iter())
+            .copied()
+            .zip(rhs.coords.iter().copied())
             .map(|(a_elem, b_elem)| a_elem * b_elem)
             .sum()
     }
 }
 
-impl Add for Vec4 {
+impl<T: Float> Add for Vec4<T> {
     type Output = Self;
 
     fn add(self, rhs: Self) -> Self {
-        let mut coords = [0.0; 4];
-        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+        let mut coords = [T::zero(); 4];
+        for (r, a, b) in izip!(
+            coords.iter_mut(),
+            self.coords.iter().copied(),
+            rhs.coords.iter().copied()
+        ) {
             *r = a + b;
         }
         Vec4 { coords }
     }
 }
 
-impl AddAssign for Vec4 {
+impl<T: Float> AddAssign for Vec4<T> {
     fn add_assign(&mut self, rhs: Self) {
-        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter()) {
+        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
             *a += b;
         }
     }
 }
 
-impl Sub for Vec4 {
+impl<T: Float> Sub for Vec4<T> {
     type Output = Self;
 
     fn sub(self, rhs: Self) -> Self {
-        let mut coords = [0.0; 4];
-        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+        let mut coords = [T::zero(); 4];
+        for (r, a, b) in izip!(
+            coords.iter_mut(),
+            self.coords.iter().copied(),
+            rhs.coords.iter().copied()
+        ) {
             *r = a - b;
         }
         Vec4 { coords }
     }
 }
 
-impl SubAssign for Vec4 {
+impl<T: Float> SubAssign for Vec4<T> {
     fn sub_assign(&mut self, rhs: Self) {
-        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter()) {
+        for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
             *a -= b;
         }
     }
 }
 
-impl Mul<f64> for Vec4 {
+impl<T: Float> Mul<T> for Vec4<T> {
     type Output = Self;
 
-    fn mul(self, rhs: f64) -> Vec4 {
-        let mut coords = [0.0; 4];
-        for (r, a) in coords.iter_mut().zip(self.coords.iter()) {
+    fn mul(self, rhs: T) -> Vec4<T> {
+        let mut coords = [T::zero(); 4];
+        for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
             *r = a * rhs;
         }
         Vec4 { coords }
     }
 }
 
-impl MulAssign<f64> for Vec4 {
-    fn mul_assign(&mut self, rhs: f64) {
+impl<T: Float> MulAssign<T> for Vec4<T> {
+    fn mul_assign(&mut self, rhs: T) {
         for a in self.coords.iter_mut() {
             *a *= rhs;
         }

src/mesh.rs

@@ -14,7 +14,7 @@ mod wavefront_obj {
         index_tuple: &IndexTuple,
         vertex_positions: &[[f32; 3]],
         normal_positions: &[[f32; 3]],
-    ) -> (Vec3, Vec3) {
+    ) -> (Vec3<f64>, Vec3<f64>) {
         let &IndexTuple(vertex_index, _, maybe_normal_index) = index_tuple;
         (
             {
@@ -47,7 +47,7 @@ mod wavefront_obj {
     ) -> Vec<Triangle> {
         if let Some(v0_index) = polygon.iter().next() {
             let (v0_vertex, v0_normal) =
-                get_vertex_and_normal(v0_index, &vertex_positions, normal_positions);
+                get_vertex_and_normal(v0_index, vertex_positions, normal_positions);
             polygon
                 .iter()
                 .skip(1)

src/random_distributions/cosine_weighted_hemisphere.rs

@@ -16,8 +16,8 @@ impl CosineWeightedHemisphere {
     }
 }
 
-impl RandomDistribution<Vec3> for CosineWeightedHemisphere {
-    fn value(&self) -> Vec3 {
+impl RandomDistribution<Vec3<f64>> for CosineWeightedHemisphere {
+    fn value(&self) -> Vec3<f64> {
         let point_on_disc = self.unit_disc.value();
         let z = 0.0f64
             .max(
@@ -27,7 +27,7 @@ impl RandomDistribution<Vec3> for CosineWeightedHemisphere {
         Vec3::new(point_on_disc.x(), point_on_disc.y(), z)
     }
 
-    fn pdf(&self, v: Vec3) -> f64 {
+    fn pdf(&self, v: Vec3<f64>) -> f64 {
         (v.x() * v.x() + v.y() * v.y()).sqrt() / PI
     }
 }

src/random_distributions/sky_light_pdf.rs

@@ -24,8 +24,8 @@ impl Default for SkyLightPdf {
     }
 }
 
-impl RandomDistribution<Vec3> for SkyLightPdf {
-    fn value(&self) -> Vec3 {
+impl RandomDistribution<Vec3<f64>> for SkyLightPdf {
+    fn value(&self) -> Vec3<f64> {
         let mut rng = thread_rng();
         let phi = rng.sample::<f64, _>(Open01) * 2.0 * PI;
         let z = self.z_distribution.value();
@@ -33,7 +33,7 @@ impl RandomDistribution<Vec3> for SkyLightPdf {
         Vec3::new(r * phi.cos(), r * phi.sin(), z)
     }
 
-    fn pdf(&self, value: Vec3) -> f64 {
+    fn pdf(&self, value: Vec3<f64>) -> f64 {
         let z = value.z();
         if z < 0.0 {
             0.0

src/random_distributions/uniform_hemisphere.rs

@@ -16,8 +16,8 @@ impl UniformHemisphere {
     }
 }
 
-impl RandomDistribution<Vec3> for UniformHemisphere {
-    fn value(&self) -> Vec3 {
+impl RandomDistribution<Vec3<f64>> for UniformHemisphere {
+    fn value(&self) -> Vec3<f64> {
         let mut rng = thread_rng();
         let mut result = Vec3::new(
             2.0 * rng.sample::<f64, _>(Open01) - 1.0,
@@ -34,7 +34,7 @@ impl RandomDistribution<Vec3> for UniformHemisphere {
         result.normalize()
     }
 
-    fn pdf(&self, _: Vec3) -> f64 {
+    fn pdf(&self, _: Vec3<f64>) -> f64 {
         1.0 / (2.0 * PI)
     }
 }

src/random_distributions/uniform_square.rs

@@ -7,24 +7,24 @@ use super::RandomDistribution;
 
 #[derive(Debug)]
 pub struct UniformSquare {
-    corner: Vec2,
+    corner: Vec2<f64>,
     size: f64,
 }
 
 impl UniformSquare {
-    pub fn new(corner: Vec2, size: f64) -> UniformSquare {
+    pub fn new(corner: Vec2<f64>, size: f64) -> UniformSquare {
         UniformSquare { corner, size }
     }
 }
 
-impl RandomDistribution<Vec2> for UniformSquare {
-    fn value(&self) -> Vec2 {
+impl RandomDistribution<Vec2<f64>> for UniformSquare {
+    fn value(&self) -> Vec2<f64> {
         let mut rng = thread_rng();
         self.corner
             + Vec2::new(rng.sample::<f64, _>(Open01), rng.sample::<f64, _>(Open01)) * self.size
     }
 
-    fn pdf(&self, _value: Vec2) -> f64 {
+    fn pdf(&self, _value: Vec2<f64>) -> f64 {
         1.0 / (self.size * self.size)
     }
 }

src/random_distributions/unit_disc.rs

@@ -24,8 +24,8 @@ impl UnitDisc {
     }
 }
 
-impl RandomDistribution<Vec2> for UnitDisc {
-    fn value(&self) -> Vec2 {
+impl RandomDistribution<Vec2<f64>> for UnitDisc {
+    fn value(&self) -> Vec2<f64> {
         let offset = self.square_distribution.value();
         if offset.x() == 0.0 && offset.y() == 0.0 {
             offset
@@ -39,7 +39,7 @@ impl RandomDistribution<Vec2> for UnitDisc {
         }
     }
 
-    fn pdf(&self, _: Vec2) -> f64 {
+    fn pdf(&self, _: Vec2<f64>) -> f64 {
         1.0 / PI
     }
 }

src/raycasting/axis_aligned_bounding_box.rs

@@ -48,7 +48,7 @@ mod tests {
         interval.contains_value(wrap_value_in_interval(v, interval))
     }
 
-    fn wrap_point_into_bounding_box(point: Vec3, bounds: &BoundingBox) -> Vec3 {
+    fn wrap_point_into_bounding_box(point: Vec3<f64>, bounds: &BoundingBox) -> Vec3<f64> {
         let mut coords = [0.0; 3];
         for i in 0..3 {
             coords[i] = wrap_value_in_interval(point[i], bounds.bounds[i]);
@@ -58,10 +58,10 @@ mod tests {
 
     #[quickcheck]
     fn correctly_detects_intersections(
-        ray_origin: Vec3,
-        corner1: Vec3,
-        corner2: Vec3,
-        random_point: Vec3,
+        ray_origin: Vec3<f64>,
+        corner1: Vec3<f64>,
+        corner2: Vec3<f64>,
+        random_point: Vec3<f64>,
     ) -> bool {
         let bounds = BoundingBox::from_corners(corner1, corner2);
         let point_in_bounds = wrap_point_into_bounding_box(random_point, &bounds);
@@ -71,10 +71,10 @@ mod tests {
 
     #[quickcheck]
     fn intersect_always_true_when_ray_origin_is_inside_bounds(
-        ray_origin: Vec3,
-        corner1: Vec3,
-        corner2: Vec3,
-        random_point: Vec3,
+        ray_origin: Vec3<f64>,
+        corner1: Vec3<f64>,
+        corner2: Vec3<f64>,
+        random_point: Vec3<f64>,
     ) -> TestResult {
         let bounds = BoundingBox::from_corners(corner1, corner2);
         let ray_origin = wrap_point_into_bounding_box(ray_origin, &bounds);
@@ -84,10 +84,10 @@ mod tests {
 
     #[quickcheck]
     fn no_intersection_when_behind_ray(
-        ray_origin: Vec3,
-        corner1: Vec3,
-        corner2: Vec3,
-        random_point: Vec3,
+        ray_origin: Vec3<f64>,
+        corner1: Vec3<f64>,
+        corner2: Vec3<f64>,
+        random_point: Vec3<f64>,
     ) -> TestResult {
         let bounds = BoundingBox::from_corners(corner1, corner2);
         if bounds.contains_point(ray_origin) {

src/raycasting/bounding_volume_hierarchy.rs

@@ -27,7 +27,7 @@ pub enum BoundingVolumeHierarchy {
     },
 }
 
-fn centre(bounds: &BoundingBox) -> Vec3 {
+fn centre(bounds: &BoundingBox) -> Vec3<f64> {
     Vec3::new(
         (bounds.bounds[0].get_min() + bounds.bounds[0].get_max()) / 2.00,
         (bounds.bounds[1].get_min() + bounds.bounds[1].get_max()) / 2.0,

src/raycasting/mod.rs

@@ -28,17 +28,17 @@ pub mod vec_aggregate;
 #[derive(Clone, Debug)]
 pub struct Ray {
     /// The start point of the ray
-    pub origin: Vec3,
+    pub origin: Vec3<f64>,
 
     /// The direction the ray goes in.
     ///
     /// This vector should always be kept normalized
-    pub direction: Vec3,
+    pub direction: Vec3<f64>,
 }
 
 impl Ray {
     /// Create a new ray
-    pub fn new(origin: Vec3, direction: Vec3) -> Ray {
+    pub fn new(origin: Vec3<f64>, direction: Vec3<f64>) -> Ray {
         Ray {
             origin,
             direction: direction.normalize(),
@@ -46,7 +46,7 @@ impl Ray {
     }
 
     /// Return the point on the ray that is `t` units from the start
-    pub fn point_at(&self, t: f64) -> Vec3 {
+    pub fn point_at(&self, t: f64) -> Vec3<f64> {
         self.origin + self.direction * t
     }
 
@@ -70,26 +70,26 @@ pub struct IntersectionInfo {
     pub distance: f64,
 
     /// The intersection point
-    pub location: Vec3,
+    pub location: Vec3<f64>,
 
     /// The surface normal at the intersection point
-    pub normal: Vec3,
+    pub normal: Vec3<f64>,
 
     /// The surface tangent at the intersection point
     ///
     /// Which surface tangent direction returned is dependent on the [Primitive](Primitive)
     /// but should generally be smooth over any given surface
-    pub tangent: Vec3,
+    pub tangent: Vec3<f64>,
 
     /// Another surface tangent, perpendicular to `tangent`
     ///
     /// The cross product or `normal` and `tangent`
-    pub cotangent: Vec3,
+    pub cotangent: Vec3<f64>,
 
     /// The direction from the intersection point back towards the ray
     ///
     /// Equal to `-ray.direction`
-    pub retro: Vec3,
+    pub retro: Vec3<f64>,
 
     /// The [Material](crate::materials::Material) which describes the optical
     /// properties of the intersected surface
@@ -146,8 +146,8 @@ mod tests {
     use quickcheck::{Arbitrary, Gen};
     impl Arbitrary for Ray {
         fn arbitrary<G: Gen>(g: &mut G) -> Ray {
-            let origin = <Vec3 as Arbitrary>::arbitrary(g);
-            let direction = <Vec3 as Arbitrary>::arbitrary(g);
+            let origin = <Vec3<f64> as Arbitrary>::arbitrary(g);
+            let direction = <Vec3<f64> as Arbitrary>::arbitrary(g);
             return Ray::new(origin, direction);
         }
     }

src/raycasting/plane.rs

@@ -7,15 +7,15 @@ use std::sync::Arc;
 
 #[derive(Clone)]
 pub struct Plane {
-    normal: Vec3,
-    tangent: Vec3,
-    cotangent: Vec3,
+    normal: Vec3<f64>,
+    tangent: Vec3<f64>,
+    cotangent: Vec3<f64>,
     distance_from_origin: f64,
     material: Arc<dyn Material>,
 }
 
 impl Plane {
-    pub fn new(normal: Vec3, distance_from_origin: f64, material: Arc<dyn Material>) -> Plane {
+    pub fn new(normal: Vec3<f64>, distance_from_origin: f64, material: Arc<dyn Material>) -> Plane {
         let normal = normal.normalize();
         let mut axis_closest_to_tangent = Vec3::zeros();
         axis_closest_to_tangent[normal.smallest_coord()] = 1.0;
@@ -77,22 +77,22 @@ impl Intersect for Plane {
 impl HasBoundingBox for Plane {
     fn bounding_box(&self) -> BoundingBox {
         let p0 = self.normal * self.distance_from_origin;
-        let f = |v: Vec3| {
+        let f = |v: Vec3<f64>| {
             Vec3::new(
                 if v.x() == 0.0 {
                     0.0
                 } else {
-                    std::f64::INFINITY
+                    f64::INFINITY
                 },
                 if v.y() == 0.0 {
                     0.0
                 } else {
-                    std::f64::INFINITY
+                    f64::INFINITY
                 },
                 if v.z() == 0.0 {
                     0.0
                 } else {
-                    std::f64::INFINITY
+                    f64::INFINITY
                 },
             )
         };

src/raycasting/sphere.rs

@@ -7,13 +7,13 @@ use std::sync::Arc;
 
 #[derive(Clone, Debug)]
 pub struct Sphere {
-    centre: Vec3,
+    centre: Vec3<f64>,
     radius: f64,
     material: Arc<dyn Material>,
 }
 
 impl Sphere {
-    pub fn new(centre: Vec3, radius: f64, material: Arc<dyn Material>) -> Sphere {
+    pub fn new(centre: Vec3<f64>, radius: f64, material: Arc<dyn Material>) -> Sphere {
         Sphere {
             centre,
             radius,
@@ -163,8 +163,8 @@ mod tests {
 
     #[quickcheck]
     fn ray_intersects_sphere_centre_at_correct_distance(
-        ray_origin: Vec3,
-        sphere_centre: Vec3,
+        ray_origin: Vec3<f64>,
+        sphere_centre: Vec3<f64>,
         radius: f64,
     ) -> TestResult {
         if radius <= 0.0 || radius + 0.000001 >= (ray_origin - sphere_centre).norm() {
@@ -184,7 +184,7 @@ mod tests {
     }
 
     #[quickcheck]
-    fn all_points_on_sphere_are_in_bounding_box(sphere_centre: Vec3, radius_vector: Vec3) -> bool {
+    fn all_points_on_sphere_are_in_bounding_box(sphere_centre: Vec3<f64>, radius_vector: Vec3<f64>) -> bool {
         let target_sphere = Sphere::new(
             sphere_centre,
             radius_vector.norm(),
@@ -196,9 +196,9 @@ mod tests {
 
     /*#[quickcheck]
     fn translation_moves_centre(
-        sphere_centre: Vec3,
+        sphere_centre: Vec3<f64>,
         radius: f64,
-        translation_vector: Vec3,
+        translation_vector: Vec3<f64>,
     ) -> TestResult {
         if radius <= 0.0 {
             return TestResult::discard();
@@ -217,9 +217,9 @@ mod tests {
 
     #[quickcheck]
     fn translation_does_not_change_radius(
-        sphere_centre: Vec3,
+        sphere_centre: Vec3<f64>,
         radius: f64,
-        translation_vector: Vec3,
+        translation_vector: Vec3<f64>,
     ) -> TestResult {
         if radius <= 0.0 {
             return TestResult::discard();
@@ -238,9 +238,9 @@ mod tests {
 
     #[quickcheck]
     fn rotation_about_centre_does_not_move_centre(
-        sphere_centre: Vec3,
+        sphere_centre: Vec3<f64>,
         radius: f64,
-        rotation_vector: Vec3,
+        rotation_vector: Vec3<f64>,
     ) -> TestResult {
         if radius <= 0.0 {
             return TestResult::discard();

src/raycasting/triangle.rs

@@ -7,8 +7,8 @@ use std::sync::Arc;
 
 #[derive(Debug, Clone)]
 pub struct Triangle {
-    pub vertices: [Vec3; 3],
-    pub normals: [Vec3; 3],
+    pub vertices: [Vec3<f64>; 3],
+    pub normals: [Vec3<f64>; 3],
     pub material: Arc<dyn Material>,
 }
 
@@ -38,7 +38,7 @@ impl Intersect for Triangle {
         let indices = indices_with_index_of_largest_element_last(&ray.direction);
         let permuted_ray_direction = permute_vector_elements(&ray.direction, &indices);
         let shear_slopes = calculate_shear_to_z_axis(&permuted_ray_direction);
-        let transformed_vertices: Vec<Vec3> = self
+        let transformed_vertices: Vec<Vec3<f64>> = self
             .vertices
             .iter()
             .map(|elem| {
@@ -70,7 +70,7 @@ impl Intersect for Triangle {
                 .map(|(&barycentric_coord, vertex)| vertex * barycentric_coord)
                 .fold(Vec3::zeros(), |a, e| a + e);
             let distance = (ray.origin - location).norm();
-            let normal: Vec3 = barycentric_coordinates
+            let normal: Vec3<f64> = barycentric_coordinates
                 .coords
                 .iter()
                 .zip(self.normals.iter())
@@ -105,7 +105,8 @@ impl HasBoundingBox for Triangle {
 
 impl Primitive for Triangle {}
 
-fn indices_with_index_of_largest_element_last(v: &Vec3) -> [usize; 3] {
+fn indices_with_index_of_largest_element_last(v: &Vec3<f64>) -> [usize; 3] {
+    #[allow(clippy::collapsible_else_if)]
     if v.x() > v.y() {
         if v.z() > v.x() {
             [0, 1, 2]
@@ -122,27 +123,27 @@ fn indices_with_index_of_largest_element_last(v: &Vec3) -> [usize; 3] {
 }
 
 fn is_valid_permutation(indices: &[usize; 3]) -> bool {
-    (0..2).all(|i: usize| indices.iter().any(|&j| j == i))
+    (0..2).all(|i: usize| indices.contains(&i))
 }
 
-fn permute_vector_elements(v: &Vec3, indices: &[usize; 3]) -> Vec3 {
+fn permute_vector_elements(v: &Vec3<f64>, indices: &[usize; 3]) -> Vec3<f64> {
     debug_assert!(is_valid_permutation(indices));
     Vec3::new(v[indices[0]], v[indices[1]], v[indices[2]])
 }
 
-fn calculate_shear_to_z_axis(v: &Vec3) -> Vec2 {
+fn calculate_shear_to_z_axis(v: &Vec3<f64>) -> Vec2<f64> {
     Vec2::new(-v.x() / v.z(), -v.y() / v.z())
 }
 
-fn apply_shear_to_z_axis(v: &Vec3, s: &Vec2) -> Vec3 {
+fn apply_shear_to_z_axis(v: &Vec3<f64>, s: &Vec2<f64>) -> Vec3<f64> {
     Vec3::new(v.x() + s.x() * v.z(), v.y() + s.y() * v.z(), v.z())
 }
 
-fn signed_edge_function(a: &Vec3, b: &Vec3) -> f64 {
+fn signed_edge_function(a: &Vec3<f64>, b: &Vec3<f64>) -> f64 {
     a.x() * b.y() - b.x() * a.y()
 }
 
-fn signed_edge_functions(vertices: &[Vec3]) -> Vec3 {
+fn signed_edge_functions(vertices: &[Vec3<f64>]) -> Vec3<f64> {
     // Iterate over the inputs in such a way that each output element is calculated
     // from the twoother elements of the input. ( (y,z) -> x, (z,x) -> y, (x,y) -> z )
     let coords: Vec<_> = vertices
@@ -157,7 +158,7 @@ fn signed_edge_functions(vertices: &[Vec3]) -> Vec3 {
     Vec3::new(coords[0], coords[1], coords[2])
 }
 
-fn barycentric_coordinates_from_signed_edge_functions(e: Vec3) -> Vec3 {
+fn barycentric_coordinates_from_signed_edge_functions(e: Vec3<f64>) -> Vec3<f64> {
     e * (1.0 / e.coords.iter().fold(0.0, |a, &b| a + b))
 }
 
@@ -251,43 +252,43 @@ mod tests {
         use quickcheck_macros::quickcheck;
 
         #[quickcheck]
-        fn result_is_valid_permutation(v: Vec3) -> bool {
+        fn result_is_valid_permutation(v: Vec3<f64>) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
             is_valid_permutation(&indices)
         }
 
         #[quickcheck]
-        fn result_includes_x(v: Vec3) -> bool {
+        fn result_includes_x(v: Vec3<f64>) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
             indices.iter().any(|&i| i == 0)
         }
 
         #[quickcheck]
-        fn result_includes_y(v: Vec3) -> bool {
+        fn result_includes_y(v: Vec3<f64>) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
             indices.iter().any(|&i| i == 1)
         }
 
         #[quickcheck]
-        fn result_includes_z(v: Vec3) -> bool {
+        fn result_includes_z(v: Vec3<f64>) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
             indices.iter().any(|&i| i == 2)
         }
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_x(v: Vec3) -> bool {
+        fn last_index_is_greater_than_or_equal_to_x(v: Vec3<f64>) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
             v[indices[2]] >= v.x()
         }
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_y(v: Vec3) -> bool {
+        fn last_index_is_greater_than_or_equal_to_y(v: Vec3<f64>) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
             v[indices[2]] >= v.y()
         }
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_z(v: Vec3) -> bool {
+        fn last_index_is_greater_than_or_equal_to_z(v: Vec3<f64>) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
             v[indices[2]] >= v.z()
         }
@@ -298,19 +299,19 @@ mod tests {
         use quickcheck_macros::quickcheck;
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_x(v: Vec3) -> bool {
+        fn last_index_is_greater_than_or_equal_to_x(v: Vec3<f64>) -> bool {
             let p = permute_vector_elements(&v, &indices_with_index_of_largest_element_last(&v));
             p.z() >= v.x()
         }
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_y(v: Vec3) -> bool {
+        fn last_index_is_greater_than_or_equal_to_y(v: Vec3<f64>) -> bool {
             let p = permute_vector_elements(&v, &indices_with_index_of_largest_element_last(&v));
             p.z() >= v.y()
         }
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_z(v: Vec3) -> bool {
+        fn last_index_is_greater_than_or_equal_to_z(v: Vec3<f64>) -> bool {
             let p = permute_vector_elements(&v, &indices_with_index_of_largest_element_last(&v));
             p.z() >= v.z()
         }
@@ -321,19 +322,19 @@ mod tests {
         use quickcheck_macros::quickcheck;
 
         #[quickcheck]
-        fn shear_to_z_axis_makes_x_zero(v: Vec3) -> bool {
+        fn shear_to_z_axis_makes_x_zero(v: Vec3<f64>) -> bool {
             let s = calculate_shear_to_z_axis(&v);
             apply_shear_to_z_axis(&v, &s).x().abs() < 0.00001
         }
 
         #[quickcheck]
-        fn shear_to_z_axis_makes_y_zero(v: Vec3) -> bool {
+        fn shear_to_z_axis_makes_y_zero(v: Vec3<f64>) -> bool {
             let s = calculate_shear_to_z_axis(&v);
             apply_shear_to_z_axis(&v, &s).y().abs() < 0.00001
         }
 
         #[quickcheck]
-        fn shear_to_z_axis_leaves_z_unchanged(v: Vec3) -> bool {
+        fn shear_to_z_axis_leaves_z_unchanged(v: Vec3<f64>) -> bool {
             let s = calculate_shear_to_z_axis(&v);
             apply_shear_to_z_axis(&v, &s).z() == v.z()
         }
@@ -345,7 +346,7 @@ mod tests {
         use quickcheck_macros::quickcheck;
 
         #[quickcheck]
-        fn sign_of_signed_edge_function_matches_winding(a: Vec3, b: Vec3) -> TestResult {
+        fn sign_of_signed_edge_function_matches_winding(a: Vec3<f64>, b: Vec3<f64>) -> TestResult {
             let a_2d = Vec2::new(a.x(), a.y());
             let b_2d = Vec2::new(b.x(), b.y());
             let c_2d = Vec2::new(0.0, 0.0);
@@ -361,9 +362,9 @@ mod tests {
 
         #[quickcheck]
         fn signed_edge_functions_has_same_result_as_signed_edge_function(
-            a: Vec3,
-            b: Vec3,
-            c: Vec3,
+            a: Vec3<f64>,
+            b: Vec3<f64>,
+            c: Vec3<f64>,
         ) -> bool {
             let es = signed_edge_functions(&vec![a, b, c]);
             es[0] == signed_edge_function(&b, &c)
@@ -372,7 +373,7 @@ mod tests {
         }
 
         #[quickcheck]
-        fn barycentric_coordinates_sum_to_one(a: Vec3, b: Vec3, c: Vec3) -> bool {
+        fn barycentric_coordinates_sum_to_one(a: Vec3<f64>, b: Vec3<f64>, c: Vec3<f64>) -> bool {
             let barycentric_coordinates =
                 barycentric_coordinates_from_signed_edge_functions(signed_edge_functions(&vec![
                     a, b, c,
@@ -496,15 +497,15 @@ mod tests {
         }
 
         fn intersect_with_centroid_and_test_result<
-            F: Fn(Option<IntersectionInfo>, Vec3) -> bool,
+            F: Fn(Option<IntersectionInfo>, Vec3<f64>) -> bool,
         >(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
             test: F,
         ) -> TestResult {
-            let centroid: Vec3 = [vertex0, vertex1, vertex2]
+            let centroid: Vec3<f64> = [vertex0, vertex1, vertex2]
                 .iter()
                 .fold(Vec3::new(0.0, 0.0, 0.0), |acc, &elem| acc + elem)
                 * (1.0 / 3.0);
@@ -530,12 +531,12 @@ mod tests {
 
         #[quickcheck]
         fn intersection_with_centroid_hits(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
         ) -> TestResult {
-            let centroid: Vec3 = [vertex0, vertex1, vertex2]
+            let centroid: Vec3<f64> = [vertex0, vertex1, vertex2]
                 .iter()
                 .fold(Vec3::new(0.0, 0.0, 0.0), |acc, &elem| acc + elem)
                 * (1.0 / 3.0);
@@ -565,10 +566,10 @@ mod tests {
 
         #[quickcheck]
         fn intersection_with_centroid_hits_centroid(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
         ) -> TestResult {
             intersect_with_centroid_and_test_result(
                 vertex0,
@@ -587,10 +588,10 @@ mod tests {
 
         #[quickcheck]
         fn intersection_with_centroid_hits_at_expected_distance(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
         ) -> TestResult {
             intersect_with_centroid_and_test_result(
                 vertex0,
@@ -609,10 +610,10 @@ mod tests {
 
         #[quickcheck]
         fn intersection_with_centroid_has_expected_normal(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
         ) -> TestResult {
             intersect_with_centroid_and_test_result(
                 vertex0,
@@ -633,10 +634,10 @@ mod tests {
 
         #[quickcheck]
         fn intersection_with_centroid_has_expected_retro(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
         ) -> TestResult {
             intersect_with_centroid_and_test_result(
                 vertex0,
@@ -672,12 +673,12 @@ mod tests {
         }
 
         fn intersect_with_barycentric_and_test_result<
-            F: Fn(Option<IntersectionInfo>, Vec3) -> bool,
+            F: Fn(Option<IntersectionInfo>, Vec3<f64>) -> bool,
         >(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
             barycentric_coords: BarycentricCoords,
             test: F,
         ) -> TestResult {
@@ -706,10 +707,10 @@ mod tests {
 
         #[quickcheck]
         fn point_with_arbitrary_barycentric_coords_hits(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
             barycentric_coords: BarycentricCoords,
         ) -> TestResult {
             intersect_with_barycentric_and_test_result(
@@ -730,10 +731,10 @@ mod tests {
 
         #[quickcheck]
         fn point_with_arbitrary_barycentric_coords_has_expected_normal(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
             barycentric_coords: BarycentricCoords,
         ) -> TestResult {
             intersect_with_barycentric_and_test_result(
@@ -756,10 +757,10 @@ mod tests {
 
         #[quickcheck]
         fn point_with_arbitrary_barycentric_coords_has_expected_distance(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
             barycentric_coords: BarycentricCoords,
         ) -> TestResult {
             intersect_with_barycentric_and_test_result(
@@ -781,10 +782,10 @@ mod tests {
 
         #[quickcheck]
         fn point_with_arbitrary_barycentric_coords_has_expected_retro(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
             barycentric_coords: BarycentricCoords,
         ) -> TestResult {
             intersect_with_barycentric_and_test_result(
@@ -806,11 +807,11 @@ mod tests {
 
         #[quickcheck]
         fn intersection_fails_when_ray_outside_first_edge(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
-            uv: Vec2,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
+            uv: Vec2<f64>,
         ) -> bool {
             let uv_origin = Vec3::from(vertex0);
             let u_axis = (vertex1 - vertex0).normalize();
@@ -834,11 +835,11 @@ mod tests {
 
         #[quickcheck]
         fn intersection_fails_when_ray_outside_second_edge(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
-            uv: Vec2,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
+            uv: Vec2<f64>,
         ) -> bool {
             let uv_origin = Vec3::from(vertex0);
             let u_axis = (vertex2 - vertex1).normalize();
@@ -862,11 +863,11 @@ mod tests {
 
         #[quickcheck]
         fn intersection_fails_when_ray_outside_third_edge(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
-            uv: Vec2,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
+            uv: Vec2<f64>,
         ) -> bool {
             let uv_origin = Vec3::from(vertex0);
             let u_axis = (vertex0 - vertex2).normalize();
@@ -890,10 +891,10 @@ mod tests {
 
         #[quickcheck]
         fn intersection_fails_when_triangle_is_behind_ray(
-            vertex0: Vec3,
-            vertex1: Vec3,
-            vertex2: Vec3,
-            ray_origin: Vec3,
+            vertex0: Vec3<f64>,
+            vertex1: Vec3<f64>,
+            vertex2: Vec3<f64>,
+            ray_origin: Vec3<f64>,
             barycentric_coords: BarycentricCoords,
         ) -> bool {
             let point_behind_ray = vertex0 * barycentric_coords.alpha

src/scene.rs

@@ -3,6 +3,6 @@ use crate::math::Vec3;
 use crate::raycasting::Aggregate;
 
 pub struct Scene {
-    pub camera_location: Vec3,
+    pub camera_location: Vec3<f64>,
     pub objects: Vec<Box<dyn Aggregate>>,
 }

src/util/algebra_utils.rs

@@ -1,6 +1,6 @@
 use crate::math::{Mat3, Vec3};
 
-pub fn try_change_of_basis_matrix(x: &Vec3, y: &Vec3, z: &Vec3) -> Option<Mat3> {
+pub fn try_change_of_basis_matrix(x: &Vec3<f64>, y: &Vec3<f64>, z: &Vec3<f64>) -> Option<Mat3<f64>> {
     Some(Mat3::from_rows(x, y, z))
 }
 
@@ -15,15 +15,15 @@ mod tests {
 
         #[test]
         fn produces_isentity_when_passed_axes() {
-            let target: Mat3 =
+            let target: Mat3<f64> =
                 try_change_of_basis_matrix(&Vec3::unit_x(), &Vec3::unit_y(), &Vec3::unit_z())
                     .unwrap();
             assert!(target == Mat3::identity())
         }
 
         #[quickcheck]
-        fn swap_xy_does_not_change_z(v: Vec3) {
-            let target: Mat3 =
+        fn swap_xy_does_not_change_z(v: Vec3<f64>) {
+            let target: Mat3<f64> =
                 try_change_of_basis_matrix(&Vec3::unit_y(), &Vec3::unit_x(), &Vec3::unit_z())
                     .unwrap();
             let v2 = target * v;
@@ -31,8 +31,8 @@ mod tests {
         }
 
         #[quickcheck]
-        fn swap_xy_copies_y_to_x(v: Vec3) {
-            let target: Mat3 =
+        fn swap_xy_copies_y_to_x(v: Vec3<f64>) {
+            let target: Mat3<f64> =
                 try_change_of_basis_matrix(&Vec3::unit_y(), &Vec3::unit_x(), &Vec3::unit_z())
                     .unwrap();
             let v2 = target * v;
@@ -40,8 +40,8 @@ mod tests {
         }
 
         #[quickcheck]
-        fn swap_xy_copies_x_to_y(v: Vec3) {
-            let target: Mat3 =
+        fn swap_xy_copies_x_to_y(v: Vec3<f64>) {
+            let target: Mat3<f64> =
                 try_change_of_basis_matrix(&Vec3::unit_y(), &Vec3::unit_x(), &Vec3::unit_z())
                     .unwrap();
             let v2 = target * v;

src/util/axis_aligned_bounding_box.rs

@@ -9,7 +9,7 @@ pub struct BoundingBox {
 }
 
 impl BoundingBox {
-    pub fn from_corners(a: Vec3, b: Vec3) -> Self {
+    pub fn from_corners(a: Vec3<f64>, b: Vec3<f64>) -> Self {
         let mut result = BoundingBox {
             bounds: [Interval::infinite(); 3],
         };
@@ -27,7 +27,7 @@ impl BoundingBox {
         }
     }
 
-    pub fn from_point(p: Vec3) -> Self {
+    pub fn from_point(p: Vec3<f64>) -> Self {
         BoundingBox {
             bounds: [
                 Interval::degenerate(p.x()),
@@ -39,14 +39,14 @@ impl BoundingBox {
 
     pub fn from_points<'a, I>(points: I) -> Self
     where
-        I: IntoIterator<Item = &'a Vec3>,
+        I: IntoIterator<Item = &'a Vec3<f64>>,
     {
         points
             .into_iter()
             .fold(BoundingBox::empty(), |acc, p| acc.expand_to_point(p))
     }
 
-    pub fn expand_to_point(&self, p: &Vec3) -> Self {
+    pub fn expand_to_point(&self, p: &Vec3<f64>) -> Self {
         BoundingBox {
             bounds: [
                 self.bounds[0].expand_to_value(p.x()),
@@ -56,7 +56,7 @@ impl BoundingBox {
         }
     }
 
-    pub fn contains_point(&self, p: Vec3) -> bool {
+    pub fn contains_point(&self, p: Vec3<f64>) -> bool {
         self.bounds
             .iter()
             .zip(p.coords.iter())
@@ -123,7 +123,7 @@ mod tests {
         let corner_110 = Vec3::new(1.0, 1.0, 0.0);
         let corner_111 = Vec3::new(1.0, 1.0, 1.0);
 
-        let test_inputs: Vec<(Vec3, Vec3)> = vec![
+        let test_inputs: Vec<(Vec3<f64>, Vec3<f64>)> = vec![
             (corner_000, corner_111),
             (corner_001, corner_110),
             (corner_010, corner_101),
@@ -143,7 +143,7 @@ mod tests {
     }
 
     #[quickcheck]
-    fn union_with_self_yields_self(a: Vec3, b: Vec3) -> bool {
+    fn union_with_self_yields_self(a: Vec3<f64>, b: Vec3<f64>) -> bool {
         let target = BoundingBox::from_corners(a, b);
         let result = target.union(&target);
         target
@@ -154,7 +154,7 @@ mod tests {
     }
 
     #[quickcheck]
-    fn union_yields_full_ranges(a: Vec3, b: Vec3, c: Vec3, d: Vec3) -> bool {
+    fn union_yields_full_ranges(a: Vec3<f64>, b: Vec3<f64>, c: Vec3<f64>, d: Vec3<f64>) -> bool {
         let target1 = BoundingBox::from_corners(a, b);
         let target2 = BoundingBox::from_corners(c, d);
         let result = target1.union(&target2);
@@ -172,19 +172,19 @@ mod tests {
     }
 
     #[quickcheck]
-    fn empty_box_contains_no_points(p: Vec3) -> bool {
+    fn empty_box_contains_no_points(p: Vec3<f64>) -> bool {
         let target = BoundingBox::empty();
         !target.contains_point(p)
     }
 
     #[quickcheck]
     fn from_points_produces_box_that_contains_only_points_bounded_by_inputs_on_all_axes(
-        p: Vec3,
-        a: Vec3,
-        b: Vec3,
-        c: Vec3,
-        d: Vec3,
-        e: Vec3,
+        p: Vec3<f64>,
+        a: Vec3<f64>,
+        b: Vec3<f64>,
+        c: Vec3<f64>,
+        d: Vec3<f64>,
+        e: Vec3<f64>,
     ) -> bool {
         let points = vec![a, b, c, d, e];
         let target = BoundingBox::from_points(&points);

src/util/interval.rs

@@ -15,15 +15,15 @@ impl Interval {
 
     pub fn empty() -> Self {
         Interval {
-            min: std::f64::INFINITY,
-            max: std::f64::NEG_INFINITY,
+            min: f64::INFINITY,
+            max: f64::NEG_INFINITY,
         }
     }
 
     pub fn infinite() -> Self {
         Interval {
-            min: std::f64::NEG_INFINITY,
-            max: std::f64::INFINITY,
+            min: f64::NEG_INFINITY,
+            max: f64::INFINITY,
         }
     }
 

src/util/morton.rs

@@ -9,7 +9,7 @@ fn spread_bits(v: u32) -> u32 {
     result
 }
 
-pub fn morton_order_value_3d(p: Vec3) -> u32 {
+pub fn morton_order_value_3d(p: Vec3<f64>) -> u32 {
     let x = p.x().normalized_to_u32(10);
     let y = p.y().normalized_to_u32(10);
     let z = p.z().normalized_to_u32(10);

src/util/normalizer.rs

@@ -47,7 +47,7 @@ impl Point3Normalizer {
         normalizer
     }
 
-    pub fn normalize(&self, point: Vec3) -> Vec3 {
+    pub fn normalize(&self, point: Vec3<f64>) -> Vec3<f64> {
         let mut result = Vec3::new(0.0, 0.0, 0.0);
         for (value_out, &value_in, normalizer) in izip!(
             result.coords.iter_mut(),
@@ -59,7 +59,7 @@ impl Point3Normalizer {
         result
     }
 
-    pub fn normalize_and_clamp(&self, point: Vec3) -> Vec3 {
+    pub fn normalize_and_clamp(&self, point: Vec3<f64>) -> Vec3<f64> {
         let mut result = Vec3::new(0.0, 0.0, 0.0);
         for (value_out, &value_in, normalizer) in izip!(
             result.coords.iter_mut(),
@@ -145,7 +145,11 @@ mod test {
     }
 
     #[quickcheck]
-    fn normalize_point3_is_the_same_as_normalize_each_dimension(a: Vec3, b: Vec3, c: Vec3) -> bool {
+    fn normalize_point3_is_the_same_as_normalize_each_dimension(
+        a: Vec3<f64>,
+        b: Vec3<f64>,
+        c: Vec3<f64>,
+    ) -> bool {
         let x_normalizer = RealNormalizer::new(Interval::new(a.x().min(b.x()), a.x().max(b.x())));
         let y_normalizer = RealNormalizer::new(Interval::new(a.y().min(b.y()), a.y().max(b.y())));
         let z_normalizer = RealNormalizer::new(Interval::new(a.z().min(b.z()), a.z().max(b.z())));
@@ -158,9 +162,9 @@ mod test {
 
     #[quickcheck]
     fn normalize_and_clamp_point3_is_the_same_as_normalize_and_clamp_each_dimension(
-        a: Vec3,
-        b: Vec3,
-        c: Vec3,
+        a: Vec3<f64>,
+        b: Vec3<f64>,
+        c: Vec3<f64>,
     ) -> bool {
         let x_normalizer = RealNormalizer::new(Interval::new(a.x().min(b.x()), a.x().max(b.x())));
         let y_normalizer = RealNormalizer::new(Interval::new(a.y().min(b.y()), a.y().max(b.y())));

src/util/polyhedra.rs

@@ -7,8 +7,8 @@ use crate::raycasting::{Primitive, Triangle};
 use std::sync::Arc;
 
 pub fn triangulate_polygon(
-    vertices: &[Vec3],
-    normal: &Vec3,
+    vertices: &[Vec3<f64>],
+    normal: &Vec3<f64>,
     material: Arc<dyn Material>,
 ) -> Vec<Arc<dyn Primitive>> {
     assert!(vertices.len() >= 3);
@@ -25,7 +25,7 @@ pub fn triangulate_polygon(
 }
 
 pub fn generate_dodecahedron(
-    centre: Vec3,
+    centre: Vec3<f64>,
     size: f64,
     material: Arc<dyn Material>,
 ) -> Vec<Arc<dyn Primitive>> {

src/util/tile_iterator.rs

@@ -27,7 +27,7 @@ pub struct TileIterator {
 impl TileIterator {
     pub fn new(total_width: usize, total_height: usize, tile_size: usize) -> TileIterator {
         // If tile_size*2 is greater than usize::max_value(), increment would overflow
-        assert!(tile_size > 0 && tile_size * 2 < usize::max_value());
+        assert!(tile_size > 0 && tile_size * 2 < usize::MAX);
         TileIterator {
             tile_size,
             total_width,