diff --git a/Cargo.toml b/Cargo.toml
index 6fdd605a..24a552a5 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -55,17 +55,33 @@ retrofire-core = { version = "0.4.0", path = "core" }
 retrofire-front = { version = "0.4.0", path = "front" }
 retrofire-geom = { version = "0.4.0", path = "geom" }
 
+[dev-dependencies]
+divan = "0.1.21"
+
 [profile.release]
 opt-level = 2
-codegen-units = 1
+codegen-units = 4
 lto = "thin"
 debug = 1
 
 [profile.bench]
 opt-level = 3
 codegen-units = 1
-lto = "fat"
+lto = "thin"
 
 [profile.dev]
 opt-level = 1
 split-debuginfo = "unpacked"
+
+
+[[bench]]
+name = "fill"
+harness = false
+
+[[bench]]
+name = "clip"
+harness = false
+
+[[bench]]
+name = "isect"
+harness = false
diff --git a/README.md b/README.md
index 20aee6c3..2e264210 100644
--- a/README.md
+++ b/README.md
@@ -54,13 +54,13 @@ for custom allocators is planned in order to make `alloc` optional as well.
 * Type-tagged affine and linear transforms and projections
 * Perspective-correct texture mapping
 * Triangle mesh data structure and a library of shapes
+* Cubic Bézier, Hermite, Catmull–Rom, and B-splines
+* Simple random number generation and distributions
 * Simple text rendering with bitmap fonts
 * Fully customizable rasterization stage
 * Collecting rendering performance data
 * Reading and writing pnm image files
 * Reading and writing Wavefront .obj files
-* Cubic Bezier curves and splines
-* Simple random number generation and distributions
 * Minifb, SDL2, and Wasm frontends
 * Forever emoji-free README and docs
 * Forever LLM-free code
diff --git a/benches/clip.rs b/benches/clip.rs
new file mode 100644
index 00000000..c2450ff3
--- /dev/null
+++ b/benches/clip.rs
@@ -0,0 +1,91 @@
+//! Triangle clipping benchmarks.
+
+use core::{array, iter::repeat_with};
+
+use divan::{Bencher, counter::ItemsCount};
+
+use retrofire_core::{
+    geom::{Tri, vertex},
+    math::rand::{DEFAULT_RNG, DefaultRng, Distrib},
+    math::{orthographic, pt3},
+    render::clip::{ClipVert, view_frustum},
+};
+
+//#[global_allocator]
+//static ALLOC: AllocProfiler = AllocProfiler::system();
+
+#[divan::bench(args = [1, 10, 100, 1000, 10_000])]
+fn clip_mixed(b: Bencher, n: usize) {
+    let rng = &mut DefaultRng::default();
+    let pts = pt3(-10.0, -10.0, -10.0)..pt3(10.0, 10.0, 10.0);
+    let proj = orthographic(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0));
+
+    b.with_inputs(|| {
+        repeat_with(|| {
+            let vs = array::from_fn(|_| {
+                ClipVert::new(vertex(proj.apply(&pts.sample(rng)), ()))
+            });
+            Tri(vs)
+        })
+        .take(n)
+        .collect::<Vec<_>>()
+    })
+    .input_counter(|tris| ItemsCount::of_iter(tris))
+    .bench_local_values(|tris| {
+        let mut out = Vec::new();
+        view_frustum::clip(tris.as_slice(), &mut out);
+        out
+    })
+}
+
+#[divan::bench(args = [1, 10, 100, 1000, 10_000])]
+fn clip_all_inside(b: Bencher, n: usize) {
+    let rng = &mut DefaultRng::default();
+    let pts = pt3(-1.0, -1.0, -1.0)..pt3(1.0, 1.0, 1.0);
+    let proj = orthographic(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0));
+
+    b.with_inputs(|| {
+        repeat_with(|| {
+            let vs = array::from_fn(|_| {
+                ClipVert::new(vertex(proj.apply(&pts.sample(rng)), ()))
+            });
+            Tri(vs)
+        })
+        .take(n)
+        .collect::<Vec<_>>()
+    })
+    .input_counter(|tris| ItemsCount::of_iter(tris))
+    .bench_local_values(|tris| {
+        let mut out = Vec::new();
+        view_frustum::clip(tris.as_slice(), &mut out);
+        out
+    })
+}
+
+#[divan::bench(args = [1, 10, 100, 1000, 10_000])]
+fn clip_all_outside(b: Bencher, n: usize) {
+    let mut rng = DEFAULT_RNG;
+    let pts = pt3(2.0, -10.0, -10.0)..pt3(10.0, 10.0, 10.0);
+    let proj = orthographic(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0));
+
+    b.with_inputs(|| {
+        repeat_with(|| {
+            let vs = ([pts.start; 3]..[pts.end; 3])
+                .sample(&mut rng)
+                .map(|pt| ClipVert::new(vertex(proj.apply(&pt), ())));
+            Tri(vs)
+        })
+        .take(n)
+        .collect::<Vec<_>>()
+    })
+    .input_counter(|tris| ItemsCount::of_iter(tris))
+    .bench_local_values(|tris| {
+        let mut out = Vec::with_capacity(tris.len());
+        view_frustum::clip(tris.as_slice(), &mut out);
+        out
+    })
+}
+
+fn main() {
+    divan::main()
+}
diff --git a/benches/fill.rs b/benches/fill.rs
new file mode 100644
index 00000000..a6f38133
--- /dev/null
+++ b/benches/fill.rs
@@ -0,0 +1,101 @@
+//! Fillrate benchmarks.
+
+use core::iter::zip;
+
+use divan::{Bencher, counter::ItemsCount};
+
+use retrofire_core::{
+    geom::{Tri, vertex},
+    math::{Color3, Color3f, color::gray, pt3, rgb},
+    render::{
+        Texture, raster::ScreenPt, raster::tri_fill, tex::SamplerRepeatPot, uv,
+    },
+    util::{buf::Buf2, pnm::save_ppm},
+};
+
+const SIZES: [f32; 5] = [4.0, 16.0, 64.0, 256.0, 1024.0];
+
+const VERTS: [ScreenPt; 3] =
+    [pt3(0.1, 0.1, 0.0), pt3(0.9, 0.3, 0.5), pt3(0.4, 0.9, 1.0)];
+
+#[divan::bench(args = SIZES)]
+fn flat(b: Bencher, sz: f32) {
+    let mut buf: Buf2<Color3> = Buf2::new((1024, 1024));
+
+    b.with_inputs(|| VERTS.map(|p| vertex(p * sz, ())))
+        .input_counter(move |vs| ItemsCount::new(Tri(*vs).area() as usize))
+        .bench_local_values(|vs| {
+            tri_fill(vs, |sl| {
+                buf[sl.y][sl.xs].fill(gray(0xCC));
+            });
+        });
+
+    save_ppm("benches_fill_flat.ppm", buf).unwrap();
+}
+#[divan::bench(args = SIZES)]
+fn gouraud(b: Bencher, sz: f32) {
+    let mut buf: Buf2<Color3f> = Buf2::new((1024, 1024));
+
+    b.with_inputs(|| {
+        [
+            vertex(VERTS[0] * sz, rgb(0.9, 0.1, 0.0)),
+            vertex(VERTS[1] * sz, rgb(0.1, 0.8, 0.1)),
+            vertex(VERTS[2] * sz, rgb(0.2, 0.3, 1.0)),
+        ]
+    })
+    .input_counter(move |vs| ItemsCount::new(Tri(*vs).area() as usize))
+    .bench_local_values(|vs| {
+        tri_fill(vs, |sl| {
+            let y = sl.y;
+            let xs = sl.xs.clone();
+            let span = &mut buf[y][xs];
+
+            for ((_, col), pix) in zip(sl.vs, span) {
+                *pix = col;
+            }
+        });
+    });
+
+    let buf = Buf2::new_from(
+        (1024, 1024),
+        buf.data().into_iter().map(|c| c.to_color3()),
+    );
+    save_ppm("benches_fill_color.ppm", buf).unwrap();
+}
+
+#[divan::bench(args = SIZES)]
+fn texture(b: Bencher, sz: f32) {
+    let mut buf: Buf2<Color3> = Buf2::new((1024, 1024));
+
+    let tex = Texture::from(Buf2::<Color3>::new_from(
+        (2, 2),
+        [gray(0xFF), gray(0x33), gray(0x33), gray(0xFF)],
+    ));
+    let sampler = SamplerRepeatPot::new(&tex);
+
+    b.with_inputs(|| {
+        [
+            vertex(VERTS[0] * sz, uv(0.0, 0.0)),
+            vertex(VERTS[1] * sz, uv(4.0, 0.0)),
+            vertex(VERTS[2] * sz, uv(0.0, 4.0)),
+        ]
+    })
+    .input_counter(move |vs| ItemsCount::new(Tri(*vs).area() as usize))
+    .bench_local_values(|vs| {
+        tri_fill(vs, |sl| {
+            let y = sl.y;
+            let xs = sl.xs.clone();
+            let span = &mut buf[y][xs];
+
+            for ((_, uv), pix) in zip(sl.vs, span) {
+                *pix = sampler.sample(&tex, uv);
+            }
+        });
+    });
+
+    save_ppm("benches_fill_tex.ppm", buf).unwrap();
+}
+
+fn main() {
+    divan::main()
+}
diff --git a/benches/isect.rs b/benches/isect.rs
new file mode 100644
index 00000000..2775dc42
--- /dev/null
+++ b/benches/isect.rs
@@ -0,0 +1,151 @@
+//! Intersection testing benchmarks.
+
+use core::hint::black_box;
+
+use divan::{Bencher, counter::ItemsCount};
+
+use retrofire::core::{
+    geom::{Ray, Sphere},
+    math::rand::*,
+    math::{Point3, degs, pt3, spherical},
+    render::scene::BBox,
+};
+use retrofire::geom::Intersect;
+
+#[divan::bench]
+fn ray_bbox_hit(b: Bencher) {
+    let mut rng = DefaultRng::default();
+    let bbox = BBox::<()>(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0));
+
+    b.with_inputs(|| {
+        let v = 100.0 * UnitSphere.sample(&mut rng);
+        Ray(v.to_pt(), -v * (0.0..10.0).sample(&mut rng))
+    })
+    .counter(ItemsCount::new(1usize))
+    .bench_local_values(|ray| {
+        assert!(ray.intersect(&black_box(bbox)).is_some())
+    });
+}
+#[divan::bench]
+fn ray_bbox_hit_2(b: Bencher) {
+    let mut rng = DefaultRng::default();
+    let bbox = BBox::<()>(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0));
+
+    let min = spherical(0.0, degs(-180.0), degs(-90.0));
+    let max = spherical(10.0, degs(180.0), degs(-45.0));
+    b.with_inputs(|| {
+        let v = (min..max).sample(&mut rng);
+        Ray(pt3(0.0, 2.0, 0.0), v.to_cart())
+    })
+    .counter(ItemsCount::new(1usize))
+    .bench_local_values(|ray| {
+        assert!(ray.intersect(&black_box(bbox)).is_some())
+    });
+}
+#[divan::bench]
+fn ray_bbox_inside(b: Bencher) {
+    let mut rng = DefaultRng::default();
+    let bbox = BBox::<()>(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0));
+
+    b.with_inputs(|| {
+        let pt = PointsInUnitBall.sample(&mut rng);
+        let dir = VectorsInUnitBall.sample(&mut rng);
+        Ray(pt, dir)
+    })
+    .counter(ItemsCount::new(1usize))
+    .bench_local_values(|ray| {
+        assert!(ray.intersect(&black_box(bbox)).is_some())
+    });
+}
+
+#[divan::bench]
+fn ray_bbox_miss(b: Bencher) {
+    let mut rng = DefaultRng::default();
+    let bbox = BBox::<()>(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0));
+
+    b.with_inputs(|| {
+        let v = (spherical(0.0, degs(-180.0), degs(-45.0))
+            ..spherical(10.0, degs(180.0), degs(90.0)))
+            .sample(&mut rng);
+
+        Ray(pt3(0.0, 3.0, 0.0), v.to_cart())
+    })
+    .counter(ItemsCount::new(1usize))
+    .bench_local_values(|ray| {
+        assert!(ray.intersect(&black_box(bbox)).is_none())
+    });
+}
+
+#[divan::bench]
+fn ray_bbox_mixed(b: Bencher) {
+    let mut rng = DefaultRng::default();
+    let (p, q) = (pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0));
+    let bbox = BBox::<()>(p, q);
+
+    b.with_inputs(|| {
+        // Approximately one third of the rays hits the box
+        let orig = (p..q).sample(&mut rng);
+        let dir = (p..q).sample(&mut rng);
+        Ray(2.0 * orig, 100.0 * dir.to_vec())
+    })
+    .counter(ItemsCount::new(1usize))
+    .bench_local_values(|ray| ray.intersect(&black_box(bbox)));
+}
+
+#[divan::bench]
+fn ray_sphere_miss(b: Bencher) {
+    let mut rng = DefaultRng::default();
+
+    let sphere = Sphere(<Point3>::origin(), 1.0);
+
+    let min = spherical(0.0, degs(-180.0), degs(-45.0));
+    let max = spherical(10.0, degs(180.0), degs(90.0));
+    b.with_inputs(|| {
+        let v = (min..max).sample(&mut rng);
+        Ray(pt3(0.0, 3.0, 0.0), v.to_cart())
+    })
+    .counter(ItemsCount::new(1usize))
+    .bench_local_values(|ray| {
+        let ip = ray.intersect(&black_box(sphere));
+        assert!(ip.is_none());
+        ip
+    });
+}
+
+#[divan::bench]
+fn ray_sphere_hit(b: Bencher) {
+    let mut rng = DefaultRng::default();
+
+    let sphere = Sphere(<Point3>::origin(), 1.0);
+
+    b.with_inputs(|| {
+        let v = (spherical(0.0, degs(-180.0), degs(-90.0))
+            ..spherical(10.0, degs(180.0), degs(-45.0)))
+            .sample(&mut rng);
+        Ray(pt3(0.0, 2.0f32.sqrt(), 0.0), v.to_cart())
+    })
+    .counter(ItemsCount::new(1usize))
+    .bench_local_values(|ray| {
+        let ip = ray.intersect(&black_box(sphere));
+        assert!(ip.is_some());
+        ip
+    });
+}
+
+#[divan::bench]
+fn ray_sphere_mixed(b: Bencher) {
+    let mut rng = DefaultRng::default();
+
+    let sphere = Sphere(<Point3>::origin(), 1.0);
+
+    b.with_inputs(|| {
+        let v = VectorsInUnitBall.sample(&mut rng);
+        Ray(pt3(0.0, 2.0, 0.0), v)
+    })
+    .counter(ItemsCount::new(1usize))
+    .bench_local_values(|ray| ray.intersect(&black_box(sphere)));
+}
+
+fn main() {
+    divan::main();
+}
diff --git a/core/examples/hello_tri.rs b/core/examples/hello_tri.rs
index 20feab14..3e33b97e 100644
--- a/core/examples/hello_tri.rs
+++ b/core/examples/hello_tri.rs
@@ -49,7 +49,7 @@ fn main() {
     if cfg!(feature = "fp") {
         assert_eq!(center_pixel, rgba(151, 128, 187, 255));
     } else {
-        assert_eq!(center_pixel, rgba(114, 102, 127, 255));
+        assert_eq!(center_pixel, rgba(114, 102, 128, 255));
     }
     #[cfg(feature = "std")]
     {
diff --git a/core/src/geom/prim.rs b/core/src/geom/prim.rs
index dfaa0212..5cd48129 100644
--- a/core/src/geom/prim.rs
+++ b/core/src/geom/prim.rs
@@ -297,7 +297,6 @@ impl<A, B> Tri<Vertex3<A, B>> {
     /// );
     /// assert_eq!(tri.area(), 6.0);
     /// ```
-    #[cfg(feature = "fp")]
     pub fn area(&self) -> f32 {
         let [t, u] = self.tangents();
         t.cross(&u).len() / 2.0
@@ -495,11 +494,7 @@ impl<B> Plane3<B> {
     /// ```
     #[inline]
     pub fn signed_dist(&self, pt: Point3<B>) -> f32 {
-        use crate::math::float::*;
-        let len_sqr = self.abc().len_sqr();
-        // TODO use to_homog once committed
-        let pt = [pt.x(), pt.y(), pt.z(), 1.0].into();
-        self.0.dot(&pt) * f32::recip_sqrt(len_sqr)
+        self.dot(pt) / self.abc().len()
     }
 
     /// Returns whether a point is in the half-space that the normal of `self`
@@ -516,10 +511,15 @@ impl<B> Plane3<B> {
     /// assert!(<Plane3>::XY.is_inside(pt));
     /// ```
     // TODO "plane.is_inside(point)" reads wrong
-    #[cfg(feature = "fp")]
     #[inline]
     pub fn is_inside(&self, pt: Point3<B>) -> bool {
-        self.signed_dist(pt) <= 0.0
+        self.dot(pt) <= 0.0
+    }
+
+    fn dot(&self, pt: Point3<B>) -> f32 {
+        // TODO add to_homog method
+        let [x, y, z] = pt.0;
+        self.0.dot(&[x, y, z, 1.0].into())
     }
 
     /// Returns an orthonormal affine basis on `self`.
@@ -543,13 +543,14 @@ impl<B> Plane3<B> {
         let up = self.abc();
 
         let right: Vec3<B> =
-            if up.x().abs() < up.y().abs() && up.x().abs() < up.z().abs() {
+            if up.x().abs() <= up.y().abs() && up.x().abs() <= up.z().abs() {
                 Vec3::X
             } else {
                 Vec3::Z
             };
         let fwd = right.cross(&up).normalize();
-        let right = up.normalize().cross(&fwd);
+        let up = up.normalize();
+        let right = up.cross(&fwd);
 
         let origin = self.offset() * up;
 
@@ -630,7 +631,6 @@ impl<const N: usize, B> Polyline<Point<[f32; N], Real<N, B>>> {
     ///
     /// assert_eq!(pline.len(), 4.0);
     /// ```
-    #[cfg(feature = "fp")]
     pub fn len(&self) -> f32 {
         self.len_by(Point::distance)
     }
@@ -777,9 +777,11 @@ impl<B> Default for Sphere<B> {
 
 #[cfg(test)]
 mod tests {
-    use crate::assert_approx_eq;
-    use crate::math::*;
     use alloc::vec;
+    use core::f32::consts::FRAC_1_SQRT_2;
+
+    use crate::math::*;
+    use crate::{assert_approx_eq, mat};
 
     use super::*;
 
@@ -821,7 +823,6 @@ mod tests {
         let tri = tri(pt2(-1.0, 0.0), pt2(2.0, 0.0), pt2(2.0, 1.0));
         assert_eq!(tri.area(), 1.5);
     }
-    #[cfg(feature = "fp")]
     #[test]
     fn triangle_area_3() {
         // base = 3, height = 2
@@ -880,7 +881,6 @@ mod tests {
         assert_approx_eq!(p.normal(), vec3(0.0, 0.0, 1.0));
         assert_approx_eq!(p.offset(), -3.0);
     }
-    #[cfg(feature = "fp")]
     #[test]
     fn plane_is_point_inside_xz() {
         let p = <Plane3>::from_point_and_normal(pt3(1.0, 2.0, 3.0), Vec3::Y);
@@ -894,7 +894,6 @@ mod tests {
         assert!(!p.is_inside(pt3(0.0, 3.0, 0.0)));
         assert!(!p.is_inside(pt3(1.0, 3.0, 3.0)));
     }
-    #[cfg(feature = "fp")]
     #[test]
     fn plane_is_point_inside_neg_xz() {
         let p = <Plane3>::from_point_and_normal(pt3(1.0, 2.0, 3.0), -Vec3::Y);
@@ -908,7 +907,6 @@ mod tests {
         assert!(p.is_inside(pt3(0.0, 3.0, 0.0)));
         assert!(p.is_inside(pt3(1.0, 3.0, 3.0)));
     }
-    #[cfg(feature = "fp")]
     #[test]
     fn plane_is_point_inside_diagonal() {
         let p = <Plane3>::from_point_and_normal(pt3(0.0, 1.0, 0.0), splat(1.0));
@@ -939,6 +937,23 @@ mod tests {
         );
     }
 
+    #[test]
+    fn plane_basis() {
+        let p = <Plane3>::new(0.0, 2.0, 2.0, 4.0);
+
+        let m = p.basis::<Plane3>();
+
+        assert_approx_eq!(
+            m,
+            mat![
+                1.0, 0.0, 0.0, 0.0;
+                0.0, FRAC_1_SQRT_2, -FRAC_1_SQRT_2, 1.0;
+                0.0, FRAC_1_SQRT_2, FRAC_1_SQRT_2, 1.0;
+                0.0, 0.0, 0.0, 1.0;
+            ]
+        );
+    }
+
     #[test]
     fn polyline_eval_f32() {
         let pl = Polyline(vec![0.0, 1.0, -0.5]);
diff --git a/core/src/lib.rs b/core/src/lib.rs
index f7c6d37e..11d2a072 100644
--- a/core/src/lib.rs
+++ b/core/src/lib.rs
@@ -19,9 +19,9 @@
 //! # Crate features
 //!
 //! * `std`:
-//!   Makes available items requiring I/O, timekeeping, or any floating-point
-//!   functions not included in `core`. In particular this means trigonometric
-//!   and transcendental functions.
+//!   Enabled by default. Makes available items requiring I/O, timekeeping, or
+//!   any floating-point functions not included in `core`. In particular, this
+//!   means trigonometric and transcendental functions.
 //!
 //!   If this feature is disabled, the crate only depends on `alloc`.
 //!
@@ -33,7 +33,10 @@
 //!   Provides fast approximate implementations of floating-point functions
 //!   via the [micromath](https://crates.io/crates/micromath) crate.
 //!
-//! All features are disabled by default.
+//! If none of the above features is enabled, fallback implementations of
+//! a critical subset of floating-point functions are used. However, all APIs
+//! whose implementation relies on trigonometric or transcendental functions
+//! are disabled.
 //!
 //! # Example
 //!
diff --git a/core/src/math/angle.rs b/core/src/math/angle.rs
index 9ccc1a24..83c30016 100644
--- a/core/src/math/angle.rs
+++ b/core/src/math/angle.rs
@@ -2,7 +2,7 @@
 
 use core::{
     f32::consts::{PI, TAU},
-    fmt::{self, Debug, Display},
+    fmt::{self, Debug, Display, Formatter},
     marker::PhantomData,
     ops::{Add, Div, Mul, Neg, Rem, Sub},
     ops::{AddAssign, DivAssign, MulAssign, SubAssign},
@@ -26,11 +26,11 @@ use crate::math::{Vec2, Vec3, float::f32, vec2, vec3};
 pub struct Angle(f32);
 
 /// Tag type for a polar coordinate space
-#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)]
+#[derive(Copy, Clone, Default, Eq, PartialEq)]
 pub struct Polar<B>(PhantomData<B>);
 
 /// Tag type for a spherical coordinate space.
-#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)]
+#[derive(Copy, Clone, Default, Eq, PartialEq)]
 pub struct Spherical<B>(PhantomData<B>);
 
 /// A polar coordinate vector, with radius and azimuth components.
@@ -199,6 +199,22 @@ impl Angle {
     pub fn clamp(self, min: Self, max: Self) -> Self {
         Self(self.0.clamp(min.0, max.0))
     }
+
+    /// Returns `self` "wrapped around" to the range `min..max`.
+    ///
+    /// # Examples
+    /// ```
+    /// use retrofire_core::assert_approx_eq;
+    /// use retrofire_core::math::{degs, turns};
+    ///
+    /// // 400 (mod 360) = 40
+    /// assert_approx_eq!(degs(400.0).wrap(turns(0.0), turns(1.0)), degs(40.0))
+    /// ```
+    #[must_use]
+    pub fn wrap(self, min: Self, max: Self) -> Self {
+        use super::float::f32;
+        Self(min.0 + f32::rem_euclid(self.0 - min.0, max.0 - min.0))
+    }
 }
 
 #[cfg(feature = "fp")]
@@ -247,21 +263,6 @@ impl Angle {
     pub fn tan(self) -> f32 {
         f32::tan(self.0)
     }
-
-    /// Returns `self` "wrapped around" to the range `min..max`.
-    ///
-    /// # Examples
-    /// ```
-    /// use retrofire_core::assert_approx_eq;
-    /// use retrofire_core::math::{degs, turns};
-    ///
-    /// // 400 (mod 360) = 40
-    /// assert_approx_eq!(degs(400.0).wrap(turns(0.0), turns(1.0)), degs(40.0))
-    /// ```
-    #[must_use]
-    pub fn wrap(self, min: Self, max: Self) -> Self {
-        Self(min.0 + f32::rem_euclid(self.0 - min.0, max.0 - min.0))
-    }
 }
 
 impl<B> PolarVec<B> {
@@ -330,17 +331,34 @@ impl<B> SphericalVec<B> {
     /// Returns `self` converted to the equivalent Cartesian 3-vector.
     ///
     /// # Examples
-    /// TODO examples
+    /// ```
+    /// use retrofire_core::assert_approx_eq;
+    /// use retrofire_core::math::{degs, spherical, vec3, SphericalVec};
+    ///
+    /// let mut v = spherical::<()>(1.0, degs(0.0), degs(0.0));
+    /// assert_approx_eq!(v.to_cart(), vec3(1.0, 0.0, 0.0));
+    ///
+    /// v = spherical(2.0, degs(90.0), degs(0.0));
+    /// assert_approx_eq!(v.to_cart(), vec3(0.0, 0.0, -2.0));
+    ///
+    /// v = spherical(3.0, degs(0.0), degs(90.0));
+    /// assert_approx_eq!(v.to_cart(), vec3(0.0, 3.0, 0.0));
+    /// ```
     #[cfg(feature = "fp")]
     pub fn to_cart(&self) -> Vec3<B> {
+        // First about z by alt, then about y by az:
+        //
+        // ( caz  0 saz )   ( calt -salt  0 )   ( r )
+        // (   0  1   0 ) · ( salt  calt  0 ) · ( 0 )
+        // (-saz  0 caz )   (   0      0  1 )   ( 0 )
+        //
+        //   ( caz  0 saz )   ( calt·r )   ( caz·calt·r )
+        // = (   0  1   0 ) · ( salt·r ) = (     salt·r )
+        //   (-saz  0 caz )   (      0 )   (-saz·calt·r )
+
         let (sin_alt, cos_alt) = self.alt().sin_cos();
         let (sin_az, cos_az) = self.az().sin_cos();
-
-        let x = cos_az * cos_alt;
-        let z = sin_az * cos_alt;
-        let y = sin_alt;
-
-        self.r() * vec3(x, y, z)
+        self.r() * vec3(cos_az * cos_alt, sin_alt, -sin_az * cos_alt)
     }
 }
 
@@ -392,12 +410,11 @@ impl<B> Vec2<B> {
 impl<B> Vec3<B> {
     /// Converts `self` into the equivalent spherical coordinate vector.
     ///
-    /// The `r` component of the result equals `self.len()`.
-    ///
-    /// The `az` component is the angle between `self` and the xy-plane in the
-    /// range (-180°, 180°] such that positive `z` maps to positive `az`.
-    ///
-    /// The `alt` component is the angle between `self` and the xz-plane in the
+    /// Returns a vector (r, az, alt) such that:
+    /// * `r` equals `self.len()`
+    /// * `az`is the angle between `self` and the xy-plane in the range
+    ///   (-180°, 180°] such that positive `z` maps to *negative* `az`, and
+    /// * `alt` is the angle between `self` and the xz-plane in the
     /// range [-90°, 90°] such that positive `y` maps to positive `alt`.
     ///
     /// # Examples
@@ -406,23 +423,23 @@ impl<B> Vec3<B> {
     ///
     /// // The positive x-axis lies at zero azimuth and altitude
     /// assert_eq!(
-    ///     vec3(2.0, 0.0, 0.0).to_spherical(),
-    ///     spherical::<()>(2.0, degs(0.0), degs(0.0))
+    ///     vec3(1.0, 0.0, 0.0).to_spherical(),
+    ///     spherical::<()>(1.0, degs(0.0), degs(0.0))
     /// );
     /// // The positive y-axis lies at 90° altitude
     /// assert_eq!(
     ///     vec3(0.0, 2.0, 0.0).to_spherical(),
     ///     spherical::<()>(2.0, degs(0.0), degs(90.0))
     /// );
-    /// // The positive z axis lies at 90° azimuth
+    /// // The positive z-axis lies at *-90°* azimuth
     /// assert_eq!(
-    ///     vec3(0.0, 0.0, 2.0).to_spherical(),
-    ///     spherical::<()>(2.0, degs(90.0), degs(0.0))
+    ///     vec3(0.0, 0.0, 3.0).to_spherical(),
+    ///     spherical::<()>(3.0, degs(-90.0), degs(0.0))
     /// );
     /// ```
     pub fn to_spherical(&self) -> SphericalVec<B> {
         let [x, y, z] = self.0;
-        let az = atan2(z, x);
+        let az = atan2(-z, x);
         let alt = atan2(y, f32::sqrt(x * x + z * z));
         let r = self.len();
         spherical(r, az, alt)
@@ -477,6 +494,17 @@ impl ZDiv for Angle {}
 // Foreign trait impls
 //
 
+impl<B> Debug for Polar<B> {
+    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
+        f.write_str("Pol")
+    }
+}
+impl<B> Debug for Spherical<B> {
+    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
+        f.write_str("Sph")
+    }
+}
+
 impl<B> Default for SphericalVec<B> {
     fn default() -> Self {
         Self::new([1.0, 0.0, 0.0])
@@ -618,6 +646,7 @@ impl<B> From<Vec3<B>> for SphericalVec<B> {
 #[allow(unused, nonstandard_style)]
 mod tests {
     use core::f32::consts::{PI, TAU};
+    use std::eprintln;
 
     use crate::{
         assert_approx_eq,
@@ -626,8 +655,12 @@ mod tests {
 
     use super::*;
 
-    const vec2: fn(f32, f32) -> Vec2 = math::vec2;
-    const vec3: fn(f32, f32, f32) -> Vec3 = math::vec3;
+    const fn vec2(x: f32, y: f32) -> Vec2 {
+        math::vec2(x, y)
+    }
+    const fn vec3(x: f32, y: f32, z: f32) -> Vec3 {
+        math::vec3(x, y, z)
+    }
 
     #[test]
     fn rads_to_degs() {
@@ -710,7 +743,6 @@ mod tests {
         assert_approx_eq!(atan2(-1.0, 1.0), degs(-45.0));
     }
 
-    #[cfg(feature = "fp")]
     #[test]
     fn wrapping() {
         use crate::assert_approx_eq;
@@ -788,68 +820,46 @@ mod tests {
         assert_eq!(vec2(0.0, -4.0).to_polar(), polar(4.0, degs(-90.0)));
     }
 
-    #[cfg(feature = "fp")]
-    #[test]
-    fn spherical_to_cartesian() {
-        let spherical = spherical::<()>;
-        assert_eq!(
-            spherical(0.0, degs(0.0), degs(0.0)).to_cart(),
-            vec3(0.0, 0.0, 0.0)
-        );
-        assert_eq!(
-            spherical(1.0, degs(0.0), degs(0.0)).to_cart(),
-            vec3(1.0, 0.0, 0.0)
-        );
-        assert_approx_eq!(
-            spherical(2.0, degs(60.0), degs(0.0)).to_cart(),
-            vec3(1.0, 0.0, SQRT_3)
-        );
-        assert_approx_eq!(
-            spherical(2.0, degs(90.0), degs(0.0)).to_cart(),
-            vec3(0.0, 0.0, 2.0)
-        );
-        assert_approx_eq!(
-            spherical(3.0, degs(123.0), degs(90.0)).to_cart(),
-            vec3(0.0, 3.0, 0.0)
-        );
+    const fn sph(r: f32, az: f32, alt: f32) -> SphericalVec {
+        spherical(r, degs(az), degs(alt))
     }
+    #[rustfmt::skip]
+    const CART_SPH: [(Vec3, SphericalVec); 10] = [
+        (vec3(   0.0,    0.0,     0.0), sph(0.0,   0.0,  0.0)),
+
+        (vec3(    1.0,   0.0,     0.0), sph(1.0,   0.0,  0.0)),
+        (vec3( SQRT_3,   0.0,    -1.0), sph(2.0,  30.0,  0.0)),
+        (vec3(    1.0,   0.0, -SQRT_3), sph(2.0,  60.0,  0.0)),
+        (vec3(    0.0,   0.0,    -2.0), sph(2.0,  90.0,  0.0)),
+        (vec3(-SQRT_3,   0.0,    -1.0), sph(2.0, 150.0,  0.0)),
+
+        // Doesn't roundtrip due to imprecision and
+        // the discontinuity from 180° to -180° :(
+        (vec3(  -3.0,    0.0,     0.0), sph(3.0, 180.0,  0.0)),
+
+        (vec3(SQRT_3,    1.0,     0.0), sph(2.0,  0.0,  30.0)),
+        (vec3(   1.0, SQRT_3,     0.0), sph(2.0,  0.0,  60.0)),
+        (vec3(   0.0,    2.0,     0.0), sph(2.0,  0.0,  90.0)),
+        (vec3(   0.0,   -3.0,     0.0), sph(3.0,  0.0, -90.0)),
+    ];
 
     #[cfg(feature = "fp")]
     #[test]
-    fn cartesian_to_spherical_zero_alt() {
-        assert_approx_eq!(
-            vec3(0.0, 0.0, 0.0).to_spherical(),
-            spherical(0.0, degs(0.0), degs(0.0))
-        );
-        assert_eq!(
-            vec3(1.0, 0.0, 0.0).to_spherical(),
-            spherical(1.0, degs(0.0), degs(0.0))
-        );
-        assert_approx_eq!(
-            vec3(1.0, SQRT_3, 0.0).to_spherical(),
-            spherical(2.0, degs(0.0), degs(60.0))
-        );
-        assert_eq!(
-            vec3(0.0, 2.0, 0.0).to_spherical(),
-            spherical(2.0, degs(0.0), degs(90.0))
-        );
+    fn spherical_to_cartesian() {
+        for (cart, sp) in CART_SPH {
+            let actual = sp.to_cart();
+            eprintln!("Testing {sp:?} -> {cart:?}");
+            assert_approx_eq!(actual, cart);
+        }
     }
 
     #[cfg(feature = "fp")]
     #[test]
     fn cartesian_to_spherical() {
-        use core::f32::consts::SQRT_2;
-        assert_approx_eq!(
-            vec3(SQRT_3, 0.0, 1.0).to_spherical(),
-            spherical(2.0, degs(30.0), degs(0.0))
-        );
-        assert_approx_eq!(
-            vec3(1.0, SQRT_2, 1.0).to_spherical(),
-            spherical(2.0, degs(45.0), degs(45.0))
-        );
-        assert_approx_eq!(
-            vec3(0.0, 0.0, 3.0).to_spherical(),
-            spherical(3.0, degs(90.0), degs(0.0))
-        );
+        for (cart, sp) in CART_SPH {
+            let actual = cart.to_spherical();
+            eprintln!("Testing {cart:?} -> {sp:?}");
+            assert_approx_eq!(actual, sp);
+        }
     }
 }
diff --git a/core/src/math/float.rs b/core/src/math/float.rs
index 73a56ed7..956dd65b 100644
--- a/core/src/math/float.rs
+++ b/core/src/math/float.rs
@@ -26,8 +26,9 @@ pub mod libm {
 
     pub use super::fallback::rem_euclid;
 
+    #[inline]
     pub fn recip_sqrt(x: f32) -> f32 {
-        powf(x, -0.5)
+        1.0 / sqrt(x)
     }
 }
 
@@ -47,7 +48,8 @@ pub mod mm {
     #[inline]
     pub fn sqrt(x: f32) -> f32 {
         let y = mm::sqrt(x);
-        // One round of Newton's method
+        // Two rounds of Newton's method
+        let y = 0.5 * (y + (x / y));
         0.5 * (y + (x / y))
     }
     /// Returns the approximate reciprocal of the square root of `x`.
@@ -61,6 +63,7 @@ pub mod mm {
     pub fn powf(x: f32, y: f32) -> f32 {
         mm::powf(x, y)
     }
+
     #[inline]
     pub fn sin(x: f32) -> f32 {
         mm::sin(x)
@@ -90,27 +93,45 @@ pub mod mm {
         }
         mm::atan2(y, x)
     }
+
+    #[inline]
+    pub fn exp(x: f32) -> f32 {
+        mm::exp(x)
+    }
+    #[inline]
+    pub fn log2(x: f32) -> f32 {
+        mm::log2(x)
+    }
 }
 
 pub mod fallback {
     /// Returns the largest integer less than or equal to `x`.
     #[inline]
     pub fn floor(x: f32) -> f32 {
-        (x as i64 - x.is_sign_negative() as i64) as f32
+        (x as i64 - (x < 0.0) as i64) as f32
     }
-    // Returns the least non-negative remainder of `x` (mod `m`).
+    /// Returns the least non-negative remainder of `x` (mod `m`).
     #[inline]
     pub fn rem_euclid(x: f32, m: f32) -> f32 {
-        x % m + (x.is_sign_negative() as u32 as f32) * m
+        let r = x % m;
+        r + if r < 0.0 { m.abs() } else { 0.0 }
     }
     /// Returns the approximate reciprocal of the square root of `x`.
     #[inline]
     pub fn recip_sqrt(x: f32) -> f32 {
+        if x < 0.0 {
+            return f32::NAN;
+        }
         // https://en.wikipedia.org/wiki/Fast_inverse_square_root
         let y = f32::from_bits(0x5f37_5a86 - (x.to_bits() >> 1));
-        // A round of Newton's method
+        // Two rounds of Newton's method
+        let y = y * (1.5 - 0.5 * x * y * y);
         y * (1.5 - 0.5 * x * y * y)
     }
+    #[inline]
+    pub fn sqrt(x: f32) -> f32 {
+        1.0 / recip_sqrt(x)
+    }
 }
 
 #[cfg(feature = "std")]
@@ -141,59 +162,134 @@ pub use fallback as f32;
 #[cfg(test)]
 #[allow(unused_imports)]
 mod tests {
+    use core::f32::consts::*;
 
+    use super::{RecipSqrt, f32, *};
     use crate::assert_approx_eq;
 
     #[cfg(feature = "libm")]
     #[test]
     fn libm_functions() {
-        use super::libm;
-        use core::f32::consts::PI;
-        assert_eq!(libm::cos(PI), -1.0);
+        assert_eq!(libm::floor(1.5), 1.0);
+        assert_eq!(libm::floor(0.99), 0.0);
+        assert_eq!(libm::floor(-0.0), 0.0);
+        assert_eq!(libm::floor(-1.1), -2.0);
+
+        assert_approx_eq!(libm::rem_euclid(1.6, 0.5), 0.1);
+        assert_approx_eq!(libm::rem_euclid(-1.6, 0.5), 0.4);
+        assert_approx_eq!(libm::rem_euclid(1.6, -0.5), 0.1);
+        assert_approx_eq!(libm::rem_euclid(-1.6, -0.5), 0.4);
+
         assert_eq!(libm::sqrt(9.0), 3.0);
+        assert_eq!(libm::sqrt(16.0), 4.0);
+        assert!(libm::sqrt(-1.0).is_nan());
+        assert_eq!(libm::recip_sqrt(9.0), 1.0 / 3.0);
+        assert_eq!(libm::recip_sqrt(0.0), f32::INFINITY);
+        assert!(libm::recip_sqrt(-1.0).is_nan());
+
+        assert_eq!(libm::powf(3.0, 2.0), 9.0);
+        assert_eq!(libm::powf(-3.0, 2.0), 9.0);
+        assert_eq!(libm::powf(3.0, -2.0), 1.0 / 9.0);
+        assert_eq!(libm::powf(-3.0, 3.0), -27.0);
+
+        assert_approx_eq!(libm::sin(FRAC_PI_6), 0.5);
+        assert_eq!(libm::cos(PI), -1.0);
+
+        assert_eq!(libm::exp(1.0), E);
+        assert_approx_eq!(libm::exp(2.0), E * E);
+        assert_eq!(libm::log2(8.0), 3.0);
+        assert!(libm::log2(-1.0).is_nan());
     }
 
     #[cfg(feature = "mm")]
     #[test]
     fn mm_functions() {
-        use core::f32::consts::*;
+        assert_eq!(mm::floor(1.5), 1.0);
+        assert_eq!(mm::floor(0.99), 0.0);
+        assert_eq!(mm::floor(-0.0), 0.0);
+        assert_eq!(mm::floor(-1.1), -2.0);
 
-        use super::f32;
+        assert_approx_eq!(mm::rem_euclid(1.6, 0.5), 0.1);
+        assert_approx_eq!(mm::rem_euclid(-1.6, 0.5), 0.4);
+        assert_approx_eq!(mm::rem_euclid(1.6, -0.5), 0.1);
+        assert_approx_eq!(mm::rem_euclid(-1.6, -0.5), 0.4);
 
-        assert_approx_eq!(f32::sin(FRAC_PI_6), 0.5);
-        assert_eq!(f32::cos(PI), -1.0);
-        assert_eq!(f32::sqrt(16.0), 4.0);
-        assert_approx_eq!(f32::sqrt(9.0), 3.0, eps = 1e-3);
+        assert_approx_eq!(mm::sqrt(9.0), 3.0);
+        assert_eq!(mm::sqrt(16.0), 4.0);
+        assert!(mm::sqrt(-1.0).is_nan());
+        assert_approx_eq!(mm::recip_sqrt(9.0), 1.0 / 3.0);
+        // mm doesn't check for zero, just gives a big number
+        assert_approx_eq!(mm::recip_sqrt(0.0), 1.9818e19);
+        // mm doesn't check for negative, panics due to sub overflow
+        //assert!(mm::recip_sqrt(-1.0).is_nan());
+
+        assert_approx_eq!(mm::powf(3.0, 2.0), 9.0);
+        assert_approx_eq!(mm::powf(-3.0, 2.0), 9.0);
+        assert_approx_eq!(mm::powf(3.0, -2.0), 1.0 / 9.0);
+        assert_approx_eq!(mm::powf(-3.0, 3.0), -27.0);
+
+        assert_approx_eq!(mm::sin(FRAC_PI_6), 0.5);
+        assert_eq!(mm::cos(PI), -1.0);
+
+        assert_eq!(mm::exp(1.0), E);
+        assert_approx_eq!(mm::exp(2.0), E * E);
+        assert_approx_eq!(mm::log2(8.0), 3.0);
+        // mm doesn't check for negative, panics due to sub overflow
+        //assert!(mm::log2(-1.0).is_nan());
     }
 
     #[cfg(feature = "std")]
     #[test]
     fn std_functions() {
-        use super::f32;
-        use core::f32::consts::PI;
-        assert_eq!(f32::cos(PI), -1.0);
+        assert_eq!(f32::floor(-0.0), 0.0);
+
+        assert_approx_eq!(f32::rem_euclid(1.6, 0.5), 0.1);
+        assert_approx_eq!(f32::rem_euclid(-1.6, 0.5), 0.4);
+        assert_approx_eq!(f32::rem_euclid(1.6, -0.5), 0.1);
+        assert_approx_eq!(f32::rem_euclid(-1.6, -0.5), 0.4);
+
         assert_eq!(f32::sqrt(9.0), 3.0);
+        assert!(f32::sqrt(-1.0).is_nan());
+        assert_eq!(f32::recip_sqrt(9.0), 1.0 / 3.0);
+        assert_eq!(f32::recip_sqrt(0.0), f32::INFINITY);
+        assert!(f32::recip_sqrt(-1.0).is_nan());
+
+        assert_eq!(f32::cos(PI), -1.0);
     }
 
     #[cfg(not(feature = "fp"))]
     #[test]
     fn fallback_functions() {
-        use super::{RecipSqrt, f32};
+        use fallback as fb;
+        assert_eq!(fb::floor(1.5), 1.0);
+        assert_eq!(fb::floor(0.99), 0.0);
+        assert_eq!(fb::floor(-0.0), 0.0);
+        assert_eq!(fb::floor(-1.1), -2.0);
 
-        assert_eq!(f32::floor(1.23), 1.0);
-        assert_eq!(f32::floor(0.0), 0.0);
-        assert_eq!(f32::floor(-1.23), -2.0);
+        assert_approx_eq!(fb::rem_euclid(1.6, 0.5), 0.1);
+        assert_approx_eq!(fb::rem_euclid(-1.6, 0.5), 0.4);
+        assert_approx_eq!(fb::rem_euclid(1.6, -0.5), 0.1);
+        assert_approx_eq!(fb::rem_euclid(-1.6, -0.5), 0.4);
 
-        assert_approx_eq!(f32::rem_euclid(1.23, 4.0), 1.23);
-        assert_approx_eq!(f32::rem_euclid(4.0, 4.0), 0.0);
-        assert_approx_eq!(f32::rem_euclid(5.67, 4.0), 1.67);
-        assert_approx_eq!(f32::rem_euclid(-1.23, 4.0), 2.77);
-    }
+        assert_approx_eq!(fb::sqrt(9.0), 3.0);
+        assert_approx_eq!(fb::sqrt(16.0), 4.0);
+        assert!(fb::sqrt(-1.0).is_nan());
+        assert_approx_eq!(fb::recip_sqrt(9.0), 1.0 / 3.0);
+        // doesn't check for infinity, just returns a big number
+        assert_approx_eq!(fb::recip_sqrt(0.0), 2.9727e19);
+        assert!(fb::recip_sqrt(-1.0).is_nan());
 
-    #[test]
-    fn recip_sqrt() {
-        use super::{RecipSqrt, f32};
-        assert_approx_eq!(f32::recip_sqrt(2.0), f32::sqrt(0.5), eps = 1e-2);
-        assert_approx_eq!(f32::recip_sqrt(9.0), 1.0 / 3.0, eps = 1e-3);
+        // assert_eq!(fb::powf(3.0, 2.0), 9.0);
+        // assert_eq!(fb::powf(-3.0, 2.0), 9.0);
+        // assert_eq!(fb::powf(3.0, -2.0), 1.0 / 9.0);
+        // assert_eq!(fb::powf(-3.0, 3.0), -27.0);
+        //
+        // assert_approx_eq!(libm::sin(FRAC_PI_6), 0.5);
+        // assert_eq!(libm::cos(PI), -1.0);
+        //
+        // assert_eq!(libm::exp(1.0), E);
+        // assert_approx_eq!(libm::exp(2.0), E * E);
+        // assert_eq!(libm::log2(8.0), 3.0);
+        // assert!(libm::log2(-1.0).is_nan());
     }
 }
diff --git a/core/src/math/grad.rs b/core/src/math/grad.rs
index 0fab3c6c..2df68f17 100644
--- a/core/src/math/grad.rs
+++ b/core/src/math/grad.rs
@@ -30,7 +30,6 @@ pub enum Kind<Pt> {
     /// |P - Q| to *t* values such that:
     /// * *t* = 0 when P = Q, and
     /// * *t* = 1 when |P - Q| >= *r*.
-    #[cfg(feature = "fp")]
     Radial(Pt, f32),
     /// TODO
     #[cfg(feature = "fp")]
@@ -61,8 +60,7 @@ impl<T: Lerp> Gradient2<T> {
     pub fn eval(&self, p: Point2) -> T {
         let t = match self.kind {
             Kind::Linear(p0, p1) => (p - p0).scalar_project(&(p1 - p0)),
-            #[cfg(feature = "fp")]
-            Kind::Radial(p0, r) => (p - p0).len() / r,
+            Kind::Radial(p0, r) => p.distance(&p0) / r,
             #[cfg(feature = "fp")]
             Kind::Conical(p0) => {
                 let angle = (p - p0).atan();
@@ -186,7 +184,6 @@ mod tests {
         assert_eq!(g.eval(pt2(2.0, -100.0)), 0.1);
     }
 
-    #[cfg(feature = "fp")]
     #[test]
     fn radial_gradient() {
         let g = Gradient2::new(
@@ -196,11 +193,11 @@ mod tests {
 
         assert_eq!(g.eval(pt2(2.0, -1.0)), 0.9); // t=0.0
 
-        assert_eq!(g.eval(pt2(2.0, -1.5)), 0.9); // t=0.25
-        assert_eq!(g.eval(pt2(1.5, -1.0)), 0.9); // t=0.25
+        assert_approx_eq!(g.eval(pt2(2.0, -1.5)), 0.9); // t=0.25
+        assert_approx_eq!(g.eval(pt2(1.5, -1.0)), 0.9); // t=0.25
 
-        assert_eq!(g.eval(pt2(2.0, -1.75)), 0.5); // t=0.375
-        assert_eq!(g.eval(pt2(1.25, -1.0)), 0.5); // t=0.375
+        assert_approx_eq!(g.eval(pt2(2.0, -1.75)), 0.5); // t=0.375
+        assert_approx_eq!(g.eval(pt2(1.25, -1.0)), 0.5); // t=0.375
 
         assert_eq!(g.eval(pt2(2.0, -3.0)), 0.1); // t=0.5
         assert_eq!(g.eval(pt2(0.0, -1.0)), 0.1); // t=0.5
diff --git a/core/src/math/mat.rs b/core/src/math/mat.rs
index 9a166604..a0d472c0 100644
--- a/core/src/math/mat.rs
+++ b/core/src/math/mat.rs
@@ -1691,7 +1691,7 @@ mod tests {
         fn orientation_no_op() {
             let m = orient_y(Y, X);
 
-            assert_eq!(m.apply(&X), X);
+            assert_approx_eq!(m.apply(&X), X);
             assert_eq!(m.apply(&X.to_pt()), X.to_pt());
 
             assert_eq!(m.apply(&Y), Y);
@@ -1705,7 +1705,7 @@ mod tests {
         fn orientation_y_to_z() {
             let m = orient_y(Z, X);
 
-            assert_eq!(m.apply(&X), X);
+            assert_approx_eq!(m.apply(&X), X);
             assert_eq!(m.apply(&X.to_pt()), X.to_pt());
 
             assert_eq!(m.apply(&Y), Z);
@@ -1719,7 +1719,7 @@ mod tests {
         fn orientation_z_to_y() {
             let m = orient_z(Y, X);
 
-            assert_eq!(m.apply(&X), X);
+            assert_approx_eq!(m.apply(&X), X);
             assert_eq!(m.apply(&X.to_pt()), X.to_pt());
 
             assert_eq!(m.apply(&Y), -Z);
diff --git a/core/src/math/param.rs b/core/src/math/param.rs
index 8941df56..4fd6d4e8 100644
--- a/core/src/math/param.rs
+++ b/core/src/math/param.rs
@@ -4,6 +4,7 @@ use super::Lerp;
 
 /// Represents a single-variable parametric curve.
 // TODO More documentation
+// TODO Associated type instead of parameter?
 pub trait Parametric<T> {
     /// Returns the value of `self` at `t`.
     ///
diff --git a/core/src/math/point.rs b/core/src/math/point.rs
index b3da545f..3dc9a513 100644
--- a/core/src/math/point.rs
+++ b/core/src/math/point.rs
@@ -109,7 +109,6 @@ impl<const N: usize, B> Point<[f32; N], Real<N, B>> {
     /// let y4 = pt2(0.0, 4.0);
     /// assert_eq!(x3.distance(&y4), 5.0);
     /// ```
-    #[cfg(feature = "fp")]
     #[inline]
     pub fn distance(&self, other: &Self) -> f32 {
         self.sub(other).len()
@@ -183,7 +182,6 @@ impl<const N: usize, B> Point<[f32; N], Real<N, B>> {
     /// assert_eq!(a.approach(&a, 1.0), a);
     /// assert_eq!(a.approach(&a, -1.0), a);
     /// ```
-    #[cfg(feature = "fp")]
     #[must_use]
     pub fn approach(&self, other: &Self, d: f32) -> Self {
         let v = *other - *self;
@@ -473,6 +471,7 @@ mod tests {
 
     mod f32 {
         use super::*;
+        use crate::assert_approx_eq;
 
         const pt2: fn(f32, f32) -> Point2 = super::pt2;
         const pt3: fn(f32, f32, f32) -> Point3 = super::pt3;
@@ -522,10 +521,12 @@ mod tests {
             );
         }
         #[test]
-        #[cfg(feature = "fp")]
         fn point_point_distance() {
-            assert_eq!(pt2(1.0, -1.0).distance(&pt2(-2.0, 3.0)), 5.0);
-            assert_eq!(pt3(1.0, -3.0, 2.0).distance(&pt3(-2.0, 3.0, 4.0)), 7.0);
+            assert_approx_eq!(pt2(1.0, -1.0).distance(&pt2(-2.0, 3.0)), 5.0);
+            assert_approx_eq!(
+                pt3(1.0, -3.0, 2.0).distance(&pt3(-2.0, 3.0, 4.0)),
+                7.0
+            );
         }
         #[test]
         fn point2_clamp() {
diff --git a/core/src/math/rand.rs b/core/src/math/rand.rs
index 3bb6a573..cb37df7c 100644
--- a/core/src/math/rand.rs
+++ b/core/src/math/rand.rs
@@ -462,7 +462,6 @@ where
     }
 }
 
-#[cfg(feature = "fp")]
 impl Distrib for UnitCircle {
     type Sample = Vec2;
 
@@ -494,7 +493,7 @@ impl Distrib for VectorsOnUnitDisk {
     /// let rng = &mut DefaultRng::default();
     ///
     /// let vec = VectorsOnUnitDisk.sample(rng);
-    /// assert!(vec.len_sqr() <= 1.0);
+    /// assert!(vec.len() <= 1.0);
     /// ```
     fn sample(&self, rng: &mut DefaultRng) -> Vec2 {
         let d = Uniform([-1.0f32; 2]..[1.0; 2]);
@@ -508,7 +507,6 @@ impl Distrib for VectorsOnUnitDisk {
     }
 }
 
-#[cfg(feature = "fp")]
 impl Distrib for UnitSphere {
     type Sample = Vec3;
 
@@ -521,7 +519,7 @@ impl Distrib for UnitSphere {
     /// let rng = &mut DefaultRng::default();
     ///
     /// let vec = UnitSphere.sample(rng);
-    /// assert_approx_eq!(vec.len_sqr(), 1.0);
+    /// assert_approx_eq!(vec.len(), 1.0);
     /// ```
     fn sample(&self, rng: &mut DefaultRng) -> Vec3 {
         let d = Uniform([-1.0; 3]..[1.0; 3]);
@@ -680,7 +678,6 @@ mod tests {
         assert_eq!(approx_100, 82);
     }
 
-    #[cfg(feature = "fp")]
     #[test]
     fn unit_circle() {
         use crate::assert_approx_eq;
@@ -696,7 +693,6 @@ mod tests {
         }
     }
 
-    #[cfg(feature = "fp")]
     #[test]
     fn unit_sphere() {
         use crate::assert_approx_eq;
diff --git a/core/src/math/spline.rs b/core/src/math/spline.rs
index ba3a07aa..7dfb5e80 100644
--- a/core/src/math/spline.rs
+++ b/core/src/math/spline.rs
@@ -3,8 +3,12 @@ use alloc::vec::Vec;
 use core::{array::from_fn, fmt::Debug, marker::PhantomData};
 
 use crate::geom::{Polyline, Ray};
+use crate::mat;
 
-use super::{Affine, Lerp, Linear, Parametric, Vector};
+use super::{
+    Affine, Lerp, Linear, Mat4, Parametric, Point, Vary, Vector, inv_lerp,
+    space::Real,
+};
 
 /// A cubic Bézier curve, defined by four control points.
 ///
@@ -63,6 +67,26 @@ pub struct BezierSpline<T>(Vec<T>);
 //       to include the correct `T::Diff: Trait` bounds
 pub struct HermiteSpline<T: Affine>(Vec<Ray<T>>, PhantomData<T::Diff>);
 
+/// A piecewise curve composed of concatenated cubic curves.
+#[derive(Debug, Clone, Eq, PartialEq)]
+pub struct CatmullRomSpline<T>(Vec<T>);
+
+/// A piecewise curve composed of concatenated cubic curves.
+#[derive(Debug, Clone, Eq, PartialEq)]
+pub struct BSpline<T>(Vec<T>);
+
+/// Euclidean (arc-length) parameterization for splines.
+///
+/// Instead of *t* ∈ [0, 1], a `Euclidean` spline is parameterized by *s*,
+/// the actual Euclidean distance travelled along the curve. This makes it easy,
+/// for example, to position objects along the spline at regular intervals, or
+/// to travel along the spline at a desired speed.
+///
+/// Because arc-length parameterization is not generally possible in closed form,
+/// this implementation uses a look-up table to map *s* values to *t* values,
+/// interpolating linearly between entries.
+pub struct Euclidean<Spl>(Spl, Vec<(f32, f32)>);
+
 /// Interpolates smoothly from 0.0 to 1.0 as `t` goes from 0.0 to 1.0.
 ///
 /// Returns 0 for all `t` <= 0 and 1 for all `t` >= 1. Has a continuous
@@ -534,6 +558,220 @@ where
     }
 }
 
+impl<T> CatmullRomSpline<T>
+where
+    T: Affine<Diff: Linear<Scalar = f32>> + Clone,
+{
+    const _CHAR_MAT: Mat4 = mat![
+         0.0,  1.0,  0.0,  0.0;
+        -0.5,  0.0,  0.5,  0.0;
+         1.0, -2.5,  2.0, -0.5;
+        -0.5,  1.5, -1.5,  0.5;
+    ];
+
+    pub fn new(pts: impl IntoIterator<Item = T>) -> Self {
+        let pts: Vec<_> = pts.into_iter().collect();
+        assert!(
+            pts.len() >= 4,
+            "a Catmull–Rom spline requires at least four points"
+        );
+        Self(pts)
+    }
+
+    /// Returns the point of `self` at the given *t* value.
+    ///
+    /// Values of *t* outside the interval [0, 1] are accepted and extrapolate
+    /// the curve beyond the control points.
+    pub fn eval(&self, t: f32) -> T {
+        let (t, [p0, p1, p2, p3]) = crb_segment(&self.0, t);
+        let [_0, t1, t2, t3] = [1.0, t, t * t, t * t * t];
+
+        let _0 = (-t1 + 2.0 * t2 - t3) / 2.0;
+        let b1 = (2.0 - 5.0 * t2 + 3.0 * t3) / 2.0;
+        let b2 = (t1 + 4.0 * t2 - 3.0 * t3) / 2.0;
+        let b3 = (-t2 + t3) / 2.0;
+
+        // b0 + b1 + b2 + b3 = 1
+        // b0 = 1 - b1 - b2 - b3
+
+        //    b0·P0 + b1·P1 + b2·P2 + b3·P3
+        // = (1 - b1 - b2 - b3)·P0 + b1·P1 + b2·P2 + b3·P3
+        // = P0 - b1·P0 - b2·P0 - b3·P0 + b1·P1 + b2·P2 + b3·P3
+        // = P0 + b1·(P1 - P0) + b2·(P2 - P0) + b3·(P3 - P0)
+
+        p0.add(&p1.sub(&p0).mul(b1))
+            .add(&p2.sub(&p0).mul(b2))
+            .add(&p3.sub(&p0).mul(b3))
+    }
+
+    /// Returns the gradient, or velocity, vector of `self` at the given *t*
+    /// value.
+    ///
+    /// Values of *t* outside the interval [0, 1] are accepted and extrapolate
+    /// the curve beyond the control points.
+    pub fn gradient(&self, t: f32) -> T::Diff {
+        let (t, [p0, p1, p2, p3]) = crb_segment(&self.0, t);
+        let [_0, _1, t2, t3] = [0.0, 1.0, 2.0 * t, 3.0 * t * t];
+
+        // ⎛ b0 ⎞        ⎛ -1 + 2·t2 -   t3 ⎞
+        // ⎜ b1 ⎟ =  1/2 ⎜    - 5·t2 + 3·t3 ⎟
+        // ⎜ b2 ⎟        ⎜  1 + 4·t2 - 3·t3 ⎟
+        // ⎝ b3 ⎠        ⎝    -   t2 +   t3 ⎠
+
+        let b1 = -2.5 * t2 + 1.5 * t3;
+        let b2 = 0.5 + 2.0 * t2 - 1.5 * t3;
+        let b3 = -0.5 * t2 + 0.5 * t3;
+
+        // b0 + b1 + b2 + b3 = 0 <=> b0 = -(b1 + b2 + b3)
+        //
+        //    b0·P0 + b1·P1 + b2·P2 + b3·P3
+        // = -(b1 + b2 + b3)·P0 + b1·P1 + b2·P2 + b3·P3
+        // = b1·P1 + b2·P2 + b3·P3 - b1·P0 - b2·P0 - b3·P0
+        // = b1·(P1 - P0) + b2·(P2 - P0) + b3·(P3 - P0)
+
+        p1.sub(&p0)
+            .mul(b1)
+            .add(&p2.sub(&p0).mul(b2))
+            .add(&p3.sub(&p0).mul(b3))
+    }
+}
+
+impl<T> BSpline<T>
+where
+    T: Affine<Diff: Linear<Scalar = f32>> + Clone,
+{
+    const _CHAR_MAT: Mat4 = {
+        const _1_6: f32 = 1.0 / 6.0;
+        const _2_3: f32 = 2.0 / 3.0;
+        mat![
+            _1_6, _2_3, _1_6,  0.0;
+            -0.5,  0.0,  0.5,  0.0;
+             0.5, -1.0,  0.5,  0.0;
+           -_1_6,  0.5, -0.5, _1_6;
+        ]
+    };
+
+    pub fn new(pts: impl IntoIterator<Item = T>) -> Self {
+        let pts: Vec<_> = pts.into_iter().collect();
+        assert!(pts.len() >= 4, "a B-spline requires at least four points");
+        Self(pts)
+    }
+
+    /// Returns the point of `self` at the given *t* value.
+    ///
+    /// Values of *t* outside the interval [0, 1] are accepted and extrapolate
+    /// the curve beyond the control points.
+    pub fn eval(&self, t: f32) -> T {
+        let (t, [p0, p1, p2, p3]) = crb_segment(&self.0, t);
+        let [_0, t1, t2, t3] = [1.0, t, t * t, t * t * t];
+
+        let _0 = (1.0 - 3.0 * t1 + 3.0 * t2 + t3) / 6.0;
+        let b1 = (4.0 - 6.0 * t2 + 3.0 * t3) / 6.0;
+        let b2 = (1.0 + 3.0 * t1 + 3.0 * t2 - 3.0 * t3) / 6.0;
+        let b3 = t3 / 6.0;
+
+        p0.add(&p1.sub(&p0).mul(b1))
+            .add(&p2.sub(&p0).mul(b2))
+            .add(&p3.sub(&p0).mul(b3))
+    }
+
+    /// Returns the gradient, or velocity, vector of `self` at the given *t*
+    /// value.
+    ///
+    /// Values of *t* outside the interval [0, 1] are accepted and extrapolate
+    /// the curve beyond the control points.
+    pub fn gradient(&self, t: f32) -> T::Diff {
+        let (t, [p0, p1, p2, p3]) = crb_segment(&self.0, t);
+        let [_0, _1, t2, t3] = [0.0, 1.0, 2.0 * t, 3.0 * t * t];
+
+        // ⎛ b0 ⎞        ⎛ -3 + 3·t2 -   t3 ⎞
+        // ⎜ b1 ⎟ = 1/6  ⎜    - 6·t2 + 3·t3 ⎟
+        // ⎜ b2 ⎟        ⎜  3 + 3·t2 - 3·t3 ⎟
+        // ⎝ b3 ⎠        ⎝               t3 ⎠
+
+        let b1 = -t2 + 0.5 * t3;
+        let b2 = 0.5 + 0.5 * t2 - 0.5 * t3;
+        let b3 = t3 / 6.0;
+
+        // b0 + b1 + b2 + b3 = 0 <=> b0 = -(b1 + b2 + b3)
+        //
+        //   b0·P0 + b1·P1 + b2·P2 + b3·P3
+        // = -(b1 + b2 + b3)·P0 + b1·P1 + b2·P2 + b3·P3
+        // = b1·P1 + b2·P2 + b3·P3 - b1·P0 - b2·P0 - b3·P0
+        // = b1·(P1 - P0) + b2·(P2 - P0) + b3·(P3 - P0)
+
+        p1.sub(&p0)
+            .mul(b1)
+            .add(&p2.sub(&p0).mul(b2))
+            .add(&p3.sub(&p0).mul(b3))
+    }
+}
+
+/// Returns the curve segment and local *t* value corresponding to
+/// the global *t* value of a Catmull–Rom or B-spline.
+fn crb_segment<T: Clone>(pts: &[T], t: f32) -> (f32, &[T; 4]) {
+    let t = 1.0 + t * (pts.len() as f32 - 3.0);
+
+    let i = (t as usize).clamp(1, pts.len() - 3);
+    let u = t - i as f32;
+    let pts = pts[i - 1..i + 3] // OK: 1 <= i < len - 3
+        .try_into()
+        .expect("3 - (-1) = 4");
+    (u, pts)
+}
+
+impl<Spl> Euclidean<Spl> {
+    /// Creates a new `Euclidean` wrapper of the given spline.
+    pub fn new<B, const N: usize>(spline: Spl) -> Self
+    where
+        Spl: Parametric<Point<[f32; N], Real<N, B>>>,
+    {
+        let mut lut = Vec::new();
+        let mut s = 0.0;
+        let mut p0 = spline.eval(0.0);
+
+        // TODO smarter sampling
+        for t in 0.0.vary_to(1.0, 256) {
+            lut.push((s, t));
+            let p1 = spline.eval(t);
+            s += p1.distance(&p0);
+            p0 = p1;
+        }
+        Self(spline, lut)
+    }
+
+    /// Returns the point of `self` at distance *s* from the start,
+    /// as measured along the curve.
+    pub fn eval<T>(&self, s: f32) -> T
+    where
+        Spl: Parametric<T>,
+    {
+        self.0.eval(self.t(s))
+    }
+
+    /// Returns the approximate arc length of the spline.
+    pub fn len(&self) -> f32 {
+        self.1.last().map_or(0.0, |m| m.0)
+    }
+
+    /// Returns the approximate *t* value for the given *s*.
+    fn t(&self, s: f32) -> f32 {
+        let lut = &self.1;
+        let i = lut.binary_search_by(|x| x.0.total_cmp(&s));
+        match i {
+            Ok(i) => lut[i].1,
+            Err(0) => lut[0].1,
+            Err(i) if i == lut.len() => lut[lut.len() - 1].1,
+            Err(i) => {
+                let (s0, t0) = lut[i - 1]; // Ok: 0 < i
+                let (s1, t1) = lut[i]; // Ok: i < lut.len()
+                // Interpolate t in [t0, t1] given s in [s0, s1]
+                t0.lerp(&t1, inv_lerp(s, s0, s1))
+            }
+        }
+    }
+}
+
 //
 // Local trait impls
 //
@@ -547,6 +785,15 @@ where
     }
 }
 
+impl<T> Parametric<T> for CubicHermite<T, T::Diff>
+where
+    T: Affine<Diff: Linear<Scalar = f32> + Clone> + Clone,
+{
+    fn eval(&self, t: f32) -> T {
+        self.eval(t)
+    }
+}
+
 impl<T> Parametric<T> for BezierSpline<T>
 where
     T: Affine<Diff: Linear<Scalar = f32> + Clone> + Clone,
@@ -558,13 +805,40 @@ where
 
 impl<T> Parametric<T> for HermiteSpline<T>
 where
-    T: Affine<Diff: Linear<Scalar = f32> + Clone> + Lerp,
+    T: Affine<Diff: Linear<Scalar = f32> + Clone> + Clone,
+{
+    fn eval(&self, t: f32) -> T {
+        self.eval(t)
+    }
+}
+
+impl<T> Parametric<T> for CatmullRomSpline<T>
+where
+    T: Affine<Diff: Linear<Scalar = f32> + Clone> + Clone,
+{
+    fn eval(&self, t: f32) -> T {
+        self.eval(t)
+    }
+}
+
+impl<T> Parametric<T> for BSpline<T>
+where
+    T: Affine<Diff: Linear<Scalar = f32> + Clone> + Clone,
 {
     fn eval(&self, t: f32) -> T {
         self.eval(t)
     }
 }
 
+impl<T, Spl> Parametric<T> for Euclidean<Spl>
+where
+    Spl: Parametric<T>,
+{
+    fn eval(&self, s: f32) -> T {
+        self.eval(s)
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use crate::assert_approx_eq;
@@ -832,4 +1106,76 @@ mod tests {
 
         assert_eq!(expected, actual);
     }
+
+    #[test]
+    fn catmull_rom_spline_point2_eval() {
+        #[rustfmt::skip]
+        let c = CatmullRomSpline::<Point2>::new([
+            pt2(-1.0, 0.0), pt2(0.0, 0.0), pt2(1.0, 0.0),
+            pt2(0.0, 1.0), pt2(1.0, 1.0), pt2(2.0, 1.0),
+        ]);
+
+        #[rustfmt::skip]
+        let expected = [
+            [33.0, -18.0], [0.0, 0.0], [0.890625, -0.0703125], [0.5, 0.5],
+            [0.109375, 1.0703125], [1.0, 1.0], [-32.0, 19.0]
+        ];
+        let actual = TEST_T_VALS.map(|t| c.eval(t).0);
+
+        assert_eq!(expected, actual);
+    }
+
+    #[test]
+    fn catmull_rom_spline_point2_gradient() {
+        #[rustfmt::skip]
+        let c = CatmullRomSpline::<Point2>::new([
+            pt2(-1.0, 0.0), pt2(0.0, 0.0), pt2(1.0, 0.0),
+            pt2(0.0, 1.0), pt2(1.0, 1.0), pt2(2.0, 1.0),
+        ]);
+
+        #[rustfmt::skip]
+        let expected = [
+            [-32.0, 16.5], [1.0, 0.0], [0.8125, 0.09375], [-1.5, 1.25],
+            [0.8125, 0.09375], [1.0, 0.0], [-32.0, 16.5]
+        ];
+        let actual = TEST_T_VALS.map(|t| c.gradient(t).0);
+
+        assert_eq!(expected, actual);
+    }
+
+    #[test]
+    fn b_spline_point2_eval() {
+        #[rustfmt::skip]
+        let b = BSpline::<Point2>::new([
+            pt2(-1.0, 0.0), pt2(0.0, 0.0), pt2(1.0, 0.0),
+            pt2(0.0, 1.0), pt2(1.0, 1.0), pt2(2.0, 1.0),
+        ]);
+
+        #[rustfmt::skip]
+        let expected = [
+            [6.0, -4.5], [0.0, 0.0], [0.609375, 0.0703125], [0.5, 0.5],
+            [0.39062497, 0.92968756], [1.0, 1.0], [-5.0, 5.5]
+        ];
+        let actual = TEST_T_VALS.map(|t| b.eval(t).0);
+
+        assert_approx_eq!(expected, actual);
+    }
+
+    #[test]
+    fn b_spline_point2_gradient() {
+        #[rustfmt::skip]
+        let b = BSpline::<Point2>::new([
+            pt2(-1.0, 0.0), pt2(0.0, 0.0), pt2(1.0, 0.0),
+            pt2(0.0, 1.0), pt2(1.0, 1.0), pt2(2.0, 1.0),
+        ]);
+
+        #[rustfmt::skip]
+        let expected = [
+            [-8.0, 4.5], [1.0, 0.0], [0.4375, 0.28125], [-0.5, 0.75],
+            [0.4375, 0.28125], [1.0, 0.0], [-8.0, 4.5]
+        ];
+        let actual = TEST_T_VALS.map(|t| b.gradient(t).0);
+
+        assert_eq!(expected, actual);
+    }
 }
diff --git a/core/src/math/vec.rs b/core/src/math/vec.rs
index 5809d980..dcc91f03 100644
--- a/core/src/math/vec.rs
+++ b/core/src/math/vec.rs
@@ -130,7 +130,6 @@ impl<R, Sp> Vector<R, Sp> {
 // TODO Many of these functions could be more generic
 impl<Sp, const N: usize> Vector<[f32; N], Sp> {
     /// Returns the length (magnitude) of `self`.
-    #[cfg(feature = "fp")]
     #[inline]
     pub fn len(&self) -> f32 {
         super::float::f32::sqrt(self.dot(self))
@@ -363,8 +362,8 @@ impl<Sc: Copy, Sp, const N: usize> Vector<[Sc; N], Sp> {
     /// ```
     #[inline]
     #[must_use]
-    pub fn map<T>(self, mut f: impl FnMut(Sc) -> T) -> Vector<[T; N], Sp> {
-        array::from_fn(|i| f(self.0[i])).into()
+    pub fn map<T>(self, f: impl FnMut(Sc) -> T) -> Vector<[T; N], Sp> {
+        self.0.map(f).into()
     }
     /// Returns a vector of the same dimension as `self` by applying `f`
     /// component-wise to `self` and `other`.
@@ -387,6 +386,70 @@ impl<Sc: Copy, Sp, const N: usize> Vector<[Sc; N], Sp> {
     ) -> Vector<[U; N], Sp> {
         array::from_fn(|i| f(self.0[i], other.0[i])).into()
     }
+
+    #[inline]
+    pub fn max(&self) -> Sc
+    where
+        Sc: PartialOrd,
+    {
+        const { assert!(N > 0, "0D vectors have no maximum") }
+        let mut max = self.0[0];
+        for c in self.0 {
+            if c > max {
+                max = c;
+            }
+        }
+        max
+    }
+
+    #[inline]
+    pub fn min(&self) -> Sc
+    where
+        Sc: PartialOrd,
+    {
+        const { assert!(N > 0, "0D vectors have no minimum") }
+        let mut min = self.0[0];
+        for c in self.0 {
+            if c < min {
+                min = c;
+            }
+        }
+        min
+    }
+
+    #[inline]
+    pub fn argmax(&self) -> usize
+    where
+        Sc: PartialOrd,
+    {
+        const { assert!(N > 0, "0D vectors have no maximum") }
+        let mut max = self.0[0];
+        let mut max_i = 0;
+        for (c, i) in zip(self.0, 0..) {
+            if c > max {
+                max = c;
+                max_i = i;
+            }
+        }
+        max_i
+    }
+
+    #[inline]
+    pub fn argmin(&self) -> usize
+    where
+        Sc: PartialOrd,
+    {
+        const { assert!(N > 0, "0D vectors have no minimum") }
+        let mut min = self.0[0];
+        let mut min_i = 0;
+        for (c, i) in zip(self.0, 0..) {
+            if c < min {
+                min = c;
+                min_i = i;
+            }
+        }
+        min_i
+    }
 }
 
 impl<Sc: Copy, B> Vector<[Sc; 2], Real<2, B>> {
@@ -438,17 +501,17 @@ impl<B> Vec2<B> {
     ///
     /// where *θ* is the (signed) angle between **a** and **b**. In particular,
     /// the result is zero if **a** and **b** are parallel (or either is zero),
-    /// positive if the angle from **a** to **b** is positive, and negative if
-    /// the angle is negative:
+    /// negative if the angle from **a** to **b** is negative (clockwise), and
+    /// positive if the angle is positive (counter-clockwise)
     ///
     /// ```text
-    ///       ^ b               ^ a
-    ///      /           ^ b   /        ^ a
-    ///     ^ a           \   /          \
-    ///    /               \ /            \
-    ///   O                 O              O-----> b
+    ///       ^ b
+    ///      /                 a ^        ^ a
+    ///     ^ a                 /          \
+    ///    /                   /            \
+    ///   O           b <-----O              O-----> b
     ///
-    ///  a⟂·b = 0        a⟂·b > 0       a⟂·b < 0
+    ///  a⟂·b = 0        a⟂·b > 0         a⟂·b < 0
     /// ```
     ///
     /// # Examples
@@ -511,15 +574,13 @@ where
     /// proportional to the area of the parallelogram formed by the vectors.
     /// Specifically, the length is given by the identity:
     ///
-    /// ```text
-    ///     |𝗮 × 𝗯| = |𝗮| |𝗯| sin 𝜽
-    /// ```
+    /// |**a** × **b**| = |**a**| |**b**| sin *θ*,
     ///
     /// where |·| denotes the length of a vector and 𝜽 equals the angle
-    /// between 𝗮 and 𝗯. Specifically, the result has unit length if 𝗮 and 𝗯
-    /// are orthogonal and |𝗮| = |𝗯| = 1. The cross product can be used to
-    /// produce an *orthonormal basis* from any two non-parallel non-zero
-    /// 3-vectors.
+    /// between **a** and **b**. Specifically, the result has unit length if
+    /// **a** and **b** are orthogonal and |**a**| = |**b**| = 1. The cross
+    /// product can be used to produce an *orthonormal basis* from any two
+    /// non-parallel non-zero 3-vectors.
     ///
     /// ```text
     ///        ^
@@ -729,7 +790,7 @@ where
     /// Note that `Self::DIM` can be less than the number of elements in `R`.
     #[inline]
     fn index(&self, i: usize) -> &Self::Output {
-        assert!(i < Self::DIM, "index {i} out of bounds ({})", Self::DIM);
+        debug_assert!(i < Self::DIM, "index {i} out of bounds ({})", Self::DIM);
         &self.0[i]
     }
 }
@@ -746,7 +807,7 @@ where
     /// Note that `Self::DIM` can be less than the number of elements in `R`.
     #[inline]
     fn index_mut(&mut self, i: usize) -> &mut Self::Output {
-        assert!(i < Self::DIM, "index {i} out of bounds ({})", Self::DIM);
+        debug_assert!(i < Self::DIM, "index {i} out of bounds ({})", Self::DIM);
         &mut self.0[i]
     }
 }
@@ -916,7 +977,6 @@ mod tests {
     mod f32 {
         use super::*;
 
-        #[cfg(feature = "fp")]
         #[test]
         fn length() {
             assert_approx_eq!(vec2(1.0, 1.0).len(), SQRT_2);
@@ -1043,8 +1103,8 @@ mod tests {
 
         #[test]
         fn perp() {
-            assert_eq!(Vec2::<()>::zero().perp(), Vec2::zero());
-            assert_eq!(Vec2::<()>::X.perp(), Vec2::Y);
+            assert_eq!(<Vec2>::zero().perp(), <Vec2>::zero());
+            assert_eq!(<Vec2>::X.perp(), Vec2::Y);
             assert_eq!(vec2(-0.2, -1.5).perp(), vec2(1.5, -0.2));
         }
 
diff --git a/core/tests/rendering.rs b/core/tests/rendering.rs
index 2a144736..314537d7 100644
--- a/core/tests/rendering.rs
+++ b/core/tests/rendering.rs
@@ -44,8 +44,8 @@ fn textured_quad() {
     assert_eq!(framebuf[255][0], rgb(0x7F, 0, 0));
     assert_eq!(framebuf[0][255], rgb(0x7F, 0, 0));
 
-    let comp = *include_bytes!("textured_quad.ppm");
-    let comp = parse_pnm(comp).expect("should be a valid ppm");
+    static COMP: &[u8] = include_bytes!("textured_quad.ppm");
+    let comp = parse_pnm(COMP.iter().copied()).expect("should be a valid ppm");
 
     assert_eq!(framebuf, comp);
 
diff --git a/demos/src/bin/bezier.rs b/demos/src/bin/bezier.rs
index 8558ffda..26823715 100644
--- a/demos/src/bin/bezier.rs
+++ b/demos/src/bin/bezier.rs
@@ -2,10 +2,12 @@ use core::ops::ControlFlow::Continue;
 
 use re::prelude::*;
 
-use re::core::geom::{Edge, Ray};
-use re::core::math::rand::{Distrib, Uniform, VectorsOnUnitDisk, Xorshift64};
-use re::core::math::spline::approximate;
-use re::core::render::raster::line;
+use re::core::{
+    geom::{Edge, Ray},
+    math::rand::{Distrib, Uniform, VectorsOnUnitDisk, Xorshift64},
+    math::spline::{Euclidean, approximate},
+    render::raster::line,
+};
 use re::front::{Frame, dims, minifb::Window};
 
 fn main() {
@@ -25,18 +27,18 @@ fn main() {
     let vel = VectorsOnUnitDisk;
 
     let mut pos_vels: Vec<(Point2, Vec2)> =
-        (pos, vel).samples(rng).take(32).collect();
+        (pos, vel).samples(rng).take(16).collect();
 
     // Disable some unneeded things
     win.ctx.color_clear = None;
     win.ctx.depth_clear = None;
 
-    win.run(|Frame { dt, buf, .. }| {
+    win.run(|Frame { t, dt, buf, .. }| {
         let buf = &mut buf.borrow_mut().color_buf.buf;
 
         // Fade out previous frame a bit
         buf.iter_mut()
-            .for_each(|c| *c = c.saturating_sub(0x08_08_02));
+            .for_each(|c| *c = c.saturating_sub(0xFF_FF_FF));
 
         let rays: Vec<Ray<_>> = pos_vels
             .chunks(2)
@@ -50,10 +52,21 @@ fn main() {
         for Edge(p0, p1) in approx.edges() {
             let vs = [p0, p1].map(|p| vertex(p.to_pt3().to(), ()));
             line(vs, |sl| {
-                buf[sl.y][sl.xs].fill(0xFF_FF_FF);
+                buf[sl.y][sl.xs].fill(0x33_00_00);
             })
         }
 
+        let euc = Euclidean::new(b);
+        for s in 0.0.vary_to(euc.len(), 32) {
+            let p = euc.eval((s + 40.0 * t.as_secs_f32()) % euc.len());
+            let x = p.x() as usize;
+            let y = p.y() as usize;
+            buf[y - 1][x] = 0xFF_FF_FF;
+            buf[y][x - 1] = 0xFF_FF_FF;
+            buf[y + 1][x] = 0xFF_FF_FF;
+            buf[y][x + 1] = 0xFF_FF_FF;
+        }
+
         let dt = dt.as_secs_f32();
         for (pos, vel) in &mut pos_vels {
             *pos = (*pos + 80.0 * *vel * dt).clamp(&min, &max);