diff --git a/cspell.json b/cspell.json
index 93512ffe1..902e3d08e 100644
--- a/cspell.json
+++ b/cspell.json
@@ -14,6 +14,7 @@
     "eisel",
     "elementwise",
     "elif",
+    "expf",
     "extrinsics",
     "frameless",
     "freelist",
@@ -27,6 +28,7 @@
     "indexmap",
     "itertools",
     "lemire",
+    "logf",
     "miri",
     "msun",
     "muls",
diff --git a/root/numeric/F32/F32.vi b/root/numeric/F32/F32.vi
index 155192b9e..73a2a1948 100644
--- a/root/numeric/F32/F32.vi
+++ b/root/numeric/F32/F32.vi
@@ -5,7 +5,12 @@ use util::{duplicate, erase};
 #[builtin = "F32"]
 pub type F32;
 
+mod common;
+
 pub mod F32 {
+  pub mod exp;
+  pub mod ln;
+
   pub const nan: F32 = 0.0 / 0.0;
   pub const inf: F32 = 1.0 / 0.0;
   pub const neg_inf: F32 = -inf;
diff --git a/root/numeric/F32/common.vi b/root/numeric/F32/common.vi
new file mode 100644
index 000000000..07a397050
--- /dev/null
+++ b/root/numeric/F32/common.vi
@@ -0,0 +1,51 @@
+
+const num_mantissa_bits: N32 = Float::num_mantissa_bits[F32];
+const max_exponent: I32 = Float::max_exponent[F32];
+const exp_mask: N32 = (1[N32] << Float::num_exponent_bits[F32]) - 1;
+
+/// normalize returns a normal number y and exponent exp
+/// satisfying x == y × 2**exp. It assumes x is finite and non-zero.
+pub fn .normalize(x: F32) -> (F32, I32) {
+  if x.to_bits() & (exp_mask << num_mantissa_bits) == 0 {
+    return (x * (1 << num_mantissa_bits) as F32, -num_mantissa_bits);
+  }
+  return (x, +0);
+}
+
+/// Calculates frac * 2 ^ exp.
+pub fn apply_exponent(frac: F32, exp: I32) -> F32 {
+  // special cases
+  if frac == 0.0 or frac == F32::inf or frac == F32::neg_inf or frac.is_nan() {
+    return frac;
+  }
+  let (frac, e) = frac.normalize();
+  exp += e;
+  let x = frac.to_bits();
+  exp += ((x >> num_mantissa_bits) as N32 & exp_mask) as I32 - max_exponent;
+  if exp < -149 {
+    // underflow
+    return if frac < 0.0 {
+      -0.0
+    } else {
+      0.0
+    };
+  }
+  if exp > +127 {
+    // overflow
+    return if frac < 0.0 {
+      F32::neg_inf
+    } else {
+      F32::inf
+    };
+  }
+  let m = 1.0;
+  if exp < -126 {
+    // denormalize
+    exp += +24
+    // 2**-24
+    m = 1.0 / (1 << 24) as F32
+  }
+  x &= !(exp_mask << num_mantissa_bits);
+  x |= (exp + max_exponent) as N32 << num_mantissa_bits;
+  return m * F32::from_bits(x);
+}
diff --git a/root/numeric/F32/exp.vi b/root/numeric/F32/exp.vi
new file mode 100644
index 000000000..a8347364b
--- /dev/null
+++ b/root/numeric/F32/exp.vi
@@ -0,0 +1,82 @@
+
+// Adapted from https://github.com/rust-lang/compiler-builtins/blob/libm-v0.2.16/libm/src/math/exp.rs
+//
+// origin: FreeBSD /usr/src/lib/msun/src/e_expf.c
+//
+// Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+
+const ln2_hi: F32 = 6.9314575195e-01;
+const ln2_lo: F32 = 1.4286067653e-06;
+const inv_ln2: F32 = 1.4426950216e+00;
+
+const P1: F32 = 1.6666625440e-1;
+const P2: F32 = -2.7667332906e-3;
+
+const overflow: F32 = 88.7228394[F32];
+const underflow: F32 = -103.972084[F32];
+
+/// Exponential, base *e* (F32)
+///
+/// Calculate the exponential of `f`, that is, *e* raised to the power `f`
+/// (where *e* is the base of the natural system of logarithms, approximately 2.71828).
+pub fn .exp(f: F32) -> F32 {
+  when {
+    f.is_nan() { f }
+    f >= overflow { F32::inf }
+    f <= underflow { 0.0 }
+    _ {
+      let hf = f.to_bits();
+      let is_negative = (hf >> 31) == 1;
+      // high word of |f|
+      hf &= 0x7fffffff;
+
+      // argument reduction
+      when {
+        hf > 0x3eb17218 {
+          // if |f| > 0.5 ln2
+          let k = when {
+            hf > 0x3f851592 {
+              if is_negative {
+                // if f < 1.5 ln2
+                (inv_ln2 * f - 0.5[F32]) as I32
+              } else {
+                // if f > 1.5 ln2
+                (inv_ln2 * f + 0.5[F32]) as I32
+              }
+            }
+            is_negative { -1 }
+            _ { +1 }
+          };
+          let kf = k as F32;
+          let hi = f - kf * ln2_hi;
+          // k*ln2hi is exact here
+          let lo = kf * ln2_lo;
+          f = hi - lo;
+          exp_poly_approximation(f, hi, lo, k)
+        }
+        hf > 0x39000000 { exp_poly_approximation(f, f, 0.0, +0) }
+        _ { 1.0 + f }
+      }
+    }
+  }
+}
+
+fn exp_poly_approximation(f: F32, hi: F32, lo: F32, k: I32) -> F32 {
+  let ff = f * f;
+  let c = f - ff * (P1 + ff * P2);
+  let y = 1.0 + (f * c / (2.0 - c) - lo + hi);
+  if k == +0 {
+    y
+  } else {
+    common::apply_exponent(y, k)
+  }
+}
diff --git a/root/numeric/F32/ln.vi b/root/numeric/F32/ln.vi
new file mode 100644
index 000000000..fb60b9e2c
--- /dev/null
+++ b/root/numeric/F32/ln.vi
@@ -0,0 +1,73 @@
+
+// Adapted from https://github.com/rust-lang/compiler-builtins/blob/libm-v0.2.16/libm/src/math/logf.rs
+//
+// origin: FreeBSD /usr/src/lib/msun/src/e_logf.c */
+//
+// Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+
+const ln2_hi: F32 = 6.9313812256e-01;
+const ln2_lo: F32 = 9.0580006145e-06;
+// |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]).
+const lg1: F32 = 0.66666662693;
+const lg2: F32 = 0.40000972152;
+const lg3: F32 = 0.28498786688;
+const lg4: F32 = 0.24279078841;
+
+/// The natural logarithm of `f` (F32).
+pub fn .ln(f: F32) -> F32 {
+  // 0x1p25f === 2 ^ 25
+  let x1p25 = F32::from_bits(0x4c000000);
+  let bits = f.to_bits();
+  let k = +0;
+
+  when {
+    bits < 0x00800000 or (bits >> 31) != 0 {
+      // f < 2**-126 (f small or negative)
+      if bits << 1 == 0 {
+        // log(+-0)=-inf
+        return -1.0 / (f * f);
+      }
+      if bits as I32 < +0 {
+        return F32::nan;
+      }
+      // subnormal number, scale up f
+      k -= +25;
+      f *= x1p25;
+      bits = f.to_bits();
+    }
+    bits >= 0x7f800000 {
+      // f is infinite or NaN
+      return f;
+    }
+    bits == 0x3f800000 {
+      // f == 1.0
+      return 0.0;
+    }
+  }
+
+  // reduce f into [sqrt(2)/2, sqrt(2)]
+  bits += (0x3f800000 - 0x3f3504f3);
+  k += ((bits >> 23) as I32) - +0x7f;
+  bits = (bits & 0x007fffff) + 0x3f3504f3;
+  f = F32::from_bits(bits);
+
+  f -= 1.0;
+  let s = f / (2.0 + f);
+  let z = s * s;
+  let w = z * z;
+  let t1 = w * (lg2 + w * lg4);
+  let t2 = z * (lg1 + w * lg3);
+  let r = t2 + t1;
+  let hfsq = 0.5 * f * f;
+  let dk = k as F32;
+  s * (hfsq + r) + dk * ln2_lo - hfsq + f + dk * ln2_hi
+}
diff --git a/root/numeric/F64/common.vi b/root/numeric/F64/common.vi
index 1af5062dd..2eae99fc2 100644
--- a/root/numeric/F64/common.vi
+++ b/root/numeric/F64/common.vi
@@ -5,15 +5,12 @@ pub const log2e: F64 = 1.44269504088896338700e+00[F64];
 
 const num_mantissa_bits: N32 = Float::num_mantissa_bits[F64];
 const max_exponent: I32 = Float::max_exponent[F64];
-
 const exp_mask: N32 = (1[N32] << Float::num_exponent_bits[F64]) - 1;
 
-// normalize returns a normal number y and exponent exp
-// satisfying x == y × 2**exp. It assumes x is finite and non-zero.
-pub fn normalize(x: F64) -> (F64, I32) {
-  // 2**-1022
-  const SmallestNormal: F64 = 2.2250738585072014e-308[F64];
-  if x.abs() < SmallestNormal {
+/// normalize returns a normal number y and exponent exp
+/// satisfying x == y × 2**exp. It assumes x is finite and non-zero.
+pub fn .normalize(x: F64) -> (F64, I32) {
+  if (x.to_bits() >> num_mantissa_bits) as N32 & exp_mask == 0 {
     return (x * (1[N64] << 52) as F64, -52);
   }
   return (x, +0);
@@ -33,7 +30,7 @@ pub fn fraction_exponent(f: F64) -> { frac: F64, exp: I32 } {
   if f == 0.0[F64] or f.is_nan() or f == F64::inf or f == F64::neg_inf {
     return { frac: f, exp: +0 };
   }
-  let (frac, exp) = normalize(f);
+  let (frac, exp) = f.normalize();
   let x = frac.to_bits();
   exp += ((x >> num_mantissa_bits) as N32 & exp_mask) as I32 - (max_exponent - +1);
   x &= !(exp_mask as N64 << num_mantissa_bits);
@@ -48,7 +45,7 @@ pub fn apply_exponent(frac: F64, exp: I32) -> F64 {
   if frac == 0.0[F64] or frac == F64::inf or frac == F64::neg_inf or frac.is_nan() {
     return frac;
   }
-  let (frac, e) = normalize(frac);
+  let (frac, e) = frac.normalize();
   exp += e;
   let x = frac.to_bits();
   exp += ((x >> num_mantissa_bits) as N32 & exp_mask) as I32 - max_exponent;
diff --git a/tests/programs/f32_ops.vi b/tests/programs/f32_ops.vi
index ce33aa161..9f613bd13 100644
--- a/tests/programs/f32_ops.vi
+++ b/tests/programs/f32_ops.vi
@@ -1,4 +1,9 @@
 
+use #root::rng::Pcg32;
+
+#[configurable]
+const max: N32 = 10_000;
+
 pub fn main(&io: &IO) {
   let interesting_floats = [
     +0.0,
@@ -30,10 +35,15 @@ pub fn main(&io: &IO) {
     1.0e10,
   ];
 
+  let rng = Pcg32::seeded("how many logs could a log calculator calculate if a log calculator could calculate logs");
+  for _ in 0..max {
+    interesting_floats.push_back(F32::from_bits(rng.gen_n32()));
+  }
+
   for f in interesting_floats {
-    let sqrt = f.sqrt();
-    io.println("{f} = {f.to_bits().to_hex()}");
-    io.println("sqrt({f}) = {sqrt} ({sqrt.to_bits().to_hex()})");
+    io.println("sqrt {f.to_bits()} {f.sqrt().to_bits()}");
+    io.println("exp {f.to_bits()} {f.exp().to_bits()}");
+    io.println("ln {f.to_bits()} {f.ln().to_bits()}");
     io.println("");
   }
 }
diff --git a/tests/programs/f64_ops.vi b/tests/programs/f64_ops.vi
index 1b799b856..539139786 100644
--- a/tests/programs/f64_ops.vi
+++ b/tests/programs/f64_ops.vi
@@ -1,4 +1,9 @@
 
+use #root::rng::Pcg32;
+
+#[configurable]
+const max: N32 = 10_000;
+
 pub fn main(&io: &IO) {
   let interesting_floats = [
     +0.0[F64],
@@ -33,14 +38,15 @@ pub fn main(&io: &IO) {
     1.0e10[F64],
   ];
 
+  let rng = Pcg32::seeded("how many logs could a log calculator calculate if a log calculator could calculate logs");
+  for _ in 0..max {
+    interesting_floats.push_back(F64::from_bits(N64(rng.gen_n32(), rng.gen_n32())));
+  }
+
   for f in interesting_floats {
-    io.println("{f} = {f.to_bits().to_hex()}");
-    let sqrt = f.sqrt();
-    io.println("sqrt({f}) = {sqrt} ({sqrt.to_bits().to_hex()})");
-    let exp = f.exp();
-    io.println("exp({f}) = {exp} ({exp.to_bits().to_hex()})");
-    let ln = f.ln();
-    io.println("ln({f}) = {ln} ({ln.to_bits().to_hex()})");
+    io.println("sqrt {f.to_bits()} {f.sqrt().to_bits()}");
+    io.println("exp {f.to_bits()} {f.exp().to_bits()}");
+    io.println("ln {f.to_bits()} {f.ln().to_bits()}");
     io.println("");
   }
 }
diff --git a/tests/snaps/programs/f32_ops/output.txt b/tests/snaps/programs/f32_ops/output.txt
deleted file mode 100644
index d1a5970a2..000000000
--- a/tests/snaps/programs/f32_ops/output.txt
+++ /dev/null
@@ -1,81 +0,0 @@
-0.0 = 0x0
-sqrt(0.0) = 0.0 (0x0)
-
--0.0 = 0x80000000
-sqrt(-0.0) = -0.0 (0x80000000)
-
--42.0 = 0xc2280000
-sqrt(-42.0) = NaN (0x7fc00000)
-
-inf = 0x7f800000
-sqrt(inf) = inf (0x7f800000)
-
--inf = 0xff800000
-sqrt(-inf) = NaN (0x7fc00000)
-
-NaN = 0x7fc00000
-sqrt(NaN) = NaN (0x7fc00000)
-
-1.40129846e-45 = 0x1
-sqrt(1.40129846e-45) = 3.74339207e-23 (0x1a3504f3)
-
-2.80259693e-45 = 0x2
-sqrt(2.80259693e-45) = 5.29395592e-23 (0x1a800000)
-
-5.87747175e-39 = 0x400000
-sqrt(5.87747175e-39) = 7.66646695e-20 (0x1fb504f3)
-
-1.17549421e-38 = 0x7fffff
-sqrt(1.17549421e-38) = 1.08420211e-19 (0x1fffffff)
-
-1.17549435e-38 = 0x800000
-sqrt(1.17549435e-38) = 1.08420217e-19 (0x20000000)
-
-1.0 = 0x3f800000
-sqrt(1.0) = 1.0 (0x3f800000)
-
-0.99999994 = 0x3f7fffff
-sqrt(0.99999994) = 0.99999994 (0x3f7fffff)
-
-1.00000012 = 0x3f800001
-sqrt(1.00000012) = 1.0 (0x3f800000)
-
-0.999999881 = 0x3f7ffffe
-sqrt(0.999999881) = 0.99999994 (0x3f7fffff)
-
-1.00000024 = 0x3f800002
-sqrt(1.00000024) = 1.00000012 (0x3f800001)
-
-0.5 = 0x3f000000
-sqrt(0.5) = 0.707106769 (0x3f3504f3)
-
-2.0 = 0x40000000
-sqrt(2.0) = 1.41421354 (0x3fb504f3)
-
-2.71828175 = 0x402df854
-sqrt(2.71828175) = 1.64872122 (0x3fd3094c)
-
-3.14159274 = 0x40490fdb
-sqrt(3.14159274) = 1.7724539 (0x3fe2dfc5)
-
-9.31322575e-10 = 0x30800000
-sqrt(9.31322575e-10) = 3.05175781e-5 (0x38000000)
-
-1.07374182e9 = 0x4e800000
-sqrt(1.07374182e9) = 32768.0 (0x47000000)
-
-2.3509887e-38 = 0x1000000
-sqrt(2.3509887e-38) = 1.53329339e-19 (0x203504f3)
-
-1.70141173e38 = 0x7effffff
-sqrt(1.70141173e38) = 1.30438176e19 (0x5f3504f3)
-
-3.40282347e38 = 0x7f7fffff
-sqrt(3.40282347e38) = 1.8446743e19 (0x5f7fffff)
-
-1.00000001e-10 = 0x2edbe6ff
-sqrt(1.00000001e-10) = 9.99999975e-6 (0x3727c5ac)
-
-1.0e10 = 0x501502f9
-sqrt(1.0e10) = 100000.0 (0x47c35000)
-
diff --git a/tests/snaps/programs/f32_ops/stats b/tests/snaps/programs/f32_ops/stats
index e4df0aae6..4657b263b 100644
--- a/tests/snaps/programs/f32_ops/stats
+++ b/tests/snaps/programs/f32_ops/stats
@@ -1,14 +1,14 @@
 
 Interactions
-  Total              983_494
-  Annihilate         441_575
-  Commute              9_350
-  Copy               135_749
-  Erase               63_894
-  Graft               95_141
-  Extrinsic          237_785
+  Total           23_617_142
+  Annihilate      10_362_109
+  Commute             20_089
+  Copy             2_411_045
+  Erase              343_879
+  Graft            2_319_560
+  Extrinsic        8_160_460
 
 Memory
-  Heap                94_480 B
-  Allocated       20_762_592 B
-  Freed           20_762_592 B
+  Heap            35_557_936 B
+  Allocated      515_474_384 B
+  Freed          515_474_384 B
diff --git a/tests/snaps/programs/f64_ops/output.txt b/tests/snaps/programs/f64_ops/output.txt
deleted file mode 100644
index 83444647c..000000000
--- a/tests/snaps/programs/f64_ops/output.txt
+++ /dev/null
@@ -1,150 +0,0 @@
-0.0 = 0x0
-sqrt(0.0) = 0.0 (0x0)
-exp(0.0) = 1.0 (0x3ff0000000000000)
-ln(0.0) = -inf (0xfff0000000000000)
-
--0.0 = 0x8000000000000000
-sqrt(-0.0) = -0.0 (0x8000000000000000)
-exp(-0.0) = 1.0 (0x3ff0000000000000)
-ln(-0.0) = -inf (0xfff0000000000000)
-
--42.0 = 0xc045000000000000
-sqrt(-42.0) = NaN (0x7ff8000000000000)
-exp(-42.0) = 5.749522264293559897e-19 (0x3c2536452ee2f75c)
-ln(-42.0) = NaN (0x7ff8000000000000)
-
-inf = 0x7ff0000000000000
-sqrt(inf) = inf (0x7ff0000000000000)
-exp(inf) = inf (0x7ff0000000000000)
-ln(inf) = inf (0x7ff0000000000000)
-
--inf = 0xfff0000000000000
-sqrt(-inf) = NaN (0x7ff8000000000000)
-exp(-inf) = 0.0 (0x0)
-ln(-inf) = NaN (0x7ff8000000000000)
-
-NaN = 0x7ff8000000000000
-sqrt(NaN) = NaN (0x7ff8000000000000)
-exp(NaN) = NaN (0x7ff8000000000000)
-ln(NaN) = NaN (0x7ff8000000000000)
-
-4.940656458412465442e-324 = 0x1
-sqrt(4.940656458412465442e-324) = 2.222758749485077483e-162 (0x1e60000000000000)
-exp(4.940656458412465442e-324) = 1.0 (0x3ff0000000000000)
-ln(4.940656458412465442e-324) = -744.4400719213812181 (0xc0874385446d71c3)
-
-9.881312916824930884e-324 = 0x2
-sqrt(9.881312916824930884e-324) = 3.143455569405257593e-162 (0x1e66a09e667f3bcd)
-exp(9.881312916824930884e-324) = 1.0 (0x3ff0000000000000)
-ln(9.881312916824930884e-324) = -743.7469247408213278 (0xc0873df9b3adb335)
-
-1.112536929253600692e-308 = 0x8000000000000
-sqrt(1.112536929253600692e-308) = 1.054768661486299963e-154 (0x1ff6a09e667f3bcd)
-exp(1.112536929253600692e-308) = 1.0 (0x3ff0000000000000)
-ln(1.112536929253600692e-308) = -709.0895657128240828 (0xc08628b76e3a7b61)
-
-2.225073858507200889e-308 = 0xfffffffffffff
-sqrt(2.225073858507200889e-308) = 1.491668146240041183e-154 (0x1fffffffffffffff)
-exp(2.225073858507200889e-308) = 1.0 (0x3ff0000000000000)
-ln(2.225073858507200889e-308) = -708.3964185322640787 (0xc086232bdd7abcd2)
-
-2.225073858507201383e-308 = 0x10000000000000
-sqrt(2.225073858507201383e-308) = 1.491668146240041349e-154 (0x2000000000000000)
-exp(2.225073858507201383e-308) = 1.0 (0x3ff0000000000000)
-ln(2.225073858507201383e-308) = -708.3964185322640787 (0xc086232bdd7abcd2)
-
-1.0 = 0x3ff0000000000000
-sqrt(1.0) = 1.0 (0x3ff0000000000000)
-exp(1.0) = 2.718281828459045091 (0x4005bf0a8b145769)
-ln(1.0) = 0.0 (0x0)
-
-0.999999999999999889 = 0x3fefffffffffffff
-sqrt(0.999999999999999889) = 0.999999999999999889 (0x3fefffffffffffff)
-exp(0.999999999999999889) = 2.718281828459044647 (0x4005bf0a8b145768)
-ln(0.999999999999999889) = -1.11022302462515654e-16 (0xbca0000000000000)
-
-1.000000000000000222 = 0x3ff0000000000001
-sqrt(1.000000000000000222) = 1.0 (0x3ff0000000000000)
-exp(1.000000000000000222) = 2.718281828459045979 (0x4005bf0a8b14576b)
-ln(1.000000000000000222) = 2.220446049250312834e-16 (0x3cafffffffffffff)
-
-0.999999999999999778 = 0x3feffffffffffffe
-sqrt(0.999999999999999778) = 0.999999999999999889 (0x3fefffffffffffff)
-exp(0.999999999999999778) = 2.718281828459044647 (0x4005bf0a8b145768)
-ln(0.999999999999999778) = -2.220446049250313574e-16 (0xbcb0000000000001)
-
-1.000000000000000444 = 0x3ff0000000000002
-sqrt(1.000000000000000444) = 1.000000000000000222 (0x3ff0000000000001)
-exp(1.000000000000000444) = 2.718281828459046423 (0x4005bf0a8b14576c)
-ln(1.000000000000000444) = 4.440892098500625176e-16 (0x3cbffffffffffffe)
-
-0.5 = 0x3fe0000000000000
-sqrt(0.5) = 0.7071067811865475727 (0x3fe6a09e667f3bcd)
-exp(0.5) = 1.648721270700128194 (0x3ffa61298e1e069c)
-ln(0.5) = -0.6931471805599452862 (0xbfe62e42fefa39ef)
-
-2.0 = 0x4000000000000000
-sqrt(2.0) = 1.414213562373095145 (0x3ff6a09e667f3bcd)
-exp(2.0) = 7.389056098930650407 (0x401d8e64b8d4ddae)
-ln(2.0) = 0.6931471805599452862 (0x3fe62e42fefa39ef)
-
-2.718281828459045091 = 0x4005bf0a8b145769
-sqrt(2.718281828459045091) = 1.648721270700128194 (0x3ffa61298e1e069c)
-exp(2.718281828459045091) = 15.15426224147926249 (0x402e4efb75e4527a)
-ln(2.718281828459045091) = 1.0 (0x3ff0000000000000)
-
-3.141592653589793116 = 0x400921fb54442d18
-sqrt(3.141592653589793116) = 1.772453850905515882 (0x3ffc5bf891b4ef6a)
-exp(3.141592653589793116) = 23.1406926327792668 (0x403724046eb09339)
-ln(3.141592653589793116) = 1.144729885849400164 (0x3ff250d048e7a1bd)
-
-9.313225746154785156e-10 = 0x3e10000000000000
-sqrt(9.313225746154785156e-10) = 3.0517578125e-5 (0x3f00000000000000)
-exp(9.313225746154785156e-10) = 1.000000000931322575 (0x3ff0000000400000)
-ln(9.313225746154785156e-10) = -20.79441541679835836 (0xc034cb5ecf0a9650)
-
-1.073741824e9 = 0x41d0000000000000
-sqrt(1.073741824e9) = 32768.0 (0x40e0000000000000)
-exp(1.073741824e9) = inf (0x7ff0000000000000)
-ln(1.073741824e9) = 20.79441541679835836 (0x4034cb5ecf0a9650)
-
-2.220446049250313081e-16 = 0x3cb0000000000000
-sqrt(2.220446049250313081e-16) = 1.490116119384765625e-8 (0x3e50000000000000)
-exp(2.220446049250313081e-16) = 1.000000000000000222 (0x3ff0000000000001)
-ln(2.220446049250313081e-16) = -36.04365338911715355 (0xc04205966f2b4f12)
-
-9.007199254740992e15 = 0x4340000000000000
-sqrt(9.007199254740992e15) = 94906265.62425155938 (0x4196a09e667f3bcd)
-exp(9.007199254740992e15) = inf (0x7ff0000000000000)
-ln(9.007199254740992e15) = 36.73680056967710073 (0x40425e4f7b2737fa)
-
-9.007199254740992e15 = 0x4340000000000000
-sqrt(9.007199254740992e15) = 94906265.62425155938 (0x4196a09e667f3bcd)
-exp(9.007199254740992e15) = inf (0x7ff0000000000000)
-ln(9.007199254740992e15) = 36.73680056967710073 (0x40425e4f7b2737fa)
-
-9.007199254740994e15 = 0x4340000000000001
-sqrt(9.007199254740994e15) = 94906265.62425155938 (0x4196a09e667f3bcd)
-exp(9.007199254740994e15) = inf (0x7ff0000000000000)
-ln(9.007199254740994e15) = 36.73680056967710073 (0x40425e4f7b2737fa)
-
-8.988465674311578541e307 = 0x7fdfffffffffffff
-sqrt(8.988465674311578541e307) = 9.480751908109175886e153 (0x5fe6a09e667f3bcc)
-exp(8.988465674311578541e307) = inf (0x7ff0000000000000)
-ln(8.988465674311578541e307) = 709.0895657128240828 (0x408628b76e3a7b61)
-
-1.797693134862315708e308 = 0x7fefffffffffffff
-sqrt(1.797693134862315708e308) = 1.340780792994259561e154 (0x5fefffffffffffff)
-exp(1.797693134862315708e308) = inf (0x7ff0000000000000)
-ln(1.797693134862315708e308) = 709.7827128933839731 (0x40862e42fefa39ef)
-
-1.000000000000000036e-10 = 0x3ddb7cdfd9d7bdbb
-sqrt(1.000000000000000036e-10) = 1.000000000000000082e-5 (0x3ee4f8b588e368f1)
-exp(1.000000000000000036e-10) = 1.000000000100000008 (0x3ff000000006df38)
-ln(1.000000000000000036e-10) = -23.02585092994045723 (0xc037069e2aa2aa5b)
-
-1.0e10 = 0x4202a05f20000000
-sqrt(1.0e10) = 100000.0 (0x40f86a0000000000)
-exp(1.0e10) = inf (0x7ff0000000000000)
-ln(1.0e10) = 23.02585092994045723 (0x4037069e2aa2aa5b)
-
diff --git a/tests/snaps/programs/f64_ops/stats b/tests/snaps/programs/f64_ops/stats
index 5ceab1558..86d027f0b 100644
--- a/tests/snaps/programs/f64_ops/stats
+++ b/tests/snaps/programs/f64_ops/stats
@@ -1,14 +1,14 @@
 
 Interactions
-  Total           16_993_196
-  Annihilate       7_748_017
-  Commute            161_350
-  Copy             2_608_277
-  Erase              611_818
-  Graft            1_425_565
-  Extrinsic        4_438_169
+  Total           79_127_100
+  Annihilate      27_206_672
+  Commute             30_089
+  Copy            14_263_643
+  Erase            4_796_515
+  Graft            5_587_975
+  Extrinsic       27_242_206
 
 Memory
-  Heap               349_248 B
-  Allocated      371_912_288 B
-  Freed          371_912_288 B
+  Heap            61_603_680 B
+  Allocated    1_634_169_136 B
+  Freed        1_634_169_136 B
diff --git a/tests/snaps/verify/f32_ops.txt b/tests/snaps/verify/f32_ops.txt
new file mode 100644
index 000000000..aeb785b3c
--- /dev/null
+++ b/tests/snaps/verify/f32_ops.txt
@@ -0,0 +1 @@
+checked 30081; 30081 correct
diff --git a/tests/snaps/verify/f64_ops.txt b/tests/snaps/verify/f64_ops.txt
new file mode 100644
index 000000000..02420113a
--- /dev/null
+++ b/tests/snaps/verify/f64_ops.txt
@@ -0,0 +1 @@
+checked 30090; 30090 correct
diff --git a/tests/snaps/verify/regex.txt b/tests/snaps/verify/regex.txt
new file mode 100644
index 000000000..6fbb038dd
--- /dev/null
+++ b/tests/snaps/verify/regex.txt
@@ -0,0 +1 @@
+
checked 10000; 669 matches
diff --git a/tests/tests.nix b/tests/tests.nix
index 3373b7a23..98a20a6cd 100644
--- a/tests/tests.nix
+++ b/tests/tests.nix
@@ -383,14 +383,16 @@ let
         else
           pkgs.runCommand "verify-${name}" { } ''
             echo verifying ${snapshot}
-            ${bin} <${snapshot}
-            touch $out
+            ${bin} <${snapshot} 2>&1 | tee $out
           '';
     in
     {
       checks."tests-verify-${name}" = check;
       verified.${snapshotName} = true;
     }
+    // pkgs.lib.optionalAttrs (!reproduce) {
+      snaps."verify/${name}.txt" = check;
+    }
   );
 
   self.apps.tests-generate-lock = flake-utils.lib.mkApp {
diff --git a/tests/verify/f32_ops.lock b/tests/verify/f32_ops.lock
new file mode 100644
index 000000000..8a9d30ad0
--- /dev/null
+++ b/tests/verify/f32_ops.lock
@@ -0,0 +1,7 @@
+# This file is automatically @generated by Cargo.
+# It is not intended for manual editing.
+version = 4
+
+[[package]]
+name = "f32_ops"
+version = "0.0.0"
diff --git a/tests/verify/f32_ops.rs b/tests/verify/f32_ops.rs
new file mode 100644
index 000000000..81d7de707
--- /dev/null
+++ b/tests/verify/f32_ops.rs
@@ -0,0 +1,45 @@
+#!/usr/bin/env -S cargo -Z script
+---
+package.edition = "2024"
+---
+
+// @snap programs/f32_ops/output.txt
+
+fn main() {
+  let mut count = 0;
+  let mut correct = 0;
+
+  for line in std::io::stdin().lines() {
+    let line = line.unwrap();
+    if line.len() == 0 {
+      continue;
+    }
+    count += 1;
+    let mut parts = line.split_whitespace();
+    let op = parts.next().unwrap();
+    let f_bits = parts.next().unwrap().parse::<u32>().unwrap();
+    let actual_bits = parts.next().unwrap().parse::<u32>().unwrap();
+    let f = f32::from_bits(f_bits);
+    let actual = f32::from_bits(actual_bits);
+    let (expected, max_ulp_diff) = match op {
+      "sqrt" => (f.sqrt(), 0),
+      "exp" => (f.exp(), 1),
+      "ln" => (f.ln(), 1),
+      _ => panic!("unknown operation: {op}"),
+    };
+    let expected_bits = expected.to_bits();
+    let ulp_diff = expected_bits.abs_diff(actual_bits);
+    if expected.is_nan() && actual.is_nan() || ulp_diff <= max_ulp_diff {
+      correct += 1;
+    } else {
+      eprintln!(
+        "incorrect {op} {f_bits:016x} ({f:?}):\n  expect: {expected_bits:016x} ({expected:?})\n  actual: {actual_bits:016x} ({actual:?})\n  ulp_diff: {ulp_diff}",
+      );
+    }
+  }
+
+  eprintln!("checked {count}; {correct} correct");
+  if correct != count {
+    std::process::exit(1);
+  }
+}
diff --git a/tests/verify/f64_ops.lock b/tests/verify/f64_ops.lock
new file mode 100644
index 000000000..cefedc4ca
--- /dev/null
+++ b/tests/verify/f64_ops.lock
@@ -0,0 +1,7 @@
+# This file is automatically @generated by Cargo.
+# It is not intended for manual editing.
+version = 4
+
+[[package]]
+name = "f64_ops"
+version = "0.0.0"
diff --git a/tests/verify/f64_ops.rs b/tests/verify/f64_ops.rs
new file mode 100644
index 000000000..ce1d9f544
--- /dev/null
+++ b/tests/verify/f64_ops.rs
@@ -0,0 +1,45 @@
+#!/usr/bin/env -S cargo -Z script
+---
+package.edition = "2024"
+---
+
+// @snap programs/f64_ops/output.txt
+
+fn main() {
+  let mut count = 0;
+  let mut correct = 0;
+
+  for line in std::io::stdin().lines() {
+    let line = line.unwrap();
+    if line.len() == 0 {
+      continue;
+    }
+    count += 1;
+    let mut parts = line.split_whitespace();
+    let op = parts.next().unwrap();
+    let f_bits = parts.next().unwrap().parse::<u64>().unwrap();
+    let actual_bits = parts.next().unwrap().parse::<u64>().unwrap();
+    let f = f64::from_bits(f_bits);
+    let actual = f64::from_bits(actual_bits);
+    let (expected, max_ulp_diff) = match op {
+      "sqrt" => (f.sqrt(), 0),
+      "exp" => (f.exp(), 1),
+      "ln" => (f.ln(), 1),
+      _ => panic!("unknown operation: {op}"),
+    };
+    let expected_bits = expected.to_bits();
+    let ulp_diff = expected_bits.abs_diff(actual_bits);
+    if expected.is_nan() && actual.is_nan() || ulp_diff <= max_ulp_diff {
+      correct += 1;
+    } else {
+      eprintln!(
+        "incorrect {op} {f_bits:016x} ({f:?}):\n  expect: {expected_bits:016x} ({expected:?})\n  actual: {actual_bits:016x} ({actual:?})\n  ulp_diff: {ulp_diff}",
+      );
+    }
+  }
+
+  eprintln!("checked {count}; {correct} correct");
+  if correct != count {
+    std::process::exit(1);
+  }
+}
diff --git a/tests/verify/regex.rs b/tests/verify/regex.rs
index 5d801a207..ee4c0f3e2 100755
--- a/tests/verify/regex.rs
+++ b/tests/verify/regex.rs
@@ -23,7 +23,7 @@ fn main() {
       re = Regex::new(&format!("^{str}$")).unwrap();
     } else {
       if re.is_match(str) != (char == 'T') {
-        println!("invalid: /{re_str}/ \"{str}\" {char}");
+        eprintln!("invalid: /{re_str}/ \"{str}\" {char}");
         invalid = true;
       }
       if char == 'T' {
@@ -32,7 +32,7 @@ fn main() {
       count += 1;
     }
   }
-  println!("\rchecked {count}; {matches} matches");
+  eprintln!("\rchecked {count}; {matches} matches");
   if invalid {
     std::process::exit(1);
   }