Fix a vector negation bug.

This helps out when I go to use the same data structure as the backing
for colors.

Cargo.lock

@@ -3331,11 +3331,18 @@ dependencies = [
  "rand_core 0.6.4",
 ]
 
+[[package]]
+name = "ray-tracer"
+version = "0.1.0"
+dependencies = [
+ "rayon",
+]
+
 [[package]]
 name = "rayon"
-version = "1.8.0"
+version = "1.10.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9c27db03db7734835b3f53954b534c91069375ce6ccaa2e065441e07d9b6cdb1"
+checksum = "b418a60154510ca1a002a752ca9714984e21e4241e804d32555251faf8b78ffa"
 dependencies = [
  "either",
  "rayon-core",
@@ -3343,9 +3350,9 @@ dependencies = [
 
 [[package]]
 name = "rayon-core"
-version = "1.12.0"
+version = "1.12.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5ce3fb6ad83f861aac485e76e1985cd109d9a3713802152be56c3b1f0e0658ed"
+checksum = "1465873a3dfdaa8ae7cb14b4383657caab0b3e8a0aa9ae8e04b044854c8dfce2"
 dependencies = [
  "crossbeam-deque",
  "crossbeam-utils",

Cargo.toml

@@ -27,5 +27,5 @@ members = [
     "sgf",
     "timezone-testing",
     "tree",
-    "visions/server",
+    "visions/server", "ray-tracer",
 ]

ray-tracer/Cargo.toml

@@ -0,0 +1,12 @@
+[package]
+name = "ray-tracer"
+version = "0.1.0"
+edition = "2021"
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+
+[dependencies]
+rayon = "1.10.0"
+
+[[bin]]
+name = "projectile"

ray-tracer/src/bin/projectile.rs

@@ -0,0 +1,53 @@
+use ray_tracer::{types::*, PPM};
+use std::{fs::File, io::Write};
+
+#[derive(Clone, Copy)]
+struct Projectile {
+    position: Point,
+    velocity: Vector,
+}
+
+struct Environment {
+    gravity: Vector,
+    wind: Vector,
+}
+
+fn tick(env: &Environment, proj: &Projectile) -> Projectile {
+    let position = proj.position + proj.velocity;
+    let velocity = proj.velocity + env.gravity + env.wind;
+    Projectile { position, velocity }
+}
+
+fn main() {
+    let mut canvas = Canvas::new(900, 550);
+
+    let start = Projectile {
+        position: Point::new(0., 1., 0.),
+        velocity: Vector::new(1., 5., 0.).normalize() * 5.,
+    };
+
+    let e = Environment {
+        gravity: Vector::new(0., -0.1, 0.),
+        wind: Vector::new(-0.01, 0., 0.),
+    };
+
+    let mut p = start;
+    while p.position.y() > 0. {
+        p = tick(&e, &p);
+
+        let x = p.position.x().round() as usize;
+        let y = p.position.y().round() as usize;
+        if x > 1 && x < 900 && y > 1 && y < 550 {
+            *canvas.pixel_mut(x, 550 - y - 1) = Color::new(1., 1., 0.);
+            *canvas.pixel_mut(x+1, 550 - y - 1) = Color::new(1., 1., 0.);
+            *canvas.pixel_mut(x-1, 550 - y - 1) = Color::new(1., 1., 0.);
+            *canvas.pixel_mut(x, 550 - y) = Color::new(1., 1., 0.);
+            *canvas.pixel_mut(x, 550 - y - 2) = Color::new(1., 1., 0.);
+        }
+    }
+
+    let ppm = PPM::from(&canvas);
+
+    let mut file = File::create("projectile.ppm").unwrap();
+    let _ = file.write(ppm.as_bytes());
+}

ray-tracer/src/bin/sphere.rs

@@ -0,0 +1,125 @@
+use ray_tracer::{types::*, PPM};
+use rayon::prelude::*;
+use std::{
+    fs::File,
+    io::Write,
+    sync::{Arc, RwLock},
+};
+
+const SIZE: usize = 1000;
+
+const EPSILON: f64 = 0.00001;
+
+fn eq_f64(l: f64, r: f64) -> bool {
+    (l - r).abs() < EPSILON
+}
+
+
+fn render(
+    camera: &Point,
+    light: &PointLight,
+    sphere: &Sphere,
+    wall_z: f64,
+    wall_size: f64,
+    filename: &str,
+) {
+    let canvas = Arc::new(RwLock::new(Canvas::new(SIZE, SIZE)));
+
+    let pixel_size = wall_size / SIZE as f64;
+    let half = wall_size / 2.;
+
+    let ray = Ray::new(camera.clone(), (Point::new(0., 0., wall_z) - camera).normalize());
+    let xs = ray.intersect(sphere);
+    /*
+    match xs.hit() {
+        Some(hit) => {
+            let point = ray.position(hit.t);
+            let normal = hit.object.normal_at(&point);
+            let eye = -ray.direction;
+            let color = hit.object.material().lighting(&light, &point, &eye, &normal);
+            println!("{:?}", color);
+        }
+        None => {}
+    }
+    */
+
+    (0..SIZE).into_par_iter().for_each(|x| {
+        (0..SIZE).into_par_iter().for_each(|y| {
+            let world_x = -half + pixel_size * x as f64;
+            let world_y = half - pixel_size * y as f64;
+            let position = Point::new(world_x, world_y, wall_z);
+
+            let ray = Ray::new(camera.clone(), (position - camera).normalize());
+            let xs = ray.intersect(sphere);
+            match xs.hit() {
+                Some(hit) => {
+                    let point = ray.position(hit.t);
+                    let normal = hit.object.normal_at(&point);
+                    let eye = -ray.direction;
+
+                    let color = hit
+                        .object
+                        .material()
+                        .lighting(&light, &point, &eye, &normal);
+                    *canvas.write().unwrap().pixel_mut(x, y) = color;
+                }
+                None => {}
+            }
+        });
+    });
+
+    let ppm = PPM::from(&*canvas.read().unwrap());
+    let mut file = File::create(filename).unwrap();
+    let _ = file.write(ppm.as_bytes());
+}
+
+fn main() {
+    let mut sphere = Sphere::default();
+
+    let mut material = Material::default();
+    material.color = Color::new(1., 0.2, 1.);
+    *sphere.material_mut() = material;
+
+    let camera = Point::new(0., 0., -5.);
+    let wall_z = 10.;
+    let wall_size = 7.;
+    let light = PointLight::new(Point::new(-5., 5., -10.), Color::new(1., 1., 1.));
+
+    render(&camera, &light, &sphere, wall_z, wall_size, "sphere_1.ppm");
+    /*
+    {
+        let light = PointLight::new(Point::new(0., 0., -10.), Color::new(1., 1., 1.));
+        println!("Light at {:0.2} {:0.2} {:0.2}", light.position.x(), light.position.y(), light.position.z());
+        render(&camera, &light, &sphere, wall_z, wall_size, "sphere_1.ppm");
+        println!();
+    }
+
+    {
+        let light = PointLight::new(Point::new(-2., 0., -10.), Color::new(1., 1., 1.));
+        println!("Light at {:0.2} {:0.2} {:0.2}", light.position.x(), light.position.y(), light.position.z());
+        render(&camera, &light, &sphere, wall_z, wall_size, "sphere_2.ppm");
+        println!();
+    }
+
+    {
+        let light = PointLight::new(Point::new(2., 0., -10.), Color::new(1., 1., 1.));
+        println!("Light at {:0.2} {:0.2} {:0.2}", light.position.x(), light.position.y(), light.position.z());
+        render(&camera, &light, &sphere, wall_z, wall_size, "sphere_3.ppm");
+        println!();
+    }
+
+    {
+        let light = PointLight::new(Point::new(0., -2., -10.), Color::new(1., 1., 1.));
+        println!("Light at {:0.2} {:0.2} {:0.2}", light.position.x(), light.position.y(), light.position.z());
+        render(&camera, &light, &sphere, wall_z, wall_size, "sphere_4.ppm");
+        println!();
+    }
+
+    {
+        let light = PointLight::new(Point::new(0., 2., -10.), Color::new(1., 1., 1.));
+        println!("Light at {:0.2} {:0.2} {:0.2}", light.position.x(), light.position.y(), light.position.z());
+        render(&camera, &light, &sphere, wall_z, wall_size, "sphere_5.ppm");
+        println!();
+    }
+    */
+}

ray-tracer/src/lib.rs

@@ -0,0 +1,5 @@
+mod ppm;
+pub mod transforms;
+pub mod types;
+
+pub use ppm::PPM;

ray-tracer/src/ppm.rs

@@ -0,0 +1,118 @@
+use crate::types::Canvas;
+
+fn color_float_to_int(val: f64) -> u8 {
+    (val * 256.).clamp(0., 255.) as u8
+}
+
+fn join_to_line_limit(data: impl IntoIterator<Item = String>) -> Vec<String> {
+    let mut lines = vec![];
+
+    let mut line = String::new();
+    for element in data {
+        if line.is_empty() {
+            line = line + &element;
+        } else if line.len() + 1 + element.len() < 70 {
+            line = line + " " + &element;
+        } else {
+            lines.push(line);
+            line = element;
+        }
+    }
+
+    if !line.is_empty() {
+        lines.push(line);
+    }
+
+    lines
+}
+
+#[derive(Debug)]
+pub struct PPM(String);
+
+impl From<&Canvas> for PPM {
+    fn from(c: &Canvas) -> Self {
+        // let v = vec![0.; c.width() * c.height() * 3];
+        let header = format!("P3\n{} {}\n255\n", c.width(), c.height());
+
+        let mut data = vec![];
+        for y in 0..c.height() {
+            let mut row = vec![];
+            for x in 0..c.width() {
+                let pixel = c.pixel(x, y);
+                row.push(color_float_to_int(pixel.red()).to_string());
+                row.push(color_float_to_int(pixel.green()).to_string());
+                row.push(color_float_to_int(pixel.blue()).to_string());
+            }
+            let mut lines = join_to_line_limit(row);
+            data.append(&mut lines);
+        }
+
+        let data = data.join("\n");
+
+        Self(format!("{header}{data}\n"))
+    }
+}
+
+impl std::ops::Deref for PPM {
+    type Target = String;
+    fn deref(&self) -> &String {
+        &self.0
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::types::Color;
+
+    #[test]
+    fn construct_ppm_header() {
+        let c = Canvas::new(5, 3);
+        let ppm = PPM::from(&c);
+
+        assert!(ppm.starts_with(
+            "P3\n\
+                5 3\n\
+                255\n"
+        ));
+    }
+
+    #[test]
+    fn construct_full_ppm_data() {
+        let mut c = Canvas::new(5, 3);
+        *c.pixel_mut(0, 0) = Color::new(1.5, 0., 0.);
+        *c.pixel_mut(2, 1) = Color::new(0., 0.5, 0.);
+        *c.pixel_mut(4, 2) = Color::new(-0.5, 0., 1.);
+
+        let expected = "P3\n\
+                        5 3\n\
+                        255\n\
+                        255 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n\
+                        0 0 0 0 0 0 0 128 0 0 0 0 0 0 0\n\
+                        0 0 0 0 0 0 0 0 0 0 0 0 0 0 255\n";
+
+        let ppm = PPM::from(&c);
+        assert_eq!(*ppm, expected);
+    }
+
+    #[test]
+    fn ppm_line_length_limit() {
+        let mut c = Canvas::new(10, 2);
+        for y in 0..2 {
+            for x in 0..10 {
+                *c.pixel_mut(x, y) = Color::new(1., 0.8, 0.6);
+            }
+        }
+
+        let expected = "P3\n\
+                        10 2\n\
+                        255\n\
+                        255 204 153 255 204 153 255 204 153 255 204 153 255 204 153 255 204\n\
+     	                153 255 204 153 255 204 153 255 204 153 255 204 153\n\
+     	                255 204 153 255 204 153 255 204 153 255 204 153 255 204 153 255 204\n\
+     	                153 255 204 153 255 204 153 255 204 153 255 204 153\n";
+
+        let ppm = PPM::from(&c);
+        assert_eq!(*ppm, expected);
+    }
+}

ray-tracer/src/transforms.rs

@@ -0,0 +1,219 @@
+use crate::types::Matrix;
+
+pub fn translation(x: f64, y: f64, z: f64) -> Matrix {
+    Matrix::from([
+        [1., 0., 0., x],
+        [0., 1., 0., y],
+        [0., 0., 1., z],
+        [0., 0., 0., 1.],
+    ])
+}
+
+pub fn scaling(x: f64, y: f64, z: f64) -> Matrix {
+    Matrix::from([
+        [x, 0., 0., 0.],
+        [0., y, 0., 0.],
+        [0., 0., z, 0.],
+        [0., 0., 0., 1.],
+    ])
+}
+
+pub fn rotation_x(r: f64) -> Matrix {
+    Matrix::from([
+        [1., 0., 0., 0.],
+        [0., r.cos(), -r.sin(), 0.],
+        [0., r.sin(), r.cos(), 0.],
+        [0., 0., 0., 1.],
+    ])
+}
+
+pub fn rotation_y(r: f64) -> Matrix {
+    Matrix::from([
+        [r.cos(), 0., r.sin(), 0.],
+        [0., 1., 0., 0.],
+        [-r.sin(), 0., r.cos(), 0.],
+        [0., 0., 0., 1.],
+    ])
+}
+
+pub fn rotation_z(r: f64) -> Matrix {
+    Matrix::from([
+        [r.cos(), -r.sin(), 0., 0.],
+        [r.sin(), r.cos(), 0., 0.],
+        [0., 0., 1., 0.],
+        [0., 0., 0., 1.],
+    ])
+}
+
+pub fn shearing(
+    x_by_y: f64,
+    x_by_z: f64,
+    y_by_x: f64,
+    y_by_z: f64,
+    z_by_x: f64,
+    z_by_y: f64,
+) -> Matrix {
+    Matrix::from([
+        [1., x_by_y, x_by_z, 0.],
+        [y_by_x, 1., y_by_z, 0.],
+        [z_by_x, z_by_y, 1., 0.],
+        [0., 0., 0., 1.],
+    ])
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::types::{Point, Vector};
+    use std::f64::consts::PI;
+
+    #[test]
+    fn multiply_by_translation_matrix() {
+        let tr = translation(5., -3., 2.);
+        let p = Point::new(-3., 4., 5.);
+        assert_eq!(&tr * p, Point::new(2., 1., 7.));
+
+        let inv = tr.inverse();
+        assert_eq!(inv * p, Point::new(-8., 7., 3.));
+    }
+
+    #[test]
+    fn translation_does_not_affect_vectors() {
+        let tr = translation(5., -3., 2.);
+        let v = Vector::new(-3., 4., 5.);
+        assert_eq!(tr * v, v);
+    }
+
+    #[test]
+    fn scaling_matrix_applied_to_point() {
+        let tr = scaling(2., 3., 4.);
+        let p = Point::new(-4., 6., 8.);
+        let v = Vector::new(-4., 6., 8.);
+
+        assert_eq!(&tr * p, Point::new(-8., 18., 32.));
+        assert_eq!(tr * v, Vector::new(-8., 18., 32.));
+    }
+
+    #[test]
+    fn reflection_is_scaling_by_negative_value() {
+        let tr = scaling(-1., 1., 1.);
+        let p = Point::new(2., 3., 4.);
+        assert_eq!(tr * p, Point::new(-2., 3., 4.));
+    }
+
+    #[test]
+    fn rotate_point_around_x() {
+        let p = Point::new(0., 1., 0.);
+        let half_quarter = rotation_x(PI / 4.);
+        let full_quarter = rotation_x(PI / 2.);
+        assert_eq!(
+            half_quarter * p,
+            Point::new(0., 2_f64.sqrt() / 2., 2_f64.sqrt() / 2.)
+        );
+        assert_eq!(full_quarter * p, Point::new(0., 0., 1.));
+    }
+
+    #[test]
+    fn inverse_rotates_opposite_direction() {
+        let p = Point::new(0., 1., 0.);
+        let half_quarter = rotation_x(PI / 4.);
+        let inv = half_quarter.inverse();
+        assert_eq!(
+            inv * p,
+            Point::new(0., 2_f64.sqrt() / 2., -2_f64.sqrt() / 2.)
+        );
+    }
+
+    #[test]
+    fn rotate_around_y() {
+        let p = Point::new(0., 0., 1.);
+        let half_quarter = rotation_y(PI / 4.);
+        let full_quarter = rotation_y(PI / 2.);
+        assert_eq!(
+            half_quarter * p,
+            Point::new(2_f64.sqrt() / 2., 0., 2_f64.sqrt() / 2.)
+        );
+        assert_eq!(full_quarter * p, Point::new(1., 0., 0.));
+    }
+
+    #[test]
+    fn rotate_around_z() {
+        let p = Point::new(0., 1., 0.);
+        let half_quarter = rotation_z(PI / 4.);
+        let full_quarter = rotation_z(PI / 2.);
+        assert_eq!(
+            half_quarter * p,
+            Point::new(-2_f64.sqrt() / 2., 2_f64.sqrt() / 2., 0.)
+        );
+        assert_eq!(full_quarter * p, Point::new(-1., 0., 0.));
+    }
+
+    #[test]
+    fn shear_moves_x_proportionally_to_y() {
+        let transform = shearing(1., 0., 0., 0., 0., 0.);
+        let p = Point::new(2., 3., 4.);
+        assert_eq!(transform * p, Point::new(5., 3., 4.));
+    }
+
+    #[test]
+    fn shear_moves_x_proportionally_to_z() {
+        let transform = shearing(0., 1., 0., 0., 0., 0.);
+        let p = Point::new(2., 3., 4.);
+        assert_eq!(transform * p, Point::new(6., 3., 4.));
+    }
+
+    #[test]
+    fn shear_moves_y_proportionally_to_x() {
+        let transform = shearing(0., 0., 1., 0., 0., 0.);
+        let p = Point::new(2., 3., 4.);
+        assert_eq!(transform * p, Point::new(2., 5., 4.));
+    }
+
+    #[test]
+    fn shear_moves_y_proportionally_to_z() {
+        let transform = shearing(0., 0., 0., 1., 0., 0.);
+        let p = Point::new(2., 3., 4.);
+        assert_eq!(transform * p, Point::new(2., 7., 4.));
+    }
+
+    #[test]
+    fn shear_moves_z_proportionally_to_x() {
+        let transform = shearing(0., 0., 0., 0., 1., 0.);
+        let p = Point::new(2., 3., 4.);
+        assert_eq!(transform * p, Point::new(2., 3., 6.));
+    }
+
+    #[test]
+    fn shear_moves_z_proportionally_to_y() {
+        let transform = shearing(0., 0., 0., 0., 0., 1.);
+        let p = Point::new(2., 3., 4.);
+        assert_eq!(transform * p, Point::new(2., 3., 7.));
+    }
+
+    #[test]
+    fn individual_transformations_are_applied_in_sequence() {
+        let p = Point::new(1., 0., 1.);
+        let rotation = rotation_x(PI / 2.);
+        let scale = scaling(5., 5., 5.);
+        let translate = translation(10., 5., 7.);
+
+        let p2 = rotation * p;
+        assert_eq!(p2, Point::new(1., -1., 0.));
+
+        let p3 = scale * p2;
+        assert_eq!(p3, Point::new(5., -5., 0.));
+
+        let p4 = translate * p3;
+        assert_eq!(p4, Point::new(15., 0., 7.));
+    }
+
+    #[test]
+    fn chained_translations() {
+        let p = Point::new(1., 0., 1.);
+        let rotation = rotation_x(PI / 2.);
+        let scale = scaling(5., 5., 5.);
+        let translate = translation(10., 5., 7.);
+        let tr = translate * scale * rotation;
+        assert_eq!(tr * p, Point::new(15., 0., 7.));
+    }
+}

ray-tracer/src/types/canvas.rs

@@ -0,0 +1,59 @@
+use crate::types::Color;
+
+pub struct Canvas {
+    width: usize,
+    height: usize,
+    pixels: Vec<Color>,
+}
+
+impl Canvas {
+    pub fn new(width: usize, height: usize) -> Self {
+        Self {
+            width,
+            height,
+            pixels: vec![Color::default(); width * height],
+        }
+    }
+
+    #[inline]
+    pub fn width(&self) -> usize {
+        self.width
+    }
+
+    #[inline]
+    pub fn height(&self) -> usize {
+        self.height
+    }
+
+    pub fn pixel(&self, x: usize, y: usize) -> &Color {
+        &self.pixels[self.addr(y, x)]
+    }
+
+    pub fn pixel_mut(&mut self, x: usize, y: usize) -> &mut Color {
+        let addr = self.addr(y, x);
+        &mut self.pixels[addr]
+    }
+
+    #[inline]
+    fn addr(&self, y: usize, x: usize) -> usize {
+        let val = y * self.width() + x;
+        if val >= self.pixels.len() {
+            eprintln!("[{val}] out of range: [{x}, {y}]");
+        }
+        val
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn addresses() {
+        let c = Canvas::new(5, 3);
+
+        assert_eq!(c.addr(0, 0), 0);
+        assert_eq!(c.addr(1, 0), 5);
+        assert_eq!(c.addr(2, 0), 10);
+    }
+}

ray-tracer/src/types/color.rs

@@ -0,0 +1,86 @@
+use crate::types::Tuple;
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Color(Tuple);
+
+impl Color {
+    pub fn new(red: f64, green: f64, blue: f64) -> Self {
+        Self(Tuple(red, green, blue, 0.))
+    }
+
+    #[inline]
+    pub fn red(&self) -> f64 {
+        self.0.0
+    }
+
+    #[inline]
+    pub fn green(&self) -> f64 {
+        self.0.1
+    }
+
+    #[inline]
+    pub fn blue(&self) -> f64 {
+        self.0.2
+    }
+}
+
+impl Default for Color {
+    fn default() -> Self {
+        Self::new(0., 0., 0.)
+    }
+}
+
+impl From<Tuple> for Color {
+    fn from(tuple: Tuple) -> Self {
+        assert_eq!(tuple.3, 0.0);
+        Self(tuple)
+    }
+}
+
+impl std::ops::Deref for Color {
+    type Target = Tuple;
+    fn deref(&self) -> &Tuple {
+        &self.0
+    }
+}
+
+impl std::ops::Add for Color {
+    type Output = Color;
+    fn add(self, r: Self) -> Self {
+        Color::from(self.0 + r.0)
+    }
+}
+
+impl std::ops::Sub for Color {
+    type Output = Color;
+    fn sub(self, r: Self) -> Self {
+        Color::from(self.0 - r.0)
+    }
+}
+
+impl std::ops::Neg for Color {
+    type Output = Color;
+    fn neg(self) -> Self::Output {
+        let mut t = -self.0;
+        t.0 = 0.;
+        Color::from(t)
+    }
+}
+
+impl std::ops::Mul for Color {
+    type Output = Color;
+    fn mul(self, r: Color) -> Self {
+        let red = self.red() * r.red();
+        let green = self.green() * r.green();
+        let blue = self.blue() * r.blue();
+        Self::new(red, green, blue)
+    }
+}
+
+impl std::ops::Mul<f64> for Color {
+    type Output = Color;
+    fn mul(self, r: f64) -> Self {
+        Color::from(self.0 * r)
+    }
+}
+

ray-tracer/src/types/intersections.rs

@@ -0,0 +1,39 @@
+use std::cmp::Ordering;
+use super::Sphere;
+
+#[derive(Clone, Debug, PartialEq)]
+pub struct Intersection<'a> {
+    pub t: f64,
+    pub object: &'a Sphere,
+}
+
+pub struct Intersections<'a>(Vec<Intersection<'a>>);
+
+impl<'a> Intersections<'a> {
+    pub fn is_empty(&'a self) -> bool {
+        self.0.is_empty()
+    }
+
+    pub fn len(&'a self) -> usize {
+        self.0.len()
+    }
+
+    pub fn hit(&'a self) -> Option<&Intersection<'a>> {
+        self.0.iter().find(|i| i.t >= 0.)
+    }
+}
+
+impl<'a> std::ops::Index<usize> for Intersections<'a> {
+    type Output = Intersection<'a>;
+    fn index(&self, idx: usize) -> &Intersection<'a> {
+        &self.0[idx]
+    }
+}
+
+impl<'a> From<Vec<Intersection<'a>>> for Intersections<'a> {
+    fn from(mut v: Vec<Intersection<'a>>) -> Self {
+        v.sort_by(|l, r| l.t.partial_cmp(&r.t).unwrap_or(Ordering::Equal));
+        Self(v)
+    }
+}
+

ray-tracer/src/types/light.rs

@@ -0,0 +1,12 @@
+use super::{Color, Point};
+
+pub struct PointLight {
+    pub position: Point,
+    pub intensity: Color,
+}
+
+impl PointLight {
+    pub fn new(position: Point, intensity: Color) -> Self {
+        Self { position, intensity }
+    }
+}

ray-tracer/src/types/material.rs

@@ -0,0 +1,152 @@
+use super::{eq_f64, Color, Point, PointLight, Vector};
+
+#[derive(Debug, PartialEq)]
+pub struct Material {
+    pub color: Color,
+    pub ambient: f64,
+    pub diffuse: f64,
+    pub specular: f64,
+    pub shininess: f64,
+}
+
+impl Default for Material {
+    fn default() -> Self {
+        Self {
+            color: Color::new(1., 1., 1.),
+            ambient: 0.1,
+            diffuse: 0.9,
+            specular: 0.9,
+            shininess: 200.,
+        }
+    }
+}
+
+impl Material {
+    pub fn lighting(
+        &self,
+        light: &PointLight,
+        position: &Point,
+        eye: &Vector,
+        normal_v: &Vector,
+    ) -> Color {
+        // let debug = eq_f64(position.x(), -0.09344)
+        //     || eq_f64(position.x(), 0.09344) && eq_f64(position.y(), 0.)
+        //     || eq_f64(position.y(), -0.09344)
+        //     || eq_f64(position.y(), 0.09344) && eq_f64(position.x(), 0.);
+        // if debug {
+        //     println!(
+        //         "Debugging {:0.5} {:0.5} {:0.5}",
+        //         position.x(),
+        //         position.y(),
+        //         position.z()
+        //     );
+        // }
+        let effective_color = self.color * light.intensity;
+
+        let light_v = (light.position - position).normalize();
+        let ambient = effective_color * self.ambient;
+
+        let light_dot_normal = light_v.dot(&normal_v);
+        let (diffuse, specular) = if light_dot_normal < 0. {
+            (Color::new(0., 0., 0.), Color::new(0., 0., 0.))
+        } else {
+            let diffuse = effective_color * self.diffuse * light_dot_normal;
+
+            let reflect_v = (-light_v).reflect(&normal_v);
+            let reflect_dot_eye = reflect_v.dot(&eye);
+
+            let specular = if reflect_dot_eye <= 0. {
+                Color::new(0., 0., 0.)
+            } else {
+                let factor = reflect_dot_eye.powf(self.shininess);
+                light.intensity * self.specular * factor
+            };
+
+            (diffuse, specular)
+        };
+
+        ambient + diffuse + specular
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::types::{Color, Point, PointLight, Vector};
+
+    #[test]
+    fn lighting_with_eye_between_light_and_surface() {
+        let m = Material::default();
+        let position = Point::default();
+
+        let eyev = Vector::new(0., 0., -1.);
+        let normalv = Vector::new(0., 0., -1.);
+        let light = PointLight::new(Point::new(0., 0., -10.), Color::new(1., 1., 1.));
+
+        assert_eq!(
+            m.lighting(&light, &position, &eyev, &normalv),
+            Color::new(1.9, 1.9, 1.9)
+        );
+    }
+
+    #[test]
+    fn lighting_with_eye_between_light_and_surface_eye_offset_45() {
+        let m = Material::default();
+        let position = Point::default();
+
+        let eyev = Vector::new(0., 2_f64.sqrt() / 2., -2_f64.sqrt() / 2.);
+        let normalv = Vector::new(0., 0., -1.);
+        let light = PointLight::new(Point::new(0., 0., -10.), Color::new(1., 1., 1.));
+
+        assert_eq!(
+            m.lighting(&light, &position, &eyev, &normalv),
+            Color::new(1., 1., 1.)
+        );
+    }
+
+    #[test]
+    fn lighting_with_eye_between_light_and_surface_light_offset_45() {
+        let m = Material::default();
+        let position = Point::default();
+
+        let eyev = Vector::new(0., 0., -1.);
+        let normalv = Vector::new(0., 0., -1.);
+        let light = PointLight::new(Point::new(0., 10., -10.), Color::new(1., 1., 1.));
+
+        assert_eq!(
+            m.lighting(&light, &position, &eyev, &normalv),
+            Color::new(0.7364, 0.7364, 0.7364)
+        );
+    }
+
+    #[test]
+    fn lighting_eye_in_path_of_reflection() {
+        let m = Material::default();
+        let position = Point::default();
+
+        // let eyev = Vector::new(0., -2_f64.sqrt() / 2., -2_f64.sqrt() / 2.);
+        let eyev = Vector::new(0., -10., -10.).normalize();
+        let normalv = Vector::new(0., 0., -1.);
+        let light = PointLight::new(Point::new(0., 10., -10.), Color::new(1., 1., 1.));
+
+        assert_eq!(
+            m.lighting(&light, &position, &eyev, &normalv),
+            Color::new(1.6364, 1.6364, 1.6364)
+        );
+    }
+
+    #[test]
+    fn lighting_with_light_behind_surface() {
+        let m = Material::default();
+        let position = Point::default();
+
+        let eyev = Vector::new(0., 0., -1.);
+        let normalv = Vector::new(0., 0., -1.);
+        let light = PointLight::new(Point::new(0., 0., 10.), Color::new(1., 1., 1.));
+
+        assert_eq!(
+            m.lighting(&light, &position, &eyev, &normalv),
+            Color::new(0.1, 0.1, 0.1)
+        );
+    }
+}

ray-tracer/src/types/matrix.rs

@@ -0,0 +1,544 @@
+use std::cmp::Ordering;
+
+use crate::types::{eq_f64, Point, Tuple, Vector};
+
+#[derive(Clone, Debug)]
+pub struct Matrix {
+    size: usize,
+    values: Vec<f64>,
+}
+
+impl Matrix {
+    pub fn new(size: usize) -> Self {
+        Self {
+            size,
+            values: vec![0.; size * size],
+        }
+    }
+
+    pub fn identity() -> Self {
+        Self::from([
+            [1., 0., 0., 0.],
+            [0., 1., 0., 0.],
+            [0., 0., 1., 0.],
+            [0., 0., 0., 1.],
+        ])
+    }
+
+    pub fn cell(&self, row: usize, column: usize) -> f64 {
+        self.values[self.addr(row, column)]
+    }
+
+    pub fn cell_mut(&mut self, row: usize, column: usize) -> &mut f64 {
+        let addr = self.addr(row, column);
+        &mut self.values[addr]
+    }
+
+    pub fn transpose(&self) -> Self {
+        let mut m = Self::new(self.size);
+
+        for row in 0..self.size {
+            for column in 0..self.size {
+                *m.cell_mut(row, column) = self.cell(column, row);
+            }
+        }
+
+        m
+    }
+
+    pub fn is_invertible(&self) -> bool {
+        !eq_f64(self.determinant(), 0.)
+    }
+
+    pub fn inverse(&self) -> Self {
+        let mut m = Matrix::new(self.size);
+        let determinant = self.determinant();
+        for row in 0..self.size {
+            for column in 0..self.size {
+                *m.cell_mut(row, column) = self.cofactor(row, column);
+            }
+        }
+
+        let mut m = m.transpose();
+        for row in 0..self.size {
+            for column in 0..self.size {
+                let value = m.cell(row, column);
+                *m.cell_mut(row, column) = value / determinant;
+            }
+        }
+
+        m
+    }
+
+    fn determinant(&self) -> f64 {
+        // TODO: optimization may not be necessary, but this can be optimized by memoizing
+        // submatrices and cofactors.
+        if self.size == 2 {
+            self.cell(0, 0) * self.cell(1, 1) - self.cell(0, 1) * self.cell(1, 0)
+        } else if self.size == 3 {
+            self.cell(0, 0) * self.cofactor(0, 0)
+                + self.cell(0, 1) * self.cofactor(0, 1)
+                + self.cell(0, 2) * self.cofactor(0, 2)
+        } else if self.size == 4 {
+            self.cell(0, 0) * self.cofactor(0, 0)
+                + self.cell(0, 1) * self.cofactor(0, 1)
+                + self.cell(0, 2) * self.cofactor(0, 2)
+                + self.cell(0, 3) * self.cofactor(0, 3)
+        } else {
+            0.
+        }
+    }
+
+    fn submatrix(&self, row: usize, column: usize) -> Self {
+        let mut m = Self::new(self.size - 1);
+
+        for r in 0..self.size {
+            for c in 0..self.size {
+                let dest_r = match r.cmp(&row) {
+                    Ordering::Less => r,
+                    Ordering::Greater => r - 1,
+                    Ordering::Equal => continue,
+                };
+                let dest_c = match c.cmp(&column) {
+                    Ordering::Less => c,
+                    Ordering::Greater => c - 1,
+                    Ordering::Equal => continue,
+                };
+                *m.cell_mut(dest_r, dest_c) = self.cell(r, c);
+            }
+        }
+
+        m
+    }
+
+    fn minor(&self, row: usize, column: usize) -> f64 {
+        self.submatrix(row, column).determinant()
+    }
+
+    fn cofactor(&self, row: usize, column: usize) -> f64 {
+        if (row + column) % 2 == 0 {
+            self.minor(row, column)
+        } else {
+            -self.minor(row, column)
+        }
+    }
+
+    #[inline]
+    fn addr(&self, row: usize, column: usize) -> usize {
+        row * self.size + column
+    }
+}
+
+impl From<[[f64; 2]; 2]> for Matrix {
+    fn from(s: [[f64; 2]; 2]) -> Self {
+        Self {
+            size: 2,
+            values: s.concat(),
+        }
+    }
+}
+
+impl From<[[f64; 3]; 3]> for Matrix {
+    fn from(s: [[f64; 3]; 3]) -> Self {
+        Self {
+            size: 3,
+            values: s.concat(),
+        }
+    }
+}
+
+impl From<[[f64; 4]; 4]> for Matrix {
+    fn from(s: [[f64; 4]; 4]) -> Self {
+        Self {
+            size: 4,
+            values: s.concat(),
+        }
+    }
+}
+
+impl PartialEq for Matrix {
+    fn eq(&self, rside: &Matrix) -> bool {
+        if self.size != rside.size {
+            return false;
+        };
+
+        self.values
+            .iter()
+            .zip(rside.values.iter())
+            .all(|(l, r)| eq_f64(*l, *r))
+    }
+}
+
+impl std::ops::Mul<&Matrix> for &Matrix {
+    type Output = Matrix;
+    fn mul(self, rside: &Matrix) -> Matrix {
+        assert_eq!(self.size, 4);
+        assert_eq!(rside.size, 4);
+
+        let mut m = Matrix::new(self.size);
+        for row in 0..4 {
+            for column in 0..4 {
+                *m.cell_mut(row, column) = self.cell(row, 0) * rside.cell(0, column)
+                    + self.cell(row, 1) * rside.cell(1, column)
+                    + self.cell(row, 2) * rside.cell(2, column)
+                    + self.cell(row, 3) * rside.cell(3, column);
+            }
+        }
+
+        m
+    }
+}
+
+impl std::ops::Mul<Matrix> for &Matrix {
+    type Output = Matrix;
+    fn mul(self, rside: Matrix) -> Matrix {
+        self * &rside
+    }
+}
+
+impl std::ops::Mul<&Matrix> for Matrix {
+    type Output = Matrix;
+    fn mul(self, rside: &Matrix) -> Matrix {
+        &self * rside
+    }
+}
+
+impl std::ops::Mul for Matrix {
+    type Output = Matrix;
+    fn mul(self, rside: Matrix) -> Matrix {
+        &self * &rside
+    }
+}
+
+impl std::ops::Mul<Tuple> for &Matrix {
+    type Output = Tuple;
+    fn mul(self, rside: Tuple) -> Tuple {
+        assert_eq!(self.size, 4);
+
+        let mut t = [0.; 4];
+
+        for (row, item) in t.iter_mut().enumerate() {
+            *item = self.cell(row, 0) * rside.0
+                + self.cell(row, 1) * rside.1
+                + self.cell(row, 2) * rside.2
+                + self.cell(row, 3) * rside.3;
+        }
+
+        Tuple::from(t)
+    }
+}
+
+impl std::ops::Mul<Tuple> for Matrix {
+    type Output = Tuple;
+    fn mul(self, rside: Tuple) -> Tuple {
+        &self * rside
+    }
+}
+
+impl std::ops::Mul<&Point> for &Matrix {
+    type Output = Point;
+    fn mul(self, rside: &Point) -> Point {
+        let t: Tuple = **rside;
+        Point::from(self * t)
+    }
+}
+
+impl std::ops::Mul<Point> for &Matrix {
+    type Output = Point;
+    fn mul(self, rside: Point) -> Point {
+        self * &rside
+    }
+}
+
+impl std::ops::Mul<&Point> for Matrix {
+    type Output = Point;
+    fn mul(self, rside: &Point) -> Point {
+        &self * rside
+    }
+}
+
+impl std::ops::Mul<Point> for Matrix {
+    type Output = Point;
+    fn mul(self, rside: Point) -> Point {
+        &self * rside
+    }
+}
+
+impl std::ops::Mul<Vector> for Matrix {
+    type Output = Vector;
+    fn mul(self, rside: Vector) -> Vector {
+        Vector::from(self * *rside)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn it_constructs_a_matrix() {
+        let m = Matrix::from([[-3., 5.], [1., -2.]]);
+        assert_eq!(m.cell(0, 0), -3.);
+        assert_eq!(m.cell(0, 1), 5.);
+        assert_eq!(m.cell(1, 0), 1.);
+        assert_eq!(m.cell(1, 1), -2.);
+
+        let m = Matrix::from([[-3., 5., 0.], [1., -2., -7.], [0., 1., 1.]]);
+        assert_eq!(m.cell(0, 0), -3.);
+        assert_eq!(m.cell(1, 1), -2.);
+        assert_eq!(m.cell(2, 2), 1.);
+
+        let m = Matrix::from([
+            [1., 2., 3., 4.],
+            [5.5, 6.5, 7.5, 8.5],
+            [9., 10., 11., 12.],
+            [13.5, 14.5, 15.5, 16.5],
+        ]);
+        assert_eq!(m.cell(0, 0), 1.);
+        assert_eq!(m.cell(0, 3), 4.);
+        assert_eq!(m.cell(1, 0), 5.5);
+        assert_eq!(m.cell(1, 2), 7.5);
+        assert_eq!(m.cell(2, 2), 11.);
+        assert_eq!(m.cell(3, 0), 13.5);
+        assert_eq!(m.cell(3, 2), 15.5);
+    }
+
+    #[test]
+    fn it_compares_two_matrices() {
+        let a = Matrix::from([
+            [1., 2., 3., 4.],
+            [5., 6., 7., 8.],
+            [9., 8., 7., 6.],
+            [5., 4., 3., 2.],
+        ]);
+        let b = Matrix::from([
+            [1., 2., 3., 4.],
+            [5., 6., 7., 8.],
+            [9., 8., 7., 6.],
+            [5., 4., 3., 2.],
+        ]);
+        let c = Matrix::from([
+            [2., 3., 4., 5.],
+            [6., 7., 8., 9.],
+            [8., 7., 6., 5.],
+            [4., 3., 2., 1.],
+        ]);
+
+        assert_eq!(a, b);
+        assert_ne!(a, c);
+    }
+
+    #[test]
+    fn multiply_two_matrices() {
+        let a = Matrix::from([
+            [1., 2., 3., 4.],
+            [5., 6., 7., 8.],
+            [9., 8., 7., 6.],
+            [5., 4., 3., 2.],
+        ]);
+
+        let b = Matrix::from([
+            [-2., 1., 2., 3.],
+            [3., 2., 1., -1.],
+            [4., 3., 6., 5.],
+            [1., 2., 7., 8.],
+        ]);
+
+        let expected = Matrix::from([
+            [20., 22., 50., 48.],
+            [44., 54., 114., 108.],
+            [40., 58., 110., 102.],
+            [16., 26., 46., 42.],
+        ]);
+
+        assert_eq!(a * b, expected);
+    }
+
+    #[test]
+    fn multiply_matrix_by_tuple() {
+        let a = Matrix::from([
+            [1., 2., 3., 4.],
+            [2., 4., 4., 2.],
+            [8., 6., 4., 1.],
+            [0., 0., 0., 1.],
+        ]);
+        let b = Tuple(1., 2., 3., 1.);
+        assert_eq!(a * b, Tuple(18., 24., 33., 1.));
+    }
+
+    #[test]
+    fn identity_matrix() {
+        let a = Matrix::from([
+            [0., 1., 2., 4.],
+            [1., 2., 4., 8.],
+            [2., 4., 8., 16.],
+            [4., 8., 16., 32.],
+        ]);
+
+        assert_eq!(&a * Matrix::identity(), a);
+    }
+
+    #[test]
+    fn transpose_a_matrix() {
+        let a = Matrix::from([
+            [0., 9., 3., 0.],
+            [9., 8., 0., 8.],
+            [1., 8., 5., 3.],
+            [0., 0., 5., 8.],
+        ]);
+
+        let expected = Matrix::from([
+            [0., 9., 1., 0.],
+            [9., 8., 8., 0.],
+            [3., 0., 5., 5.],
+            [0., 8., 3., 8.],
+        ]);
+
+        assert_eq!(a.transpose(), expected);
+    }
+
+    #[test]
+    fn calculates_2x2_determinant() {
+        let m = Matrix::from([[1., 5.], [-3., 2.]]);
+        assert_eq!(m.determinant(), 17.);
+    }
+
+    #[test]
+    fn calculates_3x3_determinant() {
+        let m = Matrix::from([[1., 2., 6.], [-5., 8., -4.], [2., 6., 4.]]);
+        assert_eq!(m.cofactor(0, 0), 56.);
+        assert_eq!(m.cofactor(0, 1), 12.);
+        assert_eq!(m.cofactor(0, 2), -46.);
+        assert_eq!(m.determinant(), -196.);
+    }
+
+    #[test]
+    fn calculates_4x4_determinant() {
+        let m = Matrix::from([
+            [-2., -8., 3., 5.],
+            [-3., 1., 7., 3.],
+            [1., 2., -9., 6.],
+            [-6., 7., 7., -9.],
+        ]);
+        assert_eq!(m.cofactor(0, 0), 690.);
+        assert_eq!(m.cofactor(0, 1), 447.);
+        assert_eq!(m.cofactor(0, 2), 210.);
+        assert_eq!(m.cofactor(0, 3), 51.);
+        assert_eq!(m.determinant(), -4071.);
+    }
+
+    #[test]
+    fn calculates_submatrix() {
+        let m = Matrix::from([[1., 5., 0.], [-3., 2., 7.], [0., 6., -3.]]);
+        let expected = Matrix::from([[-3., 2.], [0., 6.]]);
+        assert_eq!(m.submatrix(0, 2), expected);
+
+        let m = Matrix::from([
+            [-6., 1., 1., 6.],
+            [-8., 5., 8., 6.],
+            [-1., 0., 8., 2.],
+            [-7., 1., -1., 1.],
+        ]);
+        let expected = Matrix::from([[-6., 1., 6.], [-8., 8., 6.], [-7., -1., 1.]]);
+        assert_eq!(m.submatrix(2, 1), expected);
+    }
+
+    #[test]
+    fn calculates_minors_and_cofactors() {
+        let m = Matrix::from([[3., 5., 0.], [2., -1., -7.], [6., -1., 5.]]);
+        assert_eq!(m.submatrix(1, 0).determinant(), 25.);
+        assert_eq!(m.minor(0, 0), -12.);
+        assert_eq!(m.cofactor(0, 0), -12.);
+        assert_eq!(m.minor(1, 0), 25.);
+        assert_eq!(m.cofactor(1, 0), -25.);
+    }
+
+    #[test]
+    fn invert_4x4_matrix() {
+        let m = Matrix::from([
+            [6., 4., 4., 4.],
+            [5., 5., 7., 6.],
+            [4., -9., 3., -7.],
+            [9., 1., 7., -6.],
+        ]);
+
+        assert_eq!(m.determinant(), -2120.);
+        assert!(m.is_invertible());
+
+        let m = Matrix::from([
+            [-4., 2., -2., -3.],
+            [9., 6., 2., 6.],
+            [0., -5., 1., -5.],
+            [0., 0., 0., 0.],
+        ]);
+        assert_eq!(m.determinant(), 0.);
+        assert!(!m.is_invertible());
+
+        let m = Matrix::from([
+            [-5., 2., 6., -8.],
+            [1., -5., 1., 8.],
+            [7., 7., -6., -7.],
+            [1., -3., 7., 4.],
+        ]);
+        let expected = Matrix::from([
+            [0.21805, 0.45113, 0.24060, -0.04511],
+            [-0.80827, -1.45677, -0.44361, 0.52068],
+            [-0.07895, -0.22368, -0.05263, 0.19737],
+            [-0.52256, -0.81391, -0.30075, 0.30639],
+        ]);
+
+        assert_eq!(m.determinant(), 532.);
+        assert_eq!(m.cofactor(2, 3), -160.);
+        assert!(eq_f64(expected.cell(3, 2), -160. / 532.));
+        assert_eq!(m.cofactor(3, 2), 105.);
+        assert!(eq_f64(expected.cell(2, 3), 105. / 532.));
+
+        assert_eq!(m.inverse(), expected);
+
+        let m = Matrix::from([
+            [8., -5., 9., 2.],
+            [7., 5., 6., 1.],
+            [-6., 0., 9., 6.],
+            [-3., 0., -9., -4.],
+        ]);
+        let expected = Matrix::from([
+            [-0.15385, -0.15385, -0.28205, -0.53846],
+            [-0.07692, 0.12308, 0.02564, 0.03077],
+            [0.35897, 0.35897, 0.43590, 0.92308],
+            [-0.69231, -0.69231, -0.76923, -1.92308],
+        ]);
+        assert_eq!(m.inverse(), expected);
+
+        let m = Matrix::from([
+            [9., 3., 0., 9.],
+            [-5., -2., -6., -3.],
+            [-4., 9., 6., 4.],
+            [-7., 6., 6., 2.],
+        ]);
+        let expected = Matrix::from([
+            [-0.04074, -0.07778, 0.14444, -0.22222],
+            [-0.07778, 0.03333, 0.36667, -0.33333],
+            [-0.02901, -0.14630, -0.10926, 0.12963],
+            [0.17778, 0.06667, -0.26667, 0.33333],
+        ]);
+        assert_eq!(m.inverse(), expected);
+    }
+
+    #[test]
+    fn multiply_product_by_inverse() {
+        let a = Matrix::from([
+            [3., -9., 7., 3.],
+            [3., -8., 2., -9.],
+            [-4., 4., 4., 1.],
+            [-6., 5., -1., 1.],
+        ]);
+        let b = Matrix::from([
+            [8., 2., 2., 2.],
+            [3., -1., 7., 0.],
+            [7., 0., 5., 4.],
+            [6., -2., 0., 5.],
+        ]);
+        let c = &a * &b;
+        assert_eq!(c * b.inverse(), a);
+    }
+}

ray-tracer/src/types/mod.rs

@@ -0,0 +1,185 @@
+mod canvas;
+mod color;
+mod intersections;
+mod light;
+mod material;
+mod matrix;
+mod point;
+mod ray;
+mod sphere;
+mod tuple;
+mod vector;
+
+pub use canvas::Canvas;
+pub use color::Color;
+pub use intersections::{Intersections, Intersection};
+pub use light::PointLight;
+pub use material::Material;
+pub use matrix::Matrix;
+pub use point::Point;
+pub use ray::Ray;
+pub use sphere::Sphere;
+pub use tuple::Tuple;
+pub use vector::Vector;
+
+const EPSILON: f64 = 0.00001;
+
+fn eq_f64(l: f64, r: f64) -> bool {
+    (l - r).abs() < EPSILON
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn eq_f64_compares_values() {
+        assert!(eq_f64(1.0, 1.0));
+        assert!(eq_f64(0.9994, 0.9994));
+        assert!(eq_f64(0.9999994, 0.9999995));
+        assert!(!eq_f64(0.9995, 0.9994));
+    }
+
+    #[test]
+    fn add_two_tuples() {
+        let a = Tuple(3., -2., 5., 1.);
+        let b = Tuple(-2., 3., 1., 0.);
+        assert_eq!(a + b, Tuple(1., 1., 6., 1.));
+    }
+
+    #[test]
+    fn subtracts_two_tuples() {
+        let a = Tuple(3., 2., 1., 1.);
+        let b = Tuple(5., 6., 7., 1.);
+        assert_eq!(a - b, Tuple(-2., -4., -6., 0.));
+    }
+
+    #[test]
+    fn adds_point_and_vector() {
+        let a = Point::new(3., -2., 5.);
+        let b = Vector::new(-2., 3., 1.);
+        assert_eq!(a + b, Point::new(1., 1., 6.,));
+    }
+
+    #[test]
+    fn subtracts_two_points() {
+        let a = Point::new(3., 2., 1.);
+        let b = Point::new(5., 6., 7.);
+        assert_eq!(a - b, Vector::new(-2., -4., -6.,));
+    }
+
+    #[test]
+    fn subtracts_vector_from_point() {
+        let a = Point::new(3., 2., 1.);
+        let b = Vector::new(5., 6., 7.);
+        assert_eq!(a - b, Point::new(-2., -4., -6.));
+    }
+
+    #[test]
+    fn subtracts_two_vectors() {
+        let a = Vector::new(3., 2., 1.);
+        let b = Vector::new(5., 6., 7.);
+        assert_eq!(a - b, Vector::new(-2., -4., -6.));
+    }
+
+    #[test]
+    fn it_negates_primitives() {
+        assert_eq!(-Tuple(1., 2., 3., 4.), Tuple(-1., -2., -3., -4.),);
+    }
+
+    #[test]
+    fn multiply_tuple_by_scalar() {
+        assert_eq!(Tuple(1., -2., 3., -4.) * 3.5, Tuple(3.5, -7., 10.5, -14.));
+    }
+
+    #[test]
+    fn divide_tuple_by_scalar() {
+        assert_eq!(Tuple(1., -2., 3., -4.) / 2., Tuple(0.5, -1., 1.5, -2.));
+    }
+
+    #[test]
+    fn magnitude_of_vector() {
+        assert_eq!(Vector::new(1., 0., 0.).magnitude(), 1.);
+        assert_eq!(Vector::new(0., 1., 0.).magnitude(), 1.);
+        assert_eq!(Vector::new(0., 0., 1.).magnitude(), 1.);
+        assert_eq!(Vector::new(1., 2., 3.).magnitude(), 14_f64.sqrt());
+        assert_eq!(Vector::new(-1., -2., -3.).magnitude(), 14_f64.sqrt());
+    }
+
+    #[test]
+    fn normalize_vector() {
+        assert_eq!(Vector::new(4., 0., 0.).normalize(), Vector::new(1., 0., 0.));
+        assert_eq!(
+            Vector::new(1., 2., 3.).normalize(),
+            Vector::new(0.26726, 0.53452, 0.80178)
+        );
+        assert_eq!(Vector::new(1., 2., 3.).normalize().magnitude(), 1.);
+    }
+
+    #[test]
+    fn dot_product() {
+        assert_eq!(Vector::new(1., 2., 3.).dot(&Vector::new(2., 3., 4.)), 20.);
+    }
+
+    #[test]
+    fn cross_product() {
+        assert_eq!(
+            Vector::new(1., 2., 3.).cross(&Vector::new(2., 3., 4.)),
+            Vector::new(-1., 2., -1.)
+        );
+        assert_eq!(
+            Vector::new(2., 3., 4.).cross(&Vector::new(1., 2., 3.)),
+            Vector::new(1., -2., 1.)
+        );
+    }
+
+    #[test]
+    fn multiply_colors() {
+        let c1 = Color::new(1., 0.2, 0.4);
+        let c2 = Color::new(0.9, 1., 0.1);
+
+        assert_eq!(c1 * c2, Color::new(0.9, 0.2, 0.04));
+    }
+
+    #[test]
+    fn it_creates_a_canvas() {
+        let c = Canvas::new(10, 20);
+        assert_eq!(c.width(), 10);
+        assert_eq!(c.height(), 20);
+        for row in 0..20 {
+            for col in 0..10 {
+                assert_eq!(*c.pixel(row, col), Color::default());
+            }
+        }
+    }
+
+    #[test]
+    fn it_can_write_pixel() {
+        let mut c = Canvas::new(10, 20);
+        let red = Color::new(1., 0., 0.);
+        *c.pixel_mut(2, 3) = red;
+        assert_eq!(*c.pixel(2, 3), red);
+    }
+
+    #[test]
+    fn reflect_vector_approaching_at_45() {
+        let v = Vector::new(1., -1., 0.);
+        let n = Vector::new(0., 1., 0.);
+        assert_eq!(v.reflect(&n), Vector::new(1., 1., 0.));
+
+        let v = Vector::new(-1., -1., 0.);
+        let n = Vector::new(1., 0., 0.);
+        assert_eq!(v.reflect(&n), Vector::new(1., -1., 0.));
+
+        let v = Vector::new(1., 0., -1.);
+        let n = Vector::new(0., 0., 1.);
+        assert_eq!(v.reflect(&n), Vector::new(1., 0., 1.));
+    }
+
+    #[test]
+    fn reflect_off_slanted_surface() {
+        let v = Vector::new(0., -1., 0.);
+        let n = Vector::new(2_f64.sqrt() / 2., 2_f64.sqrt() / 2., 0.);
+        assert_eq!(v.reflect(&n), Vector::new(1., 0., 0.));
+    }
+}

ray-tracer/src/types/point.rs

@@ -0,0 +1,96 @@
+use crate::types::{Tuple, Vector};
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Point(Tuple);
+
+impl Point {
+    pub fn new(x: f64, y: f64, z: f64) -> Self {
+        Self(Tuple(x, y, z, 1.0))
+    }
+
+    #[inline]
+    pub fn x(&self) -> f64 {
+        self.0.0
+    }
+
+    #[inline]
+    pub fn y(&self) -> f64 {
+        self.0.1
+    }
+
+    #[inline]
+    pub fn z(&self) -> f64 {
+        self.0.2
+    }
+}
+
+impl Default for Point {
+    fn default() -> Self {
+        Self::new(0., 0., 0.)
+    }
+}
+
+impl From<Tuple> for Point {
+    fn from(tuple: Tuple) -> Self {
+        assert_eq!(tuple.3, 1.0);
+        Self(tuple)
+    }
+}
+
+impl std::ops::Deref for Point {
+    type Target = Tuple;
+    fn deref(&self) -> &Tuple {
+        &self.0
+    }
+}
+
+impl std::ops::Add<Vector> for Point {
+    type Output = Point;
+    fn add(self, r: Vector) -> Self::Output {
+        Point::from(self.0 + *r)
+    }
+}
+
+impl std::ops::Sub<&Point> for &Point {
+    type Output = Vector;
+    fn sub(self, r: &Point) -> Self::Output {
+        Vector::from(self.0 - r.0)
+    }
+}
+
+impl std::ops::Sub<Point> for &Point {
+    type Output = Vector;
+    fn sub(self, r: Point) -> Self::Output {
+        self - &r
+    }
+}
+
+impl std::ops::Sub<&Point> for Point {
+    type Output = Vector;
+    fn sub(self, r: &Point) -> Self::Output {
+        &self - r
+    }
+}
+
+impl std::ops::Sub<Point> for Point {
+    type Output = Vector;
+    fn sub(self, r: Point) -> Self::Output {
+        &self - &r
+    }
+}
+
+impl std::ops::Sub<Vector> for Point {
+    type Output = Point;
+    fn sub(self, r: Vector) -> Self::Output {
+        Point::from(self.0 - *r)
+    }
+}
+
+impl std::ops::Neg for Point {
+    type Output = Point;
+    fn neg(self) -> Self::Output {
+        let mut t = -self.0;
+        t.3 = 1.;
+        Point::from(t)
+    }
+}

ray-tracer/src/types/ray.rs

@@ -0,0 +1,191 @@
+use super::{Intersection, Intersections, Matrix, Point, Sphere, Vector};
+
+#[derive(Debug)]
+pub struct Ray {
+    origin: Point,
+    pub direction: Vector,
+}
+
+impl Ray {
+    pub fn new(origin: Point, direction: Vector) -> Self {
+        Self { origin, direction }
+    }
+
+    pub fn position(&self, t: f64) -> Point {
+        self.origin + self.direction * t
+    }
+
+    pub fn intersect<'a>(&self, s: &'a Sphere) -> Intersections<'a> {
+        let r2 = self.transform(s.transformation().inverse());
+        let sphere_to_ray = r2.origin - Point::new(0., 0., 0.);
+        let a = r2.direction.dot(&r2.direction);
+        let b = 2. * r2.direction.dot(&sphere_to_ray);
+        let c = sphere_to_ray.dot(&sphere_to_ray) - 1.;
+
+        let discriminant = b * b - 4. * a * c;
+        if discriminant < 0. {
+            return vec![].into();
+        }
+
+        let t1 = (-b - discriminant.sqrt()) / (2. * a);
+        let t2 = (-b + discriminant.sqrt()) / (2. * a);
+
+        vec![
+            Intersection { t: t1, object: s },
+            Intersection { t: t2, object: s },
+        ]
+        .into()
+    }
+
+    pub fn transform(&self, m: Matrix) -> Self {
+        Self {
+            origin: &m * self.origin,
+            direction: m * self.direction,
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::transforms::{scaling, translation};
+
+    #[test]
+    fn computing_point_from_distance() {
+        let r = Ray::new(Point::new(2., 3., 4.), Vector::new(1., 0., 0.));
+        assert_eq!(r.position(0.), Point::new(2., 3., 4.));
+        assert_eq!(r.position(1.), Point::new(3., 3., 4.));
+        assert_eq!(r.position(-1.), Point::new(1., 3., 4.));
+        assert_eq!(r.position(2.5), Point::new(4.5, 3., 4.));
+    }
+
+    #[test]
+    fn ray_intersects_sphere_at_two_points() {
+        let r = Ray::new(Point::new(0., 0., -5.), Vector::new(0., 0., 1.));
+        let s = Sphere::default();
+        let xs = r.intersect(&s);
+        assert_eq!(xs.len(), 2);
+        assert_eq!(xs[0].t, 4.);
+        assert_eq!(xs[1].t, 6.);
+        assert_eq!(*xs[0].object, s);
+        assert_eq!(*xs[1].object, s);
+    }
+
+    #[test]
+    fn ray_tangents_sphere_at_one_point() {
+        let r = Ray::new(Point::new(0., 1., -5.), Vector::new(0., 0., 1.));
+        let s = Sphere::default();
+        let xs = r.intersect(&s);
+        assert_eq!(xs.len(), 2);
+        assert_eq!(xs[0].t, 5.);
+        assert_eq!(xs[1].t, 5.);
+    }
+
+    #[test]
+    fn ray_misses_the_sphere() {
+        let r = Ray::new(Point::new(0., 2., -5.), Vector::new(0., 0., 1.));
+        let s = Sphere::default();
+        let xs = r.intersect(&s);
+        assert_eq!(xs.len(), 0);
+    }
+
+    #[test]
+    fn ray_originates_inside_the_sphere() {
+        let r = Ray::new(Point::new(0., 0., 0.), Vector::new(0., 0., 1.));
+        let s = Sphere::default();
+        let xs = r.intersect(&s);
+        assert_eq!(xs.len(), 2);
+        assert_eq!(xs[0].t, -1.);
+        assert_eq!(xs[1].t, 1.);
+    }
+
+    #[test]
+    fn sphere_is_behind_the_ray() {
+        let r = Ray::new(Point::new(0., 0., 5.), Vector::new(0., 0., 1.));
+        let s = Sphere::default();
+        let xs = r.intersect(&s);
+        assert_eq!(xs.len(), 2);
+        assert_eq!(xs[0].t, -6.);
+        assert_eq!(xs[1].t, -4.);
+    }
+
+    #[test]
+    fn hit_all_intersections_are_positive() {
+        let s = Sphere::default();
+        let i1 = Intersection { t: 1., object: &s };
+        let i2 = Intersection { t: 2., object: &s };
+        let xs = Intersections::from(vec![i1.clone(), i2.clone()]);
+
+        assert_eq!(xs.hit(), Some(&i1));
+    }
+
+    #[test]
+    fn hit_some_intersections_are_negative() {
+        let s = Sphere::default();
+        let i1 = Intersection { t: -1., object: &s };
+        let i2 = Intersection { t: 1., object: &s };
+        let xs = Intersections::from(vec![i1, i2.clone()]);
+
+        assert_eq!(xs.hit(), Some(&i2));
+    }
+
+    #[test]
+    fn hit_all_intersections_are_negative() {
+        let s = Sphere::default();
+        let i1 = Intersection { t: -2., object: &s };
+        let i2 = Intersection { t: -1., object: &s };
+        let xs = Intersections::from(vec![i1, i2]);
+
+        assert_eq!(xs.hit(), None);
+    }
+
+    #[test]
+    fn hit_is_always_lowest_nonnegative() {
+        let s = Sphere::default();
+        let i1 = Intersection { t: 5., object: &s };
+        let i2 = Intersection { t: 7., object: &s };
+        let i3 = Intersection { t: -3., object: &s };
+        let i4 = Intersection { t: 2., object: &s };
+        let xs = Intersections::from(vec![i1, i2, i3, i4.clone()]);
+
+        assert_eq!(xs.hit(), Some(&i4));
+    }
+
+    #[test]
+    fn translate_a_ray() {
+        let r = Ray::new(Point::new(1., 2., 3.), Vector::new(0., 1., 0.));
+        let m = translation(3., 4., 5.);
+        let r2 = r.transform(m);
+        assert_eq!(r2.origin, Point::new(4., 6., 8.));
+        assert_eq!(r2.direction, Vector::new(0., 1., 0.));
+    }
+
+    #[test]
+    fn scale_a_ray() {
+        let r = Ray::new(Point::new(1., 2., 3.), Vector::new(0., 1., 0.));
+        let m = scaling(2., 3., 4.);
+        let r2 = r.transform(m);
+        assert_eq!(r2.origin, Point::new(2., 6., 12.));
+        assert_eq!(r2.direction, Vector::new(0., 3., 0.,));
+    }
+
+    #[test]
+    fn intersect_scaled_sphere_with_ray() {
+        let r = Ray::new(Point::new(0., 0., -5.), Vector::new(0., 0., 1.));
+        let mut s = Sphere::default();
+        *s.transformation_mut() = scaling(2., 2., 2.);
+        let xs = r.intersect(&s);
+        assert_eq!(xs.len(), 2);
+        assert_eq!(xs[0].t, 3.);
+        assert_eq!(xs[1].t, 7.);
+    }
+
+    #[test]
+    fn intersect_translated_sphere_with_ray() {
+        let r = Ray::new(Point::new(0., 0., -5.), Vector::new(0., 0., 1.));
+        let mut s = Sphere::default();
+        *s.transformation_mut() = translation(5., 0., 0.);
+        let xs = r.intersect(&s);
+        assert_eq!(xs.len(), 0);
+    }
+}

ray-tracer/src/types/sphere.rs

@@ -0,0 +1,116 @@
+use super::{Matrix, Point, Vector, Material };
+
+#[derive(Debug, PartialEq)]
+pub struct Sphere {
+    origin: Point,
+    transformation: Matrix,
+    material: Material,
+}
+
+impl Sphere {
+    pub fn transformation(&self) -> &Matrix {
+        &self.transformation
+    }
+
+    pub fn transformation_mut(&mut self) -> &mut Matrix{
+        &mut self.transformation
+    }
+
+    pub fn material(&self) -> &Material {
+        &self.material
+    }
+
+    pub fn material_mut(&mut self) -> &mut Material {
+        &mut self.material
+    }
+
+    pub fn normal_at(&self, world_point: &Point) -> Vector {
+        let inverted_transform = self.transformation.inverse();
+        let object_point = &inverted_transform * world_point;
+        let object_normal = (object_point - Point::new(0., 0., 0.)).normalize();
+        let mut world_normal = inverted_transform.transpose() * *object_normal;
+        world_normal.3 = 0.;
+        Vector::from(world_normal).normalize()
+    }
+}
+
+impl Default for Sphere {
+    fn default() -> Self {
+        Self {
+            origin: Point::new(0., 0., 0.),
+            transformation: Matrix::identity(),
+            material: Default::default(),
+        }
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use crate::transforms::{rotation_z, scaling, translation};
+
+    #[test]
+    fn sphere_has_default_transformation() {
+        let s = Sphere::default();
+        assert_eq!(s.transformation(), &Matrix::identity());
+    }
+
+    #[test]
+    fn change_a_spheres_transformation() {
+        let mut s = Sphere::default();
+        let t = translation(2., 3., 4.);
+        *s.transformation_mut() = t.clone();
+        assert_eq!(*s.transformation(), t);
+    }
+
+    #[test]
+    fn normal_of_sphere_on_x() {
+        let s = Sphere::default();
+        let n = s.normal_at(&Point::new(1., 0., 0.));
+        assert_eq!(n, Vector::new(1., 0., 0.));
+    }
+
+    #[test]
+    fn normal_of_sphere_on_y() {
+        let s = Sphere::default();
+        let n = s.normal_at(&Point::new(0., 1., 0.));
+        assert_eq!(n, Vector::new(0., 1., 0.));
+    }
+
+    #[test]
+    fn normal_of_sphere_on_z() {
+        let s = Sphere::default();
+        let n = s.normal_at(&Point::new(0., 0., 1.));
+        assert_eq!(n, Vector::new(0., 0., 1.));
+    }
+
+    #[test]
+    fn normal_of_sphere_on_nonaxial() {
+        let s = Sphere::default();
+        let n = s.normal_at(&Point::new(3_f64.sqrt() / 3., 3_f64.sqrt() / 3., 3_f64.sqrt() / 3.));
+        assert_eq!(n, Vector::new(3_f64.sqrt() / 3., 3_f64.sqrt() / 3., 3_f64.sqrt() / 3.));
+    }
+
+    #[test]
+    fn normal_is_normalized() {
+        let s = Sphere::default();
+        let n = s.normal_at(&Point::new(3_f64.sqrt() / 3., 3_f64.sqrt() / 3., 3_f64.sqrt() / 3.));
+        assert_eq!(n.normalize(), n);
+    }
+
+    #[test]
+    fn compute_normal_of_translated_sphere() {
+        let mut s = Sphere::default();
+        *s.transformation_mut() = translation(0., 1., 0.);
+        let n = s.normal_at(&Point::new(0., 1.70711, -0.70711));
+        assert_eq!(n, Vector::new(0., 0.70711, -0.70711));
+    }
+
+    #[test]
+    fn compute_normal_of_transformed_sphere() {
+        let mut s = Sphere::default();
+        *s.transformation_mut() = scaling(1., 0.5, 1.) * rotation_z(std::f64::consts::PI / 5.);
+        let n = s.normal_at(&Point::new(0., 2_f64.sqrt() / 2., -2_f64.sqrt() / 2.));
+        assert_eq!(n, Vector::new(0., 0.97014, -0.24254));
+    }
+}

ray-tracer/src/types/tuple.rs

@@ -0,0 +1,110 @@
+use crate::types::eq_f64;
+
+#[derive(Debug, Clone, Copy)]
+pub struct Tuple(
+    pub f64, // x or red
+    pub f64, // y or green
+    pub f64, // z or blue
+    pub f64, // w, the flag which
+             // indicates point vs vec, or alpha
+);
+
+impl Tuple {
+    pub fn dot(&self, r: &Tuple) -> f64 {
+        self.0 * r.0 + self.1 * r.1 + self.2 * r.2 + self.3 * r.3
+    }
+}
+
+impl From<[f64; 4]> for Tuple {
+    fn from(source: [f64; 4]) -> Self {
+        Self(source[0], source[1], source[2], source[3])
+    }
+}
+
+impl PartialEq for Tuple {
+    fn eq(&self, r: &Tuple) -> bool {
+        eq_f64(self.0, r.0) && eq_f64(self.1, r.1) && eq_f64(self.2, r.2) && eq_f64(self.3, r.3)
+    }
+}
+
+impl std::ops::Add<&Tuple> for &Tuple {
+    type Output = Tuple;
+    fn add(self, r: &Tuple) -> Self::Output {
+        Tuple(self.0 + r.0, self.1 + r.1, self.2 + r.2, self.3 + r.3)
+    }
+}
+
+impl std::ops::Add<Tuple> for Tuple {
+    type Output = Tuple;
+    fn add(self, r: Tuple) -> Self::Output {
+        &self + &r
+    }
+}
+
+impl std::ops::Sub<&Tuple> for &Tuple {
+    type Output = Tuple;
+    fn sub(self, r: &Tuple) -> Self::Output {
+        Tuple(self.0 - r.0, self.1 - r.1, self.2 - r.2, self.3 - r.3)
+    }
+}
+
+impl std::ops::Sub<Tuple> for Tuple {
+    type Output = Tuple;
+    fn sub(self, r: Tuple) -> Self::Output {
+        &self - &r
+    }
+}
+
+impl std::ops::Neg for &Tuple {
+    type Output = Tuple;
+    fn neg(self) -> Self::Output {
+        Tuple(-self.0, -self.1, -self.2, -self.3)
+    }
+}
+
+impl std::ops::Neg for Tuple {
+    type Output = Tuple;
+    fn neg(self) -> Self::Output {
+        -&self
+    }
+}
+
+impl std::ops::Mul<f64> for &Tuple {
+    type Output = Tuple;
+    fn mul(self, scalar: f64) -> Self::Output {
+        Tuple(
+            self.0 * scalar,
+            self.1 * scalar,
+            self.2 * scalar,
+            self.3 * scalar,
+        )
+    }
+}
+
+impl std::ops::Mul<f64> for Tuple {
+    type Output = Tuple;
+    #[allow(clippy::op_ref)]
+    fn mul(self, scalar: f64) -> Self::Output {
+        &self * scalar
+    }
+}
+
+impl std::ops::Div<f64> for &Tuple {
+    type Output = Tuple;
+    fn div(self, scalar: f64) -> Self::Output {
+        Tuple(
+            self.0 / scalar,
+            self.1 / scalar,
+            self.2 / scalar,
+            self.3 / scalar,
+        )
+    }
+}
+
+impl std::ops::Div<f64> for Tuple {
+    type Output = Tuple;
+    #[allow(clippy::op_ref)]
+    fn div(self, scalar: f64) -> Self::Output {
+        &self / scalar
+    }
+}

ray-tracer/src/types/vector.rs

@@ -0,0 +1,125 @@
+use crate::types::Tuple;
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Vector(Tuple);
+
+impl Vector {
+    pub fn new(x: f64, y: f64, z: f64) -> Self {
+        Self(Tuple(x, y, z, 0.0))
+    }
+
+    #[inline]
+    pub fn x(&self) -> f64 {
+        self.0.0
+    }
+
+    #[inline]
+    pub fn y(&self) -> f64 {
+        self.0.1
+    }
+
+    #[inline]
+    pub fn z(&self) -> f64 {
+        self.0.2
+    }
+
+    pub fn magnitude(&self) -> f64 {
+        (self.x() * self.x() + self.y() * self.y() + self.z() * self.z()).sqrt()
+    }
+
+    pub fn normalize(&self) -> Self {
+        let mag = self.magnitude();
+        Self::new(self.x() / mag, self.y() / mag, self.z() / mag)
+    }
+
+    pub fn dot(&self, r: &Vector) -> f64 {
+        self.0.dot(r)
+    }
+
+    pub fn cross(&self, r: &Vector) -> Self {
+        let x = self.y() * r.z() - self.z() * r.y();
+        let y = self.z() * r.x() - self.x() * r.z();
+        let z = self.x() * r.y() - self.y() * r.x();
+        Self::new(x, y, z)
+    }
+
+    pub fn reflect(&self, n: &Vector) -> Self {
+        self - (n * 2. * self.dot(n))
+    }
+}
+
+impl Default for Vector {
+    fn default() -> Self {
+        Self::new(0., 0., 0.)
+    }
+}
+
+impl From<Tuple> for Vector {
+    fn from(tuple: Tuple) -> Self {
+        assert_eq!(tuple.3, 0.0);
+        Self(tuple)
+    }
+}
+
+impl std::ops::Deref for Vector {
+    type Target = Tuple;
+    fn deref(&self) -> &Tuple {
+        &self.0
+    }
+}
+
+impl std::ops::Add for Vector {
+    type Output = Vector;
+    fn add(self, r: Self) -> Self::Output {
+        Vector::from(self.0 + r.0)
+    }
+}
+
+impl std::ops::Sub<&Vector> for &Vector {
+    type Output = Vector;
+    fn sub(self, r: &Vector) -> Self::Output {
+        Vector::from(self.0 - r.0)
+    }
+}
+
+impl std::ops::Sub<Vector> for &Vector {
+    type Output = Vector;
+    fn sub(self, r: Vector) -> Self::Output {
+        self - &r
+    }
+}
+
+impl std::ops::Sub<&Vector> for Vector {
+    type Output = Vector;
+    fn sub(self, r: &Vector) -> Self::Output {
+        &self - r
+    }
+}
+
+impl std::ops::Sub for Vector {
+    type Output = Vector;
+    fn sub(self, r: Vector) -> Self::Output {
+        &self - &r
+    }
+}
+
+impl std::ops::Neg for Vector {
+    type Output = Vector;
+    fn neg(self) -> Self::Output {
+        Vector::from(-self.0)
+    }
+}
+
+impl std::ops::Mul<f64> for &Vector {
+    type Output = Vector;
+    fn mul(self, r: f64) -> Self::Output {
+        Vector::from(self.0 * r)
+    }
+}
+
+impl std::ops::Mul<f64> for Vector {
+    type Output = Vector;
+    fn mul(self, r: f64) -> Self {
+        Vector::from(self.0 * r)
+    }
+}