diff --git a/assignment-1c/Cargo.lock b/assignment-1c/Cargo.lock
index e8df05c..f62ab7a 100644
--- a/assignment-1c/Cargo.lock
+++ b/assignment-1c/Cargo.lock
@@ -40,6 +40,7 @@ name = "assignment-1c"
 version = "0.1.0"
 dependencies = [
  "anyhow",
+ "approx",
  "base64",
  "clap",
  "derivative",
@@ -50,6 +51,8 @@ dependencies = [
  "ordered-float",
  "rand",
  "rayon",
+ "rug",
+ "simba",
  "tracing",
  "tracing-subscriber",
 ]
@@ -60,6 +63,12 @@ version = "1.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa"
 
+[[package]]
+name = "az"
+version = "1.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7b7e4c2464d97fe331d41de9d5db0def0a96f4d823b8b32a2efd503578988973"
+
 [[package]]
 name = "backtrace"
 version = "0.3.67"
@@ -253,6 +262,16 @@ version = "0.27.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "ad0a93d233ebf96623465aad4046a8d3aa4da22d4f4beba5388838c8a434bbb4"
 
+[[package]]
+name = "gmp-mpfr-sys"
+version = "1.5.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a32868092c26fb25bb33c5ca7d8a937647979dfaa12f1f4f464beb57d726662c"
+dependencies = [
+ "libc",
+ "windows-sys",
+]
+
 [[package]]
 name = "heck"
 version = "0.4.1"
@@ -639,6 +658,19 @@ dependencies = [
  "num_cpus",
 ]
 
+[[package]]
+name = "rug"
+version = "1.19.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a465f6576b9f0844bd35749197576d68e3db169120532c4e0f868ecccad3d449"
+dependencies = [
+ "az",
+ "gmp-mpfr-sys",
+ "libc",
+ "num-integer",
+ "num-traits",
+]
+
 [[package]]
 name = "rustc-demangle"
 version = "0.1.21"
diff --git a/assignment-1c/Cargo.toml b/assignment-1c/Cargo.toml
index 2930694..aea7329 100644
--- a/assignment-1c/Cargo.toml
+++ b/assignment-1c/Cargo.toml
@@ -20,6 +20,7 @@ path = "src/main.rs"
 
 [dependencies]
 anyhow = { version = "1.0.68", features = ["backtrace"] }
+approx = "0.5.1"
 base64 = "0.21.0"
 clap = { version = "4.1.4", features = ["cargo", "derive"] }
 derivative = "2.2.0"
@@ -30,5 +31,7 @@ num = { version = "0.4.0", features = ["serde"] }
 ordered-float = "3.4.0"
 rand = "0.8.5"
 rayon = "1.6.1"
+rug = { version = "1.19.1", features = ["num-traits"] }
+simba = "0.8.0"
 tracing = "0.1.37"
 tracing-subscriber = "0.3.16"
diff --git a/assignment-1c/raytracer1b b/assignment-1c/raytracer1b
deleted file mode 100755
index d653239..000000000
Binary files a/assignment-1c/raytracer1b and /dev/null differ
diff --git a/assignment-1c/src/float.rs b/assignment-1c/src/float.rs
new file mode 100644
index 000000000..d6ebd4f
--- /dev/null
+++ b/assignment-1c/src/float.rs
@@ -0,0 +1,351 @@
+use std::{
+  fmt,
+  ops::{
+    Add, AddAssign, Deref, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign,
+    Sub, SubAssign,
+  },
+};
+
+use nalgebra::{ComplexField, Field, SimdValue, RealField};
+use num::{traits::Pow, FromPrimitive, Num, One, Zero, Signed};
+use rug::float::{ParseFloatError, Special};
+use simba::scalar::SubsetOf;
+
+#[derive(Debug, Clone, PartialEq)]
+pub struct Float(rug::Float);
+
+impl From<rug::Float> for Float {
+  fn from(f: rug::Float) -> Self { Float(f) }
+}
+
+impl PartialOrd for Float {
+  fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
+    self.0.as_ord().partial_cmp(other.0.as_ord())
+  }
+}
+
+impl fmt::Display for Float {
+  fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { self.0.fmt(f) }
+}
+
+/*
+impl Deref for Float {
+  type Target = rug::Float;
+  fn deref(&self) -> &Self::Target {
+    &self.0
+  }
+}
+*/
+
+macro_rules! impl_ops {
+  ($op_name:ty,
+   $fn_name:ident) => {
+    impl $op_name for Float {
+      type Output = Float;
+      fn $fn_name(self, rhs: Self) -> Self::Output {
+        Float(self.0.$fn_name(rhs.0))
+      }
+    }
+  };
+}
+
+impl_ops!(Add, add);
+impl_ops!(Sub, sub);
+impl_ops!(Mul, mul);
+impl_ops!(Rem, rem);
+impl_ops!(Div, div);
+impl_ops!(Pow<Float>, pow);
+
+macro_rules! impl_assign_ops {
+  ($op_name:ident,
+   $fn_name:ident) => {
+    impl $op_name for Float {
+      fn $fn_name(&mut self, rhs: Self) { self.0.$fn_name(rhs.0); }
+    }
+  };
+}
+
+impl_assign_ops!(AddAssign, add_assign);
+impl_assign_ops!(SubAssign, sub_assign);
+impl_assign_ops!(MulAssign, mul_assign);
+impl_assign_ops!(RemAssign, rem_assign);
+impl_assign_ops!(DivAssign, div_assign);
+
+impl SimdValue for Float {
+  type Element = Float;
+  type SimdBool = bool;
+
+  #[inline(always)]
+  fn lanes() -> usize { 1 }
+
+  #[inline(always)]
+  fn splat(val: Self::Element) -> Self { val }
+
+  #[inline(always)]
+  fn extract(&self, _: usize) -> Self::Element { *self }
+
+  #[inline(always)]
+  unsafe fn extract_unchecked(&self, _: usize) -> Self::Element { *self }
+
+  #[inline(always)]
+  fn replace(&mut self, _: usize, val: Self::Element) { *self = val }
+
+  #[inline(always)]
+  unsafe fn replace_unchecked(&mut self, _: usize, val: Self::Element) {
+    *self = val
+  }
+
+  #[inline(always)]
+  fn select(self, cond: Self::SimdBool, other: Self) -> Self {
+    if cond {
+      self
+    } else {
+      other
+    }
+  }
+}
+
+impl Zero for Float {
+  fn zero() -> Self { Float(rug::Float::with_val(53, Special::Zero)) }
+
+  fn is_zero(&self) -> bool { self.0.is_zero() }
+}
+
+impl One for Float {
+  fn one() -> Self { Float(rug::Float::with_val(53, 1)) }
+}
+
+impl Neg for Float {
+  type Output = Float;
+  fn neg(self) -> Self::Output { Float(self.0.neg()) }
+}
+
+impl Num for Float {
+  type FromStrRadixErr = ParseFloatError;
+
+  fn from_str_radix(
+    str: &str,
+    radix: u32,
+  ) -> Result<Self, Self::FromStrRadixErr> {
+    let incomplete = rug::Float::parse_radix(str, radix as i32)?;
+    Ok(Float(rug::Float::with_val(53, incomplete)))
+  }
+}
+
+impl Field for Float {}
+
+impl SubsetOf<Float> for Float {
+  fn to_superset(&self) -> Float { self.clone() }
+  fn from_superset_unchecked(element: &Float) -> Self { element.clone() }
+  fn is_in_subset(element: &Float) -> bool { true }
+}
+
+impl SubsetOf<Float> for f64 {
+  fn to_superset(&self) -> Float { Float(rug::Float::with_val(53, self)) }
+  fn from_superset_unchecked(element: &Float) -> Self { element.0.to_f64() }
+  fn is_in_subset(element: &Float) -> bool { true }
+}
+
+impl FromPrimitive for Float {
+  fn from_i64(n: i64) -> Option<Self> { todo!() }
+
+  fn from_u64(n: u64) -> Option<Self> { todo!() }
+}
+
+impl Signed for Float {}
+
+impl RealField for Float {
+    fn is_sign_positive(&self) -> bool {
+        todo!()
+    }
+
+    fn is_sign_negative(&self) -> bool {
+        todo!()
+    }
+
+    fn copysign(self, sign: Self) -> Self {
+        todo!()
+    }
+
+    fn max(self, other: Self) -> Self {
+        todo!()
+    }
+
+    fn min(self, other: Self) -> Self {
+        todo!()
+    }
+
+    fn clamp(self, min: Self, max: Self) -> Self {
+        todo!()
+    }
+
+    fn atan2(self, other: Self) -> Self {
+        todo!()
+    }
+
+    fn min_value() -> Option<Self> {
+        todo!()
+    }
+
+    fn max_value() -> Option<Self> {
+        todo!()
+    }
+
+    fn pi() -> Self {
+        todo!()
+    }
+
+    fn two_pi() -> Self {
+        todo!()
+    }
+
+    fn frac_pi_2() -> Self {
+        todo!()
+    }
+
+    fn frac_pi_3() -> Self {
+        todo!()
+    }
+
+    fn frac_pi_4() -> Self {
+        todo!()
+    }
+
+    fn frac_pi_6() -> Self {
+        todo!()
+    }
+
+    fn frac_pi_8() -> Self {
+        todo!()
+    }
+
+    fn frac_1_pi() -> Self {
+        todo!()
+    }
+
+    fn frac_2_pi() -> Self {
+        todo!()
+    }
+
+    fn frac_2_sqrt_pi() -> Self {
+        todo!()
+    }
+
+    fn e() -> Self {
+        todo!()
+    }
+
+    fn log2_e() -> Self {
+        todo!()
+    }
+
+    fn log10_e() -> Self {
+        todo!()
+    }
+
+    fn ln_2() -> Self {
+        todo!()
+    }
+
+    fn ln_10() -> Self {
+        todo!()
+    }
+}
+
+impl ComplexField for Float {
+  type RealField = Float;
+
+  fn from_real(re: Self::RealField) -> Self { re }
+
+  fn real(self) -> Self::RealField { todo!() }
+
+  fn imaginary(self) -> Self::RealField { todo!() }
+
+  fn modulus(self) -> Self::RealField { todo!() }
+
+  fn modulus_squared(self) -> Self::RealField { todo!() }
+
+  fn argument(self) -> Self::RealField { todo!() }
+
+  fn norm1(self) -> Self::RealField { todo!() }
+
+  fn scale(self, factor: Self::RealField) -> Self { todo!() }
+
+  fn unscale(self, factor: Self::RealField) -> Self { todo!() }
+
+  fn floor(self) -> Self { todo!() }
+
+  fn ceil(self) -> Self { todo!() }
+
+  fn round(self) -> Self { todo!() }
+
+  fn trunc(self) -> Self { todo!() }
+
+  fn fract(self) -> Self { todo!() }
+
+  fn mul_add(self, a: Self, b: Self) -> Self { todo!() }
+
+  fn abs(self) -> Self::RealField { todo!() }
+
+  fn hypot(self, other: Self) -> Self::RealField { todo!() }
+
+  fn recip(self) -> Self { todo!() }
+
+  fn conjugate(self) -> Self { todo!() }
+
+  fn sin(self) -> Self { todo!() }
+
+  fn cos(self) -> Self { todo!() }
+
+  fn sin_cos(self) -> (Self, Self) { todo!() }
+
+  fn tan(self) -> Self { todo!() }
+
+  fn asin(self) -> Self { todo!() }
+
+  fn acos(self) -> Self { todo!() }
+
+  fn atan(self) -> Self { todo!() }
+
+  fn sinh(self) -> Self { todo!() }
+
+  fn cosh(self) -> Self { todo!() }
+
+  fn tanh(self) -> Self { todo!() }
+
+  fn asinh(self) -> Self { todo!() }
+
+  fn acosh(self) -> Self { todo!() }
+
+  fn atanh(self) -> Self { todo!() }
+
+  fn log(self, base: Self::RealField) -> Self { todo!() }
+
+  fn log2(self) -> Self { todo!() }
+
+  fn log10(self) -> Self { todo!() }
+
+  fn ln(self) -> Self { todo!() }
+
+  fn ln_1p(self) -> Self { todo!() }
+
+  fn sqrt(self) -> Self { Float(self.0.sqrt()) }
+
+  fn exp(self) -> Self { todo!() }
+
+  fn exp2(self) -> Self { todo!() }
+
+  fn exp_m1(self) -> Self { todo!() }
+
+  fn powi(self, n: i32) -> Self { Float(self.0.pow(n)) }
+
+  fn powf(self, n: Self::RealField) -> Self { todo!() }
+
+  fn powc(self, n: Self) -> Self { todo!() }
+
+  fn cbrt(self) -> Self { todo!() }
+
+  fn is_finite(&self) -> bool { todo!() }
+
+  fn try_sqrt(self) -> Option<Self> { todo!() }
+}
diff --git a/assignment-1c/src/image.rs b/assignment-1c/src/image.rs
index c9d09ef..e1e90fa 100644
--- a/assignment-1c/src/image.rs
+++ b/assignment-1c/src/image.rs
@@ -9,8 +9,10 @@ use generator::{done, Gn};
 use itertools::Itertools;
 use nalgebra::Vector3;
 
+use crate::float::Float;
+
 /// A pixel color represented by a red, green, and blue value in the range 0-1.
-pub type Color = Vector3<f64>;
+pub type Color = Vector3<Float>;
 
 /// A representation of an image
 #[derive(Derivative)]
@@ -89,9 +91,9 @@ impl Image {
         None => bail!("Not enough elements"),
       };
 
-      let r = r? as f64 / max_value as f64;
-      let g = g? as f64 / max_value as f64;
-      let b = b? as f64 / max_value as f64;
+      let r = r? as Float / max_value as Float;
+      let g = g? as Float / max_value as Float;
+      let b = b? as Float / max_value as Float;
 
       let color = Color::new(r, g, b);
       data.push(color);
diff --git a/assignment-1c/src/lib.rs b/assignment-1c/src/lib.rs
index 2960125..ea6928b 100644
--- a/assignment-1c/src/lib.rs
+++ b/assignment-1c/src/lib.rs
@@ -1,5 +1,7 @@
 use nalgebra::{Vector2, Vector3};
 
+use crate::float::Float;
+
 #[macro_use]
 extern crate anyhow;
 #[macro_use]
@@ -7,11 +9,12 @@ extern crate derivative;
 #[macro_use]
 extern crate tracing;
 
+pub mod float;
 pub mod image;
 pub mod ray;
 pub mod scene;
 pub mod utils;
 
-pub type Point2 = Vector2<f64>;
-pub type Point = Vector3<f64>;
-pub type Vector = Vector3<f64>;
+pub type Point2 = Vector2<Float>;
+pub type Point = Vector3<Float>;
+pub type Vector = Vector3<Float>;
diff --git a/assignment-1c/src/ray.rs b/assignment-1c/src/ray.rs
index 5ae85e6..05c7fe5 100644
--- a/assignment-1c/src/ray.rs
+++ b/assignment-1c/src/ray.rs
@@ -1,3 +1,4 @@
+use crate::float::Float;
 use crate::{Point, Vector};
 
 /// A normalized parametric Ray of the form (origin + direction * time)
@@ -21,7 +22,7 @@ impl Ray {
   }
 
   /// Evaluate the ray at a certain point in time, yielding a point
-  pub fn eval(&self, time: f64) -> Point {
+  pub fn eval(&self, time: Float) -> Point {
     self.origin + self.direction * time
   }
 }
diff --git a/assignment-1c/src/scene/cylinder.rs b/assignment-1c/src/scene/cylinder.rs
index 18dce57..995030d 100644
--- a/assignment-1c/src/scene/cylinder.rs
+++ b/assignment-1c/src/scene/cylinder.rs
@@ -3,6 +3,7 @@ use anyhow::Result;
 use ordered_float::NotNan;
 
 use crate::utils::compute_rotation_matrix;
+use crate::float::Float;
 use crate::Vector;
 use crate::{ray::Ray, Point};
 
@@ -13,8 +14,8 @@ use super::Scene;
 pub struct Cylinder {
   pub center: Point,
   pub direction: Vector,
-  pub radius: f64,
-  pub length: f64,
+  pub radius: Float,
+  pub length: Float,
 }
 
 impl Cylinder {
diff --git a/assignment-1c/src/scene/data.rs b/assignment-1c/src/scene/data.rs
index 577b43c..213103d 100644
--- a/assignment-1c/src/scene/data.rs
+++ b/assignment-1c/src/scene/data.rs
@@ -2,6 +2,7 @@ use std::fmt::Debug;
 
 use crate::image::Color;
 use crate::utils::cross;
+use crate::float::Float;
 use crate::Point;
 
 use super::Scene;
@@ -19,10 +20,10 @@ pub struct Material {
   pub diffuse_color: Point,
   pub specular_color: Point,
 
-  pub k_a: f64,
-  pub k_d: f64,
-  pub k_s: f64,
-  pub exponent: f64,
+  pub k_a: Float,
+  pub k_d: Float,
+  pub k_s: Float,
+  pub exponent: Float,
 }
 
 #[derive(Debug)]
@@ -69,17 +70,17 @@ pub struct DepthCueing {
 
   /// Proportion of the color influenced by the depth tint when the distance is
   /// maxed (caps at 1.0)
-  pub a_max: f64,
+  pub a_max: Float,
 
   /// Proportion of the color influenced by the depth tint when the distance is
   /// at the minimum (caps at 1.0)
-  pub a_min: f64,
+  pub a_min: Float,
 
   /// The max distance that should be affected by the depth tint
-  pub dist_max: f64,
+  pub dist_max: Float,
 
   /// The min distance that should be affected by the depth tint
-  pub dist_min: f64,
+  pub dist_min: Float,
 }
 
 /// A default implementation here needs to simulate what would happen if there
@@ -101,9 +102,9 @@ impl Default for DepthCueing {
 /// Light attenuation dropoff coefficients
 #[derive(Debug)]
 pub struct Attenuation {
-  pub c1: f64,
-  pub c2: f64,
-  pub c3: f64,
+  pub c1: Float,
+  pub c2: Float,
+  pub c3: Float,
 }
 
 /// A default implementation here needs to simulate what would happen if there
@@ -122,7 +123,7 @@ impl Default for Attenuation {
 
 impl Scene {
   /// Determine the boundaries of the viewing window in world coordinates
-  pub fn compute_viewing_window(&self, distance: f64) -> Rect {
+  pub fn compute_viewing_window(&self, distance: Float) -> Rect {
     // Compute viewing directions
     let u = cross(self.view_dir, self.up_dir).normalize();
     let v = cross(u, self.view_dir).normalize();
@@ -140,7 +141,7 @@ impl Scene {
       w_over_2d * 2.0 * distance
     };
 
-    let aspect_ratio = self.image_width as f64 / self.image_height as f64;
+    let aspect_ratio = self.image_width as Float / self.image_height as Float;
     let viewing_height = viewing_width / aspect_ratio;
 
     // Compute viewing window corners
@@ -161,20 +162,20 @@ impl Scene {
   /// current scene
   pub fn pixel_translation_function(
     &self,
-    distance: f64,
+    distance: Float,
   ) -> impl Fn(usize, usize) -> Point {
     let view_window = self.compute_viewing_window(distance);
 
     let dx = view_window.upper_right - view_window.upper_left;
-    let pixel_base_x = dx / self.image_width as f64;
+    let pixel_base_x = dx / self.image_width as Float;
 
     let dy = view_window.lower_left - view_window.upper_left;
-    let pixel_base_y = dy / self.image_height as f64;
+    let pixel_base_y = dy / self.image_height as Float;
 
     // The final function to be returned
     move |px: usize, py: usize| {
-      let x_component = pixel_base_x * px as f64;
-      let y_component = pixel_base_y * py as f64;
+      let x_component = pixel_base_x * px as Float;
+      let y_component = pixel_base_y * py as Float;
 
       // Without adding this, we would be getting the top-left of the pixel's
       // rectangle. We want the center, so add half of the pixel size as
diff --git a/assignment-1c/src/scene/input_file.rs b/assignment-1c/src/scene/input_file.rs
index 1b192cf..38ac4e6 100644
--- a/assignment-1c/src/scene/input_file.rs
+++ b/assignment-1c/src/scene/input_file.rs
@@ -6,6 +6,7 @@ use anyhow::{Context, Result};
 use itertools::Itertools;
 use nalgebra::Vector3;
 
+use crate::float::Float;
 use crate::{
   image::{Color, Image},
   scene::{
@@ -93,23 +94,23 @@ impl Scene {
           }
         }
 
-        "eye" => scene.eye_pos = r!(Vector3<f64>),
-        "viewdir" => scene.view_dir = r!(Vector3<f64>),
-        "updir" => scene.up_dir = r!(Vector3<f64>),
+        "eye" => scene.eye_pos = r!(Vector3<Float>),
+        "viewdir" => scene.view_dir = r!(Vector3<Float>),
+        "updir" => scene.up_dir = r!(Vector3<Float>),
 
-        "hfov" => scene.hfov = r!(f64),
+        "hfov" => scene.hfov = r!(Float),
         "bkgcolor" => scene.bkg_color = r!(Color),
 
         // light x y z w r g b
         // attlight x y z w r g b c1 c2 c3
         "light" | "attlight" => {
-          let vec3 = r!(Vector3<f64>);
+          let vec3 = r!(Vector3<Float>);
           let w = r!(usize);
           let color = r!(Color);
 
           let attenuation = match keyword == "attlight" {
             true => {
-              let c = r!(Vector3<f64>);
+              let c = r!(Vector3<Float>);
               Some(Attenuation {
                 c1: c.x,
                 c2: c.y,
@@ -135,10 +136,10 @@ impl Scene {
         // depthcueing dcr dcg dcb amax amin distmax distmin
         "depthcueing" => {
           let color = r!(Color);
-          let a_max = r!(f64);
-          let a_min = r!(f64);
-          let dist_max = r!(f64);
-          let dist_min = r!(f64);
+          let a_max = r!(Float);
+          let a_min = r!(Float);
+          let dist_max = r!(Float);
+          let dist_min = r!(Float);
 
           scene.depth_cueing = DepthCueing {
             color,
@@ -153,10 +154,10 @@ impl Scene {
         "mtlcolor" => {
           let diffuse_color = r!(Color);
           let specular_color = r!(Color);
-          let k_a = r!(f64);
-          let k_d = r!(f64);
-          let k_s = r!(f64);
-          let exponent = r!(f64);
+          let k_a = r!(Float);
+          let k_d = r!(Float);
+          let k_s = r!(Float);
+          let exponent = r!(Float);
 
           let material = Material {
             diffuse_color,
@@ -174,7 +175,7 @@ impl Scene {
 
         "sphere" => {
           let center = r!(Point);
-          let radius = r!(f64);
+          let radius = r!(Float);
 
           scene.objects.push(Object {
             kind: ObjectKind::Sphere(Sphere { center, radius }),
@@ -186,8 +187,8 @@ impl Scene {
         "cylinder" => {
           let center = r!(Point);
           let direction = r!(Vector);
-          let radius = r!(f64);
-          let length = r!(f64);
+          let radius = r!(Float);
+          let length = r!(Float);
 
           scene.objects.push(Object {
             kind: ObjectKind::Cylinder(Cylinder {
@@ -329,10 +330,10 @@ macro_rules! impl_construct {
   };
 }
 
-impl_construct!(f64);
+impl_construct!(Float);
 impl_construct!(usize);
 
-impl Construct for Vector3<f64> {
+impl Construct for Vector3<Float> {
   type Args = ();
 
   fn construct<'a, I>(it: &mut I, _: Self::Args) -> Result<Self>
@@ -344,9 +345,9 @@ impl Construct for Vector3<f64> {
       None => bail!("Expected 3 values"),
     };
 
-    let x: f64 = x.parse()?;
-    let y: f64 = y.parse()?;
-    let z: f64 = z.parse()?;
+    let x: Float = x.parse()?;
+    let y: Float = y.parse()?;
+    let z: Float = z.parse()?;
 
     Ok(Vector3::new(x, y, z))
   }
diff --git a/assignment-1c/src/scene/sphere.rs b/assignment-1c/src/scene/sphere.rs
index 81f25a4..532d292 100644
--- a/assignment-1c/src/scene/sphere.rs
+++ b/assignment-1c/src/scene/sphere.rs
@@ -1,9 +1,10 @@
-use std::f64::consts::PI;
+use std::float::consts::PI;
 
 use anyhow::Result;
 use ordered_float::NotNan;
 
 use crate::{ray::Ray, utils::min_f64, Point};
+use crate::float::Float;
 
 use super::illumination::IntersectionContext;
 use super::Scene;
@@ -11,7 +12,7 @@ use super::Scene;
 #[derive(Debug)]
 pub struct Sphere {
   pub center: Point,
-  pub radius: f64,
+  pub radius: Float,
 }
 
 impl Sphere {
@@ -78,7 +79,7 @@ impl Sphere {
   pub fn get_texture_coord(
     &self,
     ctx: &IntersectionContext,
-  ) -> Result<(f64, f64)> {
+  ) -> Result<(Float, Float)> {
     // Reverse engineer the angles from the coordinate of the intersection
     let cosp = (ctx.point.z - self.center.z) / self.radius;
     let phi = cosp.acos();
diff --git a/assignment-1c/src/scene/triangle.rs b/assignment-1c/src/scene/triangle.rs
index b5c4c96..235c0ae 100644
--- a/assignment-1c/src/scene/triangle.rs
+++ b/assignment-1c/src/scene/triangle.rs
@@ -4,6 +4,7 @@ use anyhow::Result;
 use nalgebra::{Matrix2, Vector2, Vector3};
 use ordered_float::NotNan;
 
+use crate::float::Float;
 use crate::ray::Ray;
 use crate::utils::{cross, dot};
 use crate::{Point, Vector};
@@ -109,7 +110,7 @@ impl Triangle {
     &self,
     scene: &Scene,
     ctx: &IntersectionContext,
-  ) -> Result<(f64, f64)> {
+  ) -> Result<(Float, Float)> {
     let texture_coordinates = match self.textures {
       Some(v) => v,
       None => {
@@ -131,7 +132,7 @@ impl Triangle {
     Ok((u, v))
   }
 
-  fn convert_point(&self, scene: &Scene, point: Point) -> Vector2<f64> {
+  fn convert_point(&self, scene: &Scene, point: Point) -> Vector2<Float> {
     let (p0, e1, e2) = self.basis_vectors(scene);
 
     let ep = point - p0;
@@ -161,8 +162,8 @@ impl Triangle {
   fn compute_barycentric_coordinates(
     &self,
     scene: &Scene,
-    p: Vector2<f64>,
-  ) -> Result<(f64, f64, f64)> {
+    p: Vector2<Float>,
+  ) -> Result<(Float, Float, Float)> {
     let (_, e1, e2) = self.basis_vectors(scene);
     let d = Matrix2::new(dot(e1, e1), dot(e1, e2), dot(e2, e1), dot(e2, e2));
 
diff --git a/assignment-1c/src/utils.rs b/assignment-1c/src/utils.rs
index 68b211e..0822160 100644
--- a/assignment-1c/src/utils.rs
+++ b/assignment-1c/src/utils.rs
@@ -1,38 +1,38 @@
 use anyhow::Result;
-use nalgebra::{Matrix3, Vector3};
-use ordered_float::NotNan;
+use nalgebra::{ComplexField, Matrix3, Vector3};
+use num::{Zero, One};
+use num::traits::Pow;
 
-use crate::{Vector};
+use crate::float::Float;
+use crate::Vector;
 
+/*
 /// Finds the minimum of an iterator of f64s, ignoring any NaN values
 #[inline]
-pub fn min_f64<I>(i: I) -> Option<f64>
+pub fn min_f64<I>(i: I) -> Option<Float>
 where
-  I: Iterator<Item = f64>,
+  I: Iterator<Item = Float>,
 {
-  i.filter_map(|i| NotNan::new(i).ok())
-    .min()
-    .map(|i| i.into_inner())
+  i.map(|v| v.as_ord()).min().map(|i| i.as_float().clone())
 }
 
 /// Finds the minimum of an iterator of f64s using the given predicate, ignoring
 /// any NaN values
 #[inline]
-pub fn min_f64_by_key<I, F>(i: I, f: F) -> Option<f64>
+pub fn min_f64_by_key<I, F>(i: I, f: F) -> Option<Float>
 where
-  I: Iterator<Item = f64>,
-  F: FnMut(&NotNan<f64>),
+  I: Iterator<Item = Float>,
+  F: FnMut(&OrdFloat),
 {
-  i.filter_map(|i| NotNan::new(i).ok())
+  i.map(|v| v.as_ord().clone())
     .min_by_key(f)
-    .map(|i| i.into_inner())
+    .map(|i| i.as_float().clone())
 }
+*/
 
 /// Dot-product between two 3D vectors.
 #[inline]
-pub fn dot(a: Vector, b: Vector) -> f64 {
-  a.x * b.x + a.y * b.y + a.z * b.z
-}
+pub fn dot(a: Vector, b: Vector) -> Float { a.x * b.x + a.y * b.y + a.z * b.z }
 
 /// Cross-product between two 3D vectors.
 #[inline]
@@ -43,25 +43,42 @@ pub fn cross(a: Vector, b: Vector) -> Vector {
   Vector::new(x, y, z)
 }
 
+/// L2-norm
+#[inline]
+pub fn norm(a: Vector) -> Float {
+  let sum: Float = a.x.pow(2) + a.y.pow(2) + a.z.pow(2);
+  sum.sqrt().into()
+}
+
+/// Normalize a vector
+#[inline]
+pub fn normalize(a: Vector) -> Vector {
+  let norm = norm(a);
+  Vector::new(a.x / norm, a.y / norm, a.z / norm)
+}
+
 /// Calculate the rotation matrix between the 2 given vectors
 ///
 /// Based on the method given [here][1].
 ///
 /// [1]: https://math.stackexchange.com/a/897677
 pub fn compute_rotation_matrix(
-  a: Vector3<f64>,
-  b: Vector3<f64>,
-) -> Result<Matrix3<f64>> {
+  a: Vector3<Float>,
+  b: Vector3<Float>,
+) -> Result<Matrix3<Float>> {
   // Special case: if a and b are in the same direction, just return the
   // identity matrix.
-  if a.normalize() == b.normalize() {
+  if normalize(a) == normalize(b) {
     return Ok(Matrix3::identity());
   }
 
   let cos_t = dot(a, b);
-  let sin_t = cross(a, b).norm();
+  let sin_t = norm(cross(a, b));
 
-  let g = Matrix3::new(cos_t, -sin_t, 0.0, sin_t, cos_t, 0.0, 0.0, 0.0, 1.0);
+  let zero = Float::zero();
+  let one = Float::one();
+  let g =
+    Matrix3::new(cos_t, -sin_t, zero, sin_t, cos_t, zero, zero, zero, one);
 
   // New basis vectors
   let u = a;
diff --git a/flake.nix b/flake.nix
index ea8d474..2063ad7 100644
--- a/flake.nix
+++ b/flake.nix
@@ -26,6 +26,8 @@
             unzip
             zip
 
+            m4
+
             (python310.withPackages (p: with p; [ ipython numpy scipy sympy ]))
           ]) ++ (with toolchain; [
             cargo
diff --git a/rustfmt.toml b/rustfmt.toml
index a3be6a8..22f4704 100644
--- a/rustfmt.toml
+++ b/rustfmt.toml
@@ -1,3 +1,4 @@
 max_width = 80
 tab_spaces = 2
 wrap_comments = true
+fn_single_line = true