diff --git a/bend2/base.bend b/bend2/base.bend
index faf80db73..c02e4436f 100644
--- a/bend2/base.bend
+++ b/bend2/base.bend
@@ -1570,6 +1570,148 @@ law F32.atan2:
   for b: F32
   F32
 
+law F32.fma:
+  for a: F32
+  for b: F32
+  for c: F32
+  F32
+
+law F32.cbrt:
+  for a: F32
+  F32
+
+law F32.hypot_native:
+  for a: F32
+  for b: F32
+  F32
+
+law F32.exp2:
+  for a: F32
+  F32
+
+law F32.expm1:
+  for a: F32
+  F32
+
+law F32.exp2m1:
+  for a: F32
+  F32
+
+law F32.exp10m1:
+  for a: F32
+  F32
+
+law F32.log1p:
+  for a: F32
+  F32
+
+law F32.log2p1:
+  for a: F32
+  F32
+
+law F32.log10p1:
+  for a: F32
+  F32
+
+law F32.asinh:
+  for a: F32
+  F32
+
+law F32.acosh:
+  for a: F32
+  F32
+
+law F32.atanh:
+  for a: F32
+  F32
+
+law F32.scalbn_pos:
+  for a: F32
+  for n: U32
+  F32
+
+law F32.scalbn_neg:
+  for a: F32
+  for n: U32
+  F32
+
+law F32.logb:
+  for a: F32
+  F32
+
+law F32.roundeven:
+  for a: F32
+  F32
+
+law F32.remainder:
+  for a: F32
+  for b: F32
+  F32
+
+law F32.remquo_code:
+  for a: F32
+  for b: F32
+  U32
+
+law F32.erf:
+  for a: F32
+  F32
+
+law F32.erfc:
+  for a: F32
+  F32
+
+law F32.tgamma:
+  for a: F32
+  F32
+
+law F32.lgamma:
+  for a: F32
+  F32
+
+law F32.sinpi:
+  for a: F32
+  F32
+
+law F32.cospi:
+  for a: F32
+  F32
+
+law F32.tanpi:
+  for a: F32
+  F32
+
+law F32.asinpi:
+  for a: F32
+  F32
+
+law F32.acospi:
+  for a: F32
+  F32
+
+law F32.atanpi:
+  for a: F32
+  F32
+
+law F32.atan2pi:
+  for y: F32
+  for x: F32
+  F32
+
+law F32.rsqrt:
+  for a: F32
+  F32
+
+law F32.pown_pos:
+  for a: F32
+  for n: U32
+  F32
+
+law F32.pown_neg:
+  for a: F32
+  for n: U32
+  F32
+
 law F32.neg:
   for a: F32
   F32
@@ -1703,8 +1845,8 @@ def F32.lerp(+a: F32, b: F32, t: F32) -> F32:
 def F32.square(+a: F32) -> F32:
   F32.mul(a, a)
 
-def F32.hypot(+x: F32, +y: F32) -> F32:
-  F32.sqrt(F32.add(F32.mul(x, x), F32.mul(y, y)))
+def F32.hypot(x: F32, y: F32) -> F32:
+  F32.hypot_native(x, y)
 
 def F32.round(a: F32) -> F32:
   F32.floor(F32.add(a, 0.5))
diff --git a/bend2/comp.ts b/bend2/comp.ts
index 8576a3c74..29b7bafbe 100644
--- a/bend2/comp.ts
+++ b/bend2/comp.ts
@@ -5,6 +5,7 @@
 import * as fs from "node:fs";
 
 import * as Bend from "./bend.ts";
+import * as MathRT from "./math.ts";
 
 // Comp
 // ====
@@ -224,42 +225,7 @@ const OPERATIONS: Record<string, Intr> = Object.setPrototypeOf({
     C:  "((u64)(u32)($0))",
     JS: "Number($0 & 0xFFFFFFFFn)",
   },
-  ...tpl_ops("f32_", "add:+ sub:- mul:* div:/",
-    "f32_rewrap(f32_unbox($0) $o f32_unbox($1))", "Math.fround($0 $o $1)"),
-  f32_neg: {
-    C:  "f32_rewrap(-f32_unbox($0))",
-    JS: "(-$0)",
-  },
-  ...tpl_ops("f32_", CMPS, "((u64)(f32_unbox($0) $o f32_unbox($1)))",
-    "($0 $o $1)"),
-  ...tpl_ops("f32_", "sqrt exp log log2 log10 sin cos tan asin acos atan"
-    + " sinh cosh tanh floor ceil trunc abs:fabs:abs",
-    "f32_rewrap((f32)$o(f32_unbox($0)))", "Math.fround(Math.$o($0))"),
-  ...tpl_ops("f32_", "pow atan2",
-    "f32_rewrap((f32)$o(f32_unbox($0), f32_unbox($1)))",
-    "Math.fround(Math.$o($0, $1))"),
-  f32_mod: {
-    C:  "f32_rewrap((f32)fmod(f32_unbox($0), f32_unbox($1)))",
-    JS: "Math.fround($0 % $1)",
-  },
-  f32_to_u32: {
-    C:  "f32_to_u32($0)",
-    JS: "($0 >= 1 && $0 < 4294967296 ? Math.floor($0) : 0)",
-  },
-  f32_bits: {
-    C:  "$0",
-    JS: "f32_bits($0)",
-  },
-  f32_show: {
-    C:    "f32_show(e, $0)",
-    call: true,
-    JS:   "f32_show($0)",
-  },
-  f32_read: {
-    C:    "f32_read(e, $0)",
-    call: true,
-    JS:   "f32_read($0)",
-  },
+  ...MathRT.OPERATIONS,
   nat_add: {
     C:  "nat_chk(e, $0 + $1)",
     JS: "nat_chk($0 + $1)",
@@ -413,11 +379,7 @@ const OPTIMIZED: Record<Bend.Name, Native> = Object.setPrototypeOf({
 // Native
 // ------
 
-// sin, cos and tan are fast:: (cheap, the same pixels); the rest precise::
-const SHIMS = "sqrt exp log log2 log10 sin cos tan pow fmod".split(" ")
-  .map((n) => "#define " + n.padEnd(5) + (["sin", "cos", "tan"].includes(n)
-    ? " fast::" : " precise::") + n).join("\n")
-  + "\n#define atan2 atan2_c99";
+const SHIMS = MathRT.SHIMS;
 
 const NATIVE = {
   C: String.raw`
@@ -437,20 +399,9 @@ ${SHIMS}
 #define U32_QUO(a, b) \
   ((a) / 2 / (b) * 2 + ((a) - (a) / 2 / (b) * 2 * (b) >= (b)))
 
-INLINE f32 f32_unbox(u64 x) {
-  union { u32 u; f32 f; } p = { (u32)x };
-  return p.f;
-}
-
-INLINE u64 f32_rewrap(f32 x) {
-  union { f32 f; u32 u; } p = { x };
-  return p.u;
-}
-
-INLINE U32 f32_to_u32(U32 a) {
-  f32 v = f32_unbox(a);
-  return v >= 0.0f && v < 4294967296.0f ? (u32)v : 0;
-}
+${MathRT.C_PRE}
+${MathRT.C}
+${MathRT.C_POST}
 
 INLINE Nat nat_chk(Env e, Nat n) {
   if (n > NAT_IMM) {
@@ -464,64 +415,10 @@ INLINE Nat nat_mul(Env e, Nat a, Nat b) {
   return nat_chk(e, b != 0 && a > NAT_IMM / b ? NAT_IMM + 1 : a * b);
 }
 
-#if DEVICE
-
-#define f32_show(e, x) (err_post(e.mem, ERR_FIDS), 0)
-#define f32_read(e, s) (err_post(e.mem, ERR_FIDS), 0)
-
-#else
-
-static Term f32_show(Env e, Term x);
-static Term f32_read(Env e, Term s);
-
-#endif
+${MathRT.C_F32_IO_DECL}
 `.slice(1),
   IO: String.raw`
-static int f32_text(char* buf, f32 v) {
-  int n = 0;
-  int p = 0;
-  if (v != v) {
-    return sprintf(buf, "nan");
-  }
-  for (; p < 9; p += 1) {
-    n = snprintf(buf, 40, "%.*e", p, (double)v);
-    if (strtof(buf, NULL) == v) {
-      break;
-    }
-  }
-  char* ep = strchr(buf, 'e');
-  if (ep == NULL) {
-    return n;
-  }
-  int ex = atoi(ep + 1);
-  if (ex >= 21 || ex <= -7) {
-    n = (int)(ep - buf) + sprintf(ep, "e%c%d", ex < 0 ? '-' : '+', abs(ex));
-  } else if (ex <= p) {
-    n = snprintf(buf, 40, "%.*f", p - ex, (double)v);
-  } else {
-    int s = *buf == '-';
-    memmove(buf + s + 1, buf + s + 2, p);
-    memset(buf + s + 1 + p, '0', ex - p);
-    n = s + 1 + ex;
-  }
-  return n;
-}
-
-static Term f32_show(Env e, Term x) {
-  char buf[40];
-  return io_str(e, buf, f32_text(buf, f32_unbox(x)));
-}
-
-static Term f32_read(Env e, Term s) {
-  u64 n = 0;
-  char* text = io_cstr(e, s, &n);
-  char* end;
-  f32 v = strtof(text, &end);
-  Term out = n > 0 && (u64)(end - text) == n && strpbrk(text, "xX(") == NULL
-    ? io_box(e, CID_SOME, f32_rewrap(v)) : term_pak(CID_NONE, 0);
-  free(text);
-  return out;
-}
+${MathRT.C_IO}
 `.slice(1),
   JS: String.raw`
 function word_to_u32(w) {
@@ -558,34 +455,11 @@ function nat_chk(n) {
   return n;
 }
 
-function f32_show(x) {
-  if (x !== x) {
-    return "nan";
-  }
-  if (!Number.isFinite(x) || Object.is(x, -0)) {
-    return x < 0 ? "-inf"
-      : x === 0 ? "-0" : "inf";
-  }
-  let s = "x";
-  for (let p = 1; p <= 9 && Math.fround(Number(s)) !== x; p += 1) {
-    s = String(Number(x.toExponential(p - 1)));
-  }
-  return s;
-}
+${MathRT.JS_PRE}
 
-function f32_bits(x) {
-  return new Uint32Array(new Float32Array([x]).buffer)[0];
-}
+${MathRT.JS}
 
-function f32_from_bits(u) {
-  return new Float32Array(new Uint32Array([u]).buffer)[0];
-}
-
-function f32_read(s) {
-  const re = /^\s*[+-]?((\d+\.?\d*|\.\d+)(e[+-]?\d+)?|inf(inity)?|nan)$/i;
-  const v = Number(s.replace(/inf\w*/i, "Infinity"));
-  return re.test(s) ? {$: "Some", value: Math.fround(v)} : {$: "None"};
-}
+${MathRT.JS_POST}
 
 function char_new(code) {
   if (code > 0x10FFFF || (code >= 0xD800 && code <= 0xDFFF)) {
diff --git a/bend2/math.bend b/bend2/math.bend
new file mode 100644
index 000000000..1b6d6a878
--- /dev/null
+++ b/bend2/math.bend
@@ -0,0 +1,719 @@
+# Math: a C math.h-style F32 library for Bend 2.
+#
+# Bend 2 has one native floating-point format, F32.  This module therefore
+# exposes one spelling per math.h function instead of C's float/double/long
+# double suffix families.  Functions which return a C signed integer use S32,
+# a small sign/magnitude type, because Bend 2 has no native signed integer.
+# Consequently llround/llrint, scalbln, and the integer arguments of the C23
+# pown/rootn extras are intentionally S32 adaptations rather than pretending
+# Bend has C long/long long storage that it does not have.
+#
+# Bend also has no C floating-point environment.  nearbyint/rint use the
+# language's deterministic round-to-nearest-even operation and
+# math_errhandling() reports 0; errno/fenv exceptions are not emulated.
+#
+# The M_* definitions mirror the traditional math.h constants.  Bend has no
+# value macros, so each constant is a nullary def: M_PI(), M_SQRT2(), ...
+
+import Base
+
+# Types
+# =====
+
+type S32 is Data:
+  SNeg{mag: U32}
+  SPos{mag: U32}
+
+type FPClass is Data:
+  FP_NAN{}
+  FP_INFINITE{}
+  FP_ZERO{}
+  FP_SUBNORMAL{}
+  FP_NORMAL{}
+
+# Bit construction
+# ================
+
+def from_bits(x: U32) -> F32:
+  match x:
+    case U32{w}:
+      F32{w}
+
+# Constants
+# =========
+
+# Traditional BSD/glibc math constants, rounded to Bend's F32.
+def M_E() -> F32:
+  2.7182818284590452354
+
+def M_LOG2E() -> F32:
+  1.4426950408889634074
+
+def M_LOG10E() -> F32:
+  0.43429448190325182765
+
+def M_LN2() -> F32:
+  0.69314718055994530942
+
+def M_LN10() -> F32:
+  2.3025850929940456840
+
+def M_PI() -> F32:
+  3.14159265358979323846
+
+def M_PI_2() -> F32:
+  1.57079632679489661923
+
+def M_PI_4() -> F32:
+  0.78539816339744830962
+
+def M_1_PI() -> F32:
+  0.31830988618379067154
+
+def M_2_PI() -> F32:
+  0.63661977236758134308
+
+def M_2_SQRTPI() -> F32:
+  1.12837916709551257390
+
+def M_SQRT2() -> F32:
+  1.41421356237309504880
+
+def M_SQRT1_2() -> F32:
+  0.70710678118654752440
+
+# IEEE-754 binary32 constants.  These are made from bits so they do not depend
+# on the source parser accepting spellings such as inf or nan.
+def INFINITY() -> F32:
+  from_bits(2139095040)
+
+def HUGE_VAL() -> F32:
+  INFINITY()
+
+def HUGE_VALF() -> F32:
+  INFINITY()
+
+def NAN() -> F32:
+  from_bits(2143289344)
+
+# math_errhandling has no errno/fenv analogue in Bend; keep the C bit values
+# available while reporting that neither mechanism is provided.
+def MATH_ERRNO() -> U32:
+  1
+
+def MATH_ERREXCEPT() -> U32:
+  2
+
+def math_errhandling() -> U32:
+  0
+
+# Signed integer bridge
+# =====================
+
+def S32.zero() -> S32:
+  SPos{0}
+
+def S32.neg.if(mag: U32, zero: Bool) -> S32:
+  match zero:
+    case True{}:
+      SPos{0}
+    case False{}:
+      SNeg{mag}
+
+def S32.neg(+mag: U32) -> S32:
+  S32.neg.if(mag, U32.is_eq(mag, 0))
+
+def S32.pos(mag: U32) -> S32:
+  SPos{mag}
+
+def S32.to_f32(x: S32) -> F32:
+  match x:
+    case SNeg{mag}:
+      F32.neg(U32.to_f32(mag))
+    case SPos{mag}:
+      U32.to_f32(mag)
+
+def S32.from_f32(+x: F32) -> S32:
+  Bool.pick(S32, F32.is_lt(x, 0.0),
+    SNeg{F32.to_u32(F32.abs(x))}, SPos{F32.to_u32(x)})
+
+def S32.is_neg(x: S32) -> Bool:
+  match x:
+    case SNeg{mag}:
+      U32.is_ne(mag, 0)
+    case SPos{mag}:
+      False{}
+
+def S32.mag(x: S32) -> U32:
+  match x:
+    case SNeg{mag}:
+      mag
+    case SPos{mag}:
+      mag
+
+def S32.is_odd(x: S32) -> Bool:
+  U32.is_ne(U32.and(S32.mag(x), 1), 0)
+
+def S32.negate.pos(mag: U32, zero: Bool) -> S32:
+  match zero:
+    case True{}:
+      SPos{0}
+    case False{}:
+      SNeg{mag}
+
+def S32.negate(x: S32) -> S32:
+  match x:
+    case SNeg{mag}:
+      SPos{mag}
+    case SPos{+mag}:
+      S32.negate.pos(mag, U32.is_eq(mag, 0))
+
+def S32.show(x: S32) -> String:
+  match x:
+    case SNeg{mag}:
+      "-" ++ U32.show(mag)
+    case SPos{mag}:
+      U32.show(mag)
+
+# C uses implementation-defined integer sentinels for ilogb(0) and
+# ilogb(NaN).  These names make the Bend choice explicit.
+def FP_ILOGB0() -> S32:
+  SNeg{2147483648}
+
+def FP_ILOGBNAN() -> S32:
+  SNeg{2147483648}
+
+# Bit helpers and classification
+# ==============================
+
+def bits(x: F32) -> U32:
+  F32.bits(x)
+
+def signbit(x: F32) -> Bool:
+  U32.is_ne(U32.and(F32.bits(x), 2147483648), 0)
+
+def isnan(x: F32) -> Bool:
+  +b = F32.bits(x)
+  e = U32.and(b, 2139095040)
+  m = U32.and(b, 8388607)
+  Bool.and(U32.is_eq(e, 2139095040), U32.is_ne(m, 0))
+
+def isinf(x: F32) -> Bool:
+  U32.is_eq(U32.and(F32.bits(x), 2147483647), 2139095040)
+
+def isfinite(x: F32) -> Bool:
+  U32.is_ne(U32.and(F32.bits(x), 2139095040), 2139095040)
+
+def isnormal(x: F32) -> Bool:
+  +e = U32.and(F32.bits(x), 2139095040)
+  Bool.and(U32.is_ne(e, 0), U32.is_ne(e, 2139095040))
+
+def fpclassify(+x: F32) -> FPClass:
+  +b = F32.bits(x)
+  +e = U32.and(b, 2139095040)
+  m = U32.and(b, 8388607)
+  n = Bool.and(U32.is_eq(e, 2139095040), U32.is_ne(m, 0))
+  i = U32.is_eq(U32.and(b, 2147483647), 2139095040)
+  z = U32.is_eq(U32.and(b, 2147483647), 0)
+  sub = U32.is_eq(e, 0)
+  Bool.pick(FPClass, n, FP_NAN{},
+    Bool.pick(FPClass, i, FP_INFINITE{},
+    Bool.pick(FPClass, z, FP_ZERO{},
+    Bool.pick(FPClass, sub, FP_SUBNORMAL{}, FP_NORMAL{}))))
+
+def isgreater(x: F32, y: F32) -> Bool:
+  F32.is_gt(x, y)
+
+def isgreaterequal(x: F32, y: F32) -> Bool:
+  F32.is_ge(x, y)
+
+def isless(x: F32, y: F32) -> Bool:
+  F32.is_lt(x, y)
+
+def islessequal(x: F32, y: F32) -> Bool:
+  F32.is_le(x, y)
+
+def islessgreater(+x: F32, +y: F32) -> Bool:
+  Bool.or(F32.is_lt(x, y), F32.is_gt(x, y))
+
+def isunordered(+x: F32, +y: F32) -> Bool:
+  Bool.or(isnan(x), isnan(y))
+
+# Elementary operations
+# =====================
+
+def fabs(x: F32) -> F32:
+  F32.abs(x)
+
+def copysign(+x: F32, y: F32) -> F32:
+  xb = U32.and(F32.bits(x), 2147483647)
+  yb = U32.and(F32.bits(y), 2147483648)
+  from_bits(U32.or(xb, yb))
+
+def fmin(+x: F32, +y: F32) -> F32:
+  both_zero = Bool.and(F32.is_eq(x, 0.0), F32.is_eq(y, 0.0))
+  z = Bool.pick(F32, Bool.or(signbit(x), signbit(y)),
+    from_bits(2147483648), 0.0)
+  core = Bool.pick(F32, both_zero, z, F32.min(x, y))
+  yn = Bool.pick(F32, isnan(y), x, core)
+  Bool.pick(F32, isnan(x), y, yn)
+
+def fmax(+x: F32, +y: F32) -> F32:
+  both_zero = Bool.and(F32.is_eq(x, 0.0), F32.is_eq(y, 0.0))
+  z = Bool.pick(F32, Bool.and(signbit(x), signbit(y)),
+    from_bits(2147483648), 0.0)
+  core = Bool.pick(F32, both_zero, z, F32.max(x, y))
+  yn = Bool.pick(F32, isnan(y), x, core)
+  Bool.pick(F32, isnan(x), y, yn)
+
+def fdim(+x: F32, +y: F32) -> F32:
+  Bool.pick(F32, Bool.or(isnan(x), isnan(y)), NAN(),
+    Bool.pick(F32, F32.is_gt(x, y), F32.sub(x, y), 0.0))
+
+# Fused multiply-add through the F32 compiler primitive.
+def fma(x: F32, y: F32, z: F32) -> F32:
+  F32.fma(x, y, z)
+
+def sqrt(x: F32) -> F32:
+  F32.sqrt(x)
+
+def cbrt(x: F32) -> F32:
+  F32.cbrt(x)
+
+def hypot(x: F32, y: F32) -> F32:
+  F32.hypot(x, y)
+
+def pow(x: F32, y: F32) -> F32:
+  F32.pow(x, y)
+
+# Exponential and logarithmic functions
+# =====================================
+
+def exp(x: F32) -> F32:
+  F32.exp(x)
+
+def exp2(x: F32) -> F32:
+  F32.exp2(x)
+
+def exp10(x: F32) -> F32:
+  F32.pow(10.0, x)
+
+def expm1(x: F32) -> F32:
+  F32.expm1(x)
+
+def log(x: F32) -> F32:
+  F32.log(x)
+
+def log2(x: F32) -> F32:
+  F32.log2(x)
+
+def log10(x: F32) -> F32:
+  F32.log10(x)
+
+def log1p(x: F32) -> F32:
+  F32.log1p(x)
+
+# C23-style convenience spellings.  The p1/m1 forms use backend operations
+# so the critical product/division is evaluated wider than F32 before the
+# single final binary32 rounding.
+def exp2m1(x: F32) -> F32:
+  F32.exp2m1(x)
+
+def exp10m1(x: F32) -> F32:
+  F32.exp10m1(x)
+
+def logp1(x: F32) -> F32:
+  F32.log1p(x)
+
+def log2p1(x: F32) -> F32:
+  F32.log2p1(x)
+
+def log10p1(x: F32) -> F32:
+  F32.log10p1(x)
+
+# Trigonometric functions
+# =======================
+
+def sin(x: F32) -> F32:
+  F32.sin(x)
+
+def cos(x: F32) -> F32:
+  F32.cos(x)
+
+def tan(x: F32) -> F32:
+  F32.tan(x)
+
+def asin(x: F32) -> F32:
+  F32.asin(x)
+
+def acos(x: F32) -> F32:
+  F32.acos(x)
+
+def atan(x: F32) -> F32:
+  F32.atan(x)
+
+def atan2(y: F32, x: F32) -> F32:
+  F32.atan2(y, x)
+
+def sincos(+x: F32) -> F32 & F32:
+  (F32.sin(x), F32.cos(x))
+
+# pi-scaled forms use backend primitives with exact integer/half-integer
+# reduction instead of multiplying a large F32 by an approximation of pi.
+def sinpi(x: F32) -> F32:
+  F32.sinpi(x)
+
+def cospi(x: F32) -> F32:
+  F32.cospi(x)
+
+def tanpi(x: F32) -> F32:
+  F32.tanpi(x)
+
+def asinpi(x: F32) -> F32:
+  F32.asinpi(x)
+
+def acospi(x: F32) -> F32:
+  F32.acospi(x)
+
+def atanpi(x: F32) -> F32:
+  F32.atanpi(x)
+
+def atan2pi(y: F32, x: F32) -> F32:
+  F32.atan2pi(y, x)
+
+# Hyperbolic functions
+# ====================
+
+def sinh(x: F32) -> F32:
+  F32.sinh(x)
+
+def cosh(x: F32) -> F32:
+  F32.cosh(x)
+
+def tanh(x: F32) -> F32:
+  F32.tanh(x)
+
+def asinh(x: F32) -> F32:
+  F32.asinh(x)
+
+def acosh(x: F32) -> F32:
+  F32.acosh(x)
+
+def atanh(x: F32) -> F32:
+  F32.atanh(x)
+
+# Decomposition and scaling
+# =========================
+
+def scalbn(x: F32, n: S32) -> F32:
+  match n:
+    case SNeg{mag}:
+      F32.scalbn_neg(x, mag)
+    case SPos{mag}:
+      F32.scalbn_pos(x, mag)
+
+def scalbln(x: F32, n: S32) -> F32:
+  scalbn(x, n)
+
+def ldexp(x: F32, n: S32) -> F32:
+  scalbn(x, n)
+
+def frexp.regular(+x: F32) -> F32 & S32:
+  +e = S32.from_f32(F32.add(F32.logb(F32.abs(x)), 1.0))
+  (scalbn(x, S32.negate(e)), e)
+
+def frexp.if(x: F32, special: Bool) -> F32 & S32:
+  match special:
+    case True{}:
+      (x, S32.zero())
+    case False{}:
+      frexp.regular(x)
+
+def frexp(+x: F32) -> F32 & S32:
+  frexp.if(x, Bool.or(F32.is_eq(x, 0.0), Bool.or(isnan(x), isinf(x))))
+
+def modf.finite(+x: F32) -> F32 & F32:
+  +i = F32.trunc(x)
+  +f = F32.sub(x, i)
+  frac = Bool.pick(F32, F32.is_eq(f, 0.0), copysign(0.0, x), f)
+  (frac, i)
+
+def modf.inf(+x: F32, inf: Bool) -> F32 & F32:
+  match inf:
+    case True{}:
+      (copysign(0.0, x), x)
+    case False{}:
+      modf.finite(x)
+
+def modf.nan(+x: F32, nanv: Bool) -> F32 & F32:
+  match nanv:
+    case True{}:
+      (x, x)
+    case False{}:
+      modf.inf(x, isinf(x))
+
+def modf(+x: F32) -> F32 & F32:
+  modf.nan(x, isnan(x))
+
+def logb(x: F32) -> F32:
+  F32.logb(x)
+
+def ilogb.zero(x: F32, zero: Bool) -> S32:
+  match zero:
+    case True{}:
+      FP_ILOGB0()
+    case False{}:
+      S32.from_f32(F32.logb(x))
+
+def ilogb.inf(+x: F32, inf: Bool) -> S32:
+  match inf:
+    case True{}:
+      SPos{2147483647}
+    case False{}:
+      ilogb.zero(x, F32.is_eq(x, 0.0))
+
+def ilogb.nan(+x: F32, nanv: Bool) -> S32:
+  match nanv:
+    case True{}:
+      FP_ILOGBNAN()
+    case False{}:
+      ilogb.inf(x, isinf(x))
+
+def ilogb(+x: F32) -> S32:
+  ilogb.nan(x, isnan(x))
+
+def ceil(x: F32) -> F32:
+  F32.ceil(x)
+
+def floor(x: F32) -> F32:
+  F32.floor(x)
+
+def trunc(x: F32) -> F32:
+  F32.trunc(x)
+
+def round.small(+x: F32, neg: Bool) -> F32:
+  match neg:
+    case True{}:
+      F32.ceil(F32.sub(x, 0.5))
+    case False{}:
+      F32.floor(F32.add(x, 0.5))
+
+def round.large(+x: F32, large: Bool) -> F32:
+  match large:
+    case True{}:
+      x
+    case False{}:
+      round.small(x, signbit(x))
+
+def round(+x: F32) -> F32:
+  # Every binary32 value with |x| >= 2^23 is already integral.  Returning it
+  # directly avoids the x +/- 0.5 rounding bug at large odd integers.
+  round.large(x, F32.is_ge(F32.abs(x), 8388608.0))
+
+def roundeven(x: F32) -> F32:
+  F32.roundeven(x)
+
+def nearbyint(x: F32) -> F32:
+  roundeven(x)
+
+def rint(x: F32) -> F32:
+  roundeven(x)
+
+def lround(x: F32) -> S32:
+  S32.from_f32(round(x))
+
+def llround(x: F32) -> S32:
+  S32.from_f32(round(x))
+
+def lrint(x: F32) -> S32:
+  S32.from_f32(rint(x))
+
+def llrint(x: F32) -> S32:
+  S32.from_f32(rint(x))
+
+# Remainders
+# ==========
+
+def fmod(x: F32, y: F32) -> F32:
+  F32.mod(x, y)
+
+def remainder(x: F32, y: F32) -> F32:
+  F32.remainder(x, y)
+
+# C only guarantees the sign and at least the low three quotient bits.  We
+# return exactly those three bits in S32, avoiding overflow for huge x/y.
+def remquo.code(code: U32, neg: Bool) -> S32:
+  match neg:
+    case True{}:
+      S32.neg(U32.and(code, 7))
+    case False{}:
+      SPos{U32.and(code, 7)}
+
+def remquo(+x: F32, +y: F32) -> F32 & S32:
+  +r = F32.remainder(x, y)
+  +q = F32.remquo_code(x, y)
+  (r, remquo.code(q, U32.is_ne(U32.and(q, 2147483648), 0)))
+
+# nextafter / nexttoward
+# =====================
+
+def nextafter.zero(y: F32) -> F32:
+  s = U32.and(F32.bits(y), 2147483648)
+  from_bits(U32.or(s, 1))
+
+def nextafter.step(+x: F32, +y: F32) -> F32:
+  +b = F32.bits(x)
+  sx = signbit(x)
+  up = F32.is_lt(x, y)
+  down_bits = from_bits(U32.sub(b, 1))
+  up_bits = from_bits(U32.add(b, 1))
+  Bool.pick(F32, Bool.not(Bool.xor(sx, up)), down_bits, up_bits)
+
+def nextafter(+x: F32, +y: F32) -> F32:
+  stepped = nextafter.step(x, y)
+  zeroed = Bool.pick(F32, F32.is_eq(x, 0.0), nextafter.zero(y), stepped)
+  equaled = Bool.pick(F32, F32.is_eq(x, y), y, zeroed)
+  Bool.pick(F32, Bool.or(isnan(x), isnan(y)), F32.add(x, y), equaled)
+
+def nexttoward(x: F32, y: F32) -> F32:
+  nextafter(x, y)
+
+# Error and gamma functions
+# =========================
+
+def erf(x: F32) -> F32:
+  F32.erf(x)
+
+def erfc(x: F32) -> F32:
+  F32.erfc(x)
+
+def tgamma(x: F32) -> F32:
+  F32.tgamma(x)
+
+def lgamma(x: F32) -> F32:
+  F32.lgamma(x)
+
+# C's nan(tagp); Bend does not expose NaN payload parsing, so tag is accepted
+# for API familiarity and a canonical quiet NaN is returned.
+def nan(tag: String) -> F32:
+  NAN()
+
+# Integer-power helpers (C23 family)
+# ==================================
+
+def pown(x: F32, n: S32) -> F32:
+  match n:
+    case SNeg{mag}:
+      F32.pown_neg(x, mag)
+    case SPos{mag}:
+      F32.pown_pos(x, mag)
+
+def powr.zero_y(y: F32, zero: Bool) -> F32:
+  match zero:
+    case True{}:
+      NAN()
+    case False{}:
+      Bool.pick(F32, F32.is_lt(y, 0.0), INFINITY(), 0.0)
+
+def powr.xzero(+x: F32, +y: F32, zero: Bool) -> F32:
+  match zero:
+    case True{}:
+      powr.zero_y(y, F32.is_eq(y, 0.0))
+    case False{}:
+      F32.pow(x, y)
+
+def powr.one_inf(+x: F32, +y: F32, bad: Bool) -> F32:
+  match bad:
+    case True{}:
+      NAN()
+    case False{}:
+      powr.xzero(x, y, F32.is_eq(x, 0.0))
+
+def powr.inf_zero(+x: F32, +y: F32, bad: Bool) -> F32:
+  match bad:
+    case True{}:
+      NAN()
+    case False{}:
+      powr.one_inf(x, y, Bool.and(F32.is_eq(x, 1.0), isinf(y)))
+
+def powr.neg(+x: F32, +y: F32, neg: Bool) -> F32:
+  match neg:
+    case True{}:
+      NAN()
+    case False{}:
+      powr.inf_zero(x, y, Bool.and(isinf(x), F32.is_eq(y, 0.0)))
+
+def powr.nan(+x: F32, +y: F32, nanv: Bool) -> F32:
+  match nanv:
+    case True{}:
+      NAN()
+    case False{}:
+      powr.neg(x, y, F32.is_lt(x, 0.0))
+
+def powr.eval(+x: F32, +y: F32) -> F32:
+  # Unlike pow, powr has deliberately stricter zero/one/infinity semantics.
+  powr.nan(x, y, Bool.or(isnan(x), isnan(y)))
+
+# C23 spells this API as powr(y, x) and computes x^y.  Keep that unusual
+# argument order so code ported from <math.h> does not silently change meaning.
+def powr(+y: F32, +x: F32) -> F32:
+  powr.eval(x, y)
+
+def rootn.neg(x: F32, nf: F32, odd: Bool) -> F32:
+  match odd:
+    case True{}:
+      F32.neg(F32.pow(F32.abs(x), F32.div(1.0, nf)))
+    case False{}:
+      NAN()
+
+def rootn.sign(+x: F32, +n: S32, nf: F32, neg: Bool) -> F32:
+  match neg:
+    case True{}:
+      rootn.neg(x, nf, S32.is_odd(n))
+    case False{}:
+      F32.pow(x, F32.div(1.0, nf))
+
+def rootn.nonzero(+x: F32, +n: S32) -> F32:
+  nf = S32.to_f32(n)
+  rootn.sign(x, n, nf, signbit(x))
+
+def rootn.zero.neg(x: F32, odd: Bool) -> F32:
+  match odd:
+    case True{}:
+      copysign(INFINITY(), x)
+    case False{}:
+      INFINITY()
+
+def rootn.zero.pos(x: F32, odd: Bool) -> F32:
+  match odd:
+    case True{}:
+      x
+    case False{}:
+      0.0
+
+def rootn.zero.sign(+x: F32, +n: S32, negn: Bool) -> F32:
+  match negn:
+    case True{}:
+      rootn.zero.neg(x, S32.is_odd(n))
+    case False{}:
+      rootn.zero.pos(x, S32.is_odd(n))
+
+def rootn.xzero(+x: F32, +n: S32, zero: Bool) -> F32:
+  match zero:
+    case True{}:
+      rootn.zero.sign(x, n, S32.is_neg(n))
+    case False{}:
+      rootn.nonzero(x, n)
+
+def rootn.nzero(+x: F32, +n: S32, zero: Bool) -> F32:
+  match zero:
+    case True{}:
+      NAN()
+    case False{}:
+      rootn.xzero(x, n, F32.is_eq(x, 0.0))
+
+def rootn(+x: F32, +n: S32) -> F32:
+  rootn.nzero(x, n, U32.is_eq(S32.mag(n), 0))
+
+def rsqrt(x: F32) -> F32:
+  F32.rsqrt(x)
diff --git a/bend2/math.ts b/bend2/math.ts
new file mode 100644
index 000000000..cd913e97f
--- /dev/null
+++ b/bend2/math.ts
@@ -0,0 +1,806 @@
+// Runtime support for bend2/math.bend. Kept separate from comp.ts so the
+// compiler core remains within its repository size gate.
+
+const CMPS = "is_eq:==:=== is_ne:!=:!== is_lt:< is_le:<= is_gt:> is_ge:>=";
+
+function tpl_ops(pre: string, names: string, C: string, JS: string) {
+  const out: Record<string, { C: string; JS: string; call?: boolean }> = {};
+  for (const p of names.split(" ")) {
+    const [k, o = k, jo = o] = p.split(":");
+    out[pre + k] = { C: C.replaceAll("$o", o), JS: JS.replaceAll("$o", jo) };
+  }
+  return out;
+}
+
+export const OPERATIONS = {
+  ...tpl_ops("f32_", "add:+ sub:- mul:* div:/",
+    "f32_rewrap(f32_unbox($0) $o f32_unbox($1))", "Math.fround($0 $o $1)"),
+  f32_neg: {
+    C:  "f32_rewrap(-f32_unbox($0))",
+    JS: "(-$0)",
+  },
+  ...tpl_ops("f32_", CMPS, "((u64)(f32_unbox($0) $o f32_unbox($1)))",
+    "($0 $o $1)"),
+  ...tpl_ops("f32_", "sqrt exp log log2 log10 sin cos tan asin acos atan"
+    + " sinh cosh tanh floor ceil trunc abs:fabs:abs",
+    "f32_rewrap((f32)$o(f32_unbox($0)))", "Math.fround(Math.$o($0))"),
+  f32_pow: {
+    C:  "f32_rewrap((f32)pow(f32_unbox($0), f32_unbox($1)))",
+    JS: "f32_pow($0, $1)",
+  },
+  f32_atan2: {
+    C:  "f32_rewrap((f32)atan2(f32_unbox($0), f32_unbox($1)))",
+    JS: "Math.fround(Math.atan2($0, $1))",
+  },
+  f32_fma: {
+    C:  "f32_rewrap(fmaf(f32_unbox($0), f32_unbox($1), f32_unbox($2)))",
+    JS: "f32_fma($0, $1, $2)",
+  },
+  f32_cbrt: {
+    C:  "f32_rewrap((f32)cbrt(f32_unbox($0)))",
+    JS: "Math.fround(Math.cbrt($0))",
+  },
+  f32_hypot_native: {
+    C:  "f32_rewrap((f32)hypot(f32_unbox($0), f32_unbox($1)))",
+    JS: "Math.fround(Math.hypot($0, $1))",
+  },
+  f32_exp2: {
+    C:  "f32_rewrap((f32)exp2(f32_unbox($0)))",
+    JS: "Math.fround(2 ** $0)",
+  },
+  f32_expm1: {
+    C:  "f32_rewrap((f32)expm1(f32_unbox($0)))",
+    JS: "Math.fround(Math.expm1($0))",
+  },
+  f32_exp2m1: {
+    C:  "f32_rewrap(bend_exp2m1f(f32_unbox($0)))",
+    JS: "Math.fround(Math.expm1($0 * Math.LN2))",
+  },
+  f32_exp10m1: {
+    C:  "f32_rewrap(bend_exp10m1f(f32_unbox($0)))",
+    JS: "Math.fround(Math.expm1($0 * 2.30258509299404568402))",
+  },
+  f32_log1p: {
+    C:  "f32_rewrap((f32)log1p(f32_unbox($0)))",
+    JS: "Math.fround(Math.log1p($0))",
+  },
+  f32_log2p1: {
+    C:  "f32_rewrap(bend_log2p1f(f32_unbox($0)))",
+    JS: "Math.fround(Math.log1p($0) * Math.LOG2E)",
+  },
+  f32_log10p1: {
+    C:  "f32_rewrap(bend_log10p1f(f32_unbox($0)))",
+    JS: "Math.fround(Math.log1p($0) * Math.LOG10E)",
+  },
+  f32_asinh: {
+    C:  "f32_rewrap((f32)asinh(f32_unbox($0)))",
+    JS: "Math.fround(Math.asinh($0))",
+  },
+  f32_acosh: {
+    C:  "f32_rewrap((f32)acosh(f32_unbox($0)))",
+    JS: "Math.fround(Math.acosh($0))",
+  },
+  f32_atanh: {
+    C:  "f32_rewrap((f32)atanh(f32_unbox($0)))",
+    JS: "Math.fround(Math.atanh($0))",
+  },
+  f32_scalbn_pos: {
+    C:  "f32_rewrap(scalbnf(f32_unbox($0), (int)((u32)$1 > 512 ? 512 : (u32)$1)))",
+    JS: "f32_scalbn($0, $1, false)",
+  },
+  f32_scalbn_neg: {
+    C:  "f32_rewrap(scalbnf(f32_unbox($0), -(int)((u32)$1 > 512 ? 512 : (u32)$1)))",
+    JS: "f32_scalbn($0, $1, true)",
+  },
+  f32_logb: {
+    C:  "f32_rewrap(logbf(f32_unbox($0)))",
+    JS: "f32_logb($0)",
+  },
+  f32_roundeven: {
+    C:  "f32_rewrap(bend_roundevenf(f32_unbox($0)))",
+    JS: "f32_roundeven($0)",
+  },
+  f32_remainder: {
+    C:  "f32_rewrap(remainderf(f32_unbox($0), f32_unbox($1)))",
+    JS: "f32_remainder($0, $1)",
+  },
+  f32_remquo_code: {
+    C:  "bend_remquo_code(f32_unbox($0), f32_unbox($1))",
+    JS: "f32_remquo_code($0, $1)",
+  },
+  f32_erf: {
+    C:  "f32_rewrap(erff(f32_unbox($0)))",
+    JS: "f32_erf($0)",
+  },
+  f32_erfc: {
+    C:  "f32_rewrap(erfcf(f32_unbox($0)))",
+    JS: "f32_erfc($0)",
+  },
+  f32_tgamma: {
+    C:  "f32_rewrap(tgammaf(f32_unbox($0)))",
+    JS: "f32_tgamma($0)",
+  },
+  f32_lgamma: {
+    C:  "f32_rewrap(lgammaf(f32_unbox($0)))",
+    JS: "f32_lgamma($0)",
+  },
+  f32_sinpi: {
+    C:  "f32_rewrap(bend_sinpif(f32_unbox($0)))",
+    JS: "f32_sinpi($0)",
+  },
+  f32_cospi: {
+    C:  "f32_rewrap(bend_cospif(f32_unbox($0)))",
+    JS: "f32_cospi($0)",
+  },
+  f32_tanpi: {
+    C:  "f32_rewrap(bend_tanpif(f32_unbox($0)))",
+    JS: "f32_tanpi($0)",
+  },
+  f32_asinpi: {
+    C:  "f32_rewrap(bend_asinpif(f32_unbox($0)))",
+    JS: "Math.fround(Math.asin($0) / Math.PI)",
+  },
+  f32_acospi: {
+    C:  "f32_rewrap(bend_acospif(f32_unbox($0)))",
+    JS: "Math.fround(Math.acos($0) / Math.PI)",
+  },
+  f32_atanpi: {
+    C:  "f32_rewrap(bend_atanpif(f32_unbox($0)))",
+    JS: "Math.fround(Math.atan($0) / Math.PI)",
+  },
+  f32_atan2pi: {
+    C:  "f32_rewrap(bend_atan2pif(f32_unbox($0), f32_unbox($1)))",
+    JS: "Math.fround(Math.atan2($0, $1) / Math.PI)",
+  },
+  f32_rsqrt: {
+    C:  "f32_rewrap(bend_rsqrtf(f32_unbox($0)))",
+    JS: "Math.fround(1 / Math.sqrt($0))",
+  },
+  f32_pown_pos: {
+    C:  "f32_rewrap(bend_pownf(f32_unbox($0), (u32)$1, 0))",
+    JS: "Math.fround(Math.pow($0, ($1 >>> 0)))",
+  },
+  f32_pown_neg: {
+    C:  "f32_rewrap(bend_pownf(f32_unbox($0), (u32)$1, 1))",
+    JS: "Math.fround(Math.pow($0, -($1 >>> 0)))",
+  },
+  f32_mod: {
+    C:  "f32_rewrap((f32)fmod(f32_unbox($0), f32_unbox($1)))",
+    JS: "Math.fround($0 % $1)",
+  },
+  f32_to_u32: {
+    C:  "f32_to_u32($0)",
+    JS: "($0 >= 1 && $0 < 4294967296 ? Math.floor($0) : 0)",
+  },
+  f32_bits: {
+    C:  "$0",
+    JS: "f32_bits($0)",
+  },
+  f32_show: {
+    C:    "f32_show(e, $0)",
+    call: true,
+    JS:   "f32_show($0)",
+  },
+  f32_read: {
+    C:    "f32_read(e, $0)",
+    call: true,
+    JS:   "f32_read($0)",
+  },
+};
+
+export const SHIMS = (() => {
+  const names = ("sqrt exp log log2 log10 sin cos tan asin acos atan sinh cosh"
+    + " tanh floor ceil trunc fabs pow fmod fma cbrt hypot exp2 expm1 log1p"
+    + " asinh acosh atanh logb rint remainder erf erfc tgamma lgamma").split(" ");
+  const target = (n: string) =>
+    (["sin", "cos", "tan"].includes(n) ? " fast::" : " precise::") + n;
+  return names.map((n) => "#define " + n.padEnd(11) + target(n))
+    .concat(names.map((n) => "#define " + (n + "f").padEnd(11) + target(n)))
+    .join("\n") + "\n#define scalbnf    precise::ldexp"
+    + "\n#define atan2      atan2_c99"
+    + "\n#define atan2f     atan2_c99";
+})();
+
+export const C_PRE = String.raw`
+INLINE f32 f32_unbox(u64 x) {
+  union { u32 u; f32 f; } p = { (u32)x };
+  return p.f;
+}
+
+INLINE u64 f32_rewrap(f32 x) {
+  union { f32 f; u32 u; } p = { x };
+  return p.u;
+}
+
+`.slice(1);
+
+export const C_POST = String.raw`
+INLINE U32 f32_to_u32(U32 a) {
+  f32 v = f32_unbox(a);
+  return v >= 0.0f && v < 4294967296.0f ? (u32)v : 0;
+}
+
+`.slice(1);
+
+export const C_F32_IO_DECL = String.raw`
+#if DEVICE
+
+#define f32_show(e, x) (err_post(e.mem, ERR_FIDS), 0)
+#define f32_read(e, s) (err_post(e.mem, ERR_FIDS), 0)
+
+#else
+
+static Term f32_show(Env e, Term x);
+static Term f32_read(Env e, Term s);
+
+#endif
+`.slice(1);
+
+export const C_IO = String.raw`
+static int f32_text(char* buf, f32 v) {
+  int n = 0;
+  int p = 0;
+  if (v != v) {
+    return sprintf(buf, "nan");
+  }
+  for (; p < 9; p += 1) {
+    n = snprintf(buf, 40, "%.*e", p, (double)v);
+    if (strtof(buf, NULL) == v) {
+      break;
+    }
+  }
+  char* ep = strchr(buf, 'e');
+  if (ep == NULL) {
+    return n;
+  }
+  int ex = atoi(ep + 1);
+  if (ex >= 21 || ex <= -7) {
+    n = (int)(ep - buf) + sprintf(ep, "e%c%d", ex < 0 ? '-' : '+', abs(ex));
+  } else if (ex <= p) {
+    n = snprintf(buf, 40, "%.*f", p - ex, (double)v);
+  } else {
+    int s = *buf == '-';
+    memmove(buf + s + 1, buf + s + 2, p);
+    memset(buf + s + 1 + p, '0', ex - p);
+    n = s + 1 + ex;
+  }
+  return n;
+}
+
+static Term f32_show(Env e, Term x) {
+  char buf[40];
+  return io_str(e, buf, f32_text(buf, f32_unbox(x)));
+}
+
+static Term f32_read(Env e, Term s) {
+  u64 n = 0;
+  char* text = io_cstr(e, s, &n);
+  char* end;
+  f32 v = strtof(text, &end);
+  Term out = n > 0 && (u64)(end - text) == n && strpbrk(text, "xX(") == NULL
+    ? io_box(e, CID_SOME, f32_rewrap(v)) : term_pak(CID_NONE, 0);
+  free(text);
+  return out;
+}
+`.slice(1);
+
+export const JS_PRE = String.raw`
+function f32_show(x) {
+  if (x !== x) {
+    return "nan";
+  }
+  if (!Number.isFinite(x) || Object.is(x, -0)) {
+    return x < 0 ? "-inf"
+      : x === 0 ? "-0" : "inf";
+  }
+  let s = "x";
+  for (let p = 1; p <= 9 && Math.fround(Number(s)) !== x; p += 1) {
+    s = String(Number(x.toExponential(p - 1)));
+  }
+  return s;
+}
+
+function f32_bits(x) {
+  return new Uint32Array(new Float32Array([x]).buffer)[0];
+}
+
+function f32_from_bits(u) {
+  return new Float32Array(new Uint32Array([u]).buffer)[0];
+}
+
+`.slice(1);
+
+export const JS_POST = String.raw`
+function f32_read(s) {
+  const re = /^\s*[+-]?((\d+\.?\d*|\.\d+)(e[+-]?\d+)?|inf(inity)?|nan)$/i;
+  const v = Number(s.replace(/inf\w*/i, "Infinity"));
+  return re.test(s) ? {$: "Some", value: Math.fround(v)} : {$: "None"};
+}
+
+`.slice(1);
+
+export const C = String.raw`
+INLINE f32 bend_pownf(f32 x, u32 n, int neg_exp) {
+#ifdef __METAL_VERSION__
+  if (n == 0) return 1.0f;
+  const int neg_result = signbit(x) && (n & 1u);
+  const f32 a = fabs(x);
+  const f32 ehi = (f32)(n & 0xffff0000u);
+  const f32 elo = (f32)(n & 0x0000ffffu);
+  const f32 hi = precise::pow(a, neg_exp ? -ehi : ehi);
+  const f32 lo = precise::pow(a, neg_exp ? -elo : elo);
+  const f32 mag = hi * lo;
+  return neg_result ? -mag : mag;
+#else
+  const double e = neg_exp ? -(double)n : (double)n;
+  return (f32)pow((double)x, e);
+#endif
+}
+
+INLINE f32 bend_asinpif(f32 x) {
+#ifdef __METAL_VERSION__
+  return precise::asin(x) * 0.31830988618379067154f;
+#else
+  return (f32)(asin((double)x) * 0x1.45f306dc9c883p-2);
+#endif
+}
+
+INLINE f32 bend_acospif(f32 x) {
+#ifdef __METAL_VERSION__
+  return precise::acos(x) * 0.31830988618379067154f;
+#else
+  return (f32)(acos((double)x) * 0x1.45f306dc9c883p-2);
+#endif
+}
+
+INLINE f32 bend_atanpif(f32 x) {
+#ifdef __METAL_VERSION__
+  return precise::atan(x) * 0.31830988618379067154f;
+#else
+  return (f32)(atan((double)x) * 0x1.45f306dc9c883p-2);
+#endif
+}
+
+INLINE f32 bend_atan2pif(f32 y, f32 x) {
+#ifdef __METAL_VERSION__
+  return atan2_c99(y, x) * 0.31830988618379067154f;
+#else
+  return (f32)(atan2((double)y, (double)x) * 0x1.45f306dc9c883p-2);
+#endif
+}
+
+INLINE f32 bend_rsqrtf(f32 x) {
+#ifdef __METAL_VERSION__
+  return 1.0f / precise::sqrt(x);
+#else
+  return (f32)(1.0 / sqrt((double)x));
+#endif
+}
+
+INLINE f32 bend_exp2m1f(f32 x) {
+#ifdef __METAL_VERSION__
+  return precise::expm1(x * 0.69314718055994530942f);
+#else
+  return (f32)expm1((double)x * 0x1.62e42fefa39efp-1);
+#endif
+}
+
+INLINE f32 bend_exp10m1f(f32 x) {
+#ifdef __METAL_VERSION__
+  return precise::expm1(x * 2.30258509299404568402f);
+#else
+  return (f32)expm1((double)x * 0x1.26bb1bbb55516p+1);
+#endif
+}
+
+INLINE f32 bend_log2p1f(f32 x) {
+#ifdef __METAL_VERSION__
+  return precise::log1p(x) * 1.44269504088896340736f;
+#else
+  return (f32)(log1p((double)x) * 0x1.71547652b82fep+0);
+#endif
+}
+
+INLINE f32 bend_log10p1f(f32 x) {
+#ifdef __METAL_VERSION__
+  return precise::log1p(x) * 0.43429448190325182765f;
+#else
+  return (f32)(log1p((double)x) * 0x1.bcb7b1526e50ep-2);
+#endif
+}
+
+INLINE f32 bend_sinpif(f32 x) {
+#ifdef __METAL_VERSION__
+  return precise::sinpi(x);
+#else
+  if (!isfinite(x)) return NAN;
+  if (x == 0.0f) return x;
+  f32 a = fabsf(x);
+  if (a >= 0x1p24f) return copysignf(0.0f, x);
+  f32 r = fmodf(a, 2.0f);
+  double y;
+  if (r == 0.0f || r == 1.0f) y = 0.0;
+  else if (r == 0.5f) y = 1.0;
+  else if (r == 1.5f) y = -1.0;
+  else if (r < 0.5f) y = sin(0x1.921fb54442d18p+1 * (double)r);
+  else if (r < 1.0f) y = sin(0x1.921fb54442d18p+1 * (double)(1.0f - r));
+  else if (r < 1.5f) y = -sin(0x1.921fb54442d18p+1 * (double)(r - 1.0f));
+  else y = -sin(0x1.921fb54442d18p+1 * (double)(2.0f - r));
+  return (f32)(signbit(x) ? -y : y);
+#endif
+}
+
+INLINE f32 bend_cospif(f32 x) {
+#ifdef __METAL_VERSION__
+  return precise::cospi(x);
+#else
+  if (!isfinite(x)) return NAN;
+  f32 a = fabsf(x);
+  if (a >= 0x1p24f) return 1.0f;
+  f32 r = fmodf(a, 2.0f);
+  if (r == 0.0f) return 1.0f;
+  if (r == 0.5f || r == 1.5f) return 0.0f;
+  if (r == 1.0f) return -1.0f;
+  if (r < 0.5f) return (f32)cos(0x1.921fb54442d18p+1 * (double)r);
+  if (r < 1.0f) return (f32)-cos(0x1.921fb54442d18p+1 * (double)(1.0f - r));
+  if (r < 1.5f) return (f32)-cos(0x1.921fb54442d18p+1 * (double)(r - 1.0f));
+  return (f32)cos(0x1.921fb54442d18p+1 * (double)(2.0f - r));
+#endif
+}
+
+INLINE f32 bend_tanpif(f32 x) {
+#ifdef __METAL_VERSION__
+  return precise::tanpi(x);
+#else
+  if (!isfinite(x)) return NAN;
+  if (x == 0.0f) return x;
+  const int neg = signbit(x) != 0;
+  f32 a = fabsf(x);
+  if (a >= 0x1p24f) return copysignf(0.0f, x); // all such F32 integers are even
+  f32 p = fmodf(a, 2.0f);
+  f32 r = fmodf(a, 1.0f);
+  if (r == 0.0f) {
+    const int odd = p == 1.0f;
+    return (odd != neg) ? -0.0f : 0.0f;
+  }
+  if (r == 0.5f) {
+    const int floor_odd = p == 1.5f;
+    return (floor_odd != neg) ? -INFINITY : INFINITY;
+  }
+  double y = r < 0.5f
+    ? tan(0x1.921fb54442d18p+1 * (double)r)
+    : -tan(0x1.921fb54442d18p+1 * (double)(1.0f - r));
+  return (f32)(neg ? -y : y);
+#endif
+}
+
+INLINE f32 bend_roundevenf(f32 x) {
+#ifdef __METAL_VERSION__
+  return precise::rint(x);
+#else
+  if (!isfinite(x) || fabsf(x) >= 0x1p23f) return x;
+  f32 f = floorf(x);
+  f32 d = x - f;
+  f32 r;
+  if (d < 0.5f) r = f;
+  else if (d > 0.5f) r = f + 1.0f;
+  else r = fmodf(fabsf(f), 2.0f) == 0.0f ? f : f + 1.0f;
+  return r == 0.0f ? copysignf(0.0f, x) : r;
+#endif
+}
+
+INLINE U32 bend_remquo_code(f32 x, f32 y) {
+  if (!isfinite(x) || isnan(y) || y == 0.0f) return 0;
+  if (!isfinite(y)) return 0;
+  int q = 0;
+  f32 r;
+#ifdef __METAL_VERSION__
+  r = precise::remquo(x, y, q);
+#else
+  r = remquof(x, y, &q);
+#endif
+  u32 uq = (u32)q;
+  u32 mag = q < 0 ? 0u - uq : uq;
+  bool nonzero = q != 0 || r != x;
+  bool neg = signbit(x) != signbit(y);
+  return (mag & 7u) | (neg && nonzero ? 0x80000000u : 0u);
+}
+
+`.slice(1);
+
+export const JS = String.raw`
+function f32_decomp(x) {
+  const u = f32_bits(x);
+  const s = (u >>> 31) === 0 ? 1n : -1n;
+  const e = (u >>> 23) & 255;
+  const f = u & 0x7FFFFF;
+  if (e === 255) return null;
+  return e === 0
+    ? {s, m: BigInt(f), e: -149}
+    : {s, m: BigInt(0x800000 | f), e: e - 150};
+}
+
+function f32_round_shr(a, n) {
+  if (n <= 0) return a << BigInt(-n);
+  const sh = BigInt(n);
+  const q = a >> sh;
+  const r = a - (q << sh);
+  const h = 1n << (sh - 1n);
+  return r > h || (r === h && (q & 1n) !== 0n) ? q + 1n : q;
+}
+
+function f32_from_exact(n, e, zeroNeg = false) {
+  if (n === 0n) return zeroNeg ? -0 : 0;
+  const neg = n < 0n;
+  let a = neg ? -n : n;
+  let len = a.toString(2).length;
+  let top = e + len - 1;
+  if (top > 127) return neg ? -Infinity : Infinity;
+  let bits;
+  if (top >= -126) {
+    let q = f32_round_shr(a, len - 24);
+    if (q === 0x1000000n) {
+      q >>= 1n;
+      top += 1;
+      if (top > 127) return neg ? -Infinity : Infinity;
+    }
+    const exp = top + 127;
+    bits = (exp << 23) | Number(q - 0x800000n);
+  } else {
+    const q = f32_round_shr(a, -(e + 149));
+    if (q === 0n) return neg ? -0 : 0;
+    bits = q >= 0x800000n ? 0x00800000 : Number(q);
+  }
+  if (neg) bits |= 0x80000000;
+  return f32_from_bits(bits >>> 0);
+}
+
+function f32_fma(a, b, c) {
+  if (Number.isNaN(a) || Number.isNaN(b) || Number.isNaN(c)) return NaN;
+  const az = a === 0;
+  const bz = b === 0;
+  if ((!Number.isFinite(a) && bz) || (!Number.isFinite(b) && az)) return NaN;
+  if (!Number.isFinite(a) || !Number.isFinite(b)) {
+    const pn = (a < 0 || Object.is(a, -0)) !== (b < 0 || Object.is(b, -0));
+    if (!Number.isFinite(c) && ((c < 0) !== pn)) return NaN;
+    return pn ? -Infinity : Infinity;
+  }
+  if (!Number.isFinite(c)) return c;
+  const pa = f32_decomp(a), pb = f32_decomp(b), pc = f32_decomp(c);
+  let np = pa.s * pb.s * pa.m * pb.m;
+  const ep = pa.e + pb.e;
+  let nc = pc.s * pc.m;
+  if (np === 0n && nc === 0n) {
+    const an = a < 0 || Object.is(a, -0);
+    const bn = b < 0 || Object.is(b, -0);
+    const pn = an !== bn;
+    const cn = c < 0 || Object.is(c, -0);
+    return pn === cn && pn ? -0 : 0;
+  }
+  if (np === 0n) return c;
+  const e = Math.min(ep, pc.e);
+  np <<= BigInt(ep - e);
+  nc <<= BigInt(pc.e - e);
+  return f32_from_exact(np + nc, e);
+}
+
+function f32_pow(x, y) {
+  // ECMAScript intentionally differs from IEC 60559/C pow at a few exact
+  // special values.  Normalize those cases before using the host Math.pow.
+  if (y === 0 || x === 1) return 1;
+  if (x === -1 && (y === Infinity || y === -Infinity)) return 1;
+  return Math.fround(Math.pow(x, y));
+}
+
+function f32_scalbn(x, mag, neg) {
+  const n = Math.min(mag >>> 0, 512);
+  return Math.fround(x * 2 ** (neg ? -n : n));
+}
+
+function f32_logb(x) {
+  const u = f32_bits(x) & 0x7FFFFFFF;
+  const e = u >>> 23;
+  const f = u & 0x7FFFFF;
+  if (e === 255) return f === 0 ? Infinity : NaN;
+  if (e !== 0) return Math.fround(e - 127);
+  if (f === 0) return -Infinity;
+  return Math.fround(-149 + (31 - Math.clz32(f)));
+}
+
+function f32_roundeven(x) {
+  if (!Number.isFinite(x) || Math.abs(x) >= 0x800000) return x;
+  const f = Math.floor(x);
+  const d = x - f;
+  let r = d < 0.5 ? f : d > 0.5 ? f + 1 : (f % 2 === 0 ? f : f + 1);
+  if (r === 0) r = Object.is(x, -0) || x < 0 ? -0 : 0;
+  return Math.fround(r);
+}
+
+function f32_remainder(x, y) {
+  if (Number.isNaN(x) || Number.isNaN(y) || !Number.isFinite(x) || y === 0) {
+    return NaN;
+  }
+  if (!Number.isFinite(y)) return x;
+  const px = f32_decomp(x), py = f32_decomp(y);
+  if (px.m === 0n) return x;
+  const d = px.e - py.e;
+  const num = d >= 0 ? px.m << BigInt(d) : px.m;
+  const den = d >= 0 ? py.m : py.m << BigInt(-d);
+  let q = num / den;
+  const rem = num % den;
+  const twice = rem << 1n;
+  if (twice > den || (twice === den && (q & 1n) !== 0n)) q += 1n;
+  q *= px.s * py.s;
+  const e = Math.min(px.e, py.e);
+  const xi = px.s * (px.m << BigInt(px.e - e));
+  const yi = py.s * (py.m << BigInt(py.e - e));
+  const ri = xi - q * yi;
+  return f32_from_exact(ri, e, px.s < 0n);
+}
+
+function f32_remquo_code(x, y) {
+  if (Number.isNaN(x) || Number.isNaN(y) || !Number.isFinite(x) || y === 0
+      || !Number.isFinite(y)) {
+    return 0;
+  }
+  const px = f32_decomp(x), py = f32_decomp(y);
+  if (px.m === 0n) return 0;
+  const d = px.e - py.e;
+  const num = d >= 0 ? px.m << BigInt(d) : px.m;
+  const den = d >= 0 ? py.m : py.m << BigInt(-d);
+  let q = num / den;
+  const rem = num % den;
+  const twice = rem << 1n;
+  if (twice > den || (twice === den && (q & 1n) !== 0n)) q += 1n;
+  const neg = px.s !== py.s;
+  const mag = Number(q & 7n) >>> 0;
+  return (mag | (neg && q !== 0n ? 0x80000000 : 0)) >>> 0;
+}
+
+function f64_erf_series(x) {
+  const x2 = x * x;
+  let term = x;
+  let sum = x;
+  for (let n = 1; n < 48; ++n) {
+    term *= -x2 / n;
+    const add = term / (2 * n + 1);
+    sum += add;
+    if (Math.abs(add) < 1e-18 * Math.max(1, Math.abs(sum))) break;
+  }
+  return 1.1283791670955125739 * sum;
+}
+
+// Q(1/2,z^2) = erfc(z), evaluated by a bounded continued fraction.  Keeping
+// the iteration in binary64 avoids the severe tail cancellation of 1-erf(z).
+function f64_erfc_pos(z) {
+  if (z < 1.5) return 1 - f64_erf_series(z);
+  if (z > 27) return 0;
+  const x = z * z;
+  const a = 0.5;
+  let b = x + 1 - a;
+  let c = 1e300;
+  let d = 1 / b;
+  let h = d;
+  for (let i = 1; i <= 100; ++i) {
+    const an = -i * (i - a);
+    b += 2;
+    d = an * d + b;
+    if (Math.abs(d) < 1e-300) d = 1e-300;
+    c = b + an / c;
+    if (Math.abs(c) < 1e-300) c = 1e-300;
+    d = 1 / d;
+    const del = d * c;
+    h *= del;
+    if (Math.abs(del - 1) < 2e-16) break;
+  }
+  return Math.exp(-x + a * Math.log(x) - 0.57236494292470008707) * h;
+}
+
+function f32_erf(x) {
+  if (Number.isNaN(x)) return NaN;
+  if (x === Infinity) return 1;
+  if (x === -Infinity) return -1;
+  if (x === 0) return x;
+  const z = Math.abs(x);
+  const y = z < 1.5 ? f64_erf_series(z) : 1 - f64_erfc_pos(z);
+  return Math.fround(x < 0 ? -y : y);
+}
+
+function f32_erfc(x) {
+  if (Number.isNaN(x)) return NaN;
+  if (x === Infinity) return 0;
+  if (x === -Infinity) return 2;
+  if (x === 0) return 1;
+  const q = f64_erfc_pos(Math.abs(x));
+  return Math.fround(x < 0 ? 2 - q : q);
+}
+
+const F32_GAMMA_C = [0.99999999999980993, 676.5203681218851,
+  -1259.1392167224028, 771.32342877765313, -176.61502916214059,
+  12.507343278686905, -0.13857109526572012, 9.984369578019572e-6,
+  1.5056327351493116e-7];
+
+function f64_sinpi(x) {
+  if (!Number.isFinite(x)) return NaN;
+  if (x === 0) return x;
+  const neg = x < 0 || Object.is(x, -0);
+  const a = Math.abs(x);
+  if (a >= 16777216) return neg ? -0 : 0;
+  const r = a % 2;
+  let y;
+  if (r === 0 || r === 1) y = 0;
+  else if (r === 0.5) y = 1;
+  else if (r === 1.5) y = -1;
+  else if (r < 0.5) y = Math.sin(Math.PI * r);
+  else if (r < 1) y = Math.sin(Math.PI * (1 - r));
+  else if (r < 1.5) y = -Math.sin(Math.PI * (r - 1));
+  else y = -Math.sin(Math.PI * (2 - r));
+  return neg ? -y : y;
+}
+
+function f64_lgamma(x) {
+  if (Number.isNaN(x)) return NaN;
+  if (!Number.isFinite(x)) return Infinity;
+  if (x === 1 || x === 2) return 0;
+  if (x <= 0 && Number.isInteger(x)) return Infinity;
+  if (x < 0.5) {
+    return Math.log(Math.PI) - Math.log(Math.abs(f64_sinpi(x)))
+      - f64_lgamma(1 - x);
+  }
+  const z = x - 1;
+  let a = F32_GAMMA_C[0];
+  for (let i = 1; i < F32_GAMMA_C.length; i++) a += F32_GAMMA_C[i] / (z + i);
+  const t = z + 7.5;
+  return 0.9189385332046727 + (z + 0.5) * Math.log(t) - t + Math.log(a);
+}
+
+function f64_tgamma(x) {
+  if (Number.isNaN(x)) return NaN;
+  if (x === Infinity) return Infinity;
+  if (x === -Infinity) return NaN;
+  if (x === 0) return Object.is(x, -0) ? -Infinity : Infinity;
+  if (x < 0 && Number.isInteger(x)) return NaN;
+  if (x < 0.5) return Math.PI / (f64_sinpi(x) * f64_tgamma(1 - x));
+  return Math.exp(f64_lgamma(x));
+}
+
+function f32_tgamma(x) { return Math.fround(f64_tgamma(x)); }
+function f32_lgamma(x) { return Math.fround(f64_lgamma(x)); }
+
+function f32_sinpi(x) { return Math.fround(f64_sinpi(x)); }
+
+function f32_cospi(x) {
+  if (!Number.isFinite(x)) return NaN;
+  const a = Math.abs(x);
+  if (a >= 16777216) return 1;
+  const r = a % 2;
+  if (r === 0) return 1;
+  if (r === 0.5 || r === 1.5) return 0;
+  if (r === 1) return -1;
+  if (r < 0.5) return Math.fround(Math.cos(Math.PI * r));
+  if (r < 1) return Math.fround(-Math.cos(Math.PI * (1 - r)));
+  if (r < 1.5) return Math.fround(-Math.cos(Math.PI * (r - 1)));
+  return Math.fround(Math.cos(Math.PI * (2 - r)));
+}
+
+function f32_tanpi(x) {
+  if (!Number.isFinite(x)) return NaN;
+  if (x === 0) return x;
+  const neg = x < 0 || Object.is(x, -0);
+  const a = Math.abs(x);
+  if (a >= 16777216) return neg ? -0 : 0; // all such F32 integers are even
+  const p = a % 2;
+  const r = a % 1;
+  if (r === 0) {
+    const odd = p === 1;
+    return odd !== neg ? -0 : 0;
+  }
+  if (r === 0.5) {
+    const floorOdd = p === 1.5;
+    return floorOdd !== neg ? -Infinity : Infinity;
+  }
+  const y = r < 0.5 ? Math.tan(Math.PI * r) : -Math.tan(Math.PI * (1 - r));
+  return Math.fround(neg ? -y : y);
+}
+
+`.slice(1);
diff --git a/gates/repo.ts b/gates/repo.ts
index 352704d3a..424665dae 100644
--- a/gates/repo.ts
+++ b/gates/repo.ts
@@ -38,6 +38,8 @@ allow("WONTFIX.txt", 1500);
 allow("LICENSE", 4000);
 allow("flake.nix", 1500);
 allow("bend2/base.bend", 32000);
+allow("bend2/math.bend", 16000);
+allow("bend2/math.ts", 12000);
 allow("bend2/bend.lean", 400000);
 allow("bend2/bend.ts", 42000);
 allow("bend2/comp.ts", 65000);
diff --git a/tests/base/math_edges.bend b/tests/base/math_edges.bend
new file mode 100644
index 000000000..7d394e476
--- /dev/null
+++ b/tests/base/math_edges.bend
@@ -0,0 +1,33 @@
+import Base
+import ../../bend2/math.bend as M
+
+law main.out:
+  String
+
+def main.frexp_show(p: F32 & M.S32) -> String:
+  (m, n) = p
+  F32.show(m) ++ ":" ++ M.S32.show(n)
+
+def main.out():
+  +nz = M.from_bits(2147483648)
+  a = U32.show(F32.bits(M.nextafter(1.0, 0.0)))
+  b = U32.show(F32.bits(M.nextafter(F32.neg(1.0), 0.0)))
+  c = U32.show(F32.bits(M.nextafter(0.0, F32.neg(1.0))))
+  d = Bool.show(M.signbit(M.fmin(0.0, nz)))
+  e = Bool.show(M.signbit(M.fmax(nz, 0.0)))
+  f = F32.show(M.roundeven(2.5))
+  g = F32.show(M.roundeven(3.5))
+  h = F32.show(M.roundeven(F32.neg(2.5)))
+  i = main.frexp_show(M.frexp(12.0))
+  j = Bool.show(M.signbit(M.remainder(F32.neg(4.0), 2.0)))
+  k = M.S32.show(M.ilogb(8.0))
+  a ++ " " ++ b ++ " " ++ c ++ " " ++ d ++ " " ++ e ++ " " ++ f
+    ++ " " ++ g ++ " " ++ h ++ " " ++ i ++ " " ++ j ++ " " ++ k
+
+law main:
+  IO(Unit)
+
+def main():
+  IO.print(main.out())
+
+#|1065353215 3212836863 2147483649 True False 2 4 -2 0.75:4 True 3
diff --git a/tests/base/math_fma.bend b/tests/base/math_fma.bend
new file mode 100644
index 000000000..f2324e3e1
--- /dev/null
+++ b/tests/base/math_fma.bend
@@ -0,0 +1,13 @@
+import Base
+import ../../bend2/math.bend as M
+
+law main:
+  IO(Unit)
+
+def main():
+  a = M.from_bits(2176035992)
+  b = M.from_bits(1882851208)
+  c = M.from_bits(698933837)
+  IO.print(U32.show(F32.bits(M.fma(a, b, c))))
+
+#|2994901125
diff --git a/tests/base/math_h.bend b/tests/base/math_h.bend
new file mode 100644
index 000000000..773dbec4c
--- /dev/null
+++ b/tests/base/math_h.bend
@@ -0,0 +1,27 @@
+# math.h-style library smoke test: constants, rounding, bit stepping,
+# transcendental helpers, gamma and classification all cross a module import.
+import Base
+import ../../bend2/math.bend as M
+
+law main.out:
+  String
+
+def main.out():
+  a = F32.show(M.M_PI())
+  b = F32.show(M.M_SQRT2())
+  c = F32.show(M.round(F32.neg(1.5)))
+  d = U32.show(F32.bits(M.nextafter(1.0, 2.0)))
+  e = F32.show(M.tgamma(5.0))
+  f = F32.show(M.erf(1.0))
+  g = Bool.show(M.isfinite(M.INFINITY()))
+  h = Bool.show(M.isnan(M.NAN()))
+  a ++ " " ++ b ++ " " ++ c ++ " " ++ d ++ " " ++ e ++ " " ++ f
+    ++ " " ++ g ++ " " ++ h
+
+law main:
+  IO(Unit)
+
+def main():
+  IO.print(main.out())
+
+#|3.1415927 1.4142135 -2 1065353217 24 0.8427008 False True
diff --git a/tests/base/math_regressions.bend b/tests/base/math_regressions.bend
new file mode 100644
index 000000000..bc951681b
--- /dev/null
+++ b/tests/base/math_regressions.bend
@@ -0,0 +1,81 @@
+import Base
+import ../../bend2/math.bend as M
+
+law main.frexp_max:
+  String
+
+def main.frexp_pair(p: F32 & M.S32) -> String:
+  (m, e) = p
+  U32.show(F32.bits(m)) ++ ":" ++ M.S32.show(e)
+
+def main.frexp_max():
+  main.frexp_pair(M.frexp(M.from_bits(2139095039)))
+
+law main.modf_inf:
+  String
+
+def main.modf_pair(p: F32 & F32) -> String:
+  (f, i) = p
+  U32.show(F32.bits(f)) ++ ":" ++ Bool.show(M.isinf(i))
+
+def main.modf_inf():
+  main.modf_pair(M.modf(M.INFINITY()))
+
+law main.remquo_neg:
+  String
+
+def main.remquo_pair(p: F32 & M.S32) -> String:
+  (r, q) = p
+  U32.show(F32.bits(r)) ++ ":" ++ M.S32.show(q)
+
+def main.remquo_neg():
+  main.remquo_pair(M.remquo(F32.neg(23.5), 8.0))
+
+law main:
+  IO(Unit)
+
+def main():
+  min_norm = M.from_bits(8388608)
+  +max_fin = M.from_bits(2139095039)
+  +nz = M.from_bits(2147483648)
+  a = F32.show(M.scalbn(min_norm, M.S32.pos(128)))
+  b = F32.show(M.scalbn(max_fin, M.S32.neg(150)))
+  c = main.frexp_max()
+  d = F32.show(M.logb(max_fin))
+  e = M.S32.show(M.ilogb(max_fin))
+  f = main.modf_inf()
+  g = F32.show(M.remainder(2.0, M.INFINITY()))
+  h = U32.show(F32.bits(M.tgamma(0.0)))
+  i = U32.show(F32.bits(M.tgamma(nz)))
+  j = Bool.show(M.isnan(M.tgamma(F32.neg(M.INFINITY()))))
+  k = Bool.show(F32.is_gt(M.erfc(5.0), 0.0))
+  l = Bool.show(F32.is_ne(M.asinh(0.00000001), 0.0))
+  m = U32.show(F32.bits(M.sinpi(1.0)))
+  n = U32.show(F32.bits(M.sinpi(16777216.0)))
+  o = U32.show(F32.bits(M.pown(F32.neg(1.0), M.S32.pos(16777217))))
+  p = Bool.show(M.signbit(M.roundeven(F32.neg(0.5))))
+  q = U32.show(F32.bits(M.fma(M.from_bits(1687846912),
+    M.from_bits(2398336256), M.from_bits(2574717546))))
+  r = main.remquo_neg()
+  s = U32.show(F32.bits(M.rootn(nz, M.S32.pos(3))))
+  t = U32.show(F32.bits(M.rootn(nz, M.S32.neg(3))))
+  u = U32.show(F32.bits(M.round(8388609.0)))
+  v = U32.show(F32.bits(M.round(F32.neg(8388609.0))))
+  w = U32.show(F32.bits(M.tanpi(1.0)))
+  x = U32.show(F32.bits(M.tanpi(F32.neg(1.0))))
+  y = U32.show(F32.bits(M.tanpi(1.5)))
+  z = U32.show(F32.bits(M.tanpi(F32.neg(1.5))))
+  aa = Bool.show(M.isnan(M.powr(0.0, nz)))
+  ab = Bool.show(M.isnan(M.powr(0.0, M.INFINITY())))
+  ac = Bool.show(M.isnan(M.powr(M.INFINITY(), 1.0)))
+  ad = U32.show(F32.bits(M.powr(F32.neg(2.0), nz)))
+  ae = U32.show(F32.bits(M.powr(3.0, nz)))
+  af = U32.show(F32.bits(M.rootn(nz, M.S32.pos(2))))
+  ag = U32.show(F32.bits(M.rootn(nz, M.S32.neg(2))))
+  IO.print(a ++ " " ++ b ++ " " ++ c ++ " " ++ d ++ " " ++ e ++ " "
+    ++ f ++ " " ++ g ++ " " ++ h ++ " " ++ i ++ " " ++ j ++ " " ++ k
+    ++ " " ++ l ++ " " ++ m ++ " " ++ n ++ " " ++ o ++ " " ++ p ++ " " ++ q
+    ++ " " ++ r ++ " " ++ s ++ " " ++ t ++ " " ++ u ++ " " ++ v ++ " " ++ w
+    ++ " " ++ x ++ " " ++ y ++ " " ++ z ++ " " ++ aa ++ " " ++ ab ++ " " ++ ac
+    ++ " " ++ ad ++ " " ++ ae ++ " " ++ af ++ " " ++ ag)
+#|4 2.3841856e-7 1065353215:128 127 127 0:True 2 2139095040 4286578688 True True True 0 0 3212836864 True 3021149373 1056964608:-3 2147483648 4286578688 1258291201 3405774849 2147483648 0 4286578688 2139095040 True True True 2139095040 0 0 2139095040
diff --git a/tests/base/math_release.bend b/tests/base/math_release.bend
new file mode 100644
index 000000000..3d3f60064
--- /dev/null
+++ b/tests/base/math_release.bend
@@ -0,0 +1,28 @@
+import Base
+import ../../bend2/math.bend as M
+
+law main:
+  IO(Unit)
+
+def main():
+  +nz = M.from_bits(2147483648)
+  a = U32.show(F32.bits(M.lgamma(1.0)))
+  b = U32.show(F32.bits(M.lgamma(2.0)))
+  c = U32.show(F32.bits(M.pow(1.0, M.INFINITY())))
+  d = U32.show(F32.bits(M.pow(1.0, M.NAN())))
+  e = U32.show(F32.bits(M.pow(F32.neg(1.0), M.INFINITY())))
+  f = U32.show(F32.bits(M.fma(0.0, F32.neg(1.0), nz)))
+  g = U32.show(F32.bits(M.fma(nz, F32.neg(1.0), nz)))
+  h = U32.show(F32.bits(M.pown(2647993.0, M.S32.neg(6))))
+  i = U32.show(F32.bits(M.exp2m1(nz)))
+  j = U32.show(F32.bits(M.exp10m1(nz)))
+  k = U32.show(F32.bits(M.asinpi(nz)))
+  l = U32.show(F32.bits(M.rsqrt(nz)))
+  m = U32.show(F32.bits(M.powr(3.0, 2.0)))
+  n = Bool.show(M.isnan(M.powr(0.0, 0.0)))
+  o = U32.show(F32.bits(M.tanpi(F32.neg(1.0))))
+  p = U32.show(F32.bits(M.tanpi(1.0)))
+  IO.print(a ++ " " ++ b ++ " " ++ c ++ " " ++ d ++ " " ++ e ++ " "
+    ++ f ++ " " ++ g ++ " " ++ h ++ " " ++ i ++ " " ++ j ++ " " ++ k
+    ++ " " ++ l ++ " " ++ m ++ " " ++ n ++ " " ++ o ++ " " ++ p)
+#|0 0 1065353216 1065353216 1065353216 2147483648 0 2069991 2147483648 2147483648 2147483648 4286578688 1090519040 True 0 2147483648