floats v0.10.0: drop the max/min zero fallback and the moved sum/elementwise assembly

go.mod

@@ -4,7 +4,7 @@ go 1.27.2
 
 require (
 	github.com/go-fft/fft v0.28.11
-	github.com/go-simd/floats v0.9.0
+	github.com/go-simd/floats v0.10.0
 )
 
 require golang.org/x/sys v0.49.0 // indirect

go.sum

@@ -2,5 +2,7 @@ github.com/go-fft/fft v0.28.11 h1:ZfT78fSmPjZU2E9PgkSI/0f0E3sesYpl13GHmIFzRGc=
 github.com/go-fft/fft v0.28.11/go.mod h1:GN5V4TWCjV05ECqH1ONfjFM6Z1m2shqCdssxKGqQ41s=
 github.com/go-simd/floats v0.9.0 h1:P/sMlwYWwaqCah/cMcBRfffK1JcUxXZVj76f4AwQRN0=
 github.com/go-simd/floats v0.9.0/go.mod h1:Aq5a6oY5pUqm9Hpi+iV2nrNmaaUb0rwHvUHIerXWSvA=
+github.com/go-simd/floats v0.10.0 h1:NFMzLl2wmc1+jNE9J0W9li0O+OmtW9TE1m9XS27068k=
+github.com/go-simd/floats v0.10.0/go.mod h1:Aq5a6oY5pUqm9Hpi+iV2nrNmaaUb0rwHvUHIerXWSvA=
 golang.org/x/sys v0.49.0 h1:XbzkgYJHdqh/8m2Uu0W/dQv8nktxx4BFHp1M0gROTXA=
 golang.org/x/sys v0.49.0/go.mod h1:hNLxWAXmnKAxqDtdwIYC4bM9oQPEecfsnNMuSxOs3og=

internal/kernels/asmgen/loong64/gen.go

@@ -29,15 +29,14 @@ const (
 
 func main() {
 	f := emit.NewFile("loong64")
-	f.Add(sumKernel())
 	f.Add(dotKernel())
-	f.Add(sqrtKernel())
+	// The float64 sum, sqrt and Add/Sub/Mul/Div kernels live in go-simd/floats
+	// (lanes_loong64_gen.go there) since floats v0.10.0.
 	for _, k := range []struct{ name, vec, sca string }{
 		{"addLASX", "XVADDD", "ADDD"}, {"subLASX", "XVSUBD", "SUBD"},
 		{"mulLASX", "XVMULD", "MULD"}, {"divLASX", "XVDIVD", "DIVD"},
 	} {
-		f.Add(binKernel(k.name, k.vec, k.sca))
-		// "XVADDD" -> "VADDF", "ADDD" -> "ADDF": the LSX single-precision op.
+		//"XVADDD" -> "VADDF", "ADDD" -> "ADDF": the LSX single-precision op.
 		f.Add(bin32Kernel(k.name[:3]+"32LSX", k.vec[1:len(k.vec)-1]+"F", k.sca[:len(k.sca)-1]+"F"))
 	}
 	f.Add(gemmKernel())
@@ -74,32 +73,6 @@ func zero(b *loong64.Builder, n int) {
 	}
 }
 
-// sumKernel builds sumLASX(a *float64, n int) float64: four LASX accumulators
-// (16 doubles per iteration), folded, plus a scalar tail. A regrouping of the
-// sum, held to a tolerance like the other targets' SIMD sums.
-func sumKernel() *emit.Function {
-	sig := loong64.Layout([]string{"a", "n"}, []loong64.Type{loong64.Ptr, loong64.Int64},
-		[]string{"ret"}, []loong64.Type{loong64.Float64})
-	b := loong64.NewFunc("sumLASX", sig, 0)
-	b.LoadArg("a", "R4").LoadArg("n", "R5")
-	zero(b, 4)
-	b.Label("sloop").Raw("MOVV $16, R12").Raw("BLT R5, R12, sfold")
-	for k := 0; k < 4; k++ {
-		b.Raw("XVMOVQ %d(R4), X%d", 32*k, 4+k)
-	}
-	for k := 0; k < 4; k++ {
-		b.Raw("XVADDD X%d, X%d, X%d", 4+k, k, k)
-	}
-	b.Raw("ADDV $128, R4").Raw("ADDV $-16, R5").Raw("JMP sloop")
-	b.Label("sfold")
-	fold(b)
-	b.Label("stail").Raw("BEQ R5, sdone").
-		Raw("MOVD 0(R4), F1").Raw("ADDD F1, F0").
-		Raw("ADDV $8, R4").Raw("ADDV $-1, R5").Raw("JMP stail")
-	b.Label("sdone").StoreRet("F0", "ret").Ret()
-	return b.Func()
-}
-
 // dotKernel builds dotLASX(a, b *float64, n int) float64: four fused
 // multiply-add accumulators (xvfmadd.d), folded, plus a fused scalar tail.
 func dotKernel() *emit.Function {
@@ -126,59 +99,6 @@ func dotKernel() *emit.Function {
 	return b.Func()
 }
 
-// sqrtKernel builds sqrtLASX(dst, src *float64, n int): XVFSQRTD over 16
-// doubles per iteration, SQRTD for the tail. Both are the correctly rounded
-// IEEE square root, so the result is bit-identical to math.Sqrt.
-func sqrtKernel() *emit.Function {
-	sig := loong64.Layout([]string{"dst", "src", "n"}, []loong64.Type{loong64.Ptr, loong64.Ptr, loong64.Int64}, nil, nil)
-	b := loong64.NewFunc("sqrtLASX", sig, 0)
-	b.LoadArg("dst", "R4").LoadArg("src", "R6").LoadArg("n", "R5")
-	b.Label("qloop").Raw("MOVV $16, R12").Raw("BLT R5, R12, qtail")
-	for k := 0; k < 4; k++ {
-		b.Raw("XVMOVQ %d(R6), X%d", 32*k, k)
-	}
-	for k := 0; k < 4; k++ {
-		b.Raw("XVFSQRTD X%d, X%d", k, 4+k)
-	}
-	for k := 0; k < 4; k++ {
-		b.Raw("XVMOVQ X%d, %d(R4)", 4+k, 32*k)
-	}
-	b.Raw("ADDV $128, R4").Raw("ADDV $128, R6").Raw("ADDV $-16, R5").Raw("JMP qloop")
-	b.Label("qtail").Raw("BEQ R5, qdone").
-		Raw("MOVD 0(R6), F1").Raw("SQRTD F1, F1").Raw("MOVD F1, 0(R4)").
-		Raw("ADDV $8, R4").Raw("ADDV $8, R6").Raw("ADDV $-1, R5").Raw("JMP qtail")
-	b.Label("qdone").Ret()
-	return b.Func()
-}
-
-// binKernel builds name(dst, a, b *float64, n int): dst[i] = a[i] OP b[i], 16
-// doubles per iteration with the LASX op, a scalar tail with sca. Each lane is
-// one correctly rounded IEEE operation, bit-identical to the scalar oracle.
-func binKernel(name, vec, sca string) *emit.Function {
-	sig := loong64.Layout([]string{"dst", "a", "b", "n"},
-		[]loong64.Type{loong64.Ptr, loong64.Ptr, loong64.Ptr, loong64.Int64}, nil, nil)
-	b := loong64.NewFunc(name, sig, 0)
-	b.LoadArg("dst", "R4").LoadArg("a", "R6").LoadArg("b", "R7").LoadArg("n", "R5")
-	l := name[:3]
-	b.Label(l+"loop").Raw("MOVV $16, R12").Raw("BLT R5, R12, %stail", l)
-	for k := 0; k < 4; k++ {
-		b.Raw("XVMOVQ %d(R6), X%d", 32*k, k)
-		b.Raw("XVMOVQ %d(R7), X%d", 32*k, 4+k)
-	}
-	for k := 0; k < 4; k++ {
-		b.Raw("%s X%d, X%d, X%d", vec, 4+k, k, 8+k) // X(8+k) = X(k) OP X(4+k)
-	}
-	for k := 0; k < 4; k++ {
-		b.Raw("XVMOVQ X%d, %d(R4)", 8+k, 32*k)
-	}
-	b.Raw("ADDV $128, R4").Raw("ADDV $128, R6").Raw("ADDV $128, R7").Raw("ADDV $-16, R5").Raw("JMP %sloop", l)
-	b.Label(l+"tail").Raw("BEQ R5, %sdone", l).
-		Raw("MOVD 0(R6), F1").Raw("MOVD 0(R7), F2").Raw("%s F2, F1", sca).Raw("MOVD F1, 0(R4)").
-		Raw("ADDV $8, R4").Raw("ADDV $8, R6").Raw("ADDV $8, R7").Raw("ADDV $-1, R5").Raw("JMP %stail", l)
-	b.Label(l + "done").Ret()
-	return b.Func()
-}
-
 // lanes is a 16-byte vector constant: v in every lane of size bytes, each lane
 // little-endian.
 func lanes(size int, v uint64) []byte {

internal/kernels/asmgen/ppc64le/gen.go

@@ -31,13 +31,12 @@ const (
 
 func main() {
 	f := emit.NewFile("ppc64le")
-	f.Add(sumKernel())
 	f.Add(dotKernel())
-	f.Add(sqrtKernel())
+	// The float64 sum, sqrt and Add/Sub/Mul/Div kernels live in go-simd/floats
+	// (lanes_ppc64le_gen.go there) since floats v0.10.0.
 	for _, k := range []struct {
 		name, vec, sca string
 	}{{"addVSX", "add", "FADD"}, {"subVSX", "sub", "FSUB"}, {"mulVSX", "mul", "FMUL"}, {"divVSX", "div", "FDIV"}} {
-		f.Add(binKernel(k.name, k.vec, k.sca))
 		f.Add(bin32Kernel(k.name[:3]+"32VSX", k.vec, k.sca+"S"))
 	}
 	f.Add(gemmKernel())
@@ -85,33 +84,6 @@ func reduceSig() ppc64.Signature {
 		[]string{"ret"}, []ppc64.Type{ppc64.Float64})
 }
 
-// sumKernel builds sumVSX(a *float64, n int) float64: four VSX accumulators
-// (8 doubles per iteration), folded, plus a scalar tail. A regrouping of the
-// sum, held to a tolerance like the other targets' SIMD sums.
-func sumKernel() *emit.Function {
-	b := ppc64.NewFunc("sumVSX", reduceSig(), 0)
-	b.LoadArg("a", "R3").LoadArg("n", "R4")
-	offsets(b, false)
-	for r := 0; r < 4; r++ {
-		b.Raw("XXLXOR VS%d, VS%d, VS%d", r, r, r)
-	}
-	b.Label("sloop").Raw("CMP R4, $8").Raw("BLT sfold")
-	for k := 0; k < 4; k++ {
-		b.Raw("LXVD2X (R3)(%s), VS%d", vidx(k), 4+k)
-	}
-	for k := 0; k < 4; k++ {
-		b.XVADDDP(k, k, 4+k)
-	}
-	b.Raw("ADD $64, R3").Raw("ADD $-8, R4").Raw("BR sloop")
-	b.Label("sfold")
-	fold(b)
-	b.Label("stail").Raw("CMP R4, $0").Raw("BEQ sdone").
-		Raw("FMOVD 0(R3), F1").Raw("FADD F1, F0").
-		Raw("ADD $8, R3").Raw("ADD $-1, R4").Raw("BR stail")
-	b.Label("sdone").StoreRet("F0", "ret").Ret()
-	return b.Func()
-}
-
 // dotKernel builds dotVSX(a, b *float64, n int) float64: four fused
 // multiply-add accumulators (XVMADDADP), folded, plus a scalar tail.
 func dotKernel() *emit.Function {
@@ -141,70 +113,6 @@ func dotKernel() *emit.Function {
 	return b.Func()
 }
 
-// sqrtKernel builds sqrtVSX(dst, src *float64, n int): XVSQRTDP over 8 doubles
-// per iteration, FSQRT for the tail. Both are the correctly rounded IEEE
-// square root, so the result is bit-identical to math.Sqrt.
-func sqrtKernel() *emit.Function {
-	sig := ppc64.Layout([]string{"dst", "src", "n"}, []ppc64.Type{ppc64.Ptr, ppc64.Ptr, ppc64.Int64}, nil, nil)
-	b := ppc64.NewFunc("sqrtVSX", sig, 0)
-	b.LoadArg("dst", "R3").LoadArg("src", "R5").LoadArg("n", "R4")
-	offsets(b, false)
-	b.Label("qloop").Raw("CMP R4, $8").Raw("BLT qtail")
-	for k := 0; k < 4; k++ {
-		b.Raw("LXVD2X (R5)(%s), VS%d", vidx(k), 32+k)
-	}
-	for k := 0; k < 4; k++ {
-		b.XVSQRTDP(36+k, 32+k)
-	}
-	for k := 0; k < 4; k++ {
-		b.Raw("STXVD2X VS%d, (R3)(%s)", 36+k, vidx(k))
-	}
-	b.Raw("ADD $64, R3").Raw("ADD $64, R5").Raw("ADD $-8, R4").Raw("BR qloop")
-	b.Label("qtail").Raw("CMP R4, $0").Raw("BEQ qdone").
-		Raw("FMOVD 0(R5), F1").Raw("FSQRT F1, F1").Raw("FMOVD F1, 0(R3)").
-		Raw("ADD $8, R3").Raw("ADD $8, R5").Raw("ADD $-1, R4").Raw("BR qtail")
-	b.Label("qdone").Ret()
-	return b.Func()
-}
-
-// binKernel builds name(dst, a, b *float64, n int): dst[i] = a[i] OP b[i], 8
-// doubles per iteration with the VSX op and a scalar tail with sca. Each lane
-// is one correctly rounded IEEE operation, bit-identical to the scalar oracle.
-func binKernel(name, op, sca string) *emit.Function {
-	sig := ppc64.Layout([]string{"dst", "a", "b", "n"},
-		[]ppc64.Type{ppc64.Ptr, ppc64.Ptr, ppc64.Ptr, ppc64.Int64}, nil, nil)
-	b := ppc64.NewFunc(name, sig, 0)
-	b.LoadArg("dst", "R3").LoadArg("a", "R5").LoadArg("b", "R6").LoadArg("n", "R4")
-	offsets(b, false)
-	l := name[:3]
-	b.Label(l+"loop").Raw("CMP R4, $8").Raw("BLT %stail", l)
-	for k := 0; k < 4; k++ {
-		b.Raw("LXVD2X (R5)(%s), VS%d", vidx(k), 32+k)
-		b.Raw("LXVD2X (R6)(%s), VS%d", vidx(k), 36+k)
-	}
-	for k := 0; k < 4; k++ {
-		switch op {
-		case "add":
-			b.XVADDDP(40+k, 32+k, 36+k)
-		case "sub":
-			b.XVSUBDP(40+k, 32+k, 36+k)
-		case "mul":
-			b.XVMULDP(40+k, 32+k, 36+k)
-		case "div":
-			b.XVDIVDP(40+k, 32+k, 36+k)
-		}
-	}
-	for k := 0; k < 4; k++ {
-		b.Raw("STXVD2X VS%d, (R3)(%s)", 40+k, vidx(k))
-	}
-	b.Raw("ADD $64, R3").Raw("ADD $64, R5").Raw("ADD $64, R6").Raw("ADD $-8, R4").Raw("BR %sloop", l)
-	b.Label(l+"tail").Raw("CMP R4, $0").Raw("BEQ %sdone", l).
-		Raw("FMOVD 0(R5), F1").Raw("FMOVD 0(R6), F2").Raw("%s F2, F1", sca).Raw("FMOVD F1, 0(R3)").
-		Raw("ADD $8, R3").Raw("ADD $8, R5").Raw("ADD $8, R6").Raw("ADD $-1, R4").Raw("BR %stail", l)
-	b.Label(l + "done").Ret()
-	return b.Func()
-}
-
 // lanes is a 16-byte vector constant: v in every lane of size bytes, each lane
 // little-endian, as LXVD2X/LXVW4X load it on ppc64le (every lane equal, so
 // their lane order does not matter).

internal/kernels/asmgen/riscv64/gen.go

@@ -26,14 +26,13 @@ import (
 
 func main() {
 	f := emit.NewFile("riscv64")
-	f.Add(reduceKernel("sumRVV", false))
 	f.Add(reduceKernel("dotRVV", true))
-	f.Add(unaryKernel("sqrtRVV", "VFSQRTV V0, V16"))
+	// The float64 sum, sqrt and Add/Sub/Mul/Div kernels live in go-simd/floats
+	// (lanes_riscv64_gen.go there) since floats v0.10.0.
 	for _, k := range []struct{ name, op string }{
 		{"addRVV", "VFADDVV"}, {"subRVV", "VFSUBVV"},
 		{"mulRVV", "VFMULVV"}, {"divRVV", "VFDIVVV"},
 	} {
-		f.Add(binKernel(k.name, k.op))
 		f.Add(bin32Kernel(k.name[:3]+"32RVV", k.op))
 	}
 	f.Add(gemmKernel())
@@ -106,40 +105,6 @@ func reduceKernel(name string, dot bool) *emit.Function {
 	return b.Func()
 }
 
-// unaryKernel builds NAME(dst, src *float64, n int) for n >= 1: op reads V0 and
-// writes V16.
-func unaryKernel(name, op string) *emit.Function {
-	sig := riscv64.Layout([]string{"dst", "src", "n"},
-		[]riscv64.Type{riscv64.Ptr, riscv64.Ptr, riscv64.Int64}, nil, nil)
-	b := riscv64.NewFunc(name, sig, 0)
-	b.LoadArg("dst", "X10").LoadArg("src", "X11").LoadArg("n", "X13")
-	b.Raw("loop:").
-		Raw("VSETVLI X13, E64, M8, TA, MA, X14").
-		Raw("VLE64V (X11), V0").
-		Raw(op).
-		Raw("VSE64V V16, (X10)")
-	advance(b, "X14", "X10", "X11")
-	b.Raw("SUB X14, X13, X13").Raw("BNEZ X13, loop").Ret()
-	return b.Func()
-}
-
-// binKernel builds NAME(dst, a, b *float64, n int) for n >= 1: dst = a OP b.
-func binKernel(name, op string) *emit.Function {
-	sig := riscv64.Layout([]string{"dst", "a", "b", "n"},
-		[]riscv64.Type{riscv64.Ptr, riscv64.Ptr, riscv64.Ptr, riscv64.Int64}, nil, nil)
-	b := riscv64.NewFunc(name, sig, 0)
-	b.LoadArg("dst", "X10").LoadArg("a", "X11").LoadArg("b", "X12").LoadArg("n", "X13")
-	b.Raw("loop:").
-		Raw("VSETVLI X13, E64, M8, TA, MA, X14").
-		Raw("VLE64V (X11), V0").
-		Raw("VLE64V (X12), V8").
-		Raw("%s V8, V0, V16", op). // V16 = a OP b
-		Raw("VSE64V V16, (X10)")
-	advance(b, "X14", "X10", "X11", "X12")
-	b.Raw("SUB X14, X13, X13").Raw("BNEZ X13, loop").Ret()
-	return b.Func()
-}
-
 // extremumKernel builds NAME(a *floatN, n int) floatN for n >= 1, the
 // NaN-propagating maximum (or minimum), N = bits:
 //

internal/kernels/asmgen/s390x/gen.go

@@ -27,14 +27,13 @@ import (
 
 func main() {
 	f := emit.NewFile("s390x")
-	f.Add(reduceKernel("sumVX", false))
 	f.Add(reduceKernel("dotVX", true))
-	f.Add(sqrtKernel())
+	// The float64 sum, sqrt and Add/Sub/Mul/Div kernels live in go-simd/floats
+	// (lanes_s390x_gen.go there) since floats v0.10.0.
 	for _, k := range []struct{ name, vec, sca string }{
 		{"addVX", "VFADB", "FADD"}, {"subVX", "VFSDB", "FSUB"},
 		{"mulVX", "VFMDB", "FMUL"}, {"divVX", "VFDDB", "FDIV"},
 	} {
-		f.Add(binKernel(k.name, k.vec, k.sca))
 		f.Add(bin32Kernel(k.name[:3]+"32VXE", k.name[:3], k.sca+"S"))
 	}
 	f.Add(gemmKernel())
@@ -103,51 +102,6 @@ func reduceKernel(name string, dot bool) *emit.Function {
 	return b.Func()
 }
 
-// sqrtKernel builds sqrtVX(dst, src *float64, n int).
-func sqrtKernel() *emit.Function {
-	sig := s390x.Layout([]string{"dst", "src", "n"},
-		[]s390x.Type{s390x.Ptr, s390x.Ptr, s390x.Int64}, nil, nil)
-	b := s390x.NewFunc("sqrtVX", sig, 0)
-	b.LoadArg("dst", "R1").LoadArg("src", "R2").LoadArg("n", "R4")
-	b.Raw("loop:").Raw("CMP R4, $8").Raw("BLT tail").
-		Raw("VLM (R2), V0, V3")
-	for r := 0; r < 4; r++ {
-		b.Raw("VFSQDB V%d, V%d", r, r)
-	}
-	b.Raw("VSTM V0, V3, (R1)").
-		Raw("ADD $64, R1").Raw("ADD $64, R2").Raw("SUB $8, R4").Raw("BR loop").
-		Raw("tail:").Raw("CMP R4, $0").Raw("BEQ done").
-		Raw("FMOVD (R2), F0").Raw("FSQRT F0, F0").Raw("FMOVD F0, (R1)").
-		Raw("ADD $8, R1").Raw("ADD $8, R2").Raw("SUB $1, R4").Raw("BR tail").
-		Raw("done:").Ret()
-	return b.Func()
-}
-
-// binKernel builds NAME(dst, a, b *float64, n int): dst = a OP b. vec computes
-// Vd = Vn OP Vm from "vec Vm, Vn, Vd"; sca computes F0 = F0 OP F1 from
-// "sca F1, F0".
-func binKernel(name, vec, sca string) *emit.Function {
-	sig := s390x.Layout([]string{"dst", "a", "b", "n"},
-		[]s390x.Type{s390x.Ptr, s390x.Ptr, s390x.Ptr, s390x.Int64}, nil, nil)
-	b := s390x.NewFunc(name, sig, 0)
-	b.LoadArg("dst", "R1").LoadArg("a", "R2").LoadArg("b", "R3").LoadArg("n", "R4")
-	b.Raw("loop:").Raw("CMP R4, $8").Raw("BLT tail").
-		Raw("VLM (R2), V0, V3"). // a
-		Raw("VLM (R3), V4, V7")  // b
-	for r := 0; r < 4; r++ {
-		b.Raw("%s V%d, V%d, V%d", vec, 4+r, r, r) // V_r = a OP b
-	}
-	b.Raw("VSTM V0, V3, (R1)").
-		Raw("ADD $64, R1").Raw("ADD $64, R2").Raw("ADD $64, R3").Raw("SUB $8, R4").Raw("BR loop").
-		Raw("tail:").Raw("CMP R4, $0").Raw("BEQ done").
-		Raw("FMOVD (R2), F0").Raw("FMOVD (R3), F1").
-		Raw(sca + " F1, F0").
-		Raw("FMOVD F0, (R1)").
-		Raw("ADD $8, R1").Raw("ADD $8, R2").Raw("ADD $8, R3").Raw("SUB $1, R4").Raw("BR tail").
-		Raw("done:").Ret()
-	return b.Func()
-}
-
 // extremum32Kernel builds NAME(a *float32, n int) float32, the maximum (or
 // minimum) of a, for n >= 16. VFMAXSB/VFMINSB and their scalar W forms come
 // with the vector-enhancements facility 1 (z14), so the Go side gates them on