feat: log_approx_f32 polynomial natural log intrinsic

CHANGELOG.md

@@ -11,6 +11,7 @@
 - `tanh_approx_f32(v: f32xN) -> f32xN` — dedicated fast vector tanh. Rational `P(x²) · x / Q(x²)` approximation in the Eigen / TensorFlow / JAX family: degree-13 numerator (odd in x), degree-6 denominator (even in x²), one fdiv per call. Clamped internally to `[-9, 9]`; max absolute error ~3e-7 across the body. Avoids `@llvm.tanh.v*f32`, which LLVM scalarizes to per-lane libm `tanhf`. Replaces the catastrophic-cancellation-prone `(exp_poly_f32(2x) - 1) / (exp_poly_f32(2x) + 1)` workaround the cookbook previously documented for tanh-GELU. Motivating consumer: Olorin's `gemma4_gelu` activation path. See `docs/superpowers/specs/2026-05-19-tanh-approx-f32-design.md`.
 - `u16x32` vector token + `lo256_u16x32(u16x32) -> u16x16` / `hi256_u16x32(u16x32) -> u16x16` lane extractors. Completes the i16/u16 symmetry — the signed pair (`lo256_i16x32` / `hi256_i16x32`) shipped in v1.12.0; PR #10 explicitly deferred the unsigned siblings because the `u16x32` token itself did not exist. Pure dispatch additions on the codegen side (typeck reuses `check_lo_extract` / `check_hi_extract`, codegen reuses width-generic `compile_lo_extract` / `compile_hi_extract`); the actual new code is in the lexer/parser layer for the type token. ARM rejection inherits from the existing >128-bit guard.
 - `wmul_u64(u32x4, u32x4) -> u64x4` — fused full-width widening multiply, completing the `wmul_u64` family alongside the existing `wmul_u64_lo` / `wmul_u64_hi` pair (v1.12.0). One call widens all four lanes; lowers to two `vpmuludq` + interleave via LLVM's `mul(zext, zext)` pattern-match. x86-only: the `u64x4` return type is 256-bit and rejected by the existing ARM >128-bit guard. ARM callers continue to use the lo/hi pair (each returning `u64x2`). Wider-input variants (`u32x8` / `u32x16` inputs) explicitly deferred — those require new lexer tokens and have no documented consumer yet.
+- `log_approx_f32(v: f32xN) -> f32xN` — natural log via Eigen/Cephes-family polynomial approximation. Bit-level decomposition (`x = m · 2^e` with `m ∈ [0.5, 1)`), √2/2 rebalance to center the polynomial range, degree-8 Horner in `(m - 1)`, Cody-Waite recombine with `e · ln(2)`. Max absolute error ~3e-6 across `(0, +∞)`. Avoids `@llvm.log.v*f32`, which LLVM scalarizes to per-lane libm `logf`. Companion to `exp_poly_f32` (v1.11.0) and `tanh_approx_f32` (v1.14.0); composes with `exp_poly_f32` to roundtrip-test pinned at ~1e-4 relative error.
 
 ## v1.13.0 — 2026-05-15 — ea bench + first aarch64 baselines + Specification umbrella
 

ROADMAP.md

@@ -28,6 +28,10 @@ Forward-looking notes. Ordered by leverage, not by effort.
 
 `wmul_u64(u32x4, u32x4) -> u64x4` widens all four lanes in a single intrinsic call, replacing the manual `wmul_u64_lo` + `wmul_u64_hi` + concat dance. Lowers to two `vpmuludq` + interleave via LLVM's `mul(zext, zext)` pattern-match. x86-only: the `u64x4` return type is 256-bit and rejected by the existing ARM >128-bit guard before the intrinsic dispatcher runs. ARM callers continue to use the lo/hi pair (each returning the NEON-fitting `u64x2`). Wider-input variants (`wmul_u64(u32x8, ...) -> u64x8`) explicitly deferred — they require new `u32x8` / `u32x16` lexer tokens and have no documented consumer yet. See `docs/superpowers/specs/2026-05-19-wmul-u64-fused-design.md`.
 
+### log_approx_f32
+
+`log_approx_f32(v: f32xN) -> f32xN`. Bit-level decomposition of `x = m · 2^e` (frexp convention, `m ∈ [0.5, 1)`), √2/2 rebalance to center the polynomial range, degree-8 Eigen-coefficient Horner in `(m - 1)`, then Cody-Waite recombine with `e · ln(2)`. Avoids `@llvm.log.v*f32`, which LLVM scalarizes to per-lane libm `logf`. Max absolute error ~3e-6 across `(0, +∞)`; matches `exp_poly_f32`'s 2⁻¹⁸ relative target. Composes cleanly with `exp_poly_f32` — pin-tested via a 4-input roundtrip kernel (`exp_poly_f32(log_approx_f32(x)) ≈ x` to ~1e-4 relative). See `docs/superpowers/specs/2026-05-19-log-approx-f32-design.md`.
+
 ## Shipped in v1.12.0 (2026-05-13)
 
 - **Deprecation-warning infrastructure** + `docs/migrations/` directory + `cargo public-api` CI gate (PR #6).
@@ -60,7 +64,7 @@ Today the language spec is spread across `docs/src/reference/*.md` (types, intri
 
 ## Future API consistency
 
-- **`log_approx_f32`, `sin_cos_approx_f32`** — polynomial approximations following the `exp_poly_f32` pattern. `tanh_approx_f32` shipped in v1.14.0; the remaining two are speculative until a real consumer asks.
+- **`sin_cos_approx_f32`** — polynomial approximation following the `exp_poly_f32` / `tanh_approx_f32` / `log_approx_f32` pattern. The remaining transcendental in the original "Future API consistency" trio after `tanh_approx_f32` and `log_approx_f32` shipped in v1.14.0. Speculative until a real consumer asks — angle range, periodicity strategy, and sin-vs-cos-vs-pair API shape all open design questions.
 - **Wider-input `wmul_u64` variants** — AVX2/AVX-512 widths (`wmul_u64(u32x8, u32x8) -> u64x8` on AVX-512, etc.). Requires `u32x8` / `u32x16` lexer tokens which don't exist yet. The fused `wmul_u64(u32x4, u32x4) -> u64x4` shipped in v1.14.0; wider widths gated on a consumer asking *and* providing the input tokens.
 
 ## Future additions

src/codegen/mod.rs

@@ -23,6 +23,8 @@ mod simd_fp16;
 #[cfg(feature = "llvm")]
 mod simd_lane;
 #[cfg(feature = "llvm")]
+mod simd_log_approx;
+#[cfg(feature = "llvm")]
 mod simd_masked;
 #[cfg(feature = "llvm")]
 mod simd_math;

src/codegen/simd.rs

@@ -36,6 +36,7 @@ impl<'ctx> CodeGenerator<'ctx> {
                 | "exp"
                 | "exp_poly_f32"
                 | "tanh_approx_f32"
+                | "log_approx_f32"
                 | "reduce_add"
                 | "reduce_add_fast"
                 | "reduce_max"
@@ -260,6 +261,7 @@ impl<'ctx> CodeGenerator<'ctx> {
             "exp" => self.compile_exp(args, function),
             "exp_poly_f32" => self.compile_exp_poly_f32(args, function),
             "tanh_approx_f32" => self.compile_tanh_approx_f32(args, function),
+            "log_approx_f32" => self.compile_log_approx_f32(args, function),
             "reduce_add" | "reduce_add_fast" | "reduce_max" | "reduce_min" => {
                 if self.call_uses_f16(args, None) {
                     return self.compile_reduce_f16(args, name, function);

src/codegen/simd_log_approx.rs

@@ -0,0 +1,281 @@
+//! Polynomial-based vector natural log for f32 vectors.
+//!
+//! Emits an Eigen / Cephes-family log_approx: bit-level decomposition of
+//! x = m · 2^e (frexp convention, m ∈ [0.5, 1)), √2/2 rebalance to center
+//! the polynomial range, degree-8 Horner in (m - 1), then Cody-Waite
+//! recombine with e · ln(2). Avoids `@llvm.log.v*f32`, which LLVM
+//! scalarizes to N sequential libm `logf` calls on every supported
+//! architecture.
+//!
+//! Defined input range: (0, +∞). For x ≤ 0, NaN, or ±∞, output is
+//! undefined. Matches `exp_poly_f32`'s "bounded input contract" style.
+//! Maximum absolute error: ~3e-6 across the defined range (compatible
+//! with `exp_poly_f32`'s 2⁻¹⁸ relative target).
+
+use inkwell::FloatPredicate;
+use inkwell::values::{BasicValueEnum, FunctionValue, VectorValue};
+
+use crate::ast::Expr;
+use crate::error::CompileError;
+
+use super::CodeGenerator;
+
+// Cody-Waite split of ln(2): LN2_HI + LN2_LO ≈ ln(2) to f32 precision.
+// LN2_HI is exact-representable in f32 (0.693359375 = 22188800 / 2^25);
+// LN2_LO carries the residual correction. Accumulating e · LN2_HI and
+// e · LN2_LO separately, with the low part added before the dominant
+// `+u` linear term, preserves precision when the integer exponent is
+// large.
+// LN2_HI is exact in f32 (binary `0.10110001 1`, fits in 9 mantissa bits)
+// but clippy's excessive-precision lint flags any literal with this many
+// digits regardless.
+#[allow(clippy::excessive_precision)]
+const LN2_HI: f32 = 0.693_359_375;
+const LN2_LO: f32 = -2.121_944_4e-4;
+
+// Eigen MathFunctionsImpl.h polynomial coefficients (MPL2).
+// Degree-8 Horner in u = m - 1, where m ∈ [√2/2, √2) after rebalance,
+// fitting the (log(1+u) - u + u²/2) / u³ tail of the Taylor expansion.
+// Literals kept at Eigen's published double precision for diff-friendliness;
+// the compiler truncates each to the nearest representable f32.
+#[allow(clippy::excessive_precision)]
+mod coeffs {
+    pub(super) const P0: f32 = 7.0376836292e-2;
+    pub(super) const P1: f32 = -1.1514610310e-1;
+    pub(super) const P2: f32 = 1.1676998740e-1;
+    pub(super) const P3: f32 = -1.2420140846e-1;
+    pub(super) const P4: f32 = 1.4249322787e-1;
+    pub(super) const P5: f32 = -1.6668057665e-1;
+    pub(super) const P6: f32 = 2.0000714765e-1;
+    pub(super) const P7: f32 = -2.4999993993e-1;
+    pub(super) const P8: f32 = 3.3333331174e-1;
+}
+use coeffs::*;
+
+// Frexp-convention masks. To extract m ∈ [0.5, 1):
+//   m_bits = (bits & MANTISSA_MASK) | EXP_HALF_BITS
+//   e_raw  = (bits >> 23) & 0xFF        (raw biased exponent)
+//   e      = e_raw - 126                (so x = m · 2^e with m ∈ [0.5, 1))
+const MANTISSA_MASK: i32 = 0x807F_FFFFu32 as i32; // sign + mantissa bits
+const EXP_HALF_BITS: i32 = 0x3F00_0000; // exponent = 126 (i.e. 2^-1)
+const EXP_BIAS: i32 = 126;
+// √2/2 — rebalance boundary. Spelled via the f32 const to avoid clippy's
+// approx_constant lint; the bit pattern matches the literal Eigen uses.
+const SQRT_HALF: f32 = std::f32::consts::FRAC_1_SQRT_2;
+
+impl<'ctx> CodeGenerator<'ctx> {
+    /// Compile `log_approx_f32(v: f32xN) -> f32xN`. Width inferred from operand.
+    /// Emits the bit-decomp + polynomial + Cody-Waite recombine directly —
+    /// never calls @llvm.log.
+    pub(super) fn compile_log_approx_f32(
+        &mut self,
+        args: &[Expr],
+        function: FunctionValue<'ctx>,
+    ) -> crate::error::Result<BasicValueEnum<'ctx>> {
+        let val = self.compile_expr(&args[0], function)?;
+        let v = match val {
+            BasicValueEnum::VectorValue(vv) => vv,
+            _ => {
+                return Err(CompileError::codegen_error(
+                    "log_approx_f32 expects f32 vector, got scalar; use log() for scalar libm-precision",
+                ));
+            }
+        };
+        let vec_ty = v.get_type();
+        let elem_ty = vec_ty.get_element_type();
+        if !elem_ty.is_float_type() || elem_ty.into_float_type() != self.context.f32_type() {
+            return Err(CompileError::codegen_error(
+                "log_approx_f32 expects f32 element type",
+            ));
+        }
+        let width = vec_ty.get_size();
+        let i32_vec_ty = self.context.i32_type().vec_type(width);
+
+        // 1. Bitcast input to integer for bit-level extraction.
+        let bits = self
+            .builder
+            .build_bit_cast(v, i32_vec_ty, "log_bits")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?
+            .into_vector_value();
+
+        // 2. e_raw = bits >> 23 (logical right shift — high bits become zero).
+        let shift_23 = self.splat_i32_const_la(23, width)?;
+        let shifted = self
+            .builder
+            .build_right_shift(bits, shift_23, false, "log_shr")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+        // Mask off sign bit: e_raw = shifted & 0xFF.
+        let mask_ff = self.splat_i32_const_la(0xFF, width)?;
+        let e_raw = self
+            .builder
+            .build_and(shifted, mask_ff, "log_e_raw")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+
+        // 3. e = e_raw - 126 (frexp bias).
+        let bias = self.splat_i32_const_la(EXP_BIAS, width)?;
+        let e_int = self
+            .builder
+            .build_int_sub(e_raw, bias, "log_e_int")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+
+        // 4. Mantissa: clear original exponent bits, set them to 126 (i.e. 2^-1).
+        //    m_bits = (bits & 0x807FFFFF) | 0x3F000000;  m = bitcast<f32>(m_bits).
+        let mantissa_mask = self.splat_i32_const_la(MANTISSA_MASK, width)?;
+        let exp_half = self.splat_i32_const_la(EXP_HALF_BITS, width)?;
+        let masked = self
+            .builder
+            .build_and(bits, mantissa_mask, "log_mant_masked")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+        let m_bits = self
+            .builder
+            .build_or(masked, exp_half, "log_m_bits")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+        let m = self
+            .builder
+            .build_bit_cast(m_bits, vec_ty, "log_m")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?
+            .into_vector_value();
+        // m ∈ [0.5, 1).
+
+        // 5. √2/2 rebalance. If m < √2/2, double m and decrement e.
+        //    Result: m ∈ [√2/2, √2), u = m - 1 ∈ [-0.293, 0.414].
+        let sqrt_half = self.splat_f32_const_la(SQRT_HALF, width)?;
+        let m_lt = self
+            .builder
+            .build_float_compare(FloatPredicate::OLT, m, sqrt_half, "log_m_lt_sqrth")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+        let m_doubled = self
+            .builder
+            .build_float_add(m, m, "log_m_doubled")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+        let m_rebalanced = self
+            .builder
+            .build_select(m_lt, m_doubled, m, "log_m_rebal")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?
+            .into_vector_value();
+        let one_i32 = self.splat_i32_const_la(1, width)?;
+        let e_decremented = self
+            .builder
+            .build_int_sub(e_int, one_i32, "log_e_dec")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+        let e_rebalanced = self
+            .builder
+            .build_select(m_lt, e_decremented, e_int, "log_e_rebal")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?
+            .into_vector_value();
+
+        // 6. Convert e (i32) to f32.
+        let e_f32 = self
+            .builder
+            .build_signed_int_to_float(e_rebalanced, vec_ty, "log_e_f32")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+
+        // 7. u = m - 1.
+        let one_f32 = self.splat_f32_const_la(1.0, width)?;
+        let u = self
+            .builder
+            .build_float_sub(m_rebalanced, one_f32, "log_u")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+
+        // 8. u² (reused in polynomial-tail correction).
+        let u2 = self
+            .builder
+            .build_float_mul(u, u, "log_u2")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+
+        // 9. Horner: P(u) = ((((((((p0·u + p1)·u + p2)·u + p3)·u + p4)·u + p5)·u + p6)·u + p7)·u + p8.
+        let p0 = self.splat_f32_const_la(P0, width)?;
+        let p1 = self.splat_f32_const_la(P1, width)?;
+        let p2 = self.splat_f32_const_la(P2, width)?;
+        let p3 = self.splat_f32_const_la(P3, width)?;
+        let p4 = self.splat_f32_const_la(P4, width)?;
+        let p5 = self.splat_f32_const_la(P5, width)?;
+        let p6 = self.splat_f32_const_la(P6, width)?;
+        let p7 = self.splat_f32_const_la(P7, width)?;
+        let p8 = self.splat_f32_const_la(P8, width)?;
+
+        let poly = self.fma_la(p0, u, p1, "log_poly1", width)?;
+        let poly = self.fma_la(poly, u, p2, "log_poly2", width)?;
+        let poly = self.fma_la(poly, u, p3, "log_poly3", width)?;
+        let poly = self.fma_la(poly, u, p4, "log_poly4", width)?;
+        let poly = self.fma_la(poly, u, p5, "log_poly5", width)?;
+        let poly = self.fma_la(poly, u, p6, "log_poly6", width)?;
+        let poly = self.fma_la(poly, u, p7, "log_poly7", width)?;
+        let poly = self.fma_la(poly, u, p8, "log_poly8", width)?;
+        // Now poly = P(u), degree 8.
+
+        // 10. y = P(u) · u · u² = P(u) · u³.
+        let y = self
+            .builder
+            .build_float_mul(poly, u, "log_pu")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+        let y = self
+            .builder
+            .build_float_mul(y, u2, "log_y_corr")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+
+        // 11. Cody-Waite recombine — small-magnitude terms first.
+        //     y = y - 0.5·u²    (the -u²/2 term of log(1+u) Taylor)
+        let neg_half = self.splat_f32_const_la(-0.5, width)?;
+        let y = self.fma_la(u2, neg_half, y, "log_y_minus_half_u2", width)?;
+
+        //     y = y + e_f32 · LN2_LO    (low part of ln(2)·e)
+        let ln2_lo = self.splat_f32_const_la(LN2_LO, width)?;
+        let y = self.fma_la(e_f32, ln2_lo, y, "log_y_plus_e_lo", width)?;
+
+        //     y = y + u    (linear term — dominant for small u)
+        let y = self
+            .builder
+            .build_float_add(y, u, "log_y_plus_u")
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?;
+
+        //     y = y + e_f32 · LN2_HI    (dominant ln(2)·e contribution last)
+        let ln2_hi = self.splat_f32_const_la(LN2_HI, width)?;
+        let y = self.fma_la(e_f32, ln2_hi, y, "log_y_plus_e_hi", width)?;
+
+        Ok(y.into())
+    }
+
+    fn splat_f32_const_la(
+        &self,
+        value: f32,
+        width: u32,
+    ) -> crate::error::Result<VectorValue<'ctx>> {
+        let scalar = self.context.f32_type().const_float(value as f64);
+        self.build_splat(BasicValueEnum::FloatValue(scalar), width)
+    }
+
+    fn splat_i32_const_la(
+        &self,
+        value: i32,
+        width: u32,
+    ) -> crate::error::Result<VectorValue<'ctx>> {
+        let scalar = self.context.i32_type().const_int(value as u64, true);
+        self.build_splat(BasicValueEnum::IntValue(scalar), width)
+    }
+
+    fn fma_la(
+        &mut self,
+        a: VectorValue<'ctx>,
+        b: VectorValue<'ctx>,
+        c: VectorValue<'ctx>,
+        name: &str,
+        width: u32,
+    ) -> crate::error::Result<VectorValue<'ctx>> {
+        let vec_ty = self.context.f32_type().vec_type(width);
+        let intrinsic_name = format!("llvm.fma.v{width}f32");
+        let fn_ty = vec_ty.fn_type(&[vec_ty.into(), vec_ty.into(), vec_ty.into()], false);
+        let intrinsic = self
+            .module
+            .get_function(&intrinsic_name)
+            .unwrap_or_else(|| self.module.add_function(&intrinsic_name, fn_ty, None));
+        let result = self
+            .builder
+            .build_call(intrinsic, &[a.into(), b.into(), c.into()], name)
+            .map_err(|e| CompileError::codegen_error(e.to_string()))?
+            .try_as_basic_value()
+            .basic()
+            .ok_or_else(|| CompileError::codegen_error("fma returned no value"))?;
+        Ok(result.into_vector_value())
+    }
+}