feat(ops)!: exact (erf) GELU and erf op (#1342, part 1)

CHANGELOG.md

@@ -2,6 +2,16 @@
 
 ## [Unreleased]
 
+### Added
+- **`TensorOps.erf` and exact GELU** (#1342): `gelu(tensor, approximate: Boolean = true)` keeps the tanh form as the default and selects the exact `0.5·x·(1 + erf(x/√2))` with `approximate = false`, the form BERT-class encoders and the Laya decision models were trained with (the two differ by up to 4.7e-4 near |x| ≈ 2.7). CPU: Abramowitz–Stegun 7.1.26 erf evaluated in double (|err| < 1.5e-7) with a JVM Vector API lane kernel; `GELU(approximate)` module, `network {}` `gelu(approximate)`, Java `TensorJavaOps.gelu(a, approximate)` / `erf(a)` and `SequentialModelBuilder.gelu(approximate)`; exact backward rule and `erfBackward`.
+- **GELU carries its form as op metadata**: the flag is recorded as `parameters["approximate"]` on the graph node (`GeluOperation`, `ErfOperation`), read by `ComputeGraphExecutor`, the backward rule and the StableHLO converter, so the compiled path cannot silently diverge from eager execution. StableHLO lowers `erf` to `chlo.erf` and exact GELU to `0.5·x·(1 + chlo.erf(x·√½))`; one `GeluLowering` builder emits both forms.
+
+### Changed
+- `TensorOps.gelu` gained the `approximate` parameter: implementors of the interface must update the override (callers are source-compatible through the default). The legacy recording decorator now records `gelu`/`erf` as graph operations instead of delegating them unrecorded.
+
+### Fixed
+- `DefaultCpuOps.gelu` printed a literal `${tensor.dtype}` in its unsupported-dtype message.
+
 ## [0.58.0] - 2026-10-12
 
 Headline: **input front ends and model identity (SKEEP-006).** Tokenizers leave `skainet-io-core`

docs/modules/ROOT/pages/reference/operators/generated/index.adoc

@@ -1,6 +1,6 @@
 = AI-NET Operators Reference
 
-Generated from version `0.58.0` on 2026-10-12
+Generated from version `0.58.0` on 2026-10-14
 
 == Operators by Modality
 

docs/modules/ROOT/pages/reference/operators/generated/tensorops.adoc

@@ -778,12 +778,13 @@ fun silu(tensor:Tensor): Tensor
 
 [source,kotlin]
 ----
-fun gelu(tensor:Tensor): Tensor
+fun gelu(tensor:Tensor, approximate:Boolean): Tensor
 ----
 
 === Parameters
 
 * `tensor: Tensor`
+* `approximate: Boolean`
 
 === Return Type
 
@@ -1287,6 +1288,25 @@ fun expm1(tensor:Tensor): Tensor
 
 `Tensor`
 
+== erf
+
+[.darc-none]#✖ Generated facts only (no human prose)#
+
+=== Signature
+
+[source,kotlin]
+----
+fun erf(tensor:Tensor): Tensor
+----
+
+=== Parameters
+
+* `tensor: Tensor`
+
+=== Return Type
+
+`Tensor`
+
 == sin
 
 [.darc-none]#✖ Generated facts only (no human prose)#

docs/modules/ROOT/pages/reference/ops-status-matrix.adoc

@@ -1,7 +1,7 @@
 = Operator Coverage Matrix
 :description: Cross-backend status for every operator function in SKaiNET.
 
-Generated from `operators.json` version `0.58.0` on 2026-10-12.
+Generated from `operators.json` version `0.58.0` on 2026-10-14.
 
 Rows are `Operator.function` pairs. The `Validated` column shows whether the function's documentation has been DARC-validated by a reviewer (see xref:contributing/darc-workflow.adoc[DARC workflow]). Remaining columns are backends that appear in any function's `statusByBackend` map — a missing entry means the backend makes no claim about the function (treat it as "unknown", not "not supported").
 
@@ -70,12 +70,13 @@ Rows are `Operator.function` pairs. The `Validated` column shows whether the fun
 | `TensorOps.indexSelect` | ✖ | — | — | — 
 | `TensorOps.exp` | ✖ | — | — | — 
 | `TensorOps.expm1` | ✖ | — | — | — 
+| `TensorOps.erf` | ✖ | — | — | — 
 | `TensorOps.sin` | ✖ | — | — | — 
 | `TensorOps.cos` | ✖ | — | — | — 
 | `TensorOps.scaledDotProductAttention` | ✖ | — | — | — 
 | `Similarity.cosineDistance` | ✖ | ✅ | ✅ | ✅ 
 
-| *Done* | *1 / 65* | *1 / 65* | *1 / 65* | *1 / 65* 
+| *Done* | *1 / 66* | *1 / 66* | *1 / 66* | *1 / 66* 
 |===
 
 Per-function detail including notes lives in xref:reference/operators/generated/index.adoc[Operator reference].

skainet-backends/skainet-backend-cpu/api/jvm/skainet-backend-cpu.api

@@ -300,13 +300,14 @@ public class sk/ainet/exec/tensor/ops/DefaultCpuOpsBase : sk/ainet/context/sched
 	protected final fun elementwise (Lsk/ainet/lang/tensor/Tensor;Lsk/ainet/lang/tensor/Tensor;Lkotlin/jvm/functions/Function3;)Lsk/ainet/lang/tensor/Tensor;
 	public fun elu (Lsk/ainet/lang/tensor/Tensor;F)Lsk/ainet/lang/tensor/Tensor;
 	protected fun ensureKernelProviders ()V
+	public fun erf (Lsk/ainet/lang/tensor/Tensor;)Lsk/ainet/lang/tensor/Tensor;
 	public fun exp (Lsk/ainet/lang/tensor/Tensor;)Lsk/ainet/lang/tensor/Tensor;
 	public fun expm1 (Lsk/ainet/lang/tensor/Tensor;)Lsk/ainet/lang/tensor/Tensor;
 	public fun flatten (Lsk/ainet/lang/tensor/Tensor;II)Lsk/ainet/lang/tensor/Tensor;
 	protected final fun floatResult (Lsk/ainet/lang/tensor/Shape;Lkotlin/reflect/KClass;[F[Lsk/ainet/lang/tensor/Tensor;)Lsk/ainet/lang/tensor/Tensor;
 	public fun gather (Lsk/ainet/lang/tensor/Tensor;Lsk/ainet/lang/tensor/Tensor;I)Lsk/ainet/lang/tensor/Tensor;
 	public fun ge (Lsk/ainet/lang/tensor/Tensor;F)Lsk/ainet/lang/tensor/Tensor;
-	public fun gelu (Lsk/ainet/lang/tensor/Tensor;)Lsk/ainet/lang/tensor/Tensor;
+	public fun gelu (Lsk/ainet/lang/tensor/Tensor;Z)Lsk/ainet/lang/tensor/Tensor;
 	protected final fun getDataFactory ()Lsk/ainet/lang/tensor/data/TensorDataFactory;
 	public final fun getSchedule ()Lsk/ainet/context/schedule/Schedule;
 	protected final fun gradStateFrom ([Lsk/ainet/lang/tensor/Tensor;)Lsk/ainet/lang/tensor/GradState;

skainet-backends/skainet-backend-cpu/src/commonMain/kotlin/sk/ainet/exec/tensor/ops/DefaultCpuOps.kt

@@ -2596,8 +2596,8 @@ public open class DefaultCpuOpsBase(
     }
 
     @TensorOp()
-    @InProgress("cpu", owner = "team:cpu", issue = "task-ops.md#op-gelu")
-    override fun <T : DType, V> gelu(tensor: Tensor<T, V>): Tensor<T, V> {
+    override fun <T : DType, V> gelu(tensor: Tensor<T, V>, approximate: Boolean): Tensor<T, V> {
+        if (!approximate) return geluExact(tensor)
         floatUnaryFast(tensor) { x ->
             val inner = x + 0.044715f * (x * x * x)
             0.5f * x * (1.0f + kotlin.math.tanh(0.7978845608f * inner))
@@ -2615,7 +2615,39 @@ public open class DefaultCpuOpsBase(
                     @Suppress("UNCHECKED_CAST")
                     y as V
                 }
-                else -> throw IllegalArgumentException("Unsupported dtype for gelu: ${'$'}{tensor.dtype}")
+                else -> throw IllegalArgumentException("Unsupported dtype for gelu: ${tensor.dtype}")
+            }
+        }
+        return newTensor(outData, tensor.dtype, tensor)
+    }
+
+    /** Exact GELU `0.5·x·(1 + erf(x/√2))`; erf from [ErfKernel] (|err| < 1.5e-7, evaluated in double). */
+    private fun <T : DType, V> geluExact(tensor: Tensor<T, V>): Tensor<T, V> {
+        floatUnaryFast(tensor) { x -> ErfKernel.geluExact(x) }?.let { return it }
+        val outData = dataFactory.init<T, V>(tensor.shape, tensor.dtype) { idx ->
+            when (tensor.dtype) {
+                sk.ainet.lang.types.FP32::class, sk.ainet.lang.types.FP16::class -> {
+                    val x = tensor.data.get(*idx) as Float
+                    @Suppress("UNCHECKED_CAST")
+                    ErfKernel.geluExact(x) as V
+                }
+                else -> throw IllegalArgumentException("Unsupported dtype for gelu: ${tensor.dtype}")
+            }
+        }
+        return newTensor(outData, tensor.dtype, tensor)
+    }
+
+    @TensorOp()
+    override fun <T : DType, V> erf(tensor: Tensor<T, V>): Tensor<T, V> {
+        floatUnaryFast(tensor) { x -> ErfKernel.erf(x.toDouble()).toFloat() }?.let { return it }
+        val outData = dataFactory.init<T, V>(tensor.shape, tensor.dtype) { idx ->
+            when (tensor.dtype) {
+                sk.ainet.lang.types.FP32::class, sk.ainet.lang.types.FP16::class -> {
+                    val x = tensor.data.get(*idx) as Float
+                    @Suppress("UNCHECKED_CAST")
+                    ErfKernel.erf(x.toDouble()).toFloat() as V
+                }
+                else -> throw IllegalArgumentException("Unsupported dtype for erf: ${tensor.dtype}")
             }
         }
         return newTensor(outData, tensor.dtype, tensor)

skainet-backends/skainet-backend-cpu/src/commonMain/kotlin/sk/ainet/exec/tensor/ops/ErfKernel.kt

@@ -0,0 +1,67 @@
+package sk.ainet.exec.tensor.ops
+
+import kotlin.math.abs
+import kotlin.math.exp
+import kotlin.math.sqrt
+
+/**
+ * Scalar Gauss error function and the exact GELU built on it (#1342).
+ *
+ * `erf` is Abramowitz–Stegun 7.1.26 (|err| < 1.5e-7) evaluated in double precision; the float
+ * result absorbs that error below its last bit. The expression order is the one the consumers
+ * validated against their reference goldens, so the kernel reproduces them bit-for-bit.
+ *
+ * The `*Float` variants evaluate the same polynomial in single precision with the same operation
+ * order; they serve the scalar tail of the JVM Vector API lane kernel, so lanes and tail agree.
+ */
+internal object ErfKernel {
+    private const val P = 0.3275911
+    private const val A1 = 0.254829592
+    private const val A2 = -0.284496736
+    private const val A3 = 1.421413741
+    private const val A4 = -1.453152027
+    private const val A5 = 1.061405429
+
+    const val P_F: Float = 0.3275911f
+    const val A1_F: Float = 0.254829592f
+    const val A2_F: Float = -0.284496736f
+    const val A3_F: Float = 1.421413741f
+    const val A4_F: Float = -1.453152027f
+    const val A5_F: Float = 1.061405429f
+    /** 1/√2 in single precision. */
+    const val SQRT1_2_F: Float = 0.70710678f
+
+    /** `erf(x)` with absolute error below 1.5e-7, odd, saturating at ±1. */
+    fun erf(x: Double): Double {
+        // The polynomial's coefficients sum to 1 only to ~1e-9, so pin the origin (and keep the
+        // sign of zero) instead of returning that residue with the sign of +0.
+        if (x == 0.0) return x
+        val sign = if (x < 0) -1.0 else 1.0
+        val a = abs(x)
+        val t = 1.0 / (1.0 + P * a)
+        val y = 1.0 - (((((A5 * t + A4) * t) + A3) * t + A2) * t + A1) * t * exp(-a * a)
+        return sign * y
+    }
+
+    /** `0.5 · x · (1 + erf(x / √2))` in double precision, rounded once to float. */
+    fun geluExact(x: Float): Float {
+        val xd = x.toDouble()
+        return (0.5 * xd * (1.0 + erf(xd / sqrt(2.0)))).toFloat()
+    }
+
+    /** Single-precision erf with the lane kernel's operation order (see `DefaultCpuOpsJvm.erfVector`). */
+    fun erfFloat(x: Float): Float {
+        if (x == 0f) return x
+        val a = abs(x)
+        val t = 1f / (a * P_F + 1f)
+        var poly = t * A5_F + A4_F
+        poly = poly * t + A3_F
+        poly = poly * t + A2_F
+        poly = poly * t + A1_F
+        val y = 1f - poly * t * exp(-(a * a))
+        return if (x < 0f) -y else y
+    }
+
+    /** Single-precision exact GELU matching the lane kernel. */
+    fun geluExactFloat(x: Float): Float = 0.5f * x * (1f + erfFloat(x * SQRT1_2_F))
+}

skainet-backends/skainet-backend-cpu/src/commonTest/kotlin/sk/ainet/sk/ainet/exec/tensor/ops/DefaultCpuOpsActivationsTest.kt

@@ -151,4 +151,103 @@ class DefaultCpuOpsActivationsTest {
         }
         assertTrue(threw, "Expected IllegalArgumentException for Int32 tanh")
     }
+
+    // ---- erf / exact GELU (#1342) ----------------------------------------
+
+    /** Abramowitz–Stegun 7.1.26 in double, the host reference the kernel must reproduce. */
+    private fun erfRef(x: Double): Double {
+        val sign = if (x < 0) -1.0 else 1.0
+        val a = kotlin.math.abs(x)
+        val t = 1.0 / (1.0 + 0.3275911 * a)
+        val poly = ((((1.061405429 * t - 1.453152027) * t) + 1.421413741) * t - 0.284496736) * t + 0.254829592
+        return sign * (1.0 - poly * t * kotlin.math.exp(-a * a))
+    }
+
+    private fun geluExactRef(x: Float): Float {
+        val xd = x.toDouble()
+        return (0.5 * xd * (1.0 + erfRef(xd / kotlin.math.sqrt(2.0)))).toFloat()
+    }
+
+    private fun geluTanhRef(x: Float): Float {
+        val inner = x + 0.044715f * (x * x * x)
+        return 0.5f * x * (1.0f + kotlin.math.tanh(0.7978845608f * inner))
+    }
+
+    @Test
+    fun erf_matches_reference_table() {
+        val xs = floatArrayOf(0f, 0.5f, -0.5f, 1f, -1f, 1.5f, -1.5f, 2f, -2f, 3f, -3f, 4f, -4f, 6f, -6f)
+        val table = floatArrayOf(
+            0f, 0.5204998778f, -0.5204998778f, 0.8427007929f, -0.8427007929f, 0.9661051465f, -0.9661051465f,
+            0.9953222650f, -0.9953222650f, 0.9999779095f, -0.9999779095f, 0.9999999846f, -0.9999999846f, 1f, -1f,
+        )
+        val out = cpuOps.erf(fTensor(Shape(xs.size), xs))
+        assertEquals(Shape(xs.size), out.shape)
+        assertEquals(FP32::class, out.dtype)
+        for (i in xs.indices) {
+            assertAlmostEquals(table[i], out.data[i] as Float, 2e-7f, "erf(${xs[i]}) mismatch: ${out.data[i]}")
+        }
+    }
+
+    @Test
+    fun erf_is_odd_and_saturates() {
+        val xs = FloatArray(41) { -4f + it * 0.2f }
+        val pos = cpuOps.erf(fTensor(Shape(xs.size), xs))
+        val neg = cpuOps.erf(fTensor(Shape(xs.size), FloatArray(xs.size) { -xs[it] }))
+        for (i in xs.indices) {
+            assertAlmostEquals(-(pos.data[i] as Float), neg.data[i] as Float, 0f, "erf must be odd at x=${xs[i]}")
+        }
+        val sat = cpuOps.erf(fTensor(Shape(2), floatArrayOf(6f, -6f)))
+        assertEquals(1f, sat.data[0] as Float)
+        assertEquals(-1f, sat.data[1] as Float)
+    }
+
+    @Test
+    fun gelu_exact_matches_host_formula() {
+        val xs = FloatArray(121) { -6f + it * 0.1f }
+        val out = cpuOps.gelu(fTensor(Shape(xs.size), xs), approximate = false)
+        for (i in xs.indices) {
+            // Bit-equal on the JVM; JS/Wasm evaluate Float in double precision, so allow one ulp there.
+            assertAlmostEquals(geluExactRef(xs[i]), out.data[i] as Float, 2e-7f, "exact gelu at x=${xs[i]}")
+        }
+    }
+
+    @Test
+    fun gelu_default_is_unchanged_tanh() {
+        val xs = floatArrayOf(-3f, -2f, -1f, -0.5f, 0f, 0.5f, 1f, 2f, 3f)
+        val t = fTensor(Shape(xs.size), xs)
+        val byDefault = cpuOps.gelu(t)
+        val explicit = cpuOps.gelu(t, approximate = true)
+        for (i in xs.indices) {
+            assertEquals(byDefault.data[i] as Float, explicit.data[i] as Float, "default must equal approximate=true")
+            assertAlmostEquals(geluTanhRef(xs[i]), byDefault.data[i] as Float, 1e-6f, "tanh gelu at x=${xs[i]}")
+        }
+    }
+
+    @Test
+    fun gelu_exact_vs_tanh_max_difference() {
+        val xs = FloatArray(1201) { -6f + it * 0.01f }
+        val t = fTensor(Shape(xs.size), xs)
+        val exact = cpuOps.gelu(t, approximate = false)
+        val tanh = cpuOps.gelu(t, approximate = true)
+        var maxDiff = 0f
+        var argMax = 0f
+        for (i in xs.indices) {
+            val d = kotlin.math.abs((exact.data[i] as Float) - (tanh.data[i] as Float))
+            if (d > maxDiff) { maxDiff = d; argMax = xs[i] }
+        }
+        // The tanh approximation deviates from the exact form by ≈ 4.7e-4 near |x| ≈ 2.7.
+        assertTrue(maxDiff > 4.5e-4f && maxDiff < 5.0e-4f, "max |exact - tanh| = $maxDiff at x=$argMax")
+        assertTrue(kotlin.math.abs(argMax) in 2.5f..2.9f, "arg max at x=$argMax")
+    }
+
+    @Test
+    fun erf_and_exact_gelu_unsupported_dtype_int32_throws() {
+        val input = iTensor(Shape(2), intArrayOf(1, 2))
+        @Suppress("UNCHECKED_CAST")
+        val asInt = input as sk.ainet.lang.tensor.Tensor<Int32, Int>
+        var threw = 0
+        try { cpuOps.erf(asInt) } catch (e: IllegalArgumentException) { threw++; assertTrue(e.message?.contains("Unsupported dtype for erf") == true) }
+        try { cpuOps.gelu(asInt, approximate = false) } catch (e: IllegalArgumentException) { threw++; assertTrue(e.message?.contains("Unsupported dtype for gelu: ") == true && e.message?.contains("Int32") == true, e.message) }
+        assertEquals(2, threw)
+    }
 }

skainet-backends/skainet-backend-cpu/src/jvmMain/kotlin/sk/ainet/exec/tensor/ops/DefaultCpuOpsJvm.kt

@@ -859,6 +859,39 @@ internal class DefaultCpuOpsJvm(
         return floatResult(Shape(tensor.shape.dimensions.copyOf()), tensor.dtype, out)
     }
 
+    override fun <T : DType, V> erf(tensor: Tensor<T, V>): Tensor<T, V> {
+        vectorFloatUnary(tensor, { v -> erfVector(v) }, ErfKernel::erfFloat)?.let { return it }
+        return super.erf(tensor)
+    }
+
+    override fun <T : DType, V> gelu(tensor: Tensor<T, V>, approximate: Boolean): Tensor<T, V> {
+        if (approximate) return super.gelu(tensor, true)
+        vectorFloatUnary(tensor, { v ->
+            val half = FloatVector.broadcast(floatSpecies, 0.5f)
+            val one = FloatVector.broadcast(floatSpecies, 1f)
+            half.mul(v).mul(one.add(erfVector(v.mul(ErfKernel.SQRT1_2_F))))
+        }, ErfKernel::geluExactFloat)?.let { return it }
+        return super.gelu(tensor, false)
+    }
+
+    /**
+     * Lane-wise Abramowitz–Stegun 7.1.26 erf in single precision — the same operation order as
+     * [ErfKernel.erfFloat], which serves the scalar tail, so lanes and tail agree up to the
+     * last-ulp difference between `VectorOperators.EXP` and `kotlin.math.exp`.
+     */
+    private fun erfVector(x: FloatVector): FloatVector {
+        val one = FloatVector.broadcast(floatSpecies, 1f)
+        val a = x.abs()
+        val t = one.div(a.mul(ErfKernel.P_F).add(1f))
+        var poly = t.mul(ErfKernel.A5_F).add(ErfKernel.A4_F)
+        poly = poly.mul(t).add(ErfKernel.A3_F)
+        poly = poly.mul(t).add(ErfKernel.A2_F)
+        poly = poly.mul(t).add(ErfKernel.A1_F)
+        val e = a.mul(a).neg().lanewise(VectorOperators.EXP)
+        val y = one.sub(poly.mul(t).mul(e))
+        return y.blend(y.neg(), x.lt(0f))
+    }
+
     override fun <T : DType, V> sum(tensor: Tensor<T, V>, dim: Int?): Tensor<T, V> {
         if (dim == null) {
             vectorFloatReduceAllSum<T, V>(tensor)?.let { return it }

skainet-backends/skainet-backend-cpu/src/jvmTest/kotlin/sk/ainet/exec/tensor/ops/DefaultCpuOpsJvmElementwiseTest.kt

@@ -55,6 +55,24 @@ class DefaultCpuOpsJvmElementwiseTest {
         }
     }
 
+    @Test
+    fun erf_and_exact_gelu_match_scalar_reference_with_tails() {
+        val species = jdk.incubator.vector.FloatVector.SPECIES_PREFERRED.length()
+        val lengths = listOf(1, 3, 7, species - 1, species, species + 1, species * 2 + 5, 257)
+        val scalar = DefaultCpuOps(dataFactory)
+        val jvm = DefaultCpuOpsJvm(dataFactory)
+        for (n in lengths) {
+            val a = FloatArray(n) { -8f + 16f * it / n.coerceAtLeast(2) }
+            val tA = floatTensor(Shape(n), a)
+            // erf: lane-wise single-precision polynomial vs the double-evaluated scalar kernel
+            compareTensors(scalar.erf(tA), jvm.erf(tA), eps = 1e-6f)
+            // exact gelu: |x| up to 8, so allow two ulps at that magnitude
+            compareTensors(scalar.gelu(tA, approximate = false), jvm.gelu(tA, approximate = false), eps = 2e-6f)
+            // the tanh default must stay byte-identical: the JVM tier delegates to the base kernel
+            compareTensors(scalar.gelu(tA), jvm.gelu(tA), eps = 0f)
+        }
+    }
+
     @Test
     fun vector_flag_on_and_off_produce_same_results() {
         val n = 257 // ensure tail