diff --git a/internal/kernels/asmgen/arm64/gen.go b/internal/kernels/asmgen/arm64/gen.go
index d07eb3f..6d3f989 100644
--- a/internal/kernels/asmgen/arm64/gen.go
+++ b/internal/kernels/asmgen/arm64/gen.go
@@ -66,6 +66,8 @@ func main() {
 		}
 	}
 	f.Add(blocks8Kernel())
+	f.Add(exp32CRKernel())
+	f.Add(log32CRKernel())
 	if err := os.WriteFile("sum_arm64.s", []byte(f.String()), 0o644); err != nil {
 		fmt.Fprintln(os.Stderr, err)
 		os.Exit(1)
@@ -766,3 +768,235 @@ func sumKernel() *emit.Function {
 		Ret()
 	return b.Func()
 }
+
+// simd2 emits an Advanced SIMD two-register instruction that Go's assembler
+// does not name, by its encoding (Arm ARM, "Advanced SIMD two-register
+// miscellaneous"), with the mnemonic as a comment. A wrong encoding cannot
+// pass TestExp32NEONMatchesGo, which holds every lane to exp32CR bit for bit.
+func simd2(base uint32, mnem string, d, n int) string {
+	return fmt.Sprintf("WORD $0x%08X // "+mnem, base|uint32(n)<<5|uint32(d), d, n)
+}
+
+// facge emits FACGE Vd.4S, Vn.4S, Vm.4S ("Advanced SIMD three same", U=1,
+// opcode 11101): all ones where |Vn| >= |Vm|, zero otherwise and for NaN.
+func facge(d, n, m int) string {
+	return fmt.Sprintf("WORD $0x%08X // FACGE V%d.4S, V%d.4S, V%d.4S",
+		0x6E20EC00|uint32(m)<<16|uint32(n)<<5|uint32(d), d, n, m)
+}
+
+const (
+	fcvtl  = 0x0E617800 // FCVTL Vd.2D, Vn.2S: the low two floats, widened
+	fcvtl2 = 0x4E617800 // FCVTL2 Vd.2D, Vn.4S: the high two
+	fcvtn  = 0x0E616800 // FCVTN Vd.2S, Vn.2D: rounded once (FPCR: nearest), upper half zeroed
+	fcvtn2 = 0x4E616800 // FCVTN2 Vd.4S, Vn.2D: into the upper half, lower kept
+	uxtl   = 0x2F20A400 // UXTL Vd.2D, Vn.2S (USHLL #0): the low two words, zero-extended
+	uxtl2  = 0x6F20A400 // UXTL2 Vd.2D, Vn.4S: the high two
+	sxtl   = 0x0F20A400 // SXTL Vd.2D, Vn.2S (SSHLL #0): the low two words, sign-extended
+	sxtl2  = 0x4F20A400 // SXTL2 Vd.2D, Vn.4S: the high two
+	scvtf  = 0x4E61D800 // SCVTF Vd.2D, Vn.2D: int64 to float64 (exact below 2^53)
+)
+
+// simd3 emits an Advanced SIMD three-register instruction Go's assembler
+// does not name: Vd = Vn OP Vm, with base the encoding's fixed bits.
+func simd3(base uint32, mnem string, d, n, m int) string {
+	return fmt.Sprintf("WORD $0x%08X // "+mnem, base|uint32(m)<<16|uint32(n)<<5|uint32(d), d, n, m)
+}
+
+const (
+	cmgt4s = 0x4EA03400 // CMGT Vd.4S, Vn.4S, Vm.4S: Vn > Vm as int32, all ones or zero
+	bic16b = 0x4E601C00 // BIC Vd.16B, Vn.16B, Vm.16B: Vn AND NOT Vm
+)
+
+// log32CRKernel builds log32CRNEON(dst, src *float32, n int, k
+// *logCRNEONConsts) (done int): logCR's fast path (logcr.go, CORE-MATH's
+// logf) on eight float32 per iteration, n a positive multiple of 8, the NEON
+// twin of log32CRAVX2 (asmgen/amd64/gen.go). Per quad of four floats:
+//
+//   - in 32-bit lanes, m = ux & 0x7fffff, j = (m + 2^16) >> 17 (at most 64
+//     for any ux, so the lookups stay in the 65-entry tables), e = (ux>>23) -
+//     127, x == 1, and "not a positive normal" (ux + 0x7f800000 > -0x01000001
+//     as int32);
+//   - per pair, in float64 with logCR's unfused operations in its order
+//     (VFMUL then VFADD, never FMLA): tz = (m | 1023<<23) << 29 (UXTL), e
+//     (SXTL, SCVTF), tr[j] and tl[j] through general registers, z, r, and the
+//     bounds r and r + 0x1.f06p-33 rounded once (FCVTN, FCVTN2).
+//
+// A block of 8 with a lane that is not a positive normal, or whose bounds
+// round apart (x == 1 excepted: its result is +0), is not stored; the kernel
+// returns the count it wrote and the Go wrapper finishes that block with
+// logCR. Constants live in V16-V28.
+func log32CRKernel() *emit.Function {
+	sig := arm64.Layout(
+		[]string{"dst", "src", "n", "k"},
+		[]arm64.Type{arm64.Ptr, arm64.Ptr, arm64.Int64, arm64.Ptr},
+		[]string{"done"}, []arm64.Type{arm64.Int64},
+	)
+	const tr, tl = 208, 728 // offsets of the tables in logCRNEONConsts
+	// V16 one, V17 ln2, V18 b0, V19 b1, V20 b2, V21 w (float64 pairs); V22
+	// mmask, V23 jadd, V24 bias, V25 sub, V26 maxv, V27 oneBits (int32
+	// quads); V28 expOr (int64 pair).
+	b := arm64.NewFunc("log32CRNEON", sig, 0)
+	b.LoadArg("dst", "R0").LoadArg("src", "R1").LoadArg("n", "R2").LoadArg("k", "R8")
+	b.Raw("MOVD R8, R13")
+	b.Raw("VLD1.P 64(R13), [V16.D2, V17.D2, V18.D2, V19.D2]")
+	b.Raw("VLD1.P 64(R13), [V20.D2, V21.D2, V22.D2, V23.D2]")
+	b.Raw("VLD1.P 64(R13), [V24.D2, V25.D2, V26.D2, V27.D2]")
+	b.Raw("VLD1 (R13), [V28.D2]")
+	b.Raw(fmt.Sprintf("ADD $%d, R8, R9", tr))
+	b.Raw(fmt.Sprintf("ADD $%d, R8, R10", tl))
+	b.Raw("MOVD ZR, R3")
+	// quad runs the fast path on the four floats in x, leaving the results
+	// (x == 1 lanes as +0) in res and the lanes that may be stored in ok.
+	// Integer part in V2-V6, pair temporaries V12-V15 and V29-V31, bounds
+	// in V7 (ub) and V8 (lb).
+	quad := func(x, res, ok int) {
+		b.Raw(fmt.Sprintf("VAND V22.B16, V%d.B16, V2.B16", x)) // m
+		b.Raw("VADD V23.S4, V2.S4, V3.S4")
+		b.Raw("VUSHR $17, V3.S4, V3.S4") // j
+		b.Raw(fmt.Sprintf("VUSHR $23, V%d.S4, V4.S4", x))
+		b.Raw("VSUB V24.S4, V4.S4, V4.S4")                   // e
+		b.Raw(fmt.Sprintf("VCMEQ V27.S4, V%d.S4, V5.S4", x)) // x == 1
+		b.Raw(fmt.Sprintf("VADD V25.S4, V%d.S4, V6.S4", x))
+		b.Raw(simd3(cmgt4s, "CMGT V%d.4S, V%d.4S, V%d.4S", 6, 6, 26)) // not a positive normal
+		for h := 0; h < 2; h++ {
+			var wide, narrow, sx uint32 = uxtl, fcvtn, sxtl
+			if h == 1 {
+				wide, narrow, sx = uxtl2, fcvtn2, sxtl2
+			}
+			b.Raw(simd2(wide, "UXTL V%d.2D, V%d.S", 12, 2))
+			b.Raw("VORR V28.B16, V12.B16, V12.B16")
+			b.Raw("VSHL $29, V12.D2, V12.D2") // tz, in [1,2)
+			b.Raw(simd2(sx, "SXTL V%d.2D, V%d.S", 13, 4))
+			b.Raw(simd2(scvtf, "SCVTF V%d.2D, V%d.2D", 13, 13)) // e as float64
+			b.Raw(fmt.Sprintf("VMOV V3.S[%d], R4", 2*h))
+			b.Raw(fmt.Sprintf("VMOV V3.S[%d], R5", 2*h+1))
+			b.Raw("MOVD (R9)(R4<<3), R6").Raw("MOVD (R9)(R5<<3), R7")
+			b.Raw("MOVD (R10)(R4<<3), R11").Raw("MOVD (R10)(R5<<3), R12")
+			b.Raw("VMOV R6, V14.D[0]").Raw("VMOV R7, V14.D[1]")   // tr[j]
+			b.Raw("VMOV R11, V15.D[0]").Raw("VMOV R12, V15.D[1]") // tl[j]
+			b.Raw("VFMUL V14.D2, V12.D2, V14.D2")                 // tz*tr
+			b.Raw("VFSUB V16.D2, V14.D2, V14.D2")                 // z
+			b.Raw("VFMUL V14.D2, V14.D2, V29.D2")                 // z2
+			b.Raw("VFMUL V17.D2, V13.D2, V31.D2")                 // e*ln2
+			b.Raw("VFADD V15.D2, V31.D2, V31.D2")                 // + tl
+			b.Raw("VFMUL V18.D2, V14.D2, V12.D2")
+			b.Raw("VFADD V12.D2, V31.D2, V31.D2") // + z*b0
+			b.Raw("VFMUL V20.D2, V14.D2, V12.D2")
+			b.Raw("VFADD V19.D2, V12.D2, V12.D2") // b1 + z*b2
+			b.Raw("VFMUL V12.D2, V29.D2, V12.D2")
+			b.Raw("VFADD V12.D2, V31.D2, V31.D2") // r
+			b.Raw("VFADD V21.D2, V31.D2, V30.D2") // r + 0x1.f06p-33
+			b.Raw(simd2(narrow, "FCVTN V%d.S, V%d.2D", 7, 31))
+			b.Raw(simd2(narrow, "FCVTN V%d.S, V%d.2D", 8, 30))
+		}
+		b.Raw(fmt.Sprintf("VCMEQ V8.S4, V7.S4, V%d.S4", ok))        // the bounds agree
+		b.Raw(fmt.Sprintf("VORR V5.B16, V%d.B16, V%d.B16", ok, ok)) // or x == 1
+		b.Raw(simd3(bic16b, "BIC V%d.16B, V%d.16B, V%d.16B", ok, ok, 6))
+		b.Raw(simd3(bic16b, "BIC V%d.16B, V%d.16B, V%d.16B", res, 7, 5)) // x == 1 -> +0
+	}
+	b.Raw("lloop:")
+	b.Raw("VLD1 (R1), [V0.S4, V1.S4]")
+	quad(0, 9, 11)
+	quad(1, 10, 8)
+	b.Raw("VAND V8.B16, V11.B16, V11.B16")
+	b.Raw("VMOV V11.D[0], R6")
+	b.Raw("VMOV V11.D[1], R7")
+	b.Raw("AND R6, R7, R7")
+	b.Raw("CMN $1, R7") // every lane may be stored
+	b.Raw("BNE lbail")
+	b.Raw("VST1 [V9.S4, V10.S4], (R0)")
+	b.Raw("ADD $32, R0").Raw("ADD $32, R1").Raw("ADD $8, R3")
+	b.Raw("CMP R2, R3").Raw("BLT lloop")
+	b.Raw("lbail:")
+	b.StoreRet("R3", "done")
+	b.Ret()
+	return b.Func()
+}
+
+// exp32CRKernel builds exp32CRNEON(dst, src *float32, n int, k
+// *expCRNEONConsts) (done int): expCR (expcr.go, Arm's correctly rounded
+// cr_expf, an AdvSIMD routine in the first place) on eight float32 per
+// iteration, n a positive multiple of 8. Each quad is widened to two pairs of
+// doubles (FCVTL, FCVTL2); each pair runs expCR's fused operations in its
+// order (VFMLA and VFMLS round once, as math.FMA does); the table entries
+// z&63 are fetched through two general registers, as cr_expf does, NEON
+// having no gather; the result is rounded once (FCVTN, FCVTN2).
+//
+// A block holding a lane with |x| > 512 or a NaN (FACGE against 512 is false
+// for both) is not stored: the kernel returns the count it wrote and the Go
+// wrapper finishes that block with exp32CR. Constants live in V16-V26.
+func exp32CRKernel() *emit.Function {
+	sig := arm64.Layout(
+		[]string{"dst", "src", "n", "k"},
+		[]arm64.Type{arm64.Ptr, arm64.Ptr, arm64.Int64, arm64.Ptr},
+		[]string{"done"}, []arm64.Type{arm64.Int64},
+	)
+	const table = 176 // offset of the 64 mantissas in expCRNEONConsts
+	b := arm64.NewFunc("exp32CRNEON", sig, 0)
+	b.LoadArg("dst", "R0").LoadArg("src", "R1").LoadArg("n", "R2").LoadArg("k", "R8")
+	b.Raw("MOVD R8, R10")
+	b.Raw("VLD1.P 64(R10), [V16.D2, V17.D2, V18.D2, V19.D2]") // shift, invLn2, ln2Hi, ln2Lo
+	b.Raw("VLD1.P 64(R10), [V20.D2, V21.D2, V22.D2, V23.D2]") // c0, c1, c2, c3
+	b.Raw("VLD1 (R10), [V24.D2, V25.D2, V26.D2]")             // idxMask, expMask, 512 as 4 floats
+	b.Raw(fmt.Sprintf("ADD $%d, R8, R9", table))
+	b.Raw("MOVD ZR, R3")
+	// half runs expCR on the two doubles in x, leaving the result in t[0].
+	// t holds z, n, r2, p01, p23, idx, mant; ra and rb carry the lookup.
+	half := func(x int, t [7]int, ra, rb string) {
+		z, nn, r2, p01, p23, idx, mant := t[0], t[1], t[2], t[3], t[4], t[5], t[6]
+		b.Raw(fmt.Sprintf("VMOV V16.B16, V%d.B16", z))
+		b.Raw(fmt.Sprintf("VFMLA V17.D2, V%d.D2, V%d.D2", x, z))     // z = x*invLn2 + shift
+		b.Raw(fmt.Sprintf("VFSUB V16.D2, V%d.D2, V%d.D2", z, nn))    // n = z - shift
+		b.Raw(fmt.Sprintf("VFMLS V18.D2, V%d.D2, V%d.D2", nn, x))    // r = x - n*ln2Hi
+		b.Raw(fmt.Sprintf("VFMLS V19.D2, V%d.D2, V%d.D2", nn, x))    // r -= n*ln2Lo
+		b.Raw(fmt.Sprintf("VFMUL V%d.D2, V%d.D2, V%d.D2", x, x, r2)) // r2
+		b.Raw(fmt.Sprintf("VMOV V20.B16, V%d.B16", p01))
+		b.Raw(fmt.Sprintf("VFMLA V21.D2, V%d.D2, V%d.D2", x, p01)) // p01 = r*c1 + c0
+		b.Raw(fmt.Sprintf("VMOV V22.B16, V%d.B16", p23))
+		b.Raw(fmt.Sprintf("VFMLA V23.D2, V%d.D2, V%d.D2", x, p23))       // p23 = r*c3 + c2
+		b.Raw(fmt.Sprintf("VFMLA V%d.D2, V%d.D2, V%d.D2", p23, r2, p01)) // p04 = r2*p23 + p01
+		b.Raw(fmt.Sprintf("VFMLA V%d.D2, V%d.D2, V%d.D2", p01, r2, x))   // y = r2*p04 + r
+		b.Raw(fmt.Sprintf("VAND V24.B16, V%d.B16, V%d.B16", z, idx))     // z & 63
+		b.Raw(fmt.Sprintf("VMOV V%d.D[0], %s", idx, ra))
+		b.Raw(fmt.Sprintf("VMOV V%d.D[1], %s", idx, rb))
+		b.Raw(fmt.Sprintf("MOVD (R9)(%s<<3), %s", ra, ra))
+		b.Raw(fmt.Sprintf("MOVD (R9)(%s<<3), %s", rb, rb))
+		b.Raw(fmt.Sprintf("VMOV %s, V%d.D[0]", ra, mant))
+		b.Raw(fmt.Sprintf("VMOV %s, V%d.D[1]", rb, mant))
+		b.Raw(fmt.Sprintf("VSHL $46, V%d.D2, V%d.D2", z, z))
+		b.Raw(fmt.Sprintf("VAND V25.B16, V%d.B16, V%d.B16", z, z))       // exponent and sign
+		b.Raw(fmt.Sprintf("VORR V%d.B16, V%d.B16, V%d.B16", mant, z, z)) // s
+		b.Raw(fmt.Sprintf("VFMLA V%d.D2, V%d.D2, V%d.D2", x, z, z))      // s*y + s
+	}
+	t0 := [7]int{3, 4, 5, 6, 7, 14, 15}
+	t1 := [7]int{27, 28, 29, 30, 31, 12, 13}
+	b.Raw("eloop:")
+	b.Raw("VLD1 (R1), [V0.S4, V1.S4]")
+	b.Raw(facge(12, 26, 0)) // 512 >= |x|: false for NaN
+	b.Raw(facge(13, 26, 1))
+	b.Raw("VAND V13.B16, V12.B16, V12.B16")
+	b.Raw("VMOV V12.D[0], R6")
+	b.Raw("VMOV V12.D[1], R7")
+	b.Raw("AND R6, R7, R7")
+	b.Raw("CMN $1, R7") // every lane in range
+	b.Raw("BNE ebail")
+	b.Raw(simd2(fcvtl, "FCVTL V%d.2D, V%d.2S", 2, 0))
+	b.Raw(simd2(fcvtl2, "FCVTL2 V%d.2D, V%d.4S", 9, 0))
+	b.Raw(simd2(fcvtl, "FCVTL V%d.2D, V%d.2S", 10, 1))
+	b.Raw(simd2(fcvtl2, "FCVTL2 V%d.2D, V%d.4S", 11, 1))
+	half(2, t0, "R4", "R5")
+	half(9, t1, "R6", "R7")
+	b.Raw(simd2(fcvtn, "FCVTN V%d.2S, V%d.2D", 0, 3))
+	b.Raw(simd2(fcvtn2, "FCVTN2 V%d.4S, V%d.2D", 0, 27))
+	half(10, t0, "R4", "R5")
+	half(11, t1, "R6", "R7")
+	b.Raw(simd2(fcvtn, "FCVTN V%d.2S, V%d.2D", 1, 3))
+	b.Raw(simd2(fcvtn2, "FCVTN2 V%d.4S, V%d.2D", 1, 27))
+	b.Raw("VST1 [V0.S4, V1.S4], (R0)")
+	b.Raw("ADD $32, R0").Raw("ADD $32, R1").Raw("ADD $8, R3")
+	b.Raw("CMP R2, R3").Raw("BLT eloop")
+	b.Raw("ebail:")
+	b.StoreRet("R3", "done")
+	b.Ret()
+	return b.Func()
+}
diff --git a/internal/kernels/expcr_arm64.go b/internal/kernels/expcr_arm64.go
new file mode 100644
index 0000000..8db64f5
--- /dev/null
+++ b/internal/kernels/expcr_arm64.go
@@ -0,0 +1,43 @@
+package kernels
+
+// expCRNEONConsts is the constant block exp32CRNEON reads (asmgen/arm64/
+// gen.go, exp32CRKernel) into V16-V26, then the table it indexes.
+type expCRNEONConsts struct {
+	shift, invLn2, ln2Hi, ln2Lo [2]float64 // offsets 0, 16, 32, 48
+	c0, c1, c2, c3              [2]float64 // 64, 80, 96, 112
+	idxMask, expMask            [2]uint64  // 128, 144
+	lim                         [4]float32 // 160
+	table                       [64]uint64 // 176
+}
+
+var expCRKN = func() *expCRNEONConsts {
+	k := new(expCRNEONConsts)
+	for i := range 2 {
+		k.shift[i], k.invLn2[i], k.ln2Hi[i], k.ln2Lo[i] = crShift, crInvLn2, crLn2Hi, crLn2Lo
+		k.c0[i], k.c1[i], k.c2[i], k.c3[i] = crC0, crC1, crC2, crC3
+		k.idxMask[i], k.expMask[i] = 63, 0xfff0000000000000
+	}
+	k.lim = [4]float32{512, 512, 512, 512}
+	k.table = crMant
+	return k
+}()
+
+// exp32Fast is the correctly rounded float32 exp: blocks of 8 on the NEON
+// kernel, which hands back each block holding a lane with |x| > 512 or a NaN
+// to exp32CR, as is the tail. Both run expCR's operations in its order, so
+// they agree bit for bit.
+func exp32Fast(dst, src []float32) {
+	src = src[:len(dst)]
+	i, n := 0, len(dst)&^7
+	for i < n {
+		i += exp32CRNEON(&dst[i], &src[i], n-i, expCRKN)
+		if i < n {
+			exp32CR(dst[i:i+8], src[i:i+8])
+			i += 8
+		}
+	}
+	exp32CR(dst[i:], src[i:])
+}
+
+//go:noescape
+func exp32CRNEON(dst, src *float32, n int, k *expCRNEONConsts) (done int)
diff --git a/internal/kernels/expcr_arm64_test.go b/internal/kernels/expcr_arm64_test.go
new file mode 100644
index 0000000..19dcf16
--- /dev/null
+++ b/internal/kernels/expcr_arm64_test.go
@@ -0,0 +1,74 @@
+package kernels
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+	"unsafe"
+)
+
+// TestExpCRNEONConstsLayout holds expCRNEONConsts to the offsets
+// exp32CRKernel loads.
+func TestExpCRNEONConstsLayout(t *testing.T) {
+	var k expCRNEONConsts
+	for _, c := range []struct {
+		name      string
+		got, want uintptr
+	}{
+		{"shift", unsafe.Offsetof(k.shift), 0}, {"ln2Lo", unsafe.Offsetof(k.ln2Lo), 48},
+		{"c0", unsafe.Offsetof(k.c0), 64}, {"c3", unsafe.Offsetof(k.c3), 112},
+		{"idxMask", unsafe.Offsetof(k.idxMask), 128}, {"expMask", unsafe.Offsetof(k.expMask), 144},
+		{"lim", unsafe.Offsetof(k.lim), 160}, {"table", unsafe.Offsetof(k.table), 176},
+	} {
+		if c.got != c.want {
+			t.Errorf("%s at %d, the kernel reads %d", c.name, c.got, c.want)
+		}
+	}
+}
+
+// TestExp32NEONMatchesGo holds exp32Fast (NEON blocks, Go for the blocks it
+// hands back and the tail) to exp32CR bit for bit (NaN to NaN): random bit
+// patterns (most blocks hand back), random inputs in [-110, 110] (none
+// does), the special values and the ±512 edge, every length 0..70, and in
+// place.
+func TestExp32NEONMatchesGo(t *testing.T) {
+	r := rand.New(rand.NewSource(7))
+	src := make([]float32, 0, 300000)
+	for _, v := range []float32{0, float32(math.Copysign(0, -1)), 1, -1, 88.72283, 88.72284, 89, -103.972, -103.973, -104,
+		511.9, 512, 512.1, -512, -512.1, 1e30, -1e30, float32(math.Inf(1)), float32(math.Inf(-1)),
+		float32(math.NaN()), math.Float32frombits(0xffc00000), math.SmallestNonzeroFloat32} {
+		src = append(src, v)
+	}
+	for len(src) < 100000 {
+		src = append(src, math.Float32frombits(r.Uint32()))
+	}
+	for len(src) < cap(src) {
+		src = append(src, float32(r.Float64()*220-110))
+	}
+	want := make([]float32, len(src))
+	exp32CR(want, src)
+	same := func(a, b float32) bool { return math.Float32bits(a) == math.Float32bits(b) || (a != a && b != b) }
+	got := make([]float32, len(src))
+	exp32Fast(got, src)
+	for i := range src {
+		if !same(got[i], want[i]) {
+			t.Fatalf("exp(0x%08x): NEON 0x%08x, Go 0x%08x", math.Float32bits(src[i]), math.Float32bits(got[i]), math.Float32bits(want[i]))
+		}
+	}
+	for n := 0; n <= 70; n++ {
+		d := make([]float32, n)
+		exp32Fast(d, src[:n])
+		for i := range d {
+			if !same(d[i], want[i]) {
+				t.Fatalf("length %d [%d]: 0x%08x, want 0x%08x", n, i, math.Float32bits(d[i]), math.Float32bits(want[i]))
+			}
+		}
+	}
+	in := append([]float32(nil), src[:4099]...)
+	exp32Fast(in, in)
+	for i := range in {
+		if !same(in[i], want[i]) {
+			t.Fatalf("in place [%d]: 0x%08x, want 0x%08x", i, math.Float32bits(in[i]), math.Float32bits(want[i]))
+		}
+	}
+}
diff --git a/internal/kernels/expcr_fence_arm64_test.go b/internal/kernels/expcr_fence_arm64_test.go
new file mode 100644
index 0000000..f473f39
--- /dev/null
+++ b/internal/kernels/expcr_fence_arm64_test.go
@@ -0,0 +1,28 @@
+//go:build linux || darwin
+
+package kernels
+
+import "testing"
+
+// TestExp32NEONStaysInside fences src and dst of the NEON exp kernel
+// (fences32 needs mmap/mprotect, hence the build tag of fence_test.go).
+func TestExp32NEONStaysInside(t *testing.T) {
+	for n := 8; n <= 64; n += 8 {
+		for _, s := range fences32(t, n) {
+			for _, d := range fences32(t, n) {
+				exp32CRNEON(&d[0], &s[0], n, expCRKN)
+			}
+		}
+	}
+}
+
+// TestLog32NEONStaysInside fences src and dst of the NEON log kernel.
+func TestLog32NEONStaysInside(t *testing.T) {
+	for n := 8; n <= 64; n += 8 {
+		for _, s := range fences32(t, n) {
+			for _, d := range fences32(t, n) {
+				log32CRNEON(&d[0], &s[0], n, logCRKN)
+			}
+		}
+	}
+}
diff --git a/internal/kernels/expcr_other.go b/internal/kernels/expcr_other.go
index 80b6267..eca7ad4 100644
--- a/internal/kernels/expcr_other.go
+++ b/internal/kernels/expcr_other.go
@@ -1,4 +1,4 @@
-//go:build !amd64
+//go:build !amd64 && !arm64
 
 package kernels
 
diff --git a/internal/kernels/logcr_arm64.go b/internal/kernels/logcr_arm64.go
new file mode 100644
index 0000000..31bc972
--- /dev/null
+++ b/internal/kernels/logcr_arm64.go
@@ -0,0 +1,46 @@
+package kernels
+
+// logCRNEONConsts is the constant block log32CRNEON loads into V16-V28
+// (asmgen/arm64/gen.go, log32CRKernel), then the tables it indexes.
+type logCRNEONConsts struct {
+	one, ln2, b0, b1, b2, w               [2]float64  // offsets 0, 16, 32, 48, 64, 80
+	mmask, jadd, bias, sub, maxv, oneBits [4]uint32   // 96, 112, 128, 144, 160, 176
+	expOr                                 [2]uint64   // 192
+	tr, tl                                [65]float64 // 208, 728
+}
+
+var logCRKN = func() *logCRNEONConsts {
+	k := new(logCRNEONConsts)
+	for i := range 2 {
+		k.one[i], k.ln2[i], k.w[i] = 1, 0x1.62e42fefa39efp-1, 0x1.f06p-33
+		k.b0[i], k.b1[i], k.b2[i] = logCRB[0], logCRB[1], logCRB[2]
+		k.expOr[i] = 1023 << 23
+	}
+	for i := range 4 {
+		k.mmask[i], k.jadd[i], k.bias[i] = 1<<23-1, 1<<16, 127
+		// The positive-normal test: ux + sub (wrapping) <= maxv as int32.
+		k.sub[i], k.maxv[i], k.oneBits[i] = 0x7f800000, 0xfeffffff, 0x3f800000
+	}
+	k.tr, k.tl = logCRTr, logCRTl
+	return k
+}()
+
+// log32Fast is the correctly rounded float32 log: blocks of 8 on the NEON
+// kernel, which hands back each block it cannot finish (an input that is not
+// a positive normal float32, or bounds that round apart) to logCR, as is the
+// tail. Both run logCR's operations in its order, so they agree bit for bit.
+func log32Fast(dst, src []float32) {
+	src = src[:len(dst)]
+	i, n := 0, len(dst)&^7
+	for i < n {
+		i += log32CRNEON(&dst[i], &src[i], n-i, logCRKN)
+		if i < n {
+			logCRSlice(dst[i:i+8], src[i:i+8])
+			i += 8
+		}
+	}
+	logCRSlice(dst[i:], src[i:])
+}
+
+//go:noescape
+func log32CRNEON(dst, src *float32, n int, k *logCRNEONConsts) (done int)
diff --git a/internal/kernels/logcr_arm64_test.go b/internal/kernels/logcr_arm64_test.go
new file mode 100644
index 0000000..b29d7dd
--- /dev/null
+++ b/internal/kernels/logcr_arm64_test.go
@@ -0,0 +1,73 @@
+package kernels
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+	"unsafe"
+)
+
+// TestLogCRNEONConstsLayout holds logCRNEONConsts to the offsets
+// log32CRKernel loads.
+func TestLogCRNEONConstsLayout(t *testing.T) {
+	var k logCRNEONConsts
+	for _, c := range []struct {
+		name      string
+		got, want uintptr
+	}{
+		{"one", unsafe.Offsetof(k.one), 0}, {"w", unsafe.Offsetof(k.w), 80},
+		{"mmask", unsafe.Offsetof(k.mmask), 96}, {"oneBits", unsafe.Offsetof(k.oneBits), 176},
+		{"expOr", unsafe.Offsetof(k.expOr), 192}, {"tr", unsafe.Offsetof(k.tr), 208},
+		{"tl", unsafe.Offsetof(k.tl), 728},
+	} {
+		if c.got != c.want {
+			t.Errorf("%s at %d, the kernel reads %d", c.name, c.got, c.want)
+		}
+	}
+}
+
+// TestLog32NEONMatchesGo holds log32Fast (NEON blocks, Go for the blocks it
+// hands back and the tail) to logCR bit for bit (NaN to NaN): random bit
+// patterns (most blocks hand back), random positive normals (few do), 1.0
+// runs, the special values, every length 0..70, and in place.
+func TestLog32NEONMatchesGo(t *testing.T) {
+	r := rand.New(rand.NewSource(8))
+	src := make([]float32, 0, 300000)
+	for _, v := range []float32{1, 1, 1, 1, 1, 1, 1, 1, 0, float32(math.Copysign(0, -1)), -1,
+		float32(math.Inf(1)), float32(math.Inf(-1)), float32(math.NaN()), math.SmallestNonzeroFloat32,
+		0x1p-126, math.MaxFloat32, 0.99999994, 1.0000001, 2, 0.5} {
+		src = append(src, v)
+	}
+	for len(src) < 100000 {
+		src = append(src, math.Float32frombits(r.Uint32()))
+	}
+	for len(src) < cap(src) {
+		src = append(src, math.Float32frombits(0x00800000+uint32(r.Int63n(0x7f000000))))
+	}
+	want := make([]float32, len(src))
+	logCRSlice(want, src)
+	same := func(a, b float32) bool { return math.Float32bits(a) == math.Float32bits(b) || (a != a && b != b) }
+	got := make([]float32, len(src))
+	log32Fast(got, src)
+	for i := range src {
+		if !same(got[i], want[i]) {
+			t.Fatalf("log(0x%08x): NEON 0x%08x, Go 0x%08x", math.Float32bits(src[i]), math.Float32bits(got[i]), math.Float32bits(want[i]))
+		}
+	}
+	for n := 0; n <= 70; n++ {
+		d := make([]float32, n)
+		log32Fast(d, src[:n])
+		for i := range d {
+			if !same(d[i], want[i]) {
+				t.Fatalf("length %d [%d]: 0x%08x, want 0x%08x", n, i, math.Float32bits(d[i]), math.Float32bits(want[i]))
+			}
+		}
+	}
+	in := append([]float32(nil), src[:4099]...)
+	log32Fast(in, in)
+	for i := range in {
+		if !same(in[i], want[i]) {
+			t.Fatalf("in place [%d]: 0x%08x, want 0x%08x", i, math.Float32bits(in[i]), math.Float32bits(want[i]))
+		}
+	}
+}
diff --git a/internal/kernels/sum_arm64.s b/internal/kernels/sum_arm64.s
index 1075603..1b68e21 100644
--- a/internal/kernels/sum_arm64.s
+++ b/internal/kernels/sum_arm64.s
@@ -2021,3 +2021,305 @@ TEXT ·blocks8NEON(SB), NOSPLIT, $0-24
 	VST1 [V0.S4, V1.S4], (R0)
 	RET
 
+TEXT ·exp32CRNEON(SB), NOSPLIT, $0-40
+	MOVD dst+0(FP), R0
+	MOVD src+8(FP), R1
+	MOVD n+16(FP), R2
+	MOVD k+24(FP), R8
+	MOVD R8, R10
+	VLD1.P 64(R10), [V16.D2, V17.D2, V18.D2, V19.D2]
+	VLD1.P 64(R10), [V20.D2, V21.D2, V22.D2, V23.D2]
+	VLD1 (R10), [V24.D2, V25.D2, V26.D2]
+	ADD $176, R8, R9
+	MOVD ZR, R3
+	eloop:
+	VLD1 (R1), [V0.S4, V1.S4]
+	WORD $0x6E20EF4C // FACGE V12.4S, V26.4S, V0.4S
+	WORD $0x6E21EF4D // FACGE V13.4S, V26.4S, V1.4S
+	VAND V13.B16, V12.B16, V12.B16
+	VMOV V12.D[0], R6
+	VMOV V12.D[1], R7
+	AND R6, R7, R7
+	CMN $1, R7
+	BNE ebail
+	WORD $0x0E617802 // FCVTL V2.2D, V0.2S
+	WORD $0x4E617809 // FCVTL2 V9.2D, V0.4S
+	WORD $0x0E61782A // FCVTL V10.2D, V1.2S
+	WORD $0x4E61782B // FCVTL2 V11.2D, V1.4S
+	VMOV V16.B16, V3.B16
+	VFMLA V17.D2, V2.D2, V3.D2
+	VFSUB V16.D2, V3.D2, V4.D2
+	VFMLS V18.D2, V4.D2, V2.D2
+	VFMLS V19.D2, V4.D2, V2.D2
+	VFMUL V2.D2, V2.D2, V5.D2
+	VMOV V20.B16, V6.B16
+	VFMLA V21.D2, V2.D2, V6.D2
+	VMOV V22.B16, V7.B16
+	VFMLA V23.D2, V2.D2, V7.D2
+	VFMLA V7.D2, V5.D2, V6.D2
+	VFMLA V6.D2, V5.D2, V2.D2
+	VAND V24.B16, V3.B16, V14.B16
+	VMOV V14.D[0], R4
+	VMOV V14.D[1], R5
+	MOVD (R9)(R4<<3), R4
+	MOVD (R9)(R5<<3), R5
+	VMOV R4, V15.D[0]
+	VMOV R5, V15.D[1]
+	VSHL $46, V3.D2, V3.D2
+	VAND V25.B16, V3.B16, V3.B16
+	VORR V15.B16, V3.B16, V3.B16
+	VFMLA V2.D2, V3.D2, V3.D2
+	VMOV V16.B16, V27.B16
+	VFMLA V17.D2, V9.D2, V27.D2
+	VFSUB V16.D2, V27.D2, V28.D2
+	VFMLS V18.D2, V28.D2, V9.D2
+	VFMLS V19.D2, V28.D2, V9.D2
+	VFMUL V9.D2, V9.D2, V29.D2
+	VMOV V20.B16, V30.B16
+	VFMLA V21.D2, V9.D2, V30.D2
+	VMOV V22.B16, V31.B16
+	VFMLA V23.D2, V9.D2, V31.D2
+	VFMLA V31.D2, V29.D2, V30.D2
+	VFMLA V30.D2, V29.D2, V9.D2
+	VAND V24.B16, V27.B16, V12.B16
+	VMOV V12.D[0], R6
+	VMOV V12.D[1], R7
+	MOVD (R9)(R6<<3), R6
+	MOVD (R9)(R7<<3), R7
+	VMOV R6, V13.D[0]
+	VMOV R7, V13.D[1]
+	VSHL $46, V27.D2, V27.D2
+	VAND V25.B16, V27.B16, V27.B16
+	VORR V13.B16, V27.B16, V27.B16
+	VFMLA V9.D2, V27.D2, V27.D2
+	WORD $0x0E616860 // FCVTN V0.2S, V3.2D
+	WORD $0x4E616B60 // FCVTN2 V0.4S, V27.2D
+	VMOV V16.B16, V3.B16
+	VFMLA V17.D2, V10.D2, V3.D2
+	VFSUB V16.D2, V3.D2, V4.D2
+	VFMLS V18.D2, V4.D2, V10.D2
+	VFMLS V19.D2, V4.D2, V10.D2
+	VFMUL V10.D2, V10.D2, V5.D2
+	VMOV V20.B16, V6.B16
+	VFMLA V21.D2, V10.D2, V6.D2
+	VMOV V22.B16, V7.B16
+	VFMLA V23.D2, V10.D2, V7.D2
+	VFMLA V7.D2, V5.D2, V6.D2
+	VFMLA V6.D2, V5.D2, V10.D2
+	VAND V24.B16, V3.B16, V14.B16
+	VMOV V14.D[0], R4
+	VMOV V14.D[1], R5
+	MOVD (R9)(R4<<3), R4
+	MOVD (R9)(R5<<3), R5
+	VMOV R4, V15.D[0]
+	VMOV R5, V15.D[1]
+	VSHL $46, V3.D2, V3.D2
+	VAND V25.B16, V3.B16, V3.B16
+	VORR V15.B16, V3.B16, V3.B16
+	VFMLA V10.D2, V3.D2, V3.D2
+	VMOV V16.B16, V27.B16
+	VFMLA V17.D2, V11.D2, V27.D2
+	VFSUB V16.D2, V27.D2, V28.D2
+	VFMLS V18.D2, V28.D2, V11.D2
+	VFMLS V19.D2, V28.D2, V11.D2
+	VFMUL V11.D2, V11.D2, V29.D2
+	VMOV V20.B16, V30.B16
+	VFMLA V21.D2, V11.D2, V30.D2
+	VMOV V22.B16, V31.B16
+	VFMLA V23.D2, V11.D2, V31.D2
+	VFMLA V31.D2, V29.D2, V30.D2
+	VFMLA V30.D2, V29.D2, V11.D2
+	VAND V24.B16, V27.B16, V12.B16
+	VMOV V12.D[0], R6
+	VMOV V12.D[1], R7
+	MOVD (R9)(R6<<3), R6
+	MOVD (R9)(R7<<3), R7
+	VMOV R6, V13.D[0]
+	VMOV R7, V13.D[1]
+	VSHL $46, V27.D2, V27.D2
+	VAND V25.B16, V27.B16, V27.B16
+	VORR V13.B16, V27.B16, V27.B16
+	VFMLA V11.D2, V27.D2, V27.D2
+	WORD $0x0E616861 // FCVTN V1.2S, V3.2D
+	WORD $0x4E616B61 // FCVTN2 V1.4S, V27.2D
+	VST1 [V0.S4, V1.S4], (R0)
+	ADD $32, R0
+	ADD $32, R1
+	ADD $8, R3
+	CMP R2, R3
+	BLT eloop
+	ebail:
+	MOVD R3, done+32(FP)
+	RET
+
+TEXT ·log32CRNEON(SB), NOSPLIT, $0-40
+	MOVD dst+0(FP), R0
+	MOVD src+8(FP), R1
+	MOVD n+16(FP), R2
+	MOVD k+24(FP), R8
+	MOVD R8, R13
+	VLD1.P 64(R13), [V16.D2, V17.D2, V18.D2, V19.D2]
+	VLD1.P 64(R13), [V20.D2, V21.D2, V22.D2, V23.D2]
+	VLD1.P 64(R13), [V24.D2, V25.D2, V26.D2, V27.D2]
+	VLD1 (R13), [V28.D2]
+	ADD $208, R8, R9
+	ADD $728, R8, R10
+	MOVD ZR, R3
+	lloop:
+	VLD1 (R1), [V0.S4, V1.S4]
+	VAND V22.B16, V0.B16, V2.B16
+	VADD V23.S4, V2.S4, V3.S4
+	VUSHR $17, V3.S4, V3.S4
+	VUSHR $23, V0.S4, V4.S4
+	VSUB V24.S4, V4.S4, V4.S4
+	VCMEQ V27.S4, V0.S4, V5.S4
+	VADD V25.S4, V0.S4, V6.S4
+	WORD $0x4EBA34C6 // CMGT V6.4S, V6.4S, V26.4S
+	WORD $0x2F20A44C // UXTL V12.2D, V2.S
+	VORR V28.B16, V12.B16, V12.B16
+	VSHL $29, V12.D2, V12.D2
+	WORD $0x0F20A48D // SXTL V13.2D, V4.S
+	WORD $0x4E61D9AD // SCVTF V13.2D, V13.2D
+	VMOV V3.S[0], R4
+	VMOV V3.S[1], R5
+	MOVD (R9)(R4<<3), R6
+	MOVD (R9)(R5<<3), R7
+	MOVD (R10)(R4<<3), R11
+	MOVD (R10)(R5<<3), R12
+	VMOV R6, V14.D[0]
+	VMOV R7, V14.D[1]
+	VMOV R11, V15.D[0]
+	VMOV R12, V15.D[1]
+	VFMUL V14.D2, V12.D2, V14.D2
+	VFSUB V16.D2, V14.D2, V14.D2
+	VFMUL V14.D2, V14.D2, V29.D2
+	VFMUL V17.D2, V13.D2, V31.D2
+	VFADD V15.D2, V31.D2, V31.D2
+	VFMUL V18.D2, V14.D2, V12.D2
+	VFADD V12.D2, V31.D2, V31.D2
+	VFMUL V20.D2, V14.D2, V12.D2
+	VFADD V19.D2, V12.D2, V12.D2
+	VFMUL V12.D2, V29.D2, V12.D2
+	VFADD V12.D2, V31.D2, V31.D2
+	VFADD V21.D2, V31.D2, V30.D2
+	WORD $0x0E616BE7 // FCVTN V7.S, V31.2D
+	WORD $0x0E616BC8 // FCVTN V8.S, V30.2D
+	WORD $0x6F20A44C // UXTL V12.2D, V2.S
+	VORR V28.B16, V12.B16, V12.B16
+	VSHL $29, V12.D2, V12.D2
+	WORD $0x4F20A48D // SXTL V13.2D, V4.S
+	WORD $0x4E61D9AD // SCVTF V13.2D, V13.2D
+	VMOV V3.S[2], R4
+	VMOV V3.S[3], R5
+	MOVD (R9)(R4<<3), R6
+	MOVD (R9)(R5<<3), R7
+	MOVD (R10)(R4<<3), R11
+	MOVD (R10)(R5<<3), R12
+	VMOV R6, V14.D[0]
+	VMOV R7, V14.D[1]
+	VMOV R11, V15.D[0]
+	VMOV R12, V15.D[1]
+	VFMUL V14.D2, V12.D2, V14.D2
+	VFSUB V16.D2, V14.D2, V14.D2
+	VFMUL V14.D2, V14.D2, V29.D2
+	VFMUL V17.D2, V13.D2, V31.D2
+	VFADD V15.D2, V31.D2, V31.D2
+	VFMUL V18.D2, V14.D2, V12.D2
+	VFADD V12.D2, V31.D2, V31.D2
+	VFMUL V20.D2, V14.D2, V12.D2
+	VFADD V19.D2, V12.D2, V12.D2
+	VFMUL V12.D2, V29.D2, V12.D2
+	VFADD V12.D2, V31.D2, V31.D2
+	VFADD V21.D2, V31.D2, V30.D2
+	WORD $0x4E616BE7 // FCVTN V7.S, V31.2D
+	WORD $0x4E616BC8 // FCVTN V8.S, V30.2D
+	VCMEQ V8.S4, V7.S4, V11.S4
+	VORR V5.B16, V11.B16, V11.B16
+	WORD $0x4E661D6B // BIC V11.16B, V11.16B, V6.16B
+	WORD $0x4E651CE9 // BIC V9.16B, V7.16B, V5.16B
+	VAND V22.B16, V1.B16, V2.B16
+	VADD V23.S4, V2.S4, V3.S4
+	VUSHR $17, V3.S4, V3.S4
+	VUSHR $23, V1.S4, V4.S4
+	VSUB V24.S4, V4.S4, V4.S4
+	VCMEQ V27.S4, V1.S4, V5.S4
+	VADD V25.S4, V1.S4, V6.S4
+	WORD $0x4EBA34C6 // CMGT V6.4S, V6.4S, V26.4S
+	WORD $0x2F20A44C // UXTL V12.2D, V2.S
+	VORR V28.B16, V12.B16, V12.B16
+	VSHL $29, V12.D2, V12.D2
+	WORD $0x0F20A48D // SXTL V13.2D, V4.S
+	WORD $0x4E61D9AD // SCVTF V13.2D, V13.2D
+	VMOV V3.S[0], R4
+	VMOV V3.S[1], R5
+	MOVD (R9)(R4<<3), R6
+	MOVD (R9)(R5<<3), R7
+	MOVD (R10)(R4<<3), R11
+	MOVD (R10)(R5<<3), R12
+	VMOV R6, V14.D[0]
+	VMOV R7, V14.D[1]
+	VMOV R11, V15.D[0]
+	VMOV R12, V15.D[1]
+	VFMUL V14.D2, V12.D2, V14.D2
+	VFSUB V16.D2, V14.D2, V14.D2
+	VFMUL V14.D2, V14.D2, V29.D2
+	VFMUL V17.D2, V13.D2, V31.D2
+	VFADD V15.D2, V31.D2, V31.D2
+	VFMUL V18.D2, V14.D2, V12.D2
+	VFADD V12.D2, V31.D2, V31.D2
+	VFMUL V20.D2, V14.D2, V12.D2
+	VFADD V19.D2, V12.D2, V12.D2
+	VFMUL V12.D2, V29.D2, V12.D2
+	VFADD V12.D2, V31.D2, V31.D2
+	VFADD V21.D2, V31.D2, V30.D2
+	WORD $0x0E616BE7 // FCVTN V7.S, V31.2D
+	WORD $0x0E616BC8 // FCVTN V8.S, V30.2D
+	WORD $0x6F20A44C // UXTL V12.2D, V2.S
+	VORR V28.B16, V12.B16, V12.B16
+	VSHL $29, V12.D2, V12.D2
+	WORD $0x4F20A48D // SXTL V13.2D, V4.S
+	WORD $0x4E61D9AD // SCVTF V13.2D, V13.2D
+	VMOV V3.S[2], R4
+	VMOV V3.S[3], R5
+	MOVD (R9)(R4<<3), R6
+	MOVD (R9)(R5<<3), R7
+	MOVD (R10)(R4<<3), R11
+	MOVD (R10)(R5<<3), R12
+	VMOV R6, V14.D[0]
+	VMOV R7, V14.D[1]
+	VMOV R11, V15.D[0]
+	VMOV R12, V15.D[1]
+	VFMUL V14.D2, V12.D2, V14.D2
+	VFSUB V16.D2, V14.D2, V14.D2
+	VFMUL V14.D2, V14.D2, V29.D2
+	VFMUL V17.D2, V13.D2, V31.D2
+	VFADD V15.D2, V31.D2, V31.D2
+	VFMUL V18.D2, V14.D2, V12.D2
+	VFADD V12.D2, V31.D2, V31.D2
+	VFMUL V20.D2, V14.D2, V12.D2
+	VFADD V19.D2, V12.D2, V12.D2
+	VFMUL V12.D2, V29.D2, V12.D2
+	VFADD V12.D2, V31.D2, V31.D2
+	VFADD V21.D2, V31.D2, V30.D2
+	WORD $0x4E616BE7 // FCVTN V7.S, V31.2D
+	WORD $0x4E616BC8 // FCVTN V8.S, V30.2D
+	VCMEQ V8.S4, V7.S4, V8.S4
+	VORR V5.B16, V8.B16, V8.B16
+	WORD $0x4E661D08 // BIC V8.16B, V8.16B, V6.16B
+	WORD $0x4E651CEA // BIC V10.16B, V7.16B, V5.16B
+	VAND V8.B16, V11.B16, V11.B16
+	VMOV V11.D[0], R6
+	VMOV V11.D[1], R7
+	AND R6, R7, R7
+	CMN $1, R7
+	BNE lbail
+	VST1 [V9.S4, V10.S4], (R0)
+	ADD $32, R0
+	ADD $32, R1
+	ADD $8, R3
+	CMP R2, R3
+	BLT lloop
+	lbail:
+	MOVD R3, done+32(FP)
+	RET
+