Look inside divisions in query_maybe_zeros

[mwhincup]

Fixes #166. Also fixes the case reported in SciML/ModelingToolkit.jl#5279 (see below).

Problem

_check_allow_symbolic_parameter uses query_maybe_zeros to refuse division by an expression that may be zero. maybe_zeros_descend did not descend into BSImpl.Div, and MaybeZeroQueryFn has no Div case. So any quotient was treated as safe, including one whose numerator is built from maybe_zeros entries, such as (Ixy*Iyz + Ixz*Iyy) / det(J). A product holding such a quotient was also missed, for example r * ((p + q) / s).

Change

maybe_zeros_descend now returns true for BSImpl.Div, so SU.query visits the numerator and the denominator with the existing predicate. A quotient is flagged when:

• its numerator may be zero (by the existing rules for factors, sums whose terms are all maybe-zero, powers, max/min/ifelse, ...), or
• its denominator may be zero. The quotient is then not finite, and dividing by it is just as unsafe.

There are no other changes: the predicate and the other descent rules are as before.

Effect on tearing

This only changes anything when an expression in maybe_zeros(sys) (declared, or found by discover_maybe_zeros) appears inside a division in a coefficient or a linear-SCC matrix entry. Those coefficients are no longer divided by. Instead the variable is not solved for through that equation, and an inline linear SCC uses the numeric safe_ldiv path instead of sym_lu. All existing tests in both packages give the same results.

Tests

New testsets in lib/ModelingToolkitTearing/test/runtests.jl:

1. `maybe_zeros` is respected inside divisions
   • unit checks of query_maybe_zeros on quotients: numerator maybe-zero, denominator maybe-zero, nested in products and powers, and negative cases;
   • a structural check that (p*r + q)/s * y ~ x is not solved for y when maybe_zeros = [p, q];
   • the maybe_zeros check does not look inside a division, so tearing divides by a parameter expression that is zero (NaN in initialization) #166 case: trimming a rigid body by initialisation (D(w) ~ 0 with the moments as missing parameters, explicit adj(J)/det(J), only Ixz non-zero). The initialisation LinearProblem used to hold [-Inf NaN; -Inf NaN] and is now finite.
2. Inline linear SCC does not pivot on a quotient with a maybe-zero numerator: the inline_linear_sccs: analytical sym_lu divides by (I21 - I11*I32/I31), which is 0 when only I31 is non-zero (NaN in rigid-body Euler equations) SciML/ModelingToolkit.jl#5279 Euler-equation reproducer with inline_linear_sccs = true and only I31 non-zero. On main the RHS is NaN. Now it matches J \ (-w × Jw).

Before the fix, 12 of the 18 checks in (1) and both checks in (2) fail. After it, all pass. The full suites pass for StateSelection and ModelingToolkitTearing (Julia 1.12.5, ModelingToolkit 11.45.4, SymbolicUtils 4.51.0, Symbolics 7.45.0).

Not covered

A compiled system is still only safe for the zeros known at compile time. If no products of inertia are declared or found as maybe-zero (for example, all defaults non-zero) and they are set to zero later in ODEProblem, both tearing and the analytic linear-SCC solve can still divide by zero. That is the existing maybe_zeros contract, so I've left it out of this PR. It is discussed in SciML/ModelingToolkit.jl#5279.

🤖 Generated with Claude Code