A symbolic determinant is a quotient by its pivots, so it is NaN for ordinary matrices #992

Measured on master at 4ee698da.

Matrix.Determinant is computed by Gaussian elimination and the pivots are left as literal divisions, so the expression it returns is undefined wherever a pivot vanishes — at points where the determinant itself is perfectly well defined.

2×2

((Entity.Matrix)"[[x, 1], [2, y]]".ToEntity()).Determinant
// x * (y * x + -2) / x        -- and .Simplify() is  x * y - 2 provided not x = 0

The determinant of that matrix is x * y - 2 for every x, 0 included. What comes back refuses at x = 0:

det.Substitute("x", 0).Substitute("y", 5).Evaled   // NaN      the determinant is -2
det.Substitute("x", 3).Substitute("y", 4).Evaled   // 10       correct

The same matrix built with the zero already in it is fine, so this is the symbolic path only:

((Entity.Matrix)"[[0, 1], [2, 5]]".ToEntity()).Determinant.Evaled   // -2

3×3, where it stops being an edge case

((Entity.Matrix)"[[a, b, c], [d, e, f], [g, h, i]]".ToEntity()).Determinant.Simplify()
// a * (a ^ 2 * (i * a + -c * g) * (a * e + -b * d) + -a ^ 2 * (a * f + -c * d) * (a * h + -b * g))
//   * (a * e + -b * d) / (a ^ 2 * a * a * (a * e + -b * d))

The denominator is a ^ 4 * (a * e - b * d), so the expression is undefined whenever a = 0 or the top-left 2×2 minor vanishes. Substituting ordinary integer matrices into it:

Two of four ordinary matrices. The first has a = 0; the second has a * e = b * d. Neither is degenerate as a matrix — the first is the classic singular example every linear-algebra course opens with.

Why this is a defect and not a missing feature

NaN here means "this does not exist", and the determinant does exist — AGENTS.md is explicit that confusing that with "I could not settle this" ships a wrong answer. A caller who computes a symbolic determinant and substitutes gets NaN for a matrix whose determinant is 0, which is exactly the case they were most likely checking for.

It is also not a hard fix in principle: the determinant of a matrix over a commutative ring is a polynomial in its entries, so there is a division-free algorithm (cofactor expansion for small n, Bareiss for larger) whose result needs no condition at all. The current output is not merely conditioned — it is a quotient, so even the unconditional part has to be cancelled before it reads as a polynomial.

Found while surveying SymPy parity for #717; SymPy's Matrix([[x,1],[2,y]]).det() is x*y - 2 with no condition.