Add exact multiaffine real-stability tests via Rayleigh differences

[PerAlexandersson]

Summary

Add an exact real-stability test for multiaffine polynomials with integer
coefficients, using Brändén's Rayleigh-difference characterization. Reuse the
existing randomized stability/counterexample screening where useful, but do
not let randomized sampling certify a positive result.

This should be presented beside the existing Lorentzian functionality while
keeping the mathematical distinction explicit:

homogeneous + multiaffine + nonnegative + real stable  =>  Lorentzian,

but Lorentzian does not imply real stable.

Exact criterion

For a real multiaffine polynomial f, define

Delta_ij(f) = (partial_i f)(partial_j f) - f(partial_i partial_j f).

Then f is real stable if and only if

Delta_ij(f)(x) >= 0

for every real point x and every distinct pair i,j. Because f is
multiaffine, Delta_ij is independent of x_i and x_j, so each global
nonnegativity problem has only n-2 variables.

For comparison and independent testing, failure of real stability can also be
written directly as the existential real-algebraic formula

exists x,y: y_1>0, ..., y_n>0,
            Re(f(x+i y))=0 and Im(f(x+i y))=0.

Thus a complete fallback exists by quantifier elimination over real closed
fields. Integer coefficients make the entire question exact.

Repository architecture and reuse

The public polytool crate is currently univariate and does not contain a
multivariate polynomial representation. The monorepo already has the relevant
infrastructure in:

• sym-poly/multipoly/src/multipoly.rs: sparse MultiPoly<C> and exact
  arithmetic;
• sym-poly/multipoly/src/lorentzian.rs: Lorentzian result diagnostics,
  support tests, exact Hessian/inertia code, and a currently private
  partial_derivative implementation.

Avoid creating a second multivariate representation inside polytool.
Prefer:

1. a canonical library implementation in sym-poly-multipoly;
2. a thin user-facing command/MCP surface associated with Polytool;
3. an explicit packaging decision for the standalone polytool projection
   (published dependency, feature-gated companion crate, or documented
   monorepo-only command).

Promote or generalize shared derivative helpers instead of copying the private
Lorentzian implementation. Use BigInt as the canonical coefficient type:
forming derivatives and products in Delta_ij can overflow i64 even when
the input coefficients fit in i64.

Proposed API

Names are provisional, but the library should expose concepts at roughly this
level:

pub fn rayleigh_difference(
    f: &MultiPoly<BigInt>,
    i: usize,
    j: usize,
) -> Result<MultiPoly<BigInt>, RealStabilityInputError>;

pub fn check_multiaffine_real_stability(
    f: &MultiPoly<BigInt>,
    options: &RealStabilityOptions,
) -> RealStabilityResult;

The structured result should distinguish at least:

• Stable, backed only by a complete exact method or an exactly verified
  nonnegativity certificate;
• Unstable, including the failing pair (i,j), Delta_ij, and an exact
  witness whenever available;
• NotMultiaffine, including a monomial and offending exponent;
• ResourceLimit / UnsupportedBackend, which must not be reported as
  unstable.

The zero-polynomial convention should be chosen explicitly and documented,
since the current Lorentzian checker accepts zero while definitions of stable
polynomial often exclude it.

Implementation plan

1. Exact Rayleigh layer

• Validate exponent-vector lengths, variable indices, and multiaffinity.
• Implement exact first and mixed partial derivatives for MultiPoly<BigInt>.
• Form every Delta_ij with checked size/term budgets.
• Canonically remove the known-independent variables x_i,x_j and assert that
  no residual exponent remains in those positions.
• Deduplicate symmetric pairs and short-circuit on the first certified failure,
  while optionally returning all pair diagnostics.

2. Tiered global-nonnegativity backend

Use a common trait/result type so faster special cases and a complete general
backend have identical certificate semantics.

• Zero remaining variables (n <= 2): exact constant sign.
• One remaining variable (n = 3): exact global nonnegativity using the
  existing Sturm/root-isolation infrastructure or an extracted shared exact
  univariate kernel. A Rayleigh difference here has degree at most two, but the
  implementation need not rely on that accident.
• Two remaining variables (n = 4): add a dedicated exact path only if it
  is demonstrably simpler and certifying; otherwise route to the general
  backend. Rayleigh differences are at most quadratic in each remaining
  variable.
• General case: decide whether Delta_ij(x) < 0 has a real solution using
  an exact real-algebraic/quantifier-elimination backend. First evaluate
  available native Rust libraries and a clean backend abstraction; do not make
  an optional Mathematica/Sage installation part of the core correctness
  contract.
• An exactly checked rational SOS certificate may certify nonnegativity, but
  failure to find an SOS decomposition is only inconclusive, never
  unstable.

The existing randomized path should remain a cheap prefilter. A sampled
negative value can be converted to and rechecked as an exact rational witness;
sampling that finds no violation cannot return Stable.

3. CLI and MCP surface

Add a discoverable command/tool such as:

polytool multivariate real-stability --input polynomial.json --format json

Use a sparse exact JSON representation with variable names, exponent vectors,
and coefficients accepted as decimal strings, so JavaScript/MCP clients cannot
silently round large integers. Include resource controls for variables, terms,
coefficient bits, Rayleigh-product terms, backend work, and wall time.

Machine-readable output should include:

• final status and whether it is a proof, counterexample, or inconclusive run;
• multiaffinity validation;
• pairs checked and the first failing pair;
• exact witness coordinates and exact value of Delta_ij when unstable;
• backend/certificate provenance and resource statistics;
• whether randomized screening was used (never as positive evidence).

Expose the same schema and status vocabulary through MCP. Do not return a bare
boolean when a resource limit or unavailable exact backend prevents a decision.

4. Tests and independent validation

Include at least:

• stable products of positive linear forms and elementary symmetric
  polynomials;
• e_2(x,y,z) = xy+xz+yz as a small stable case;
• 1+xyz, for which Delta_xy = -z, as a small unstable case with witness;
• rejection of x^2+y as nonmultiaffine;
• zero-polynomial convention and unused-variable behavior;
• coefficients large enough to force BigInt intermediates;
• invariance under variable permutation and multiplication by a positive
  scalar;
• examples that are Lorentzian but not real stable, preventing accidental
  conflation of the two predicates;
• agreement between the Rayleigh implementation and the direct
  real/imaginary quantifier-elimination formulation on a bounded fixture set;
• CLI/MCP schema, exact-integer serialization, resource-limit, and diagnostic
  tests.

Use published examples where possible and cross-check small fixtures against
an independent CAS/QE implementation. External CAS output is test evidence,
not the portable runtime implementation.

Acceptance criteria

• The canonical library computes every Delta_ij exactly with BigInt and
  returns useful typed diagnostics.
• All multiaffine inputs through three variables receive a complete exact
  answer without randomized evidence.
• Higher-dimensional positive answers are returned only when backed by an
  exactly verified certificate or a complete exact backend.
• A negative result includes a checkable failing pair and exact witness.
• Lorentzian and real-stability APIs/documentation clearly state implication
  direction and non-equivalence.
• CLI and MCP accept exact sparse input, expose resource controls, and preserve
  inconclusive rather than coercing it to false.
• Existing randomized checks are reused as optional counterexample screening,
  not reimplemented.
• Library, CLI, MCP, documentation, and independent-fixture tests pass in the
  monorepo and in the chosen standalone packaging configuration.

References

• Petter Brändén, Polynomials with the half-plane property and matroid
  theory, Advances in Mathematics 216 (2007), Theorem 5.6:
  https://arxiv.org/abs/math/0605678
• Petter Brändén and June Huh, Lorentzian polynomials, Annals of Mathematics
  192 (2020): https://doi.org/10.4007/annals.2020.192.3.4
• Saugata Basu, Richard Pollack, and Marie-Françoise Roy, Algorithms in Real
  Algebraic Geometry (quantifier elimination and semialgebraic feasibility):
  https://doi.org/10.1007/978-3-662-05355-3