ported/excel_numeric/x87

/// Public ISA FPATAN operation, selected by W111 black-box complex-argument
/// witnesses. No binary provenance inference is needed for this arithmetic graph.
pub(super) fn atan(x: f64) -> f64 {
    let mut result = 0.0_f64;
    let cw_core = CW_CORE;
    let mut cw_save = 0_u16;
    // SAFETY: the two pushes and FPATAN/store pops balance the x87 stack. All
    // pointers refer to live locals and the caller's control word is restored.
    unsafe {
        asm!(
            "fnstcw word ptr [{save}]",
            "fldcw word ptr [{core}]",
            "fld qword ptr [{x}]",
            "fld1",
            "fpatan",
            "fstp qword ptr [{res}]",
            "fldcw word ptr [{save}]",
            x = in(reg) &x,
            res = in(reg) &mut result,
            core = in(reg) &cw_core,
            save = in(reg) &mut cw_save,
            options(nostack, preserves_flags),
        );
    }
    result
}

/// W111 ATAN magnitude reduction selected by independent black-box counterexamples.
/// For |x| > 1, keep atan(1/|x|) and pi/2 extended until the final subtraction.
/// FPATAN takes the two operands directly, so no rounded reciprocal is introduced.
pub(super) fn atan_reduced(x: f64) -> f64 {
    let magnitude = x.abs();
    if magnitude <= 1.0 { return atan(x); }
    let mut result = 0.0_f64;
    let mut cw_save = 0_u16;
    let cw_core = CW_CORE;
    let two = 2.0_f64;
    // SAFETY: all operands are live locals; the two FPATAN operands reduce to
    // one angle, FLDPI adds one value, and subtraction/store empty the stack.
    // The caller's control word is restored before returning.
    unsafe {
        asm!(
            "fnstcw word ptr [{save}]", "fldcw word ptr [{core}]",
            "fld1", "fld qword ptr [{x}]", "fpatan",
            "fldpi", "fdiv qword ptr [{two}]", "fsubrp st(1), st(0)",
            "fstp qword ptr [{res}]", "fldcw word ptr [{save}]",
            x = in(reg) &magnitude, two = in(reg) &two, res = in(reg) &mut result,
            core = in(reg) &cw_core, save = in(reg) &mut cw_save,
            options(nostack, preserves_flags),
        );
    }
    result.copysign(x)
}

/// `RN53(RN64(a * b))` — one x87 `FMUL` at PC=64 then a `binary64` store: the
/// **double-rounded** product. Excel forms `POWER`'s `|y|·ln x` this way; a plain
/// SSE multiply is single-rounded and differs on double-rounding-window inputs.
pub(super) fn mul(a: f64, b: f64) -> f64 {
    let mut result: f64 = 0.0;
    let cw_core: u16 = CW_CORE;
    let mut cw_save: u16 = 0;
    // SAFETY: balanced x87 stack (2 push; `fmulp` and `fstp` pop both), pointer
    // I/O only, control word saved and restored.
    unsafe {
        asm!(
            "fnstcw word ptr [{save}]",
            "fldcw word ptr [{core}]",
            "fld qword ptr [{a}]",     // st0 = a
            "fld qword ptr [{b}]",     // st0 = b, st1 = a
            "fmulp st(1), st",         // st0 = RN64(a*b)
            "fstp qword ptr [{res}]",  // RN53 store
            "fldcw word ptr [{save}]",
            a = in(reg) &a,
            b = in(reg) &b,
            res = in(reg) &mut result,
            core = in(reg) &cw_core,
            save = in(reg) &mut cw_save,
        );
    }
    result
}

/// `RN53(RN64(a + b))` — one x87 `FADD` at PC=64 then a `binary64` store: the
/// double-rounded sum. Excel's legacy x87-compiled function bodies produce
/// this for every `t = a + b;` with `t` spilled to memory (W109: XNPV's
/// `1 + rate` and its running `total += term` accumulate this way).
pub(super) fn add(a: f64, b: f64) -> f64 {
    let mut result: f64 = 0.0;
    let cw_core: u16 = CW_CORE;
    let mut cw_save: u16 = 0;
    // SAFETY: balanced x87 stack (one push; `fadd` folds in memory operand,
    // `fstp` pops), pointer I/O only, control word saved and restored.
    unsafe {
        asm!(
            "fnstcw word ptr [{save}]",
            "fldcw word ptr [{core}]",
            "fld qword ptr [{a}]",     // st0 = a
            "fadd qword ptr [{b}]",    // st0 = RN64(a+b)
            "fstp qword ptr [{res}]",  // RN53 store
            "fldcw word ptr [{save}]",
            a = in(reg) &a,
            b = in(reg) &b,
            res = in(reg) &mut result,
            core = in(reg) &cw_core,
            save = in(reg) &mut cw_save,
        );
    }
    result
}

/// `RN53(RN64(a - b))` — one x87 `FSUB` at PC=64 then a `binary64` store: the
/// double-rounded difference (W109: NPER's `tf·pmt - fv·rate`).
pub(super) fn sub(a: f64, b: f64) -> f64 {
    let mut result: f64 = 0.0;
    let cw_core: u16 = CW_CORE;
    let mut cw_save: u16 = 0;
    // SAFETY: balanced x87 stack (one push; `fsub` folds in memory operand,
    // `fstp` pops), pointer I/O only, control word saved and restored.
    unsafe {
        asm!(
            "fnstcw word ptr [{save}]",
            "fldcw word ptr [{core}]",
            "fld qword ptr [{a}]",     // st0 = a
            "fsub qword ptr [{b}]",    // st0 = RN64(a-b)
            "fstp qword ptr [{res}]",  // RN53 store
            "fldcw word ptr [{save}]",
            a = in(reg) &a,
            b = in(reg) &b,
            res = in(reg) &mut result,
            core = in(reg) &cw_core,
            save = in(reg) &mut cw_save,
        );
    }
    result
}

/// `RN53(RN64(a / b))` — one x87 `FDIV` at PC=64 then a `binary64` store: the
/// double-rounded quotient (W109: XNPV's per-term `value / pow`).
pub(super) fn div(a: f64, b: f64) -> f64 {
    let mut result: f64 = 0.0;
    let cw_core: u16 = CW_CORE;
    let mut cw_save: u16 = 0;
    // SAFETY: balanced x87 stack (one push; `fdiv` folds in memory operand,
    // `fstp` pops), pointer I/O only, control word saved and restored.
    unsafe {
        asm!(
            "fnstcw word ptr [{save}]",
            "fldcw word ptr [{core}]",
            "fld qword ptr [{a}]",     // st0 = a
            "fdiv qword ptr [{b}]",    // st0 = RN64(a/b)
            "fstp qword ptr [{res}]",  // RN53 store
            "fldcw word ptr [{save}]",
            a = in(reg) &a,
            b = in(reg) &b,
            res = in(reg) &mut result,
            core = in(reg) &cw_core,
            save = in(reg) &mut cw_save,
        );
    }
    result
}

/// `RN53(RN64(1.0 / x))` — one x87 `FDIV` at PC=64 then a `binary64` store: the
/// double-rounded reciprocal Excel uses to stage `POWER(x, y)` for `y < 0`
/// (compute the positive power, store it, then reciprocate).
pub(super) fn recip(x: f64) -> f64 {
    let mut result: f64 = 0.0;
    let cw_core: u16 = CW_CORE;
    let mut cw_save: u16 = 0;
    // SAFETY: balanced x87 stack (`fld1` pushes, `fstp` pops), pointer I/O only.
    unsafe {
        asm!(
            "fnstcw word ptr [{save}]",
            "fldcw word ptr [{core}]",
            "fld1",                    // st0 = 1
            "fdiv qword ptr [{x}]",    // st0 = RN64(1/x)
            "fstp qword ptr [{res}]",  // RN53 store
            "fldcw word ptr [{save}]",
            x = in(reg) &x,
            res = in(reg) &mut result,
            core = in(reg) &cw_core,
            save = in(reg) &mut cw_save,
        );
    }
    result
}