AVX2 kernel for float64 Log: 1.55x faster than NumPy

README.md

@@ -131,9 +131,9 @@ is wrong on subnormals
   elements (2.4–5.6×), the `*Into` forms of `Add`/`Mul`/`Sqrt` (faster or at
   parity at every size), and `MatMul` against single-threaded OpenBLAS
   (3.7–5.9× from 512²); in a `Workspace`, `x + y` and `sqrt(x*y + x)` match or
-  beat NumPy at every size measured.
-- **At parity:** `SumAxis(0)`, `Inner` against 16-thread OpenBLAS, `Log` on
-  one core.
+  beat NumPy at every size measured; float64 `Log` on one core with AVX2
+  (1.55× at 16 Ki, 1.68× at 1 Mi).
+- **At parity:** `SumAxis(0)`, `Inner` against 16-thread OpenBLAS.
 - **Slower:** `MatMul` against 16-thread OpenBLAS at 512² and up (0.42× at
   512², 0.51× at 1024², measured in one run with OpenBLAS; at 256² ours is
   1.44× faster); `Dot` and mat·vec against NumPy's 16-thread BLAS (0.69× and 0.52× since v0.6.0's helper pool, measured in one run; 0.16× before); `Exp`

SECURITY.md

@@ -66,6 +66,15 @@ standard library's known vulnerabilities; `govulncheck ./...` reports none.
 
 ## Audit log
 
+**2026-11-07 (v0.36.0).** The AVX2 float64 log kernel (`logAVX2`) reads and
+writes at most `n` doubles, a positive multiple of 4 chosen by the Go
+wrapper, and returns how many it wrote; the wrapper finishes the block it
+stopped at in Go and resumes past it, so the loop always advances, and a
+block it stops at is not written (`dst == src` is safe). Its gather index is
+`2·((ix - off) >> 45 & 127)` into the 128-row table, masked for any input,
+so it cannot leave the table. It is gated on AVX2 and FMA (CPU and OS) and
+fenced with PROT_NONE on both sides of each operand.
+
 **2026-11-07 (v0.35.0).** `Workspace.alloc32` turns a slice of `(n+1)/2`
 float64 from the arena into `n` float32 with `unsafe.Slice`: the float32
 slice lies inside the float64 one (2·⌈n/2⌉·4 ≥ 4n bytes), its cap is n so an

docs/perf.md

@@ -1600,10 +1600,42 @@ NumPy). In a Workspace, float32 `Exp` is about 12% faster than NumPy and
 `Log` about 5% slower; NumPy computes both in float32, at 2.52 and 3.83 ULP.
 Allocating, Go's zeroing of each result adds 4 µs at 16 Ki.
 
-**Open defect:** float64 `Log` is 1.23× slower than NumPy at 16 Ki on Zen 3
-(85 µs against 69 µs), while float64 `Exp` is 16% faster. The float32 kernels
-for the other architectures are still to write; there `Exp` and `Log` run
-the Go transcriptions.
+#### float64 `Log` in AVX2: 1.55× faster than NumPy (2026-11-07)
+
+float64 `Log` was 1.23× slower than NumPy at 16 Ki on Zen 3 (85 µs against
+69 µs). NumPy's 4.2 ns per element is glibc's scalar log, which is
+optimized-routines' log.c, the code ours ports; the gap was our scalar Go.
+
+**`logAVX2`** (go-asmgen, AVX2 and FMA) runs `Log`'s table path on four
+doubles at a time, with `Log`'s operations in its order: the one FMA it has
+(r = z·invc − 1) as `VFMADD231PD`, every other product and sum as a separate
+`VMULPD`/`VADDPD`, which is how Go compiles them on amd64.
+- `invc` and `logc` are gathered from the 16-byte table rows with
+  `VPGATHERQQ`.
+- AVX2 has no 64-bit arithmetic shift and no int64→float64 conversion, so
+  the exponent kd = int64(tmp)>>52 is built exactly from u = tmp>>52:
+  the float64 with bits 0x4330000000000000 | (u ^ 0x800), minus 2^52 + 2048.
+- A block holding a lane that is not a positive normal, or one within 2⁻⁴
+  of 1 (`Log`'s polynomial path), is not stored; it is handed back to the Go
+  loop.
+
+`TestLogAVX2MatchesGo` holds it to the Go loop bit for bit on 300 000 inputs
+covering every exponent and table index, near 1, the special values, every
+length 0–70 and in place. Dropping the `+ kd·ln2lo` term fails it, while
+`TestLogULP` alone (0.5078 ULP worst, unchanged) would not see it. Fenced;
+the constant block's layout is tested.
+
+| float64 `Log`, one Zen 3 core, load ≤ 3.7 | 16 Ki | 1 Mi |
+|---|---:|---:|
+| v0.35.0 | 83.9 µs | 5.06 ms |
+| **AVX2** | **44.4 µs** | **2.59 ms** |
+| NumPy 2.5.4, `np.log(a)` | 68.9–69.4 µs | 4.33–4.41 ms |
+
+−47% and −49% against v0.35.0 (`abbench`, 8 interleaved rounds, p = 0.000),
+both allocating, on `BenchmarkDTypeExpLog`'s data.
+
+The float32 kernels for the other architectures are still to write; there
+`Exp` and `Log` run the Go transcriptions.
 
 ## SIMD coverage
 

internal/kernels/asmgen/amd64/gen.go

@@ -57,6 +57,7 @@ func main() {
 	f.Add(amd64.FeatureProbe("hasAVX2", amd64.AVX2))
 	f.Add(exp32CRKernel())
 	f.Add(log32CRKernel())
+	f.Add(logKernel())
 	f.Add(blocks8Kernel())
 	f.Add(bin32Kernel("add32SSE2", "ADDPS", "ADDSS"))
 	f.Add(bin32Kernel("sub32SSE2", "SUBPS", "SUBSS"))
@@ -549,6 +550,101 @@ func log32CRKernel() *emit.Function {
 	return b.Func()
 }
 
+// logKernel builds logAVX2(dst, src *float64, n int, k *logConsts, tab
+// *logEntry) (done int): the table path of Log (log.go, the port of
+// optimized-routines' log) on four doubles at a time, n a positive multiple
+// of 4, with Log's operations in its order and its single FMA
+// (r = z*invc - 1, VFMADD231PD) where it has one; every other product and
+// sum is a separate VMULPD/VADDPD, as Go compiles them on amd64.
+//
+//   - tmp = ix - off; i = (tmp>>45) & 127 indexes the 16-byte table entries
+//     (VPGATHERQQ at 2i, scale 8: invc, then logc 8 bytes on); k =
+//     int64(tmp)>>52 has no AVX2 instruction (no 64-bit arithmetic shift, no
+//     int64 to float64 conversion), so kd is built exactly from u = tmp>>52:
+//     float64 bits 0x4330000000000000 | (u ^ 0x800) is 2^52 + (u ^ 0x800),
+//     and subtracting 2^52 + 2048 leaves the sign-extended k.
+//
+// A block holding a lane that is not a positive normal (ix < 0x0010<<48 or
+// above the largest finite, as int64), or one within 2^-4 of 1 (Log's
+// polynomial path), is not stored: the kernel returns the count it wrote and
+// Log finishes that block in Go. Needs AVX2 (gathers, 64-bit compares) and
+// FMA.
+func logKernel() *emit.Function {
+	sig := amd64.Layout(
+		[]string{"dst", "src", "n", "k", "tab"},
+		[]amd64.Type{amd64.Ptr, amd64.Ptr, amd64.Int64, amd64.Ptr, amd64.Ptr},
+		[]string{"done"}, []amd64.Type{amd64.Int64},
+	)
+	const (
+		lo, hiMax, n1loM1, n1hi = 0, 32, 64, 96
+		off, idxMask, topMask   = 128, 160, 192
+		magic, magicSub, x2048  = 224, 256, 288
+		negOne, ln2hi, ln2lo    = 320, 352, 384
+		a0, a1, a2, a3, a4      = 416, 448, 480, 512, 544
+	)
+	b := amd64.NewFunc("logAVX2", sig, 0)
+	b.LoadArg("dst", "DI").LoadArg("src", "SI").LoadArg("n", "CX").LoadArg("k", "R8").LoadArg("tab", "R9")
+	b.Raw("XORQ AX, AX")
+	b.Raw("dloop:")
+	b.Raw("VMOVDQU (SI), Y0") // ix
+	// The lanes Log's loop sends to log1.
+	b.Raw(fmt.Sprintf("VPCMPGTQ %d(R8), Y0, Y1", hiMax)) // ix > largest finite (as int64)
+	b.Raw(fmt.Sprintf("VMOVDQU %d(R8), Y2", lo))
+	b.Raw("VPCMPGTQ Y0, Y2, Y2") // ix < 0x0010<<48: subnormal, zero or negative
+	b.Raw(fmt.Sprintf("VPCMPGTQ %d(R8), Y0, Y3", n1loM1))
+	b.Raw(fmt.Sprintf("VMOVDQU %d(R8), Y4", n1hi))
+	b.Raw("VPCMPGTQ Y0, Y4, Y4")
+	b.Raw("VPAND Y4, Y3, Y3") // within 2^-4 of 1
+	b.Raw("VPOR Y3, Y2, Y2")
+	b.Raw("VPOR Y1, Y2, Y2")
+	b.Raw("VPTEST Y2, Y2")
+	b.Raw("JNZ dbail")
+	b.Raw(fmt.Sprintf("VPSUBQ %d(R8), Y0, Y1", off)) // tmp
+	b.Raw("VPSRLQ $45, Y1, Y2")
+	b.Raw(fmt.Sprintf("VPAND %d(R8), Y2, Y2", idxMask))
+	b.Raw("VPADDQ Y2, Y2, Y2") // 2i: entries are two float64 wide
+	b.Raw("VPSRLQ $52, Y1, Y3")
+	b.Raw(fmt.Sprintf("VPXOR %d(R8), Y3, Y3", x2048))
+	b.Raw(fmt.Sprintf("VPOR %d(R8), Y3, Y3", magic))
+	b.Raw(fmt.Sprintf("VSUBPD %d(R8), Y3, Y3", magicSub)) // kd
+	b.Raw(fmt.Sprintf("VPAND %d(R8), Y1, Y4", topMask))
+	b.Raw("VPSUBQ Y4, Y0, Y4") // iz = ix - tmp&(0xfff<<52)
+	b.Raw("VPCMPEQQ Y5, Y5, Y5")
+	b.Raw("VPGATHERQQ Y5, (R9)(Y2*8), Y6") // invc
+	b.Raw("VPCMPEQQ Y5, Y5, Y5")
+	b.Raw("VPGATHERQQ Y5, 8(R9)(Y2*8), Y7") // logc
+	b.Raw(fmt.Sprintf("VMOVUPD %d(R8), Y8", negOne))
+	b.Raw("VFMADD231PD Y6, Y4, Y8")                    // r = z*invc - 1, one rounding
+	b.Raw(fmt.Sprintf("VMULPD %d(R8), Y3, Y9", ln2hi)) // kd*ln2hi
+	b.Raw("VADDPD Y7, Y9, Y9")                         // w = ... + logc
+	b.Raw("VADDPD Y8, Y9, Y10")                        // hi = w + r
+	b.Raw("VSUBPD Y10, Y9, Y11")                       // w - hi
+	b.Raw("VADDPD Y8, Y11, Y11")                       // + r
+	b.Raw(fmt.Sprintf("VMULPD %d(R8), Y3, Y12", ln2lo))
+	b.Raw("VADDPD Y12, Y11, Y11") // lo = w - hi + r + kd*ln2lo
+	b.Raw("VMULPD Y8, Y8, Y12")   // r2
+	b.Raw(fmt.Sprintf("VMULPD %d(R8), Y8, Y13", a2))
+	b.Raw(fmt.Sprintf("VADDPD %d(R8), Y13, Y13", a1)) // A1 + r*A2
+	b.Raw(fmt.Sprintf("VMULPD %d(R8), Y8, Y14", a4))
+	b.Raw(fmt.Sprintf("VADDPD %d(R8), Y14, Y14", a3)) // A3 + r*A4
+	b.Raw("VMULPD Y14, Y12, Y14")
+	b.Raw("VADDPD Y14, Y13, Y13") // A1 + r*A2 + r2*(A3 + r*A4)
+	b.Raw("VMULPD Y12, Y8, Y14")  // r*r2
+	b.Raw("VMULPD Y13, Y14, Y14") // r*r2*(...)
+	b.Raw(fmt.Sprintf("VMULPD %d(R8), Y12, Y13", a0))
+	b.Raw("VADDPD Y13, Y11, Y11") // lo + r2*A0
+	b.Raw("VADDPD Y14, Y11, Y11") // + r*r2*(...)
+	b.Raw("VADDPD Y10, Y11, Y11") // + hi
+	b.Raw("VMOVUPD Y11, (DI)")
+	b.Raw("ADDQ $32, SI").Raw("ADDQ $32, DI").Raw("ADDQ $4, AX")
+	b.Raw("CMPQ AX, CX").Raw("JLT dloop")
+	b.Raw("dbail:")
+	b.Raw("VZEROUPPER")
+	b.StoreRet("AX", "done")
+	b.Ret()
+	return b.Func()
+}
+
 // bin32Kernel builds NAME(dst, a, b *float32, n int): dst[i] = a[i] OP b[i],
 // 16 floats per iteration over four XMM registers with the packed op, then a
 // scalar tail with sca. Each lane is one IEEE float32 operation, so the result

internal/kernels/log.go

@@ -44,7 +44,11 @@ const (
 
 // 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 }{
+// logEntry is one row of logTab: 1/c and log(c), 16 bytes (logAVX2 gathers
+// invc at 16i and logc at 16i + 8).
+type logEntry struct{ invc, logc float64 }
+
+var logTab = [logN]logEntry{
 	{0x1.734f0c3e0de9fp+0, -0x1.7cc7f79e69000p-2},
 	{0x1.713786a2ce91fp+0, -0x1.76feec20d0000p-2},
 	{0x1.6f26008fab5a0p+0, -0x1.713e31351e000p-2},
@@ -238,10 +242,15 @@ func logReduced(ix uint64) float64 {
 	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.
+// Log writes log(src[i]) into dst[i]: on the AVX2 kernel where there is one
+// (log_amd64.go), on logScalar otherwise and for the blocks the kernel hands
+// back.
+func Log(dst, src []float64) { logFast(dst[:len(src)], src) }
+
+// logScalar is Log in Go. 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) {
+func logScalar(dst, src []float64) {
 	dst = dst[:len(src)]
 	const lo, span = 0x0010 << 48, (0x7ff0 - 0x0010) << 48
 	near1lo := math.Float64bits(1.0 - 0x1p-4)

internal/kernels/log_amd64.go

@@ -0,0 +1,54 @@
+package kernels
+
+import "math"
+
+// logConsts is the constant block logAVX2 reads (asmgen/amd64/gen.go,
+// logKernel), each value as a 4-lane vector at the offset the kernel uses.
+type logConsts struct {
+	lo, hiMax, n1loM1, n1hi [4]uint64  // offsets 0, 32, 64, 96
+	off, idxMask, topMask   [4]uint64  // 128, 160, 192
+	magic                   [4]uint64  // 224
+	magicSub                [4]float64 // 256
+	x2048                   [4]uint64  // 288
+	negOne, ln2hi, ln2lo    [4]float64 // 320, 352, 384
+	a0, a1, a2, a3, a4      [4]float64 // 416, 448, 480, 512, 544
+}
+
+var logK = func() *logConsts {
+	k := new(logConsts)
+	for i := range 4 {
+		k.lo[i], k.hiMax[i] = 0x0010<<48, 0x7fefffffffffffff
+		k.n1loM1[i] = math.Float64bits(1.0-0x1p-4) - 1
+		k.n1hi[i] = math.Float64bits(1.0 + 0x1.09p-4)
+		k.off[i], k.idxMask[i], k.topMask[i] = logOff, logN-1, 0xfff<<52
+		k.magic[i], k.magicSub[i], k.x2048[i] = 0x4330000000000000, 0x1p52+2048, 0x800
+		k.negOne[i], k.ln2hi[i], k.ln2lo[i] = -1, logLn2hi, logLn2lo
+		k.a0[i], k.a1[i], k.a2[i], k.a3[i], k.a4[i] = logA0, logA1, logA2, logA3, logA4
+	}
+	return k
+}()
+
+// useLogAVX2 selects the AVX2 kernel; a var so tests can force the Go path.
+var useLogAVX2 = hasAVX2() && useFMA
+
+// logFast runs blocks of 4 on the AVX2 kernel, which hands back each block
+// holding an input Log does not take through its table path (not a positive
+// normal, or within 2^-4 of 1) to logScalar, as is the tail. Both run the
+// same operations in the same order, so they agree bit for bit.
+func logFast(dst, src []float64) {
+	i := 0
+	if useLogAVX2 {
+		n := len(src) &^ 3
+		for i < n {
+			i += logAVX2(&dst[i], &src[i], n-i, logK, &logTab[0])
+			if i < n {
+				logScalar(dst[i:i+4], src[i:i+4])
+				i += 4
+			}
+		}
+	}
+	logScalar(dst[i:], src[i:])
+}
+
+//go:noescape
+func logAVX2(dst, src *float64, n int, k *logConsts, tab *logEntry) (done int)

internal/kernels/log_amd64_test.go

@@ -0,0 +1,84 @@
+package kernels
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+	"unsafe"
+)
+
+// TestLogConstsLayout holds logConsts to the offsets logKernel uses.
+func TestLogConstsLayout(t *testing.T) {
+	var k logConsts
+	for _, c := range []struct {
+		name      string
+		got, want uintptr
+	}{
+		{"lo", unsafe.Offsetof(k.lo), 0}, {"n1hi", unsafe.Offsetof(k.n1hi), 96},
+		{"off", unsafe.Offsetof(k.off), 128}, {"topMask", unsafe.Offsetof(k.topMask), 192},
+		{"magic", unsafe.Offsetof(k.magic), 224}, {"magicSub", unsafe.Offsetof(k.magicSub), 256},
+		{"x2048", unsafe.Offsetof(k.x2048), 288}, {"negOne", unsafe.Offsetof(k.negOne), 320},
+		{"ln2lo", unsafe.Offsetof(k.ln2lo), 384}, {"a0", unsafe.Offsetof(k.a0), 416},
+		{"a4", unsafe.Offsetof(k.a4), 544},
+	} {
+		if c.got != c.want {
+			t.Errorf("%s at %d, the kernel reads %d", c.name, c.got, c.want)
+		}
+	}
+	if unsafe.Sizeof(logEntry{}) != 16 {
+		t.Errorf("logEntry is %d bytes, the kernel gathers at 16i", unsafe.Sizeof(logEntry{}))
+	}
+}
+
+// TestLogAVX2MatchesGo holds logFast (AVX2 blocks, Go for the blocks it hands
+// back and the tail) to logScalar bit for bit (NaN to NaN): random bit
+// patterns (most blocks hand back), positive normals across the whole
+// exponent range (few do, and every k from -1022 to 1023 and every table
+// index occurs), values near 1, the special values, every length 0..70, and
+// in place.
+func TestLogAVX2MatchesGo(t *testing.T) {
+	if !useLogAVX2 {
+		t.Skip("no AVX2+FMA here")
+	}
+	r := rand.New(rand.NewSource(9))
+	src := make([]float64, 0, 300000)
+	for _, v := range []float64{1, 1, 1, 1, 0, math.Copysign(0, -1), -1, math.Inf(1), math.Inf(-1), math.NaN(),
+		5e-324, 0x1p-1022, math.MaxFloat64, 1 - 0x1p-4, 1 + 0x1.09p-4, math.Nextafter(1-0x1p-4, 0), 0.25, 24.25} {
+		src = append(src, v)
+	}
+	for len(src) < 100000 {
+		src = append(src, math.Float64frombits(r.Uint64()))
+	}
+	for len(src) < 250000 {
+		src = append(src, math.Float64frombits(0x0010<<48+uint64(r.Int63n(0x7fe0<<48))))
+	}
+	for len(src) < cap(src) {
+		src = append(src, 0.9+0.2*r.Float64())
+	}
+	want := make([]float64, len(src))
+	logScalar(want, src)
+	same := func(a, b float64) bool { return math.Float64bits(a) == math.Float64bits(b) || (a != a && b != b) }
+	got := make([]float64, len(src))
+	logFast(got, src)
+	for i := range src {
+		if !same(got[i], want[i]) {
+			t.Fatalf("log(0x%016x): AVX2 0x%016x, Go 0x%016x", math.Float64bits(src[i]), math.Float64bits(got[i]), math.Float64bits(want[i]))
+		}
+	}
+	for n := 0; n <= 70; n++ {
+		d := make([]float64, n)
+		logFast(d, src[:n])
+		for i := range d {
+			if !same(d[i], want[i]) {
+				t.Fatalf("length %d [%d]: 0x%016x, want 0x%016x", n, i, math.Float64bits(d[i]), math.Float64bits(want[i]))
+			}
+		}
+	}
+	in := append([]float64(nil), src[100000:104099]...)
+	logFast(in, in)
+	for i := range in {
+		if !same(in[i], want[100000+i]) {
+			t.Fatalf("in place [%d]: 0x%016x, want 0x%016x", i, math.Float64bits(in[i]), math.Float64bits(want[100000+i]))
+		}
+	}
+}

internal/kernels/log_fence_amd64_test.go

@@ -0,0 +1,20 @@
+//go:build linux || darwin
+
+package kernels
+
+import "testing"
+
+// TestLogAVX2StaysInside fences src and dst of the AVX2 float64 log kernel
+// (fences needs mmap/mprotect, hence the build tag of fence_test.go).
+func TestLogAVX2StaysInside(t *testing.T) {
+	if !useLogAVX2 {
+		t.Skip("no AVX2+FMA here")
+	}
+	for n := 4; n <= 64; n += 4 {
+		for _, s := range fences(t, n) {
+			for _, d := range fences(t, n) {
+				logAVX2(&d[0], &s[0], n, logK, &logTab[0])
+			}
+		}
+	}
+}

internal/kernels/log_other.go

@@ -0,0 +1,6 @@
+//go:build !amd64
+
+package kernels
+
+// logFast is Log's Go loop on these architectures (no SIMD kernel yet).
+func logFast(dst, src []float64) { logScalar(dst, src) }