diff --git a/docs/perf.md b/docs/perf.md
index b213f8c..07039e2 100644
--- a/docs/perf.md
+++ b/docs/perf.md
@@ -520,6 +520,33 @@ held to the n·ε·Σ|aᵢbᵢ| bound and exact on integer data.
 NumPy's multi-threaded dot still wins on this part (24 µs): OpenMP pins its
 threads, so each one rereads its own chunk from its own L3 slice on every
 repetition, while goroutines move between core complexes.
+## Log and Log10 (2026-10-04)
+
+`Log` went through `Map(math.Log)`: 11.6 ns per element on a Zen 3 core
+against NumPy's 4.2. It is now a port of the double-precision log of Arm's
+optimized-routines (the companion of `Exp`, also glibc's since 2.28; MIT
+licence): x = 2^k·z, z split into 128 subintervals, log(x) = k·ln2 + log(c) +
+log1p(z/c − 1) with 1/c and log(c) from a table and a degree-6 polynomial;
+inputs within 2⁻⁴ of 1 take a degree-12 polynomial. z/c − 1 is one
+`math.FMA` (exact on every target), so the second table the non-FMA variant
+needs is not used. `Log10` is log(x)·(1/ln10), the formula of `math.Log10`, on
+this log.
+
+- **Accuracy:** worst 0.508 ULP over 10 000 inputs (whole range, near 1,
+  subnormals), against a 300-bit reference (log m = 2·atanh((m−1)/(m+1)));
+  `math.Log` measured 0.727 on arm64. `TestLogTable` checks each of the 128
+  entries against the properties `log_data.c` documents (1/invc inside its
+  subinterval, logc = log(c) within the rounding of invc, 0x1.8p9 + logc
+  exact).
+- **A second Go bug it removes:** on amd64, Go 1.26's `math.Log` is wrong for
+  subnormal inputs. `math.Log(5e-324)` returns −709.09 instead of −744.44, and
+  `math.Log(1e-310)` −709.09 instead of −713.80 (NumPy: −744.44, −713.80); the
+  error measured 7·10¹³ ULP. Its assembly does not normalise them. Confirmed on
+  a real Zen 3. `Log` and `Log10` inherited it; `TestLogSubnormal` pins the fix.
+  (`Log2` already normalised through `math.Frexp` and was right.)
+
+Speed, 16 Ki elements, one core: M4 37 µs (`math.Log` 72), Zen 3 73 µs
+(`math.Log` 175, NumPy 69).
 
 ## SIMD coverage
 
diff --git a/internal/kernels/log.go b/internal/kernels/log.go
new file mode 100644
index 0000000..dcb4509
--- /dev/null
+++ b/internal/kernels/log.go
@@ -0,0 +1,296 @@
+package kernels
+
+import "math"
+
+// Log is a port of the double-precision log of Arm's optimized-routines
+// (math/log.c and math/log_data.c, by Szabolcs Nagy; also glibc's log since
+// 2.28), used under its MIT licence:
+//
+//	Copyright (c) 2018-2025, Arm Limited.
+//	SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+//
+// x = 2^k z with z in [0x1.6p-1, 0x1.6p0), split into N = 128 subintervals;
+// log(x) = k*ln2 + log(c) + log1p(z/c - 1) for c near the centre of z's
+// subinterval, with 1/c and log(c) from a table and log1p from a degree-6
+// polynomial. Inputs within 2^-4 of 1 take a separate degree-12 polynomial.
+// r = z/c - 1 is computed with one fused multiply-add (math.FMA, exact on every
+// target), the "fast fma" variant of log.c, so the second table its non-FMA
+// variant needs is not used. Arm documents 0.519 ULP (0.507 near 1);
+// TestLogULP measures this port against a 300-bit reference.
+
+const (
+	logTableBits = 7
+	logN         = 1 << logTableBits
+	logOff       = 0x3fe6000000000000
+	logLn2hi     = 0x1.62e42fefa3800p-1
+	logLn2lo     = 0x1.ef35793c76730p-45
+	logA0        = -0x1.0000000000001p-1
+	logA1        = 0x1.555555551305bp-2
+	logA2        = -0x1.fffffffeb459p-3
+	logA3        = 0x1.999b324f10111p-3
+	logA4        = -0x1.55575e506c89fp-3
+	logB0        = -0x1p-1
+	logB1        = 0x1.5555555555577p-2
+	logB2        = -0x1.ffffffffffdcbp-3
+	logB3        = 0x1.999999995dd0cp-3
+	logB4        = -0x1.55555556745a7p-3
+	logB5        = 0x1.24924a344de3p-3
+	logB6        = -0x1.fffffa4423d65p-4
+	logB7        = 0x1.c7184282ad6cap-4
+	logB8        = -0x1.999eb43b068ffp-4
+	logB9        = 0x1.78182f7afd085p-4
+	logB10       = -0x1.5521375d145cdp-4
+)
+
+// logTab holds {1/c, log(c)} per subinterval, copied from log_data.c
+// (N == 128); TestLogTable checks every entry against its defining properties.
+var logTab = [logN]struct{ invc, logc float64 }{
+	{0x1.734f0c3e0de9fp+0, -0x1.7cc7f79e69000p-2},
+	{0x1.713786a2ce91fp+0, -0x1.76feec20d0000p-2},
+	{0x1.6f26008fab5a0p+0, -0x1.713e31351e000p-2},
+	{0x1.6d1a61f138c7dp+0, -0x1.6b85b38287800p-2},
+	{0x1.6b1490bc5b4d1p+0, -0x1.65d5590807800p-2},
+	{0x1.69147332f0cbap+0, -0x1.602d076180000p-2},
+	{0x1.6719f18224223p+0, -0x1.5a8ca86909000p-2},
+	{0x1.6524f99a51ed9p+0, -0x1.54f4356035000p-2},
+	{0x1.63356aa8f24c4p+0, -0x1.4f637c36b4000p-2},
+	{0x1.614b36b9ddc14p+0, -0x1.49da7fda85000p-2},
+	{0x1.5f66452c65c4cp+0, -0x1.445923989a800p-2},
+	{0x1.5d867b5912c4fp+0, -0x1.3edf439b0b800p-2},
+	{0x1.5babccb5b90dep+0, -0x1.396ce448f7000p-2},
+	{0x1.59d61f2d91a78p+0, -0x1.3401e17bda000p-2},
+	{0x1.5805612465687p+0, -0x1.2e9e2ef468000p-2},
+	{0x1.56397cee76bd3p+0, -0x1.2941b3830e000p-2},
+	{0x1.54725e2a77f93p+0, -0x1.23ec58cda8800p-2},
+	{0x1.52aff42064583p+0, -0x1.1e9e129279000p-2},
+	{0x1.50f22dbb2bddfp+0, -0x1.1956d2b48f800p-2},
+	{0x1.4f38f4734ded7p+0, -0x1.141679ab9f800p-2},
+	{0x1.4d843cfde2840p+0, -0x1.0edd094ef9800p-2},
+	{0x1.4bd3ec078a3c8p+0, -0x1.09aa518db1000p-2},
+	{0x1.4a27fc3e0258ap+0, -0x1.047e65263b800p-2},
+	{0x1.4880524d48434p+0, -0x1.feb224586f000p-3},
+	{0x1.46dce1b192d0bp+0, -0x1.f474a7517b000p-3},
+	{0x1.453d9d3391854p+0, -0x1.ea4443d103000p-3},
+	{0x1.43a2744b4845ap+0, -0x1.e020d44e9b000p-3},
+	{0x1.420b54115f8fbp+0, -0x1.d60a22977f000p-3},
+	{0x1.40782da3ef4b1p+0, -0x1.cc00104959000p-3},
+	{0x1.3ee8f5d57fe8fp+0, -0x1.c202956891000p-3},
+	{0x1.3d5d9a00b4ce9p+0, -0x1.b81178d811000p-3},
+	{0x1.3bd60c010c12bp+0, -0x1.ae2c9ccd3d000p-3},
+	{0x1.3a5242b75dab8p+0, -0x1.a45402e129000p-3},
+	{0x1.38d22cd9fd002p+0, -0x1.9a877681df000p-3},
+	{0x1.3755bc5847a1cp+0, -0x1.90c6d69483000p-3},
+	{0x1.35dce49ad36e2p+0, -0x1.87120a645c000p-3},
+	{0x1.34679984dd440p+0, -0x1.7d68fb4143000p-3},
+	{0x1.32f5cceffcb24p+0, -0x1.73cb83c627000p-3},
+	{0x1.3187775a10d49p+0, -0x1.6a39a9b376000p-3},
+	{0x1.301c8373e3990p+0, -0x1.60b3154b7a000p-3},
+	{0x1.2eb4ebb95f841p+0, -0x1.5737d76243000p-3},
+	{0x1.2d50a0219a9d1p+0, -0x1.4dc7b8fc23000p-3},
+	{0x1.2bef9a8b7fd2ap+0, -0x1.4462c51d20000p-3},
+	{0x1.2a91c7a0c1babp+0, -0x1.3b08abc830000p-3},
+	{0x1.293726014b530p+0, -0x1.31b996b490000p-3},
+	{0x1.27dfa5757a1f5p+0, -0x1.2875490a44000p-3},
+	{0x1.268b39b1d3bbfp+0, -0x1.1f3b9f879a000p-3},
+	{0x1.2539d838ff5bdp+0, -0x1.160c8252ca000p-3},
+	{0x1.23eb7aac9083bp+0, -0x1.0ce7f57f72000p-3},
+	{0x1.22a012ba940b6p+0, -0x1.03cdc49fea000p-3},
+	{0x1.2157996cc4132p+0, -0x1.f57bdbc4b8000p-4},
+	{0x1.201201dd2fc9bp+0, -0x1.e370896404000p-4},
+	{0x1.1ecf4494d480bp+0, -0x1.d17983ef94000p-4},
+	{0x1.1d8f5528f6569p+0, -0x1.bf9674ed8a000p-4},
+	{0x1.1c52311577e7cp+0, -0x1.adc79202f6000p-4},
+	{0x1.1b17c74cb26e9p+0, -0x1.9c0c3e7288000p-4},
+	{0x1.19e010c2c1ab6p+0, -0x1.8a646b372c000p-4},
+	{0x1.18ab07bb670bdp+0, -0x1.78d01b3ac0000p-4},
+	{0x1.1778a25efbcb6p+0, -0x1.674f145380000p-4},
+	{0x1.1648d354c31dap+0, -0x1.55e0e6d878000p-4},
+	{0x1.151b990275fddp+0, -0x1.4485cdea1e000p-4},
+	{0x1.13f0ea432d24cp+0, -0x1.333d94d6aa000p-4},
+	{0x1.12c8b7210f9dap+0, -0x1.22079f8c56000p-4},
+	{0x1.11a3028ecb531p+0, -0x1.10e4698622000p-4},
+	{0x1.107fbda8434afp+0, -0x1.ffa6c6ad20000p-5},
+	{0x1.0f5ee0f4e6bb3p+0, -0x1.dda8d4a774000p-5},
+	{0x1.0e4065d2a9fcep+0, -0x1.bbcece4850000p-5},
+	{0x1.0d244632ca521p+0, -0x1.9a1894012c000p-5},
+	{0x1.0c0a77ce2981ap+0, -0x1.788583302c000p-5},
+	{0x1.0af2f83c636d1p+0, -0x1.5715e67d68000p-5},
+	{0x1.09ddb98a01339p+0, -0x1.35c8a49658000p-5},
+	{0x1.08cabaf52e7dfp+0, -0x1.149e364154000p-5},
+	{0x1.07b9f2f4e28fbp+0, -0x1.e72c082eb8000p-6},
+	{0x1.06ab58c358f19p+0, -0x1.a55f152528000p-6},
+	{0x1.059eea5ecf92cp+0, -0x1.63d62cf818000p-6},
+	{0x1.04949cdd12c90p+0, -0x1.228fb8caa0000p-6},
+	{0x1.038c6c6f0ada9p+0, -0x1.c317b20f90000p-7},
+	{0x1.02865137932a9p+0, -0x1.419355daa0000p-7},
+	{0x1.0182427ea7348p+0, -0x1.81203c2ec0000p-8},
+	{0x1.008040614b195p+0, -0x1.0040979240000p-9},
+	{0x1.fe01ff726fa1ap-1, 0x1.feff384900000p-9},
+	{0x1.fa11cc261ea74p-1, 0x1.7dc41353d0000p-7},
+	{0x1.f6310b081992ep-1, 0x1.3cea3c4c28000p-6},
+	{0x1.f25f63ceeadcdp-1, 0x1.b9fc114890000p-6},
+	{0x1.ee9c8039113e7p-1, 0x1.1b0d8ce110000p-5},
+	{0x1.eae8078cbb1abp-1, 0x1.58a5bd001c000p-5},
+	{0x1.e741aa29d0c9bp-1, 0x1.95c8340d88000p-5},
+	{0x1.e3a91830a99b5p-1, 0x1.d276aef578000p-5},
+	{0x1.e01e009609a56p-1, 0x1.07598e598c000p-4},
+	{0x1.dca01e577bb98p-1, 0x1.253f5e30d2000p-4},
+	{0x1.d92f20b7c9103p-1, 0x1.42edd8b380000p-4},
+	{0x1.d5cac66fb5ccep-1, 0x1.606598757c000p-4},
+	{0x1.d272caa5ede9dp-1, 0x1.7da76356a0000p-4},
+	{0x1.cf26e3e6b2ccdp-1, 0x1.9ab434e1c6000p-4},
+	{0x1.cbe6da2a77902p-1, 0x1.b78c7bb0d6000p-4},
+	{0x1.c8b266d37086dp-1, 0x1.d431332e72000p-4},
+	{0x1.c5894bd5d5804p-1, 0x1.f0a3171de6000p-4},
+	{0x1.c26b533bb9f8cp-1, 0x1.067152b914000p-3},
+	{0x1.bf583eeece73fp-1, 0x1.147858292b000p-3},
+	{0x1.bc4fd75db96c1p-1, 0x1.2266ecdca3000p-3},
+	{0x1.b951e0c864a28p-1, 0x1.303d7a6c55000p-3},
+	{0x1.b65e2c5ef3e2cp-1, 0x1.3dfc33c331000p-3},
+	{0x1.b374867c9888bp-1, 0x1.4ba366b7a8000p-3},
+	{0x1.b094b211d304ap-1, 0x1.5933928d1f000p-3},
+	{0x1.adbe885f2ef7ep-1, 0x1.66acd2418f000p-3},
+	{0x1.aaf1d31603da2p-1, 0x1.740f8ec669000p-3},
+	{0x1.a82e63fd358a7p-1, 0x1.815c0f51af000p-3},
+	{0x1.a5740ef09738bp-1, 0x1.8e92954f68000p-3},
+	{0x1.a2c2a90ab4b27p-1, 0x1.9bb3602f84000p-3},
+	{0x1.a01a01393f2d1p-1, 0x1.a8bed1c2c0000p-3},
+	{0x1.9d79f24db3c1bp-1, 0x1.b5b515c01d000p-3},
+	{0x1.9ae2505c7b190p-1, 0x1.c2967ccbcc000p-3},
+	{0x1.9852ef297ce2fp-1, 0x1.cf635d5486000p-3},
+	{0x1.95cbaeea44b75p-1, 0x1.dc1bd3446c000p-3},
+	{0x1.934c69de74838p-1, 0x1.e8c01b8cfe000p-3},
+	{0x1.90d4f2f6752e6p-1, 0x1.f5509c0179000p-3},
+	{0x1.8e6528effd79dp-1, 0x1.00e6c121fb800p-2},
+	{0x1.8bfce9fcc007cp-1, 0x1.071b80e93d000p-2},
+	{0x1.899c0dabec30ep-1, 0x1.0d46b9e867000p-2},
+	{0x1.87427aa2317fbp-1, 0x1.13687334bd000p-2},
+	{0x1.84f00acb39a08p-1, 0x1.1980d67234800p-2},
+	{0x1.82a49e8653e55p-1, 0x1.1f8ffe0cc8000p-2},
+	{0x1.8060195f40260p-1, 0x1.2595fd7636800p-2},
+	{0x1.7e22563e0a329p-1, 0x1.2b9300914a800p-2},
+	{0x1.7beb377dcb5adp-1, 0x1.3187210436000p-2},
+	{0x1.79baa679725c2p-1, 0x1.377266dec1800p-2},
+	{0x1.77907f2170657p-1, 0x1.3d54ffbaf3000p-2},
+	{0x1.756cadbd6130cp-1, 0x1.432eee32fe000p-2},
+}
+
+// logNear1 is log(x) for x within 2^-4 of 1, where the table path would lose
+// accuracy: a degree-12 polynomial with r*r split into exact halves. The
+// conversions to float64 keep the compiler from fusing the compensated sums.
+func logNear1(x float64) float64 {
+	if x == 1 {
+		return 0
+	}
+	r := x - 1.0
+	r2 := r * r
+	r3 := r * r2
+	y := r3 * (logB1 + r*logB2 + r2*logB3 +
+		r3*(logB4+r*logB5+r2*logB6+r3*(logB7+r*logB8+r2*logB9+r3*logB10)))
+	w := float64(r * 0x1p27)
+	rhi := float64(r+w) - w
+	rlo := r - rhi
+	w = float64(float64(rhi*rhi) * logB0) // exact: rhi has 26 bits, B0 = -0.5
+	hi := r + w
+	lo := r - hi + w
+	lo += float64(float64(logB0*rlo) * float64(rhi+r))
+	y += lo
+	y += hi
+	return y
+}
+
+// log1 is log(x) for one element, every case (log.c's log).
+func log1(x float64) float64 {
+	ix := math.Float64bits(x)
+	if ix-math.Float64bits(1.0-0x1p-4) < math.Float64bits(1.0+0x1.09p-4)-math.Float64bits(1.0-0x1p-4) {
+		return logNear1(x)
+	}
+	if top := uint32(ix >> 48); top-0x0010 >= 0x7ff0-0x0010 {
+		// x < 0x1p-1022, or Inf or NaN.
+		switch {
+		case ix*2 == 0:
+			return math.Inf(-1) // log(±0)
+		case ix == math.Float64bits(math.Inf(1)):
+			return x
+		case top&0x8000 != 0 || top&0x7ff0 == 0x7ff0:
+			return math.NaN() // x < 0, or NaN
+		}
+		// Subnormal: normalise it.
+		ix = math.Float64bits(x * 0x1p52)
+		ix -= 52 << 52
+	}
+	return logReduced(ix)
+}
+
+// logReduced is the table path for a positive normal x given by its bits.
+func logReduced(ix uint64) float64 {
+	tmp := ix - logOff
+	i := (tmp >> (52 - logTableBits)) % logN
+	k := int64(tmp) >> 52
+	iz := ix - tmp&(0xfff<<52)
+	e := &logTab[i]
+	r := math.FMA(math.Float64frombits(iz), e.invc, -1.0)
+	kd := float64(k)
+	w := kd*logLn2hi + e.logc // exact (table property)
+	hi := w + r
+	lo := w - hi + r + kd*logLn2lo
+	r2 := r * r
+	return lo + r2*logA0 + r*r2*(logA1+r*logA2+r2*(logA3+r*logA4)) + hi
+}
+
+// Log writes log(src[i]) into dst[i]. Positive normal inputs away from 1 take
+// the table path straight from the loop; everything else goes through log1.
+// TestLogLoopMatchesLog1 holds the two to the same bits.
+func Log(dst, src []float64) {
+	dst = dst[:len(src)]
+	const lo, span = 0x0010 << 48, (0x7ff0 - 0x0010) << 48
+	near1lo := math.Float64bits(1.0 - 0x1p-4)
+	near1span := math.Float64bits(1.0+0x1.09p-4) - near1lo
+	for i, x := range src {
+		ix := math.Float64bits(x)
+		if ix-lo >= span || ix-near1lo < near1span {
+			dst[i] = log1(x)
+			continue
+		}
+		// logReduced, written out: it does not inline, and a call per
+		// element is a large share of the time (as for Exp).
+		tmp := ix - logOff
+		e := &logTab[(tmp>>(52-logTableBits))%logN]
+		r := math.FMA(math.Float64frombits(ix-tmp&(0xfff<<52)), e.invc, -1.0)
+		kd := float64(int64(tmp) >> 52)
+		w := kd*logLn2hi + e.logc
+		hi := w + r
+		lo := w - hi + r + kd*logLn2lo
+		r2 := r * r
+		dst[i] = lo + r2*logA0 + r*r2*(logA1+r*logA2+r2*(logA3+r*logA4)) + hi
+	}
+}
+
+// Log10 writes log10(src[i]) into dst[i] as log(x)*(1/ln10), the formula of
+// math.Log10, on this package's log.
+func Log10(dst, src []float64) {
+	Log(dst, src)
+	for i := range dst {
+		dst[i] *= 1 / math.Ln10
+	}
+}
+
+// Log10P is Log10, split across cores above ParThreshold.
+func Log10P(dst, src []float64) {
+	n := len(dst)
+	if n < ParThreshold {
+		Log10(dst, src)
+		return
+	}
+	parallelFor(n, numWorkers(n), func(lo, hi int) { Log10(dst[lo:hi], src[lo:hi]) })
+}
+
+// LogP is Log, split across cores above ParThreshold.
+func LogP(dst, src []float64) {
+	n := len(dst)
+	if n < ParThreshold {
+		Log(dst, src)
+		return
+	}
+	parallelFor(n, numWorkers(n), func(lo, hi int) { Log(dst[lo:hi], src[lo:hi]) })
+}
diff --git a/internal/kernels/log_test.go b/internal/kernels/log_test.go
new file mode 100644
index 0000000..063ba7a
--- /dev/null
+++ b/internal/kernels/log_test.go
@@ -0,0 +1,154 @@
+package kernels
+
+import (
+	"math"
+	"math/big"
+	"math/rand"
+	"testing"
+)
+
+const logRefPrec = 300
+
+func bigF(x float64) *big.Float { return new(big.Float).SetPrec(logRefPrec).SetFloat64(x) }
+
+// refAtanhSeries is 2*atanh(t) = 2*sum t^(2n+1)/(2n+1), for |t| <= 1/3.
+func refAtanhSeries(t *big.Float) *big.Float {
+	t2 := new(big.Float).SetPrec(logRefPrec).Mul(t, t)
+	sum := new(big.Float).SetPrec(logRefPrec)
+	pow := new(big.Float).SetPrec(logRefPrec).Set(t)
+	for n := 0; n < 400; n++ {
+		sum.Add(sum, new(big.Float).SetPrec(logRefPrec).Quo(pow, big.NewFloat(float64(2*n+1))))
+		pow.Mul(pow, t2)
+	}
+	return sum.Mul(sum, big.NewFloat(2))
+}
+
+var logRefLn2 = refAtanhSeries(new(big.Float).SetPrec(logRefPrec).Quo(bigF(1), bigF(3)))
+
+// refLog is log(x) for x > 0 to 300 bits, independent of the code under test:
+// x = m * 2^e with m in [1,2), log m = 2*atanh((m-1)/(m+1)).
+func refLog(x float64) *big.Float {
+	m := bigF(x)
+	e := m.MantExp(m) - 1 // x = m' * 2^(e+1), m' in [0.5,1) -> m = 2m' in [1,2)
+	m.SetMantExp(m, 1)
+	t := new(big.Float).SetPrec(logRefPrec).Quo(
+		new(big.Float).SetPrec(logRefPrec).Sub(m, bigF(1)),
+		new(big.Float).SetPrec(logRefPrec).Add(m, bigF(1)))
+	r := refAtanhSeries(t)
+	return r.Add(r, new(big.Float).SetPrec(logRefPrec).Mul(logRefLn2, big.NewFloat(float64(e))))
+}
+
+// logULP is |got - exact| in units of the spacing of doubles at exact.
+func logULP(got float64, exact *big.Float) float64 {
+	if exact.Sign() == 0 {
+		return math.Abs(got) * math.Ldexp(1, 1074)
+	}
+	e := exact.MantExp(nil) - 1
+	d := bigF(got)
+	d.Sub(d, exact)
+	d.SetMantExp(d, -(max(e-52, -1074)))
+	f, _ := d.Float64()
+	return math.Abs(f)
+}
+
+// TestLogTable checks every table entry against the properties log_data.c
+// says it was built for: c = 1/invc lies in its subinterval of [0x1.6p-1,
+// 0x1.6p0), logc is log(c) to within the rounding of invc, and 0x1.8p9 + logc
+// is exact (so k*ln2hi + logc is computed without rounding).
+func TestLogTable(t *testing.T) {
+	for i, e := range logTab {
+		lo := math.Float64frombits(logOff + uint64(i)<<45)
+		hi := math.Float64frombits(logOff + uint64(i+1)<<45)
+		if c := 1 / e.invc; c < lo || c >= hi {
+			t.Errorf("entry %d: 1/invc = %v outside [%v, %v)", i, c, lo, hi)
+		}
+		neg := refLog(e.invc)
+		d := bigF(e.logc)
+		d.Add(d, neg) // logc - log(c) ~= logc + log(invc)
+		if f, _ := d.Float64(); math.Abs(f) > 0x1p-52 {
+			t.Errorf("entry %d: logc %v differs from -log(invc) by %v", i, e.logc, f)
+		}
+		if (0x1.8p9+e.logc)-0x1.8p9 != e.logc {
+			t.Errorf("entry %d: 0x1.8p9 + logc is not exact", i)
+		}
+	}
+}
+
+// TestLogULP measures the worst error over random inputs across the whole
+// positive range (log-uniform exponent), the region near 1 and subnormals.
+func TestLogULP(t *testing.T) {
+	rng := rand.New(rand.NewSource(3))
+	var xs []float64
+	for i := 0; i < 2500; i++ {
+		xs = append(xs,
+			math.Ldexp(1+rng.Float64(), rng.Intn(2046)-1022), // whole normal range
+			1+(rng.Float64()*2-1)*0x1.2p-4,                   // near 1 (both paths)
+			math.Ldexp(rng.Float64(), -1022),                 // subnormal
+			0.5+rng.Float64()*2,                              // the table's own range
+		)
+	}
+	worst, worstX, goWorst := 0.0, 0.0, 0.0
+	for _, x := range xs {
+		if x <= 0 {
+			continue
+		}
+		ref := refLog(x)
+		if e := logULP(log1(x), ref); e > worst {
+			worst, worstX = e, x
+		}
+		goWorst = max(goWorst, logULP(math.Log(x), ref))
+	}
+	t.Logf("log: worst %.4f ULP at x=%v over %d inputs (math.Log: %.4f)", worst, worstX, len(xs), goWorst)
+	if worst > 0.53 {
+		t.Errorf("log error %.4f ULP at x=%v exceeds 0.53", worst, worstX)
+	}
+}
+
+func TestLogSpecial(t *testing.T) {
+	for _, c := range []struct{ x, want float64 }{
+		{1, 0}, {0, math.Inf(-1)}, {math.Copysign(0, -1), math.Inf(-1)},
+		{math.Inf(1), math.Inf(1)}, {math.E, 1}, {2, math.Ln2},
+	} {
+		if got := log1(c.x); got != c.want {
+			t.Errorf("log(%v) = %v, want %v", c.x, got, c.want)
+		}
+	}
+	for _, x := range []float64{-1, math.Inf(-1), math.NaN(), -1e-310} {
+		if got := log1(x); !math.IsNaN(got) {
+			t.Errorf("log(%v) = %v, want NaN", x, got)
+		}
+	}
+	if got := log1(5e-324); got != -744.4400719213812 {
+		t.Errorf("log(5e-324) = %v", got)
+	}
+}
+
+// TestLogLoopMatchesLog1: Log's loop and log1 give the same bits everywhere.
+func TestLogLoopMatchesLog1(t *testing.T) {
+	rng := rand.New(rand.NewSource(4))
+	xs := []float64{0, -0.0, 1, 1 - 0x1p-4, 1 + 0x1.09p-4, 0x1p-1022, 5e-324, -3,
+		math.Inf(1), math.Inf(-1), math.NaN(), math.MaxFloat64}
+	for i := 0; i < 20000; i++ {
+		xs = append(xs, math.Ldexp(rng.Float64()*2-0.5, rng.Intn(2100)-1074), 1+(rng.Float64()-0.5)*0.3)
+	}
+	got := make([]float64, len(xs))
+	Log(got, xs)
+	for i, x := range xs {
+		w := log1(x)
+		if math.Float64bits(got[i]) != math.Float64bits(w) && !(math.IsNaN(w) && math.IsNaN(got[i])) {
+			t.Fatalf("Log(%v) = %v, log1 = %v", x, got[i], w)
+		}
+	}
+	src := randVec(3000, 9)
+	for i := range src {
+		src[i] = math.Abs(src[i]) * 100
+	}
+	want, par := make([]float64, len(src)), make([]float64, len(src))
+	Log(want, src)
+	withThresholds(64, 1<<14, func() { LogP(par, src) })
+	eqSlice(t, par, want)
+	LogP(par[:5], src[:5])
+	Log10(want, src)
+	withThresholds(64, 1<<14, func() { Log10P(par, src) })
+	eqSlice(t, par, want)
+}
diff --git a/ufunc.go b/ufunc.go
index af2608b..daf6de6 100644
--- a/ufunc.go
+++ b/ufunc.go
@@ -61,7 +61,12 @@ func (a *Array) Exp() *Array {
 }
 
 // Log returns the elementwise natural logarithm.
-func (a *Array) Log() *Array { return a.Map(math.Log) }
+//
+// It runs a table-driven port of Arm's optimized-routines log (0.51 ULP worst
+// case, measured), split across cores: about twice as fast as math.Log, and
+// correct for subnormal inputs on amd64, where Go 1.26's math.Log is not
+// (log(5e-324) gives -709.09 instead of -744.44).
+func (a *Array) Log() *Array { return a.unary(kernels.LogP) }
 
 // Log2 returns the elementwise base-2 logarithm, within 2 ULP of numpy's.
 func (a *Array) Log2() *Array { return a.Map(log2) }
@@ -81,8 +86,18 @@ func log2(x float64) float64 {
 	return math.Log(frac)*(1/math.Ln2) + float64(exp)
 }
 
-// Log10 returns the elementwise base-10 logarithm.
-func (a *Array) Log10() *Array { return a.Map(math.Log10) }
+// Log10 returns the elementwise base-10 logarithm, log(x)*(1/ln10) as
+// math.Log10 computes it, on the same log as Log (so subnormals are right).
+func (a *Array) Log10() *Array { return a.unary(kernels.Log10P) }
+
+// unary applies a contiguous-slice kernel elementwise into a new array.
+func (a *Array) unary(kernel func(dst, src []float64)) *Array {
+	src := a.contiguousData()
+	dst := a.alloc(len(src), false)
+	kernel(dst, src)
+	cp := append([]int(nil), a.shape...)
+	return &Array{data: dst, shape: cp, strides: rowMajorStrides(cp), ws: a.ws}
+}
 
 // Sin returns the elementwise sine (radians).
 func (a *Array) Sin() *Array { return a.Map(math.Sin) }
diff --git a/ufunc_test.go b/ufunc_test.go
index 3c73e2b..ca5263b 100644
--- a/ufunc_test.go
+++ b/ufunc_test.go
@@ -116,3 +116,17 @@ func TestExpTopOfRange(t *testing.T) {
 	}
 	wantData(t, v.Exp(), []int{3}, []float64{1.3549863193146328e+308, 1.7928227943945155e+308, 0})
 }
+
+// TestLogSubnormal: Log and Log10 are right for subnormal inputs on every
+// arch. Go 1.26's amd64 math.Log returns -709.09 for both 5e-324 and 1e-310,
+// which Log and Log10 used to inherit; numpy gives -744.44 and -713.80.
+func TestLogSubnormal(t *testing.T) {
+	x := mustArr(t, ok(FromData([]float64{5e-324, 1e-310, 1, math.E}, 4)))
+	wantData(t, x.Log(), []int{4}, []float64{-744.4400719213812, -713.8013788281542, 0, 1})
+	got := x.Log10()
+	for i, want := range []float64{-323.3062153431158, -310, 0, 0.4342944819032518} {
+		if math.Abs(got.data[i]-want) > 2e-16*math.Abs(want) {
+			t.Errorf("Log10[%d] = %v, want %v", i, got.data[i], want)
+		}
+	}
+}