symplex 0.2.0

//! 0.2 numeric back-ends: C99 code generation (`to_c_fn*`).
//!
//! Structural tests plus one end-to-end test that compiles the emitted C
//! with the system C compiler (skipped when none is available) and compares
//! results against `Ex::compile`.

use std::process::Command;

use symplex::matrix::{CodegenOptions, Precision};
use symplex::prelude::*;

fn assert_valid_c(code: &str, fn_name: &str) {
    assert!(
        code.starts_with("/* Generated by symplex. */\n#include <math.h>\n"),
        "{code}"
    );
    assert!(
        code.contains(&format!(" {fn_name}(")),
        "missing function `{fn_name}` in:\n{code}"
    );
    for (open, close) in [('{', '}'), ('(', ')')] {
        let o = code.chars().filter(|&c| c == open).count();
        let c = code.chars().filter(|&c| c == close).count();
        assert_eq!(o, c, "unbalanced {open}{close} in:\n{code}");
    }
    assert!(code.trim_end().ends_with('}'));
}

#[test]
fn basic_structure_and_cse() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let y = ctx.symbol("y");
    let s = x.sin();
    let code = (&s.powi(2) + &(&s * &y.cos()))
        .to_c_fn("f", &["x", "y"])
        .unwrap();
    assert_valid_c(&code, "f");
    assert!(code.contains("double f(double x, double y) {"), "{code}");
    assert!(code.contains("const double t0 = sin(x);"), "{code}");
    assert!(code.contains("return "), "{code}");
    // Zero-argument functions take `void`.
    let code = ctx.pi().to_c_fn("pi_val", &[]).unwrap();
    assert!(code.contains("double pi_val(void) {"), "{code}");
    assert!(code.contains("3.141592653589793"));
}

#[test]
fn math_h_functions_and_powers() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let y = ctx.symbol("y");
    let code = (x.gamma() + x.log_gamma() + x.erf() + y.erfc() + x.abs() + x.floor() + y.ceiling())
        .to_c_fn("g", &["x", "y"])
        .unwrap();
    for needle in [
        "tgamma(x)",
        "lgamma(x)",
        "erf(x)",
        "erfc(y)",
        "fabs(x)",
        "floor(x)",
        "ceil(y)",
    ] {
        assert!(code.contains(needle), "missing {needle}:\n{code}");
    }
    assert!(
        !code.contains("static inline"),
        "no helpers needed:\n{code}"
    );
    let code = (x.powi(3) + x.powi(-2) + x.powi(9) + x.sqrt() + x.pow(&ctx.rational(1, 3)))
        .to_c_fn("p", &["x"])
        .unwrap();
    assert!(code.contains("x * x * x"), "{code}");
    assert!(code.contains("1.0 / (x * x)"), "{code}");
    assert!(code.contains("pow(x, 9.0)"), "{code}");
    assert!(code.contains("sqrt(x)"), "{code}");
    assert!(code.contains("cbrt(x)"), "{code}");
    let code = (&x / &y).to_c_fn("d", &["x", "y"]).unwrap();
    assert!(code.contains("return x / y;"), "{code}");
    let code = (x.exp() - 1 + (&y + 1).ln())
        .to_c_fn("n", &["x", "y"])
        .unwrap();
    assert!(code.contains("expm1(x)"), "{code}");
    assert!(code.contains("log1p(y)"), "{code}");
    let code = y.atan2(&x).to_c_fn("a", &["x", "y"]).unwrap();
    assert!(code.contains("atan2(y, x)"), "{code}");
    let code = (x.min_with(&y) + x.max_with(&y))
        .to_c_fn("m", &["x", "y"])
        .unwrap();
    assert!(
        code.contains("fmin(x, y)") && code.contains("fmax(x, y)"),
        "{code}"
    );
}

#[test]
fn helpers_are_emitted_once_in_dependency_order() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let e =
        &(&x.harmonic() + &x.digamma()) + &(&x.bessel_j(&ctx.int(1)) + &x.bessel_y(&ctx.int(2)));
    let code = e.to_c_fn("h", &["x"]).unwrap();
    assert_valid_c(&code, "h");
    assert_eq!(
        code.matches("static inline double symplex_digamma(")
            .count(),
        1
    );
    assert_eq!(
        code.matches("static inline double symplex_harmonic(")
            .count(),
        1
    );
    assert_eq!(
        code.matches("static inline void symplex_bessel_miller(")
            .count(),
        1
    );
    let pos = |s: &str| {
        code.find(s)
            .unwrap_or_else(|| panic!("missing {s}:\n{code}"))
    };
    assert!(pos("symplex_is_int(double") < pos("symplex_digamma(double"));
    assert!(pos("symplex_digamma(double") < pos("symplex_harmonic(double"));
    assert!(pos("symplex_bessel_j_series(int") < pos("symplex_bessel_j(int"));
    assert!(pos("symplex_harmonic(double") < pos("double h(double x)"));
    assert!(code.contains("symplex_bessel_j(1, x)") && code.contains("symplex_bessel_y(2, x)"));
    // emit_runtime = false: references without definitions + full runtime API.
    let opts = CodegenOptions {
        emit_runtime: false,
        ..Default::default()
    };
    let code = x
        .lambertw()
        .to_c_fn_with_options("w", &["x"], &opts)
        .unwrap();
    assert!(code.contains("symplex_lambert_w0(x)") && !code.contains("static inline"));
    let rt = opts.c_runtime();
    assert!(rt.contains("static inline double symplex_lambert_w0(double x)"));
    assert!(rt.contains("static inline double symplex_bessel_k(int n, double x)"));
}

#[test]
fn options_precision_inline_fma_asserts() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let y = ctx.symbol("y");
    let e = &(&x * &y) + &x.ln();
    let code = e.to_c_fn("f", &["x", "y"]).unwrap();
    assert!(code.contains("fma(x, y, log(x))"), "{code}");
    let opts = CodegenOptions {
        precision: Precision::F32,
        inline: true,
        use_mul_add: false,
        checked_domain: true,
        ..Default::default()
    };
    let code = e.to_c_fn_with_options("f", &["x", "y"], &opts).unwrap();
    assert!(code.contains("#include <assert.h>"));
    assert!(
        code.contains("static inline float f(float x, float y) {"),
        "{code}"
    );
    assert!(!code.contains("fma"), "{code}");
    assert!(code.contains("(assert(x > 0.0f), logf(x))"), "{code}");
    // Helpers stay double; float call sites cast.
    let code = x
        .lambertw()
        .to_c_fn_with_options("w", &["x"], &opts)
        .unwrap();
    assert!(
        code.contains("(float)symplex_lambert_w0((double)x)"),
        "{code}"
    );
    assert!(code.contains("static inline double symplex_lambert_w0(double x)"));
    let code = x
        .bessel_k(&ctx.int(2))
        .to_c_fn_with_options("k", &["x"], &opts)
        .unwrap();
    assert!(
        code.contains("(float)symplex_bessel_k(2, (double)x)"),
        "{code}"
    );
    assert!(code.contains("assert(x > 0.0f)"), "K domain check:\n{code}");
}

#[test]
fn piecewise_and_booleans() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let y = ctx.symbol("y");
    let zero = ctx.zero();
    let pw = Ex::piecewise(&[
        (&x, &x.gt(&y).and(&x.gt(&zero))),
        (&(-&x), &x.le(&y).not()),
        (&zero, &x.le(&y).or(&x.gt(&y))),
    ]);
    let code = pw.to_c_fn("pw", &["x", "y"]).unwrap();
    assert_valid_c(&code, "pw");
    assert!(code.contains("((x > y) && (x > 0.0)) ? x :"), "{code}");
    assert!(code.contains("(!(y >= x)) ? (-x) :"), "{code}");
    assert!(code.contains("((y >= x) || (x > y)) ? 0.0 : NAN"), "{code}");
    // Boolean node in numeric position → 1.0 / 0.0.
    let code = Ex::piecewise(&[(&x, &x.gt(&zero))])
        .to_c_fn("half", &["x"])
        .unwrap();
    assert!(code.contains("return (x > 0.0) ? x : NAN;"), "{code}");
}

#[test]
fn errors_name_the_offender() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let y = ctx.symbol("y");
    match (&x + &y).to_c_fn("f", &["x"]) {
        Err(SymplexError::FreeSymbol { name }) => assert_eq!(name, "y"),
        other => panic!("{other:?}"),
    }
    match ctx.i_unit().to_c_fn("f", &[]) {
        Err(SymplexError::NotImplemented(msg)) => assert!(msg.contains("ImaginaryUnit")),
        other => panic!("{other:?}"),
    }
    let n = ctx.symbol("n");
    match x.bessel_i(&n).to_c_fn("f", &["x", "n"]) {
        Err(SymplexError::NotImplemented(msg)) => assert!(msg.contains("besseli")),
        other => panic!("{other:?}"),
    }
    match (x.sin().sin()).exp().integrate(&x).to_c_fn("f", &["x"]) {
        Err(SymplexError::NotImplemented(msg)) => assert!(msg.contains("Integral")),
        other => panic!("{other:?}"),
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// End-to-end with the system C compiler
// ═══════════════════════════════════════════════════════════════════════════

fn find_cc() -> Option<String> {
    for cc in ["cc", "clang", "gcc"] {
        if Command::new(cc)
            .arg("--version")
            .output()
            .is_ok_and(|o| o.status.success())
        {
            return Some(cc.to_string());
        }
    }
    None
}

#[test]
fn generated_c_compiles_and_matches_compile() {
    let Some(cc) = find_cc() else {
        eprintln!("no C compiler found; skipping end-to-end C test");
        return;
    };
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let y = ctx.symbol("y");
    let pw = Ex::piecewise(&[(&x, &x.gt(&y)), (&y, &x.le(&y))]);
    let cases: Vec<(&str, Ex)> = vec![
        ("poly", &x.powi(2) + &(&x * 3) + 1),
        (
            "trig_cse",
            &x.sin().powi(2) + &(&x.cos().powi(2) * &y) + &(&x.sin() * &x.cos()),
        ),
        (
            "gamma_family",
            &(&x.gamma() + &x.log_gamma()) + &x.digamma(),
        ),
        ("erf_family", &(&x.erf() * &y.erfc()) + &(&x + 1).lambertw()),
        (
            "beta_binomial",
            &x.beta(&y) + &(&y.binomial(&ctx.int(2)) / &y.factorial()),
        ),
        (
            "bessel",
            &(&x.bessel_j(&ctx.int(2)) + &x.bessel_y(&ctx.int(1)))
                + &(&x.bessel_i(&ctx.int(0)) * &x.bessel_k(&ctx.int(3))),
        ),
        (
            "orthopoly",
            &(&(&x.legendre(&ctx.int(4)) + &x.chebyshev_t(&ctx.int(3)))
                + &(&x.chebyshev_u(&ctx.int(2)) + &x.hermite(&ctx.int(3))))
                + &x.laguerre(&ctx.int(2)),
        ),
        (
            "sequences",
            &(&(&y.fibonacci() + &y.lucas()) + &y.harmonic()) + &y.factorial2(),
        ),
        (
            "pochhammer",
            &x.rising_factorial(&ctx.int(3)) + &x.falling_factorial(&ctx.int(2)),
        ),
        (
            "piecewise_elem",
            &(&(&pw + &x.min_with(&y)) + &(&x.max_with(&y) + &x.sign()))
                + &(&(&(&x - &y).heaviside() + &x.floor()) + &(&y.ceiling() + &y.atan2(&x))),
        ),
        ("numopt", &(&x.exp() - 1) + &(&y + 1).ln()),
        (
            "roots_div",
            &(&x.pow(&ctx.rational(1, 3)) + &x.sqrt()) + &(&x / &(&y + 1)),
        ),
        ("neg_gamma", (-&x).gamma() * (&ctx.zero() - &y).erf()),
        (
            "real_roots",
            &(-&x).pow(&ctx.rational(2, 5)) + &(-&x).pow(&ctx.rational(3, 5)),
        ),
    ];
    let points: &[(f64, f64)] = &[
        (0.5, 2.0),
        (1.7, 3.0),
        (3.25, 1.0),
        (7.5, 5.0),
        (0.125, 4.0),
    ];

    let opts = CodegenOptions {
        emit_runtime: false,
        ..Default::default()
    };
    let mut src = String::new();
    src.push_str("#include <stdio.h>\n");
    src.push_str(&opts.c_runtime());
    for (name, expr) in &cases {
        let code = expr
            .to_c_fn_with_options(name, &["x", "y"], &opts)
            .unwrap_or_else(|e| panic!("{name}: {e}"));
        // Drop the per-function preamble (runtime already included once).
        let body = code
            .trim_start_matches("/* Generated by symplex. */\n#include <math.h>\n")
            .to_string();
        src.push_str(&body);
        src.push('\n');
    }
    src.push_str("int main(void) {\n");
    for (name, _) in &cases {
        for (px, py) in points {
            src.push_str(&format!(
                "    printf(\"%.17g\\n\", {name}({px:?}, {py:?}));\n"
            ));
        }
    }
    src.push_str("    return 0;\n}\n");

    let dir = std::env::temp_dir().join(format!(
        "symplex_c_e2e_{}_{}",
        std::process::id(),
        std::time::SystemTime::now()
            .duration_since(std::time::UNIX_EPOCH)
            .map(|d| d.as_nanos())
            .unwrap_or(0)
    ));
    std::fs::create_dir_all(&dir).unwrap();
    let c_file = dir.join("gen.c");
    std::fs::write(&c_file, &src).unwrap();
    let exe = dir.join("gen_bin");
    let out = Command::new(&cc)
        .args(["-std=c99", "-O2", "-Wall", "-Werror", "-o"])
        .arg(&exe)
        .arg(&c_file)
        .arg("-lm")
        .output()
        .expect("run C compiler");
    assert!(
        out.status.success(),
        "generated C failed to compile:\n{}\n--- source ---\n{src}",
        String::from_utf8_lossy(&out.stderr)
    );
    let run = Command::new(&exe).output().unwrap();
    assert!(run.status.success());
    let stdout = String::from_utf8_lossy(&run.stdout);
    let values: Vec<f64> = stdout
        .lines()
        .map(|l| {
            let t = l.trim();
            match t {
                "nan" | "-nan" | "nan(ind)" => f64::NAN,
                "inf" => f64::INFINITY,
                "-inf" => f64::NEG_INFINITY,
                _ => t
                    .parse::<f64>()
                    .unwrap_or_else(|_| panic!("bad C output {t:?}")),
            }
        })
        .collect();
    let _ = std::fs::remove_dir_all(&dir);
    assert_eq!(values.len(), cases.len() * points.len());

    let mut idx = 0;
    for (name, expr) in &cases {
        let f = expr.compile(&["x", "y"]).unwrap();
        for &(px, py) in points {
            let want = f(&[px, py]);
            let got = values[idx];
            idx += 1;
            if want.is_nan() {
                assert!(got.is_nan(), "{name}({px},{py}): C={got}, vm=NaN");
                continue;
            }
            // libm's tgamma/lgamma/erf differ from the shared runtime by a few ulps.
            let err = (got - want).abs() / want.abs().max(1.0);
            assert!(
                err < 1e-12,
                "{name}({px},{py}): C={got:?}, vm={want:?}, err {err:e}"
            );
        }
    }
}