fix(math): C23 sinpi/cospi/tanpi exact at integer/half-integer args (glibc parity)

sinpi/cospi/tanpi (+ f32) were implemented naively as sin(x*PI)/cos(x*PI)/
tan(x*PI), violating the entire C23 contract for the pi-scaled trig functions.
A differential probe vs host glibc found 13 divergences:

- sinpi(1.0) returned 1.2e-16 instead of EXACT 0; sinpi(1e20) returned -0.39
  (catastrophic — argument reduction lost all precision) instead of 0.
- cospi(0.5)/cospi(1.5)/... returned ~1e-16 instead of EXACT +0.
- tanpi(0.5) returned a huge finite (~1.6e16) instead of +inf + FE_DIVBYZERO;
  tanpi(1.0) lost its -0 sign.

Reimplemented via the exact identity f(n+r) with n=round(x), r=x-n in
[-0.5,0.5] (exact by Sterbenz): sinpi(x)=(-1)^n*sin(PI*r), cospi(x)=
(-1)^n*cos(PI*r), tanpi=sinpi/cospi (the (-1)^n cancels; the half-integer
pole becomes ±1/±0 which auto-raises FE_DIVBYZERO → ±inf). Integer and
half-integer args, huge args (>=2^53 are even integers), inf (→NaN+INVALID),
and NaN are all special-cased for exactness; irrational args use the in-tree
sin/cos kernels (<=4 ULP vs glibc). acospi/asinpi/atanpi/atan2pi were already
correct (plain f(x)/PI) and are unchanged.

New gate conformance_diff_pi_trig: live differential (glibc via dlsym to
bypass fl's no_mangle interposition; hardware FP flags) over 41 f64 + 17 f32
inputs — specials exact, irrational <=4 ULP, exception flags exact. 116
math_abi unit tests pass.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: Dicklesworthstone

crates/frankenlibc-abi/src/math_abi.rs

@@ -3578,12 +3578,68 @@ pub unsafe extern "C" fn atan2pif128(x: f64, y: f64) -> f64 {
     unsafe { atan2pi(x, y) }
 }
 #[cfg_attr(not(debug_assertions), unsafe(no_mangle))]
+/// Re-raise FE_INVALID on the cold pi-function domain path (x = ±inf).
+#[inline]
+fn pi_fn_raise_invalid_f64() {
+    let _ = core::hint::black_box(core::hint::black_box(0.0_f64) / core::hint::black_box(0.0_f64));
+}
+#[inline]
+fn pi_fn_raise_invalid_f32() {
+    let _ = core::hint::black_box(core::hint::black_box(0.0_f32) / core::hint::black_box(0.0_f32));
+}
+
+// C23 pi-scaled trig: sinpi(x)=sin(pi*x), cospi(x)=cos(pi*x), tanpi(x)=tan(pi*x),
+// computed via the exact identity f(n+r) with n=round(x), r=x-n in [-0.5,0.5]
+// (exact by Sterbenz for |x|<2^53). This yields EXACT results at integer and
+// half-integer arguments (sinpi(1)=+0, cospi(0.5)=+0, tanpi(0.5)=+inf) and
+// stays correct for huge arguments where the naive sin(x*PI) loses all
+// precision. |x|>=2^53 is always an even integer: sinpi=±0, cospi=1, tanpi=±0.
 pub unsafe extern "C" fn cospi(x: f64) -> f64 {
-    unsafe { cos(x * std::f64::consts::PI) }
+    if x.is_nan() {
+        return x;
+    }
+    if x.is_infinite() {
+        pi_fn_raise_invalid_f64();
+        return f64::NAN;
+    }
+    if x.abs() >= 9007199254740992.0 {
+        return 1.0; // even integer → cos(pi*even)=+1
+    }
+    let n = x.round();
+    let r = x - n;
+    let n_odd = (n as i64) & 1 != 0;
+    if r == 0.0 {
+        return if n_odd { -1.0 } else { 1.0 };
+    }
+    if r == 0.5 || r == -0.5 {
+        return 0.0; // cos at odd multiple of pi/2 is +0
+    }
+    let c = frankenlibc_core::math::cos(r * std::f64::consts::PI);
+    if n_odd { -c } else { c }
 }
 #[cfg_attr(not(debug_assertions), unsafe(no_mangle))]
 pub unsafe extern "C" fn cospif(x: f32) -> f32 {
-    unsafe { cosf(x * std::f32::consts::PI) }
+    if x.is_nan() {
+        return x;
+    }
+    if x.is_infinite() {
+        pi_fn_raise_invalid_f32();
+        return f32::NAN;
+    }
+    if x.abs() >= 16777216.0 {
+        return 1.0;
+    }
+    let n = x.round();
+    let r = x - n;
+    let n_odd = (n as i64) & 1 != 0;
+    if r == 0.0 {
+        return if n_odd { -1.0 } else { 1.0 };
+    }
+    if r == 0.5 || r == -0.5 {
+        return 0.0;
+    }
+    let c = frankenlibc_core::math::cosf(r * std::f32::consts::PI);
+    if n_odd { -c } else { c }
 }
 #[cfg_attr(not(debug_assertions), unsafe(no_mangle))]
 pub unsafe extern "C" fn cospil(x: f64) -> f64 {
@@ -3611,11 +3667,54 @@ pub unsafe extern "C" fn cospif128(x: f64) -> f64 {
 }
 #[cfg_attr(not(debug_assertions), unsafe(no_mangle))]
 pub unsafe extern "C" fn sinpi(x: f64) -> f64 {
-    unsafe { sin(x * std::f64::consts::PI) }
+    if x.is_nan() {
+        return x;
+    }
+    if x.is_infinite() {
+        pi_fn_raise_invalid_f64();
+        return f64::NAN;
+    }
+    if x.abs() >= 9007199254740992.0 {
+        // even integer → sin(pi*even)=±0 with sign of x
+        return if x.is_sign_negative() { -0.0 } else { 0.0 };
+    }
+    let n = x.round();
+    let r = x - n;
+    let n_odd = (n as i64) & 1 != 0;
+    if r == 0.0 {
+        return if x.is_sign_negative() { -0.0 } else { 0.0 };
+    }
+    if r == 0.5 || r == -0.5 {
+        let m = if r > 0.0 { 1.0 } else { -1.0 };
+        return if n_odd { -m } else { m };
+    }
+    let s = frankenlibc_core::math::sin(r * std::f64::consts::PI);
+    if n_odd { -s } else { s }
 }
 #[cfg_attr(not(debug_assertions), unsafe(no_mangle))]
 pub unsafe extern "C" fn sinpif(x: f32) -> f32 {
-    unsafe { sinf(x * std::f32::consts::PI) }
+    if x.is_nan() {
+        return x;
+    }
+    if x.is_infinite() {
+        pi_fn_raise_invalid_f32();
+        return f32::NAN;
+    }
+    if x.abs() >= 16777216.0 {
+        return if x.is_sign_negative() { -0.0 } else { 0.0 };
+    }
+    let n = x.round();
+    let r = x - n;
+    let n_odd = (n as i64) & 1 != 0;
+    if r == 0.0 {
+        return if x.is_sign_negative() { -0.0 } else { 0.0 };
+    }
+    if r == 0.5 || r == -0.5 {
+        let m = if r > 0.0 { 1.0 } else { -1.0 };
+        return if n_odd { -m } else { m };
+    }
+    let s = frankenlibc_core::math::sinf(r * std::f32::consts::PI);
+    if n_odd { -s } else { s }
 }
 #[cfg_attr(not(debug_assertions), unsafe(no_mangle))]
 pub unsafe extern "C" fn sinpil(x: f64) -> f64 {
@@ -3643,11 +3742,32 @@ pub unsafe extern "C" fn sinpif128(x: f64) -> f64 {
 }
 #[cfg_attr(not(debug_assertions), unsafe(no_mangle))]
 pub unsafe extern "C" fn tanpi(x: f64) -> f64 {
-    unsafe { tan(x * std::f64::consts::PI) }
+    if x.is_nan() {
+        return x;
+    }
+    if x.is_infinite() {
+        pi_fn_raise_invalid_f64();
+        return f64::NAN;
+    }
+    // tanpi = sinpi/cospi: the (-1)^n factors cancel, the half-integer pole
+    // becomes ±1/±0 (auto-raising FE_DIVBYZERO → ±inf), and the integer zero
+    // becomes ±0/±1 with the correct sign.
+    let s = unsafe { sinpi(x) };
+    let c = unsafe { cospi(x) };
+    s / c
 }
 #[cfg_attr(not(debug_assertions), unsafe(no_mangle))]
 pub unsafe extern "C" fn tanpif(x: f32) -> f32 {
-    unsafe { tanf(x * std::f32::consts::PI) }
+    if x.is_nan() {
+        return x;
+    }
+    if x.is_infinite() {
+        pi_fn_raise_invalid_f32();
+        return f32::NAN;
+    }
+    let s = unsafe { sinpif(x) };
+    let c = unsafe { cospif(x) };
+    s / c
 }
 #[cfg_attr(not(debug_assertions), unsafe(no_mangle))]
 pub unsafe extern "C" fn tanpil(x: f64) -> f64 {

crates/frankenlibc-abi/tests/conformance_diff_pi_trig.rs

@@ -0,0 +1,144 @@
+//! Differential gate for the C23 pi-scaled trig functions sinpi/cospi/tanpi
+//! (+ f32) vs host glibc.
+//!
+//! These must produce EXACT results at integer and half-integer arguments
+//! (sinpi(n)=±0, cospi(n+0.5)=+0, cospi(n)=±1, tanpi(n+0.5)=±inf with
+//! FE_DIVBYZERO) and stay correct for huge arguments — properties the naive
+//! `sin(x*PI)` formulation violates. fl is called via Rust paths; glibc is
+//! reached through `dlsym` on libm.so.6 so the fn pointer bypasses fl's
+//! no_mangle interposition of the same symbol. FP exception flags are hardware
+//! (MXCSR), read directly with fetestexcept — no interposition concern.
+#![cfg(target_os = "linux")]
+#![allow(unsafe_code)]
+
+use frankenlibc_abi::math_abi as fl;
+use std::ffi::{c_char, c_int, c_void};
+
+const RTLD_NOW: c_int = 2;
+const HARD: c_int = 0x1D; // INVALID|DIVBYZERO|OVERFLOW|UNDERFLOW
+
+unsafe extern "C" {
+    fn dlopen(filename: *const c_char, flag: c_int) -> *mut c_void;
+    fn dlsym(handle: *mut c_void, symbol: *const c_char) -> *mut c_void;
+    fn feclearexcept(e: c_int) -> c_int;
+    fn fetestexcept(e: c_int) -> c_int;
+}
+
+fn libm() -> *mut c_void {
+    let h = unsafe { dlopen(c"libm.so.6".as_ptr(), RTLD_NOW) };
+    assert!(!h.is_null(), "dlopen libm failed");
+    h
+}
+fn sym(h: *mut c_void, name: &std::ffi::CStr) -> *mut c_void {
+    let p = unsafe { dlsym(h, name.as_ptr()) };
+    assert!(!p.is_null(), "missing libm symbol {name:?}");
+    p
+}
+
+fn ulp_ok_f64(a: f64, b: f64) -> bool {
+    // b is the glibc reference. Exact bit match for non-finite / zero (sign
+    // matters); <=4 ULP for finite nonzero (the math conformance contract).
+    if b.is_nan() {
+        return a.is_nan();
+    }
+    if !b.is_finite() || b == 0.0 {
+        return a.to_bits() == b.to_bits();
+    }
+    if a.is_nan() || a.is_sign_negative() != b.is_sign_negative() {
+        return false;
+    }
+    let ai = a.to_bits() as i64;
+    let bi = b.to_bits() as i64;
+    (ai - bi).unsigned_abs() <= 4
+}
+fn ulp_ok_f32(a: f32, b: f32) -> bool {
+    if b.is_nan() {
+        return a.is_nan();
+    }
+    if !b.is_finite() || b == 0.0 {
+        return a.to_bits() == b.to_bits();
+    }
+    if a.is_nan() || a.is_sign_negative() != b.is_sign_negative() {
+        return false;
+    }
+    let ai = a.to_bits() as i32;
+    let bi = b.to_bits() as i32;
+    (ai - bi).unsigned_abs() <= 4
+}
+
+const XS: &[f64] = &[
+    0.0, -0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, -0.5, -1.0, -1.5, -2.0, -2.5, -3.0, 0.25, 0.75, -0.25,
+    0.1, 0.2, 0.3, 0.4, 0.45, 0.7, 1.3, 2.7, -3.3, 4.25, 10.1, 10.25, -0.123, 100.5, 1000.25,
+    1.0 / 3.0, 1e15, 1e20, -1e20, 9007199254740992.0, 0.5 + 1e15,
+    f64::INFINITY, f64::NEG_INFINITY, f64::NAN,
+];
+const XSF: &[f32] = &[
+    0.0, -0.0, 0.5, 1.0, 1.5, 2.0, -0.5, -1.0, 0.25, 0.75, 0.1, 0.3, 0.7, 10.25, 16777216.0,
+    f32::INFINITY, f32::NAN,
+];
+
+#[test]
+fn pi_trig_matches_glibc() {
+    let h = libm();
+    let g_sinpi: extern "C" fn(f64) -> f64 = unsafe { core::mem::transmute(sym(h, c"sinpi")) };
+    let g_cospi: extern "C" fn(f64) -> f64 = unsafe { core::mem::transmute(sym(h, c"cospi")) };
+    let g_tanpi: extern "C" fn(f64) -> f64 = unsafe { core::mem::transmute(sym(h, c"tanpi")) };
+    let g_sinpif: extern "C" fn(f32) -> f32 = unsafe { core::mem::transmute(sym(h, c"sinpif")) };
+    let g_cospif: extern "C" fn(f32) -> f32 = unsafe { core::mem::transmute(sym(h, c"cospif")) };
+    let g_tanpif: extern "C" fn(f32) -> f32 = unsafe { core::mem::transmute(sym(h, c"tanpif")) };
+
+    let mut div: Vec<String> = Vec::new();
+
+    macro_rules! cmp64 {
+        ($name:literal, $flf:path, $gf:expr, $x:expr) => {{
+            let x: f64 = $x;
+            unsafe { feclearexcept(HARD) };
+            let fv = unsafe { $flf(x) };
+            let ff = unsafe { fetestexcept(HARD) } & HARD;
+            unsafe { feclearexcept(HARD) };
+            let gv = $gf(x);
+            let gf2 = unsafe { fetestexcept(HARD) } & HARD;
+            if !ulp_ok_f64(fv, gv) || ff != gf2 {
+                div.push(format!(
+                    "{}({:.6e}): fl={:016x}/flags{:#x} glibc={:016x}/flags{:#x}",
+                    $name, x, fv.to_bits(), ff, gv.to_bits(), gf2
+                ));
+            }
+        }};
+    }
+    macro_rules! cmp32 {
+        ($name:literal, $flf:path, $gf:expr, $x:expr) => {{
+            let x: f32 = $x;
+            unsafe { feclearexcept(HARD) };
+            let fv = unsafe { $flf(x) };
+            let ff = unsafe { fetestexcept(HARD) } & HARD;
+            unsafe { feclearexcept(HARD) };
+            let gv = $gf(x);
+            let gf2 = unsafe { fetestexcept(HARD) } & HARD;
+            if !ulp_ok_f32(fv, gv) || ff != gf2 {
+                div.push(format!(
+                    "{}({:.6e}): fl={:08x}/flags{:#x} glibc={:08x}/flags{:#x}",
+                    $name, x, fv.to_bits(), ff, gv.to_bits(), gf2
+                ));
+            }
+        }};
+    }
+
+    for &x in XS {
+        cmp64!("sinpi", fl::sinpi, g_sinpi, x);
+        cmp64!("cospi", fl::cospi, g_cospi, x);
+        cmp64!("tanpi", fl::tanpi, g_tanpi, x);
+    }
+    for &x in XSF {
+        cmp32!("sinpif", fl::sinpif, g_sinpif, x);
+        cmp32!("cospif", fl::cospif, g_cospif, x);
+        cmp32!("tanpif", fl::tanpif, g_tanpif, x);
+    }
+
+    assert!(
+        div.is_empty(),
+        "pi-trig divergences vs glibc ({}):\n  {}",
+        div.len(),
+        div.join("\n  ")
+    );
+}