diff --git a/CHANGELOG.md b/CHANGELOG.md index a197511..6194e78 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/ROADMAP.md b/ROADMAP.md index 0b9ab8c..8b52033 100644 --- a/ROADMAP.md +++ b/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 diff --git a/src/codegen/mod.rs b/src/codegen/mod.rs index 41a9778..e160439 100644 --- a/src/codegen/mod.rs +++ b/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; diff --git a/src/codegen/simd.rs b/src/codegen/simd.rs index edfa454..71fd3b0 100644 --- a/src/codegen/simd.rs +++ b/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); diff --git a/src/codegen/simd_log_approx.rs b/src/codegen/simd_log_approx.rs new file mode 100644 index 0000000..201cc4c --- /dev/null +++ b/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> { + 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(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> { + 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> { + 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> { + 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()) + } +} diff --git a/src/typeck/intrinsics.rs b/src/typeck/intrinsics.rs index e09d584..dbe4882 100644 --- a/src/typeck/intrinsics.rs +++ b/src/typeck/intrinsics.rs @@ -46,6 +46,7 @@ impl TypeChecker { "sqrt" | "rsqrt" | "exp" => Some(self.check_sqrt(name, args, locals, span)), "exp_poly_f32" => Some(self.check_exp_poly_f32(args, locals, span)), "tanh_approx_f32" => Some(self.check_tanh_approx_f32(args, locals, span)), + "log_approx_f32" => Some(self.check_log_approx_f32(args, locals, span)), "to_f32" | "to_f64" | "to_f16" | "to_i16" | "to_i32" | "to_i64" => { Some(self.check_conversion(name, args, locals, span)) } @@ -434,6 +435,39 @@ impl TypeChecker { } } + /// Type-check `log_approx_f32(v: f32xN) -> f32xN`. Same f32-vector-only + /// shape as `exp_poly_f32`; scalar / f64 / f16 / integer all rejected. + fn check_log_approx_f32( + &self, + args: &[Expr], + locals: &HashMap, + span: &Span, + ) -> crate::error::Result { + if args.len() != 1 { + return Err(CompileError::type_error( + format!("log_approx_f32 expects 1 argument, got {}", args.len()), + span.clone(), + )); + } + let t = self.check_expr(&args[0], locals)?; + match &t { + Type::Vector { elem, .. } if **elem == Type::F32 => Ok(t), + Type::Vector { elem, .. } => Err(CompileError::type_error( + format!("log_approx_f32 expects f32 element type, got {elem}"), + span.clone(), + )), + Type::F32 | Type::FloatLiteral => Err(CompileError::type_error( + "log_approx_f32 expects f32 vector, got scalar; use log() for scalar libm-precision" + .to_string(), + span.clone(), + )), + _ => Err(CompileError::type_error( + format!("log_approx_f32 expects float vector, got {t}"), + span.clone(), + )), + } + } + fn check_prefetch( &self, args: &[Expr], diff --git a/tests/phase14_log_approx.rs b/tests/phase14_log_approx.rs new file mode 100644 index 0000000..9a11ecb --- /dev/null +++ b/tests/phase14_log_approx.rs @@ -0,0 +1,501 @@ +#[cfg(feature = "llvm")] +mod tests { + use ea_compiler::{CompileOptions, OutputMode}; + use std::process::Command; + use tempfile::TempDir; + + /// Smoke test: log_approx_f32(splat(1.0)) should give all 0.0s. + /// Exercises the full pipeline — bit extraction gives e=0, m=1.0 + /// (after exponent reset), √2/2 rebalance doesn't fire (1.0 > √2/2), + /// u = m - 1 = 0, polynomial returns 0, recombine returns 0. + #[test] + fn test_log_approx_f32x4_at_one() { + let ea = r#" + export func k(out: *mut f32) { + let z: f32x4 = splat(1.0) + let r: f32x4 = log_approx_f32(z) + store(out, 0, r) + } + "#; + let c = r#" + #include + extern void k(float *out); + int main(void) { + float out[4] = {1, 1, 1, 1}; + k(out); + for (int i = 0; i < 4; ++i) printf("%g\n", out[i]); + return 0; + } + "#; + let dir = TempDir::new().unwrap(); + let obj = dir.path().join("k.o"); + let cpath = dir.path().join("h.c"); + let bin = dir.path().join("k_bin"); + let opts = CompileOptions { + opt_level: 3, + target_cpu: None, + extra_features: String::new(), + target_triple: None, + }; + ea_compiler::compile_with_options(ea, &obj, OutputMode::ObjectFile, &opts) + .expect("compile failed"); + std::fs::write(&cpath, c).expect("write c"); + let status = Command::new("cc") + .args([ + cpath.to_str().unwrap(), + obj.to_str().unwrap(), + "-o", + bin.to_str().unwrap(), + "-lm", + ]) + .status() + .expect("link failed"); + assert!(status.success(), "linker failed"); + let out = Command::new(&bin).output().expect("run failed"); + let stdout = String::from_utf8_lossy(&out.stdout).replace("\r\n", "\n"); + assert_eq!(stdout.trim(), "0\n0\n0\n0"); + } + + /// Compile a kernel that runs log_approx_f32 over `inputs` lane-by-lane, + /// link with C harness that calls logf, assert absolute error ≤ 3e-6. + /// + /// Absolute (not relative) error because log(x) → 0 as x → 1, making + /// relative error blow up near 1. Across the rest of the input range + /// the magnitude of log(x) is bounded modestly, so absolute is the + /// natural metric. + fn accuracy_test_impl(inputs: &[f32], vector_type: &str) { + let lanes = if vector_type == "f32x4" { 4 } else { 8 }; + + let ea = format!( + r#" + export func k(input: *f32, output: *mut f32, n: i32) {{ + let mut i: i32 = 0 + while i + {lanes} <= n {{ + let v: {vector_type} = load(input, i) + let r: {vector_type} = log_approx_f32(v) + store(output, i, r) + i = i + {lanes} + }} + }} + "# + ); + + let mut padded = inputs.to_vec(); + while !padded.len().is_multiple_of(lanes) { + padded.push(1.0); // pad with 1.0 so log_approx_f32 returns 0 on padding + } + let n = padded.len(); + let original_n = inputs.len(); + + let inputs_str = padded + .iter() + .map(|f| format!("{f:.10e}f")) + .collect::>() + .join(", "); + + let c = format!( + r#" + #include + #include + extern void k(const float *input, float *output, int n); + int main(void) {{ + float in[{n}] = {{{inputs_str}}}; + float out[{n}] = {{0}}; + k(in, out, {n}); + for (int i = 0; i < {original_n}; ++i) {{ + float ref = logf(in[i]); + float got = out[i]; + float abs_err = fabsf(got - ref); + if (abs_err > 3.0e-6f) {{ + printf("FAIL i=%d in=%g got=%g ref=%g abs=%g\n", i, in[i], got, ref, abs_err); + return 1; + }} + }} + printf("OK\n"); + return 0; + }} + "# + ); + + let dir = TempDir::new().unwrap(); + let obj = dir.path().join("k.o"); + let cpath = dir.path().join("h.c"); + let bin = dir.path().join("k_bin"); + let opts = CompileOptions { + opt_level: 3, + target_cpu: None, + extra_features: String::new(), + target_triple: None, + }; + ea_compiler::compile_with_options(&ea, &obj, OutputMode::ObjectFile, &opts) + .expect("compile failed"); + std::fs::write(&cpath, c).expect("write c"); + let status = Command::new("cc") + .args([ + cpath.to_str().unwrap(), + obj.to_str().unwrap(), + "-o", + bin.to_str().unwrap(), + "-lm", + ]) + .status() + .expect("link failed"); + assert!(status.success(), "linker failed"); + let out = Command::new(&bin).output().expect("run failed"); + let stdout = String::from_utf8_lossy(&out.stdout).replace("\r\n", "\n"); + assert_eq!( + stdout.trim(), + "OK", + "stderr: {}", + String::from_utf8_lossy(&out.stderr) + ); + } + + #[test] + fn test_log_approx_f32x4_boundary_points() { + accuracy_test_impl( + &[ + // Powers of 2 — exercise the e_int extraction across the full + // f32 exponent range + 1.0 / 1024.0, + 1.0 / 8.0, + 0.5, + std::f32::consts::FRAC_1_SQRT_2, // √2/2 — rebalance boundary + 1.0, + std::f32::consts::SQRT_2, // √2 — opposite boundary + 2.0, + 8.0, + 1024.0, + // Transcendental constants + std::f32::consts::E, + std::f32::consts::PI, + // Larger values + 100.0, + 1.0e6, + 1.0e9, + ], + "f32x4", + ); + } + + /// e^log(x) ≈ x roundtrip — combined exp_poly + log_approx accuracy guard. + /// The two intrinsics use compatible Cody-Waite ln(2) splits, so the + /// combined error should stay within ~3e-5 (each contributes ~3e-6). + #[test] + fn test_log_approx_then_exp_poly_roundtrip() { + let ea = r#" + export func k(input: *f32, output: *mut f32) { + let x: f32x4 = load(input, 0) + let l: f32x4 = log_approx_f32(x) + let r: f32x4 = exp_poly_f32(l) + store(output, 0, r) + } + "#; + let c = r#" + #include + #include + extern void k(const float *input, float *output); + int main(void) { + float in[4] = {0.5f, 2.0f, 7.5f, 42.0f}; + float out[4] = {0}; + k(in, out); + for (int i = 0; i < 4; ++i) { + float rel = fabsf(out[i] - in[i]) / in[i]; + if (rel > 1.0e-4f) { + printf("FAIL i=%d in=%g out=%g rel=%g\n", i, in[i], out[i], rel); + return 1; + } + } + printf("OK\n"); + return 0; + } + "#; + let dir = TempDir::new().unwrap(); + let obj = dir.path().join("k.o"); + let cpath = dir.path().join("h.c"); + let bin = dir.path().join("k_bin"); + let opts = CompileOptions { + opt_level: 3, + target_cpu: None, + extra_features: String::new(), + target_triple: None, + }; + ea_compiler::compile_with_options(ea, &obj, OutputMode::ObjectFile, &opts) + .expect("compile failed"); + std::fs::write(&cpath, c).expect("write c"); + let status = Command::new("cc") + .args([ + cpath.to_str().unwrap(), + obj.to_str().unwrap(), + "-o", + bin.to_str().unwrap(), + "-lm", + ]) + .status() + .expect("link failed"); + assert!(status.success(), "linker failed"); + let out = Command::new(&bin).output().expect("run failed"); + let stdout = String::from_utf8_lossy(&out.stdout).replace("\r\n", "\n"); + assert_eq!(stdout.trim(), "OK"); + } + + #[test] + fn test_log_approx_f32x4_random() { + // Sample log-uniformly over [0.01, 100] — i.e., 2^[-7, 7] roughly. + // Avoids the x→0 region where the approximation degrades. + let mut points = Vec::new(); + let mut state: u32 = 0xDEAD_BEEF; + for _ in 0..256 { + state = state.wrapping_mul(1664525).wrapping_add(1013904223); + // Map to [log(0.01), log(100)] = [-4.6, 4.6], then exp. + let u = state as f32 / u32::MAX as f32; + let logx = -4.6 + 9.2 * u; + let x = logx.exp(); + points.push(x); + } + accuracy_test_impl(&points, "f32x4"); + } + + #[cfg(target_arch = "x86_64")] + #[test] + fn test_log_approx_f32x8_boundary_points() { + accuracy_test_impl( + &[ + 1.0 / 1024.0, + 1.0 / 8.0, + 0.5, + std::f32::consts::FRAC_1_SQRT_2, + 1.0, + std::f32::consts::SQRT_2, + 2.0, + 8.0, + 1024.0, + std::f32::consts::E, + std::f32::consts::PI, + 100.0, + 1.0e6, + 1.0e9, + ], + "f32x8", + ); + } + + #[cfg(target_arch = "x86_64")] + #[test] + fn test_log_approx_f32x8_random() { + let mut points = Vec::new(); + let mut state: u32 = 0xCAFE_F00D; + for _ in 0..256 { + state = state.wrapping_mul(1664525).wrapping_add(1013904223); + let u = state as f32 / u32::MAX as f32; + let logx = -4.6 + 9.2 * u; + let x = logx.exp(); + points.push(x); + } + accuracy_test_impl(&points, "f32x8"); + } + + /// CRITICAL regression guard: a future "simplification" of + /// compile_log_approx_f32 to delegate to compile_log (which doesn't + /// exist yet but would scalarize via libm) would silently undo the + /// vectorization that motivates the intrinsic. Pins absence of @llvm.log. + #[cfg(target_arch = "x86_64")] + #[test] + fn test_log_approx_f32_does_not_emit_llvm_log() { + let src = r#" + export func k(v: f32x8) -> f32x8 { + return log_approx_f32(v) + } + "#; + let opts = CompileOptions { + opt_level: 0, + target_cpu: None, + extra_features: String::new(), + target_triple: None, + }; + let dir = TempDir::new().unwrap(); + let ir_path = dir.path().join("k.ll"); + ea_compiler::compile_with_options(src, &ir_path, OutputMode::LlvmIr, &opts) + .expect("compile failed"); + let ir = std::fs::read_to_string(&ir_path).expect("read IR"); + assert!( + !ir.contains("@llvm.log"), + "log_approx_f32 must NOT lower to @llvm.log; found scalarization-prone intrinsic:\n{ir}" + ); + } + + /// Confirm the bit-manipulation + Horner pattern is emitted for f32x8. + #[cfg(target_arch = "x86_64")] + #[test] + fn test_log_approx_f32x8_emits_expected_pattern() { + let src = r#" + export func k(v: f32x8) -> f32x8 { + return log_approx_f32(v) + } + "#; + let opts = CompileOptions { + opt_level: 0, + target_cpu: None, + extra_features: String::new(), + target_triple: None, + }; + let dir = TempDir::new().unwrap(); + let ir_path = dir.path().join("k.ll"); + ea_compiler::compile_with_options(src, &ir_path, OutputMode::LlvmIr, &opts) + .expect("compile failed"); + let ir = std::fs::read_to_string(&ir_path).expect("read IR"); + + // Horner needs ≥8 FMAs for the polynomial; plus the recombine + // (e·c1, e·c2) FMAs add ~2 more. + let fma_count = ir.matches("@llvm.fma.v8f32").count(); + assert!( + fma_count >= 8, + "expected at least 8 @llvm.fma.v8f32 calls (8 polynomial + 2 recombine); got {fma_count}\nIR:\n{ir}" + ); + // Bit manipulation: bitcast float↔int, lshr (right shift), and/or + assert!( + ir.contains("bitcast <8 x float>") && ir.contains("to <8 x i32>"), + "expected float→int bitcast for exponent extraction\nIR:\n{ir}" + ); + assert!( + ir.contains("bitcast <8 x i32>") && ir.contains("to <8 x float>"), + "expected int→float bitcast for mantissa reconstruction\nIR:\n{ir}" + ); + assert!( + ir.contains("lshr <8 x i32>"), + "expected lshr <8 x i32> for exponent extraction\nIR:\n{ir}" + ); + // sitofp to convert exponent to f32 + assert!( + ir.contains("sitofp <8 x i32>"), + "expected sitofp <8 x i32> for exponent → f32 conversion\nIR:\n{ir}" + ); + } + + #[test] + fn test_log_approx_f32x4_emits_expected_pattern() { + let src = r#" + export func k(v: f32x4) -> f32x4 { + return log_approx_f32(v) + } + "#; + let opts = CompileOptions { + opt_level: 0, + target_cpu: None, + extra_features: String::new(), + target_triple: None, + }; + let dir = TempDir::new().unwrap(); + let ir_path = dir.path().join("k.ll"); + ea_compiler::compile_with_options(src, &ir_path, OutputMode::LlvmIr, &opts) + .expect("compile failed"); + let ir = std::fs::read_to_string(&ir_path).expect("read IR"); + assert!(!ir.contains("@llvm.log"), "no @llvm.log expected:\n{ir}"); + let fma_count = ir.matches("@llvm.fma.v4f32").count(); + assert!( + fma_count >= 8, + "expected at least 8 @llvm.fma.v4f32 calls; got {fma_count}\nIR:\n{ir}" + ); + assert!( + ir.contains("sitofp <4 x i32>"), + "expected sitofp <4 x i32>:\n{ir}" + ); + } + + /// Scalar f32 should fail with helpful message pointing at log(). + #[test] + fn test_log_approx_f32_rejects_scalar_f32() { + let src = r#" + export func k(x: f32) -> f32 { + return log_approx_f32(x) + } + "#; + let opts = CompileOptions { + opt_level: 0, + target_cpu: None, + extra_features: String::new(), + target_triple: None, + }; + let dir = TempDir::new().unwrap(); + let obj = dir.path().join("k.o"); + let err = ea_compiler::compile_with_options(src, &obj, OutputMode::ObjectFile, &opts) + .expect_err("scalar f32 should fail"); + let msg = format!("{err}"); + assert!( + msg.contains("f32 vector"), + "error must mention f32 vector requirement, got: {msg}" + ); + } + + #[test] + fn test_log_approx_f32_rejects_f64_vector() { + let src = r#" + export func k(v: f64x2) -> f64x2 { + return log_approx_f32(v) + } + "#; + let opts = CompileOptions { + opt_level: 0, + target_cpu: None, + extra_features: String::new(), + target_triple: None, + }; + let dir = TempDir::new().unwrap(); + let obj = dir.path().join("k.o"); + let err = ea_compiler::compile_with_options(src, &obj, OutputMode::ObjectFile, &opts) + .expect_err("f64x2 should fail"); + let msg = format!("{err}"); + assert!( + msg.contains("f32") && msg.contains("element type"), + "error must mention f32 element type, got: {msg}" + ); + } + + #[test] + fn test_log_approx_f32_rejects_integer_vector() { + let src = r#" + export func k(v: i32x4) -> i32x4 { + return log_approx_f32(v) + } + "#; + let opts = CompileOptions { + opt_level: 0, + target_cpu: None, + extra_features: String::new(), + target_triple: None, + }; + let dir = TempDir::new().unwrap(); + let obj = dir.path().join("k.o"); + let err = ea_compiler::compile_with_options(src, &obj, OutputMode::ObjectFile, &opts) + .expect_err("i32x4 should fail"); + let msg = format!("{err}"); + assert!( + msg.contains("float") || msg.contains("f32"), + "error must mention float requirement, got: {msg}" + ); + } + + #[test] + fn test_log_approx_f32_rejects_wrong_arity() { + let src = r#" + export func k(a: f32x4, b: f32x4) -> f32x4 { + return log_approx_f32(a, b) + } + "#; + let opts = CompileOptions { + opt_level: 0, + target_cpu: None, + extra_features: String::new(), + target_triple: None, + }; + let dir = TempDir::new().unwrap(); + let obj = dir.path().join("k.o"); + let err = ea_compiler::compile_with_options(src, &obj, OutputMode::ObjectFile, &opts) + .expect_err("wrong arity should fail"); + let msg = format!("{err}"); + assert!( + msg.contains("1 argument") || msg.contains("expects 1"), + "error must mention arity, got: {msg}" + ); + } +}