From f8f1e4a83105077caf06d9b0535e1f71c458a45c Mon Sep 17 00:00:00 2001 From: Peter Lukka Date: Tue, 19 May 2026 06:16:40 +0000 Subject: [PATCH 1/3] feat(intrinsic): sin_approx_f32 + cos_approx_f32 with shared core MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two intrinsics sharing a common range-reduction + polynomial core: reduce mod π/2 via 2-piece Cody-Waite (PI_2_HI exact f32, PI_2_LO negative residual; sum reproduces π/2 to ~2e-15 in f64), then compute both sin and cos polynomials over d' ∈ [-π/4, π/4]: - sin: degree-3 in s=d², covers d through d^7 — truncation ≤3.4e-7 - cos: degree-4 in s, covers d^0 through d^8 — truncation ≤2.6e-8 The cos variant reuses the same core with q += 1 — a precision-free integer shift expressing the mathematical identity cos(x) = sin(x + π/2). Adding π/2 to v before reduction would lose bits for large |v|; shifting the integer quadrant index after reduction is exact. Final blend: swap = (k & 1), negate = (k & 2). Both polynomials are always computed (~4 wasted FMAs) for branchless SIMD execution; LLVM CSE eliminates redundant range reduction at the caller level when both intrinsics see the same input. Max abs error ~3e-6 across the documented [-1e7, 1e7] range. Avoids @llvm.sin / @llvm.cos which LLVM scalarizes to per-lane libm sinf/cosf on every supported architecture. Pair-return form deferred — Eä has no precedent for multi-return intrinsics, and the established single-return pattern composes cleanly with the rest of the transcendental family. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/codegen/mod.rs | 2 + src/codegen/simd.rs | 4 + src/codegen/simd_sin_cos_approx.rs | 296 +++++++++++++++++++++++++++++ src/typeck/intrinsics.rs | 48 +++++ 4 files changed, 350 insertions(+) create mode 100644 src/codegen/simd_sin_cos_approx.rs diff --git a/src/codegen/mod.rs b/src/codegen/mod.rs index e160439..30d9ba4 100644 --- a/src/codegen/mod.rs +++ b/src/codegen/mod.rs @@ -37,6 +37,8 @@ mod simd_pack_unsigned; #[cfg(feature = "llvm")] mod simd_saturating; #[cfg(feature = "llvm")] +mod simd_sin_cos_approx; +#[cfg(feature = "llvm")] mod simd_tanh_approx; #[cfg(feature = "llvm")] mod simd_util; diff --git a/src/codegen/simd.rs b/src/codegen/simd.rs index 71fd3b0..d21b772 100644 --- a/src/codegen/simd.rs +++ b/src/codegen/simd.rs @@ -37,6 +37,8 @@ impl<'ctx> CodeGenerator<'ctx> { | "exp_poly_f32" | "tanh_approx_f32" | "log_approx_f32" + | "sin_approx_f32" + | "cos_approx_f32" | "reduce_add" | "reduce_add_fast" | "reduce_max" @@ -262,6 +264,8 @@ impl<'ctx> CodeGenerator<'ctx> { "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), + "sin_approx_f32" => self.compile_sin_approx_f32(args, function), + "cos_approx_f32" => self.compile_cos_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_sin_cos_approx.rs b/src/codegen/simd_sin_cos_approx.rs new file mode 100644 index 0000000..dcc826f --- /dev/null +++ b/src/codegen/simd_sin_cos_approx.rs @@ -0,0 +1,296 @@ +//! Polynomial-based vector sin and cos for f32 vectors. +//! +//! Both intrinsics share a common range-reduction + polynomial core: +//! reduce input mod π/2 via a 3-piece Cody-Waite split, compute both sin +//! and cos polynomials over the reduced argument in `[-π/4, π/4]`, then +//! blend by quadrant. `cos_approx_f32` reuses the same machinery with +//! `q += 1` — a precision-free integer shift that expresses the +//! mathematical identity `cos(x) = sin(x + π/2)`. Avoids `@llvm.sin` / +//! `@llvm.cos`, which LLVM scalarizes to per-lane libm `sinf` / `cosf` +//! on every supported architecture. +//! +//! Defined input range: [-1e7, 1e7] radians. Beyond this, the 3-piece +//! Cody-Waite split loses bits in the `q · π/2` subtraction and the +//! reduced argument's precision degrades. Maximum absolute error within +//! the defined range: ~3e-6. + +use inkwell::IntPredicate; +use inkwell::values::{BasicValueEnum, FunctionValue, IntValue, VectorValue}; + +use crate::ast::Expr; +use crate::error::CompileError; + +use super::CodeGenerator; + +// 2/π — input multiplier for quadrant rounding. +const TWO_OVER_PI: f32 = std::f32::consts::FRAC_2_PI; + +// 2-piece Cody-Waite split of π/2. PI_2_HI is the f32 closest to π/2 +// (slightly larger than the true value); PI_2_LO is the negative +// residual. Sum reproduces π/2 to within ~2e-15 in f64. With FMA in +// the codegen, the effective reduction precision is ~47 bits — enough +// for inputs up to ~1e7 radians, our documented range. +// +// Earlier drafts tried a 3-piece Sleef-style split derived by halving +// Sleef's PI_A_F/B_F/C_F constants, but the values I wrote down didn't +// actually sum to π/2 (off by ~6e-8). The 2-piece Eigen-style split +// here is precision-exact and simpler. +#[allow(clippy::excessive_precision)] +const PI_2_HI: f32 = 1.5707963705062866; +#[allow(clippy::excessive_precision)] +const PI_2_LO: f32 = -4.3711388e-8; + +// Polynomial coefficients for the reduced argument d' ∈ [-π/4, π/4]. +// +// sin(d') ≈ d' + d'³·(-1/6) + d'⁵·(1/120) + d'⁷·(-1/5040) +// = d' · (1 + s·(SIN_C1 + s·(SIN_C2 + s·SIN_C3))) where s = d'² +// +// cos(d') ≈ 1 + d'²·(-1/2) + d'⁴·(1/24) + d'⁶·(-1/720) + d'⁸·(1/40320) +// = Horner(COS_C0, COS_C1, COS_C2, COS_C3, COS_C4; s) — degree 4 in s +// +// Truncation error within d' ∈ [-π/4, π/4]: +// sin: |d'|⁹/9! ≤ (π/4)⁹/362880 ≈ 3.4e-7 → well within 3e-6 target +// cos: |d'|¹⁰/10! ≤ (π/4)¹⁰/3628800 ≈ 2.6e-8 → effectively zero +// +// Coefficients are Taylor values (well within the 3e-6 target without +// extra minimax tuning at this range). +#[allow(clippy::excessive_precision)] +mod coeffs { + pub(super) const SIN_C1: f32 = -1.666_666_7e-1; // -1/6 + pub(super) const SIN_C2: f32 = 8.333_333_3e-3; // 1/120 + pub(super) const SIN_C3: f32 = -1.984_127_0e-4; // -1/5040 + + pub(super) const COS_C0: f32 = 1.0; + pub(super) const COS_C1: f32 = -5.0e-1; // -1/2 + pub(super) const COS_C2: f32 = 4.166_666_7e-2; // 1/24 + pub(super) const COS_C3: f32 = -1.388_888_9e-3; // -1/720 + pub(super) const COS_C4: f32 = 2.480_158_7e-5; // 1/40320 +} +use coeffs::*; + +impl<'ctx> CodeGenerator<'ctx> { + /// Compile `sin_approx_f32(v: f32xN) -> f32xN`. + pub(super) fn compile_sin_approx_f32( + &mut self, + args: &[Expr], + function: FunctionValue<'ctx>, + ) -> crate::error::Result> { + self.compile_sin_cos_core(args, function, false) + } + + /// Compile `cos_approx_f32(v: f32xN) -> f32xN`. + pub(super) fn compile_cos_approx_f32( + &mut self, + args: &[Expr], + function: FunctionValue<'ctx>, + ) -> crate::error::Result> { + self.compile_sin_cos_core(args, function, true) + } + + /// Shared core. `want_cos = true` adds 1 to the quadrant index `q` + /// before the blend — mathematically equivalent to `cos(x) = sin(x + π/2)` + /// but precision-free (integer shift) rather than loss-prone (adding + /// π/2 to large floating-point inputs). + fn compile_sin_cos_core( + &mut self, + args: &[Expr], + function: FunctionValue<'ctx>, + want_cos: bool, + ) -> crate::error::Result> { + let val = self.compile_expr(&args[0], function)?; + let v = match val { + BasicValueEnum::VectorValue(vv) => vv, + _ => { + let name = if want_cos { + "cos_approx_f32" + } else { + "sin_approx_f32" + }; + return Err(CompileError::codegen_error(format!( + "{name} expects f32 vector, got scalar; use {} for scalar libm-precision", + if want_cos { "cos()" } else { "sin()" } + ))); + } + }; + 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() { + let name = if want_cos { + "cos_approx_f32" + } else { + "sin_approx_f32" + }; + return Err(CompileError::codegen_error(format!( + "{name} expects f32 element type" + ))); + } + let width = vec_ty.get_size(); + let i32_vec_ty = self.context.i32_type().vec_type(width); + + // 1. Quadrant: qf = nearbyint(v · 2/π); q = fptosi(qf) + let two_over_pi = self.splat_f32_const_sc(TWO_OVER_PI, width)?; + let v_scaled = self + .builder + .build_float_mul(v, two_over_pi, "sc_v_scaled") + .map_err(|e| CompileError::codegen_error(e.to_string()))?; + + let nearbyint_name = format!("llvm.nearbyint.v{width}f32"); + let nearbyint_fn = self + .module + .get_function(&nearbyint_name) + .unwrap_or_else(|| { + let fn_ty = vec_ty.fn_type(&[vec_ty.into()], false); + self.module.add_function(&nearbyint_name, fn_ty, None) + }); + let qf = self + .builder + .build_call(nearbyint_fn, &[v_scaled.into()], "sc_qf") + .map_err(|e| CompileError::codegen_error(e.to_string()))? + .try_as_basic_value() + .basic() + .ok_or_else(|| CompileError::codegen_error("nearbyint returned no value"))? + .into_vector_value(); + let q = self + .builder + .build_float_to_signed_int(qf, i32_vec_ty, "sc_q") + .map_err(|e| CompileError::codegen_error(e.to_string()))?; + + // 2. Cody-Waite range reduction: d' = v - qf · PI_2_HI - qf · PI_2_LO + // Each fma(neg_qf, c, accumulator) preserves precision; the second + // step adds the residual that PI_2_HI can't capture. + let neg_qf = self + .builder + .build_float_neg(qf, "sc_neg_qf") + .map_err(|e| CompileError::codegen_error(e.to_string()))?; + let pi_2_hi = self.splat_f32_const_sc(PI_2_HI, width)?; + let pi_2_lo = self.splat_f32_const_sc(PI_2_LO, width)?; + let d = self.fma_sc(neg_qf, pi_2_hi, v, "sc_d1", width)?; + let d = self.fma_sc(neg_qf, pi_2_lo, d, "sc_d", width)?; + // d ∈ [-π/4, π/4] + + // 3. s = d² (reused by both polynomials) + let s = self + .builder + .build_float_mul(d, d, "sc_s") + .map_err(|e| CompileError::codegen_error(e.to_string()))?; + + // 4. sin polynomial: sin(d) = d · (1 + s·(C1 + s·(C2 + s·C3))) + // = d + d·s·u where u = C1 + s·C2 + s²·C3 + let sin_c1 = self.splat_f32_const_sc(SIN_C1, width)?; + let sin_c2 = self.splat_f32_const_sc(SIN_C2, width)?; + let sin_c3 = self.splat_f32_const_sc(SIN_C3, width)?; + let sin_u = self.fma_sc(sin_c3, s, sin_c2, "sc_sin_u1", width)?; + let sin_u = self.fma_sc(sin_u, s, sin_c1, "sc_sin_u", width)?; + let d_times_s = self + .builder + .build_float_mul(d, s, "sc_d_s") + .map_err(|e| CompileError::codegen_error(e.to_string()))?; + let sin_val = self.fma_sc(d_times_s, sin_u, d, "sc_sin_val", width)?; + + // 5. cos polynomial: cos(d) = Horner(COS_C4, COS_C3, COS_C2, COS_C1, COS_C0; s) + // = ((((C4·s + C3)·s + C2)·s + C1)·s + C0 + let cos_c0 = self.splat_f32_const_sc(COS_C0, width)?; + let cos_c1 = self.splat_f32_const_sc(COS_C1, width)?; + let cos_c2 = self.splat_f32_const_sc(COS_C2, width)?; + let cos_c3 = self.splat_f32_const_sc(COS_C3, width)?; + let cos_c4 = self.splat_f32_const_sc(COS_C4, width)?; + let cos_val = self.fma_sc(cos_c4, s, cos_c3, "sc_cos_p1", width)?; + let cos_val = self.fma_sc(cos_val, s, cos_c2, "sc_cos_p2", width)?; + let cos_val = self.fma_sc(cos_val, s, cos_c1, "sc_cos_p3", width)?; + let cos_val = self.fma_sc(cos_val, s, cos_c0, "sc_cos_val", width)?; + + // 6. Quadrant blend. + // For sin_approx: k = q. For cos_approx: k = q + 1 (precision-free phase shift). + // swap = (k & 1) != 0 → pick cos_val instead of sin_val + // negate = (k & 2) != 0 → flip result sign + let k = if want_cos { + let one = self.splat_i32_const_sc(1, width)?; + self.builder + .build_int_add(q, one, "sc_k_cos") + .map_err(|e| CompileError::codegen_error(e.to_string()))? + } else { + q + }; + + let mask_one = self.splat_i32_const_sc(1, width)?; + let mask_two = self.splat_i32_const_sc(2, width)?; + let zero = self.splat_i32_const_sc(0, width)?; + let k_and_1 = self + .builder + .build_and(k, mask_one, "sc_k_and_1") + .map_err(|e| CompileError::codegen_error(e.to_string()))?; + let swap_mask = self + .builder + .build_int_compare(IntPredicate::NE, k_and_1, zero, "sc_swap_mask") + .map_err(|e| CompileError::codegen_error(e.to_string()))?; + let k_and_2 = self + .builder + .build_and(k, mask_two, "sc_k_and_2") + .map_err(|e| CompileError::codegen_error(e.to_string()))?; + let negate_mask = self + .builder + .build_int_compare(IntPredicate::NE, k_and_2, zero, "sc_neg_mask") + .map_err(|e| CompileError::codegen_error(e.to_string()))?; + + // result = swap ? cos_val : sin_val + let blended = self + .builder + .build_select(swap_mask, cos_val, sin_val, "sc_blended") + .map_err(|e| CompileError::codegen_error(e.to_string()))? + .into_vector_value(); + // result = negate ? -blended : blended + let neg_blended = self + .builder + .build_float_neg(blended, "sc_neg_blended") + .map_err(|e| CompileError::codegen_error(e.to_string()))?; + let result = self + .builder + .build_select(negate_mask, neg_blended, blended, "sc_result") + .map_err(|e| CompileError::codegen_error(e.to_string()))?; + + Ok(result) + } + + fn splat_f32_const_sc( + &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_sc( + &self, + value: i32, + width: u32, + ) -> crate::error::Result> { + let scalar: IntValue<'ctx> = self.context.i32_type().const_int(value as u64, true); + self.build_splat(BasicValueEnum::IntValue(scalar), width) + } + + fn fma_sc( + &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 dbe4882..adac94c 100644 --- a/src/typeck/intrinsics.rs +++ b/src/typeck/intrinsics.rs @@ -47,6 +47,12 @@ impl TypeChecker { "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)), + "sin_approx_f32" => { + Some(self.check_sin_cos_approx_f32("sin_approx_f32", args, locals, span)) + } + "cos_approx_f32" => { + Some(self.check_sin_cos_approx_f32("cos_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)) } @@ -468,6 +474,48 @@ impl TypeChecker { } } + /// Type-check `sin_approx_f32(v: f32xN) -> f32xN` and + /// `cos_approx_f32(v: f32xN) -> f32xN`. Same f32-vector-only shape as + /// `exp_poly_f32` — scalar / f64 / f16 / integer all rejected. Error + /// text varies by intrinsic name to point at the right libm fallback. + fn check_sin_cos_approx_f32( + &self, + name: &str, + args: &[Expr], + locals: &HashMap, + span: &Span, + ) -> crate::error::Result { + if args.len() != 1 { + return Err(CompileError::type_error( + format!("{name} expects 1 argument, got {}", args.len()), + span.clone(), + )); + } + let libm_fallback = if name == "cos_approx_f32" { + "cos" + } else { + "sin" + }; + 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!("{name} expects f32 element type, got {elem}"), + span.clone(), + )), + Type::F32 | Type::FloatLiteral => Err(CompileError::type_error( + format!( + "{name} expects f32 vector, got scalar; use {libm_fallback}() for scalar libm-precision" + ), + span.clone(), + )), + _ => Err(CompileError::type_error( + format!("{name} expects float vector, got {t}"), + span.clone(), + )), + } + } + fn check_prefetch( &self, args: &[Expr], From 3084b6f4d33ae5aac61373bb38c0f41da8869157 Mon Sep 17 00:00:00 2001 From: Peter Lukka Date: Tue, 19 May 2026 06:16:53 +0000 Subject: [PATCH 2/3] test(intrinsic): sin/cos_approx_f32 boundary, Pythagorean, IR guards MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 19 tests covering: - at_zero (sin + cos): smoke that sin(0)=0 and cos(0)=1. - boundary_points (sin × {x4, x8} + cos × {x4, x8}): vs libm sinf/cosf across {-2π … 2π} plus moderate (10, 100) and edge (1e6) inputs within the documented [-1e7, 1e7] range. Abs error ≤ 3e-6. - random (sin × {x4, x8} + cos × x4): 256 LCG samples over [-10π, 10π] exercising range reduction over multiple quadrants. - pythagorean_identity: pins sin²(x) + cos²(x) ≈ 1 to 6e-6 across 4 inputs. A regression in either intrinsic's q-handling would surface here even if individual tolerance still passes. - does_not_emit_llvm_sin_cos: IR guard against @llvm.sin/@llvm.cos delegation (LLVM scalarizes both to libm). - emits_expected_pattern (sin × {x4, x8}): IR has ≥6 FMAs + @llvm.nearbyint + fptosi → . - Six typeck rejections (3 per intrinsic): scalar f32, f64x2, integer vector; arity for sin. An earlier draft of the Cody-Waite constants had a transcription error (3-piece Sleef-derived values that didn't sum to π/2, off by 6e-8). The boundary_points test caught it loudly — cos(1e6) showed 0.014 error vs the 3e-6 tolerance. Fix was the 2-piece Eigen-style split now in the source; lesson is to verify constants against their claimed identity, not transcribe by eye. Co-Authored-By: Claude Opus 4.7 (1M context) --- tests/phase14_sin_cos_approx.rs | 509 ++++++++++++++++++++++++++++++++ 1 file changed, 509 insertions(+) create mode 100644 tests/phase14_sin_cos_approx.rs diff --git a/tests/phase14_sin_cos_approx.rs b/tests/phase14_sin_cos_approx.rs new file mode 100644 index 0000000..4cf11c5 --- /dev/null +++ b/tests/phase14_sin_cos_approx.rs @@ -0,0 +1,509 @@ +#[cfg(feature = "llvm")] +mod tests { + use ea_compiler::{CompileOptions, OutputMode}; + use std::process::Command; + use tempfile::TempDir; + + /// Smoke test: sin_approx_f32(splat(0.0)) should give all 0.0s; the + /// range reduction yields q=0, d'=0, and the quadrant-0 branch picks + /// sin_val = d' = 0. + #[test] + fn test_sin_approx_f32x4_at_zero() { + let ea = r#" + export func k(out: *mut f32) { + let z: f32x4 = splat(0.0) + let r: f32x4 = sin_approx_f32(z) + store(out, 0, r) + } + "#; + run_smoke(ea, "0\n0\n0\n0"); + } + + /// Smoke test: cos_approx_f32(splat(0.0)) should give all 1.0s; the + /// range reduction yields q=0, the cos-shift gives k=1, swap=true → + /// pick cos_val = 1.0. + #[test] + fn test_cos_approx_f32x4_at_zero() { + let ea = r#" + export func k(out: *mut f32) { + let z: f32x4 = splat(0.0) + let r: f32x4 = cos_approx_f32(z) + store(out, 0, r) + } + "#; + run_smoke(ea, "1\n1\n1\n1"); + } + + fn run_smoke(ea: &str, expected: &str) { + let c = r#" + #include + extern void k(float *out); + int main(void) { + float out[4] = {-99, -99, -99, -99}; + 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(), expected); + } + + /// Compile a kernel that runs `intrinsic` over `inputs` lane-by-lane, + /// link with C harness that calls `ref_fn`, assert absolute error ≤ 3e-6. + /// + /// Absolute error tolerance because sin/cos pass through zero at the + /// quadrant boundaries — relative error blows up there. + fn accuracy_test_impl(inputs: &[f32], vector_type: &str, intrinsic: &str, ref_fn: &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} = {intrinsic}(v) + store(output, i, r) + i = i + {lanes} + }} + }} + "# + ); + + let mut padded = inputs.to_vec(); + while !padded.len().is_multiple_of(lanes) { + padded.push(0.0); + } + 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 = {ref_fn}(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) + ); + } + + fn boundary_points() -> Vec { + use std::f32::consts::PI; + vec![ + -2.0 * PI, + -3.0 * PI / 2.0, + -PI, + -PI / 2.0, + -PI / 4.0, + -0.1, + 0.0, + 0.1, + PI / 4.0, + PI / 2.0, + PI, + 3.0 * PI / 2.0, + 2.0 * PI, + // Moderate range + 10.0, + 100.0, + // Within the documented defined range [-1e7, 1e7] + 1.0e6, + ] + } + + #[test] + fn test_sin_approx_f32x4_boundary_points() { + accuracy_test_impl(&boundary_points(), "f32x4", "sin_approx_f32", "sinf"); + } + + #[test] + fn test_cos_approx_f32x4_boundary_points() { + accuracy_test_impl(&boundary_points(), "f32x4", "cos_approx_f32", "cosf"); + } + + /// Pythagorean identity: sin²(x) + cos²(x) ≈ 1. + /// Pin-tests that the two intrinsics share consistent range reduction + /// and quadrant logic — a regression in either's q-handling would + /// surface here even if individual tolerance still passes. + #[test] + fn test_sin_cos_pythagorean_identity() { + let ea = r#" + export func k(input: *f32, output: *mut f32) { + let x: f32x4 = load(input, 0) + let s: f32x4 = sin_approx_f32(x) + let c: f32x4 = cos_approx_f32(x) + let result: f32x4 = fma(s, s, c .* c) + store(output, 0, result) + } + "#; + let c = r#" + #include + #include + extern void k(const float *input, float *output); + int main(void) { + float in[4] = {0.1f, 1.2f, 3.4f, 5.6f}; + float out[4] = {0}; + k(in, out); + for (int i = 0; i < 4; ++i) { + float err = fabsf(out[i] - 1.0f); + if (err > 6.0e-6f) { + printf("FAIL i=%d in=%g s²+c²=%g err=%g\n", i, in[i], out[i], 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"); + } + + #[test] + fn test_sin_approx_f32x4_random() { + let mut points = Vec::new(); + let mut state: u32 = 0xDEAD_BEEF; + for _ in 0..256 { + state = state.wrapping_mul(1664525).wrapping_add(1013904223); + // Sample [-10π, 10π] — exercises range reduction over multiple + // quadrants without hitting the documented range boundary. + let u = state as f32 / u32::MAX as f32; + let x = (u - 0.5) * 20.0 * std::f32::consts::PI; + points.push(x); + } + accuracy_test_impl(&points, "f32x4", "sin_approx_f32", "sinf"); + } + + #[test] + fn test_cos_approx_f32x4_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 x = (u - 0.5) * 20.0 * std::f32::consts::PI; + points.push(x); + } + accuracy_test_impl(&points, "f32x4", "cos_approx_f32", "cosf"); + } + + #[cfg(target_arch = "x86_64")] + #[test] + fn test_sin_approx_f32x8_boundary_points() { + accuracy_test_impl(&boundary_points(), "f32x8", "sin_approx_f32", "sinf"); + } + + #[cfg(target_arch = "x86_64")] + #[test] + fn test_cos_approx_f32x8_boundary_points() { + accuracy_test_impl(&boundary_points(), "f32x8", "cos_approx_f32", "cosf"); + } + + #[cfg(target_arch = "x86_64")] + #[test] + fn test_sin_approx_f32x8_random() { + let mut points = Vec::new(); + let mut state: u32 = 0xFACE_BEEF; + for _ in 0..256 { + state = state.wrapping_mul(1664525).wrapping_add(1013904223); + let u = state as f32 / u32::MAX as f32; + let x = (u - 0.5) * 20.0 * std::f32::consts::PI; + points.push(x); + } + accuracy_test_impl(&points, "f32x8", "sin_approx_f32", "sinf"); + } + + /// CRITICAL regression guard: future "simplifications" that delegate + /// to @llvm.sin / @llvm.cos would scalarize to per-lane libm sinf/cosf + /// and silently undo the vectorization that motivates the intrinsics. + #[cfg(target_arch = "x86_64")] + #[test] + fn test_sin_cos_approx_does_not_emit_llvm_sin_cos() { + let src = r#" + export func sk(v: f32x8) -> f32x8 { return sin_approx_f32(v) } + export func ck(v: f32x8) -> f32x8 { return cos_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.sin") && !ir.contains("@llvm.cos"), + "sin/cos_approx_f32 must NOT lower to @llvm.sin/@llvm.cos:\n{ir}" + ); + } + + /// Confirm the range-reduction + polynomial pattern is emitted for f32x8. + #[cfg(target_arch = "x86_64")] + #[test] + fn test_sin_approx_f32x8_emits_expected_pattern() { + let src = r#" + export func k(v: f32x8) -> f32x8 { + return sin_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 for sin (3 FMA) + Horner for cos (4 FMA) + Cody-Waite (≥2 FMA) + let fma_count = ir.matches("@llvm.fma.v8f32").count(); + assert!( + fma_count >= 6, + "expected at least 6 @llvm.fma.v8f32 calls; got {fma_count}\nIR:\n{ir}" + ); + assert!( + ir.contains("@llvm.nearbyint"), + "expected @llvm.nearbyint for quadrant rounding:\n{ir}" + ); + assert!( + ir.contains("fptosi <8 x float>") && ir.contains("to <8 x i32>"), + "expected fptosi <8 x float> → <8 x i32> for quadrant index:\n{ir}" + ); + } + + #[test] + fn test_sin_approx_f32x4_emits_expected_pattern() { + let src = r#" + export func k(v: f32x4) -> f32x4 { + return sin_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.sin"), "no @llvm.sin expected:\n{ir}"); + let fma_count = ir.matches("@llvm.fma.v4f32").count(); + assert!( + fma_count >= 6, + "expected at least 6 @llvm.fma.v4f32 calls; got {fma_count}\nIR:\n{ir}" + ); + } + + // --- Typeck rejections (sin) --- + + #[test] + fn test_sin_approx_f32_rejects_scalar_f32() { + let src = r#" + export func k(x: f32) -> f32 { + return sin_approx_f32(x) + } + "#; + let err = compile_err(src); + assert!( + err.contains("f32 vector"), + "error must mention f32 vector requirement, got: {err}" + ); + } + + #[test] + fn test_sin_approx_f32_rejects_f64_vector() { + let src = r#" + export func k(v: f64x2) -> f64x2 { + return sin_approx_f32(v) + } + "#; + let err = compile_err(src); + assert!( + err.contains("f32") && err.contains("element type"), + "error must mention f32 element type, got: {err}" + ); + } + + #[test] + fn test_sin_approx_f32_rejects_integer_vector() { + let src = r#" + export func k(v: i32x4) -> i32x4 { + return sin_approx_f32(v) + } + "#; + let err = compile_err(src); + assert!( + err.contains("float") || err.contains("f32"), + "error must mention float requirement, got: {err}" + ); + } + + #[test] + fn test_sin_approx_f32_rejects_wrong_arity() { + let src = r#" + export func k(a: f32x4, b: f32x4) -> f32x4 { + return sin_approx_f32(a, b) + } + "#; + let err = compile_err(src); + assert!( + err.contains("1 argument") || err.contains("expects 1"), + "error must mention arity, got: {err}" + ); + } + + // --- Typeck rejections (cos) — same shape as sin --- + + #[test] + fn test_cos_approx_f32_rejects_scalar_f32() { + let src = r#" + export func k(x: f32) -> f32 { + return cos_approx_f32(x) + } + "#; + let err = compile_err(src); + assert!(err.contains("f32 vector"), "got: {err}"); + } + + #[test] + fn test_cos_approx_f32_rejects_f64_vector() { + let src = r#" + export func k(v: f64x2) -> f64x2 { + return cos_approx_f32(v) + } + "#; + let err = compile_err(src); + assert!( + err.contains("f32") && err.contains("element type"), + "got: {err}" + ); + } + + fn compile_err(src: &str) -> String { + 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("expected compile failure"); + format!("{err}") + } +} From 880a51c9ce5400e6b7ea1f30a5dae1d572a2bbaf Mon Sep 17 00:00:00 2001 From: Peter Lukka Date: Tue, 19 May 2026 06:17:01 +0000 Subject: [PATCH 3/3] docs(v1.14.0): sin/cos_approx_f32 CHANGELOG + ROADMAP closes the trio CHANGELOG entry under v1.14.0 Added. ROADMAP Shipped entry + collapses the Future API consistency section: the v1.11.0-era trio (tanh, log, sin/cos plus u16x32 and wider wmul_u64) is fully closed out by v1.14.0. The f32 transcendental approximation family is feature-complete; new entries land here as real consumers surface them. Co-Authored-By: Claude Opus 4.7 (1M context) --- CHANGELOG.md | 1 + ROADMAP.md | 7 ++++++- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6194e78..22a8c27 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ - `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. +- `sin_approx_f32(v: f32xN) -> f32xN` and `cos_approx_f32(v: f32xN) -> f32xN` — polynomial sin/cos via shared mod-π/2 range reduction (2-piece Cody-Waite + FMA), Taylor-truncated polynomials over `d' ∈ [-π/4, π/4]`, and a quadrant blend. The `cos` variant reuses the same core via `q += 1` (precision-free integer shift expressing `cos(x) = sin(x + π/2)`). Max abs error ~3e-6 across `[-1e7, 1e7]`. Closes the original "Future API consistency" trio (`tanh_approx_f32` v1.14.0, `log_approx_f32` v1.14.0, sin/cos v1.14.0); the f32 transcendental approximation family is feature-complete. ## v1.13.0 — 2026-05-15 — ea bench + first aarch64 baselines + Specification umbrella diff --git a/ROADMAP.md b/ROADMAP.md index 8b52033..dac1f20 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -32,6 +32,12 @@ Forward-looking notes. Ordered by leverage, not by effort. `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`. +### sin_approx_f32 + cos_approx_f32 + +`sin_approx_f32(v: f32xN) -> f32xN` and `cos_approx_f32(v: f32xN) -> f32xN`. Shared core: reduce mod π/2 via 2-piece Cody-Waite split with FMA-preserved precision, compute both sin and cos polynomials over the reduced argument `d' ∈ [-π/4, π/4]` (degree 3 in `s = d'²` for sin, degree 4 for cos), then quadrant-blend at the end. The `cos` variant reuses the same core with `q += 1` — a precision-free integer shift expressing the mathematical identity `cos(x) = sin(x + π/2)`. Max abs error ~3e-6 across `[-1e7, 1e7]`. Closes the original "Future API consistency" trio entry. Spec at `docs/superpowers/specs/2026-05-19-sin-cos-approx-f32-design.md`. + +The roadmap entry was titled "sin_cos_approx_f32" suggesting a pair-return, but Eä has no precedent for multi-return intrinsics. Shipping two separate intrinsics matches the established pattern (`exp_poly_f32` / `tanh_approx_f32` / `log_approx_f32`) and lets LLVM CSE handle any caller-side range-reduction redundancy. + ## Shipped in v1.12.0 (2026-05-13) - **Deprecation-warning infrastructure** + `docs/migrations/` directory + `cargo public-api` CI gate (PR #6). @@ -64,7 +70,6 @@ Today the language spec is spread across `docs/src/reference/*.md` (types, intri ## Future API consistency -- **`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