diff --git a/simplex/axpy_archsimd_amd64.go b/simplex/axpy_archsimd_amd64.go
new file mode 100644
index 0000000..6bd64e5
--- /dev/null
+++ b/simplex/axpy_archsimd_amd64.go
@@ -0,0 +1,34 @@
+//go:build goexperiment.simd && amd64
+
+package simplex
+
+import "simd/archsimd"
+
+// axpy computes dst[k] += factor*src[k], using SSE2 directly (no portable
+// dispatch stub) via simd/archsimd's Float64x2 — part of the amd64 baseline.
+func axpy(dst, src []float64, factor float64) {
+	fv := archsimd.BroadcastFloat64x2(factor)
+	n := len(dst)
+	k := 0
+	for ; k+2 <= n; k += 2 {
+		d := archsimd.LoadFloat64x2(dst[k : k+2])
+		s := archsimd.LoadFloat64x2(src[k : k+2])
+		s.MulAdd(fv, d).Store(dst[k : k+2])
+	}
+	for ; k < n; k++ {
+		dst[k] += factor * src[k]
+	}
+}
+
+// scale computes dst[k] /= divisor for all k, using SSE2 directly.
+func scale(dst []float64, divisor float64) {
+	dv := archsimd.BroadcastFloat64x2(divisor)
+	n := len(dst)
+	k := 0
+	for ; k+2 <= n; k += 2 {
+		archsimd.LoadFloat64x2(dst[k : k+2]).Div(dv).Store(dst[k : k+2])
+	}
+	for ; k < n; k++ {
+		dst[k] /= divisor
+	}
+}
diff --git a/simplex/axpy_archsimd_arm64.go b/simplex/axpy_archsimd_arm64.go
new file mode 100644
index 0000000..cb22ea6
--- /dev/null
+++ b/simplex/axpy_archsimd_arm64.go
@@ -0,0 +1,34 @@
+//go:build goexperiment.simd && arm64
+
+package simplex
+
+import "simd/archsimd"
+
+// axpy computes dst[k] += factor*src[k], using NEON directly (no portable
+// dispatch stub) via simd/archsimd's Float64x2 (128-bit, 2 lanes).
+func axpy(dst, src []float64, factor float64) {
+	fv := archsimd.BroadcastFloat64x2(factor)
+	n := len(dst)
+	k := 0
+	for ; k+2 <= n; k += 2 {
+		d := archsimd.LoadFloat64x2(dst[k : k+2])
+		s := archsimd.LoadFloat64x2(src[k : k+2])
+		s.MulAdd(fv, d).Store(dst[k : k+2])
+	}
+	for ; k < n; k++ {
+		dst[k] += factor * src[k]
+	}
+}
+
+// scale computes dst[k] /= divisor for all k, using NEON directly.
+func scale(dst []float64, divisor float64) {
+	dv := archsimd.BroadcastFloat64x2(divisor)
+	n := len(dst)
+	k := 0
+	for ; k+2 <= n; k += 2 {
+		archsimd.LoadFloat64x2(dst[k : k+2]).Div(dv).Store(dst[k : k+2])
+	}
+	for ; k < n; k++ {
+		dst[k] /= divisor
+	}
+}
diff --git a/simplex/axpy_generic.go b/simplex/axpy_generic.go
new file mode 100644
index 0000000..747e2a8
--- /dev/null
+++ b/simplex/axpy_generic.go
@@ -0,0 +1,18 @@
+//go:build !goexperiment.simd
+
+package simplex
+
+// axpy computes dst[k] += factor*src[k] for all k. Scalar fallback used by
+// default builds (Go 1.27+ with GOEXPERIMENT=simd gets a vectorized version).
+func axpy(dst, src []float64, factor float64) {
+	for k := range dst {
+		dst[k] += factor * src[k]
+	}
+}
+
+// scale computes dst[k] /= divisor for all k.
+func scale(dst []float64, divisor float64) {
+	for k := range dst {
+		dst[k] /= divisor
+	}
+}
diff --git a/simplex/axpy_simd.go b/simplex/axpy_simd.go
new file mode 100644
index 0000000..0a809ba
--- /dev/null
+++ b/simplex/axpy_simd.go
@@ -0,0 +1,37 @@
+//go:build goexperiment.simd && !arm64 && !amd64
+
+package simplex
+
+import "simd"
+
+// axpy computes dst[k] += factor*src[k], vectorized via the portable simd
+// package. arm64 and amd64 use simd/archsimd directly (axpy_archsimd_*.go).
+func axpy(dst, src []float64, factor float64) {
+	fv := simd.BroadcastFloat64s(factor)
+	width := fv.Len()
+	n := len(dst)
+	k := 0
+	for ; k+width <= n; k += width {
+		d := simd.LoadFloat64s(dst[k : k+width])
+		s := simd.LoadFloat64s(src[k : k+width])
+		s.MulAdd(fv, d).Store(dst[k : k+width])
+	}
+	for ; k < n; k++ {
+		dst[k] += factor * src[k]
+	}
+}
+
+// scale computes dst[k] /= divisor for all k, vectorized via the portable
+// simd package.
+func scale(dst []float64, divisor float64) {
+	dv := simd.BroadcastFloat64s(divisor)
+	width := dv.Len()
+	n := len(dst)
+	k := 0
+	for ; k+width <= n; k += width {
+		simd.LoadFloat64s(dst[k : k+width]).Div(dv).Store(dst[k : k+width])
+	}
+	for ; k < n; k++ {
+		dst[k] /= divisor
+	}
+}
diff --git a/simplex/simplex.go b/simplex/simplex.go
index 6314534..b39bff6 100644
--- a/simplex/simplex.go
+++ b/simplex/simplex.go
@@ -238,10 +238,7 @@ func duals(st *State, cost []float64, m int, dst []float64) {
 		if cb == 0 {
 			continue
 		}
-		row := st.binv[i*m : i*m+m]
-		for k := 0; k < m; k++ {
-			dst[k] += cb * row[k]
-		}
+		axpy(dst, st.binv[i*m:i*m+m], cb)
 	}
 }
 
@@ -523,11 +520,7 @@ func (lp *LP) pivot(st *State, q int, dir float64, a []float64, t float64, leave
 		if i == leaveRow || a[i] == 0 {
 			continue
 		}
-		factor := a[i]
-		row := st.binv[i*m : i*m+m]
-		for k := 0; k < m; k++ {
-			row[k] -= factor * rowR[k]
-		}
+		axpy(st.binv[i*m:i*m+m], rowR, -a[i])
 	}
 	st.basicOf[leaveRow] = q
 	st.status[q] = basic