Compute the determinant without dividing by its pivots (#992)

BREAKING-CHANGES.md

@@ -35,6 +35,8 @@ read first.
 | **Silent** | `limit(i, i, 0)` | unevaluated | `0` |
 | **Silent** | `integral(i, i)` | `-1/2 + C` | `i_1 ^ 2 / 2 + C` |
 | | `lambda(i, i + 1)` | `InvalidArgumentParseException` | the lambda |
+| **Silent** | the symbolic determinant of a matrix, substituted where a pivot vanishes — `[[0,1,2],[3,4,5],[6,7,8]]` through `[[a,b,c],[d,e,f],[g,h,i]]` | `NaN` | `0` |
+| **Silent** | `((Entity.Matrix)"[[x, 1], [2, y]]".ToEntity()).Determinant.Simplify()` | `x * y - 2 provided not x = 0` | `x * y - 2` |
 | **Silent** | `[1, 2] in RR`, and a matrix or a finite set against any of `BB`, `ZZ`, `QQ`, `RR`, `CC` | `True` | `False` |
 | **Silent** | `"sin(pi)".ToEntity().Differentiate(MathS.pi)`, and every derivative over `pi` or `e` | `-1`, the chain rule run over a symbol that cannot change | `0` |
 | **Silent** | `"x ^ 3".ToEntity().Differentiate(x, 2)`, and every `Differentiate(x, n)` with `n >= 1` | `(0 * x ^ 2 + 2 * x ^ 1 * 1 * 3) * 1 + 0 * 3 * x ^ 2` | `2 * x * 3` |
@@ -55,6 +57,64 @@ read first.
 | **Silent** | `integral(x, [x, y]T)` | `[[C + x ^ 2, C + x * y]]` | `[[x ^ 2 / 2 + C, x * y + C]]` |
 | **Silent** | `derivative(e ^ 2, e)`, over a named constant | `0` | `2 * e` |
 
+### The symbolic determinant is a polynomial, not a quotient by its pivots
+
+`Matrix.Determinant` was computed by Gaussian elimination, which leaves the pivots as literal
+divisions. The expression it returned was therefore undefined wherever a pivot vanishes — at points
+where the determinant itself is perfectly well defined.
+
+```csharp
+var m = (Entity.Matrix)"[[x, 1], [2, y]]".ToEntity();
+m.Determinant             // was: x * (y * x + -2) / x           now: x * y + -2
+m.Determinant.Simplify()  // was: x * y - 2 provided not x = 0   now: x * y - 2
+m.Determinant.Substitute("x", 0).Substitute("y", 5).Evaled   // was: NaN   now: -2
+```
+
+`x = 0` is an ordinary point of that matrix — its determinant there is `-2` — and `NaN` asserts the
+value **does not exist**. At 3×3 it stops being an edge case, because the denominator is
+`a ^ 4 * (a * e - b * d)` and so two conditions have to miss:
+
+```csharp
+// [[a, b, c], [d, e, f], [g, h, i]].Determinant, substituted:
+// [[0, 1, 2], [3, 4, 5], [6, 7, 8]]     was: NaN   now: 0     (the pivot a is 0)
+// [[1, 2, 3], [2, 4, 6], [1, 1, 1]]     was: NaN   now: 0     (a * e = b * d)
+// [[1, 2, 3], [4, 5, 6], [7, 8, 10]]    was: -3    now: -3
+// [[2, 1, 0], [1, 2, 1], [0, 1, 2]]     was: 4     now: 4
+```
+
+Two of those four are ordinary matrices, and the first is the singular example every linear-algebra
+course opens with. Both wrong answers were **silent**: the call succeeded and returned `NaN`, which
+is exactly the answer a caller checking for singularity was looking for.
+
+The determinant of a matrix over a commutative ring is a polynomial in its entries, so Laplace
+expansion — which never divides — needs no condition at all. It is what the property's own
+documentation and the comment above it already claimed was in use.
+
+**This is not a performance trade.** Measured on the same machine, property only, both arms built
+from source: Laplace returns a *smaller* expression at every size, and the elimination it replaces
+was the slower of the two on numeric matrices by a wide margin.
+
+| entries | Gaussian complexity | Laplace complexity | Gaussian, numeric | Laplace, numeric |
+|---|---|---|---|---|
+| 2×2 | 13 | 9 | | |
+| 3×3 | 79 | 37 | | |
+| 4×4 | 443 | 163 | | |
+| 5×5 | 2461 | 833 | | |
+| 6×6 | 13673 | 5021 | 164 ms | 126 ms |
+| 7×7 | | 35173 | | |
+| 8×8 | | | over 120 s | 358 ms |
+| 10×10 | | | over 120 s | 4478 ms |
+
+Laplace expansion is `O(n!)`, and a fully symbolic determinant has `n!` terms however it is
+computed, so that is the size of the answer rather than an overhead. The practical ceiling on a
+numeric matrix moves from 7×7 to 10×10; 11×11 does not return, where under the elimination 8×8
+already did not. A fraction-free elimination (Bareiss) would raise it further and is worth having,
+but it is a separate change and this one is not blocked on it.
+
+Measured: the whole suite passes with the fix in, and the corpus is unchanged — 116/119 with 0
+wrong, 0 error, 0 timeout, and no case's verdict or answer altered.
+[#992](https://github.com/asc-community/AngouriMath/issues/992).
+
 ### A matrix is no longer a member of every special set at once
 
 `SpecialSet.TryContains` decided membership by asking `MayContain`, which is deliberately permissive

Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs

@@ -264,7 +264,14 @@ public Entity AsScalar()
             private LazyPropertyA<Matrix> t;
 
             // We do not need to use Gaussian elimination here
-            // since we anyway get N! memory use
+            // since we anyway get N! memory use.
+            // Gaussian elimination is also not merely no cheaper, it is wrong here:
+            // it leaves the pivots as literal divisions, so the expression it returns
+            // is undefined wherever a pivot vanishes -- at points where the
+            // determinant itself is perfectly well defined. The determinant of a
+            // matrix over a commutative ring is a polynomial in its entries, and
+            // Laplace expansion never divides, so there is nothing to exclude.
+            // https://github.com/asc-community/AngouriMath/issues/992
             /// <summary>
             /// Finds the symbolical determinant via Laplace's method
             /// </summary>
@@ -273,7 +280,7 @@ public Entity AsScalar()
                 {
                     if (!@this.IsSquare) 
                         return null;
-                    return @this.InnerMatrix.DeterminantGaussianSafeDivision().InnerSimplified;
+                    return @this.InnerMatrix.DeterminantLaplace().InnerSimplified;
                 },
                 this
                 );

Sources/Tests/UnitTests/Algebra/MatrixTest.cs

@@ -5,6 +5,7 @@
 // Website: https://am.angouri.org.
 //
 
+using System.Linq;
 using AngouriMath;
 using AngouriMath.Core.Exceptions;
 using Xunit;
@@ -871,5 +872,56 @@ [Fact] public void Concat13()
                     )
                 )
             );
+
+        // The determinant of a matrix over a commutative ring is a polynomial in
+        // its entries, so it is defined wherever the entries are. Computing it by
+        // Gaussian elimination left the pivots as literal divisions, which made the
+        // returned expression undefined wherever a pivot vanishes -- so substituting
+        // an ordinary matrix into a symbolic determinant gave NaN for a determinant
+        // that exists. Laplace expansion never divides, so there is nothing to
+        // exclude. https://github.com/asc-community/AngouriMath/issues/992
+        [Theory]
+        // the 2x2 of the issue: the pivot is x, and x = 0 is an ordinary point
+        [InlineData("[[x, 1], [2, y]]", "x", 0, "y", 5, -2)]
+        [InlineData("[[x, 1], [2, y]]", "x", 3, "y", 4, 10)]
+        // a vanishing pivot in the other position
+        [InlineData("[[1, x], [y, 0]]", "x", 2, "y", 0, 0)]
+        public void SymbolicDeterminantIsDefinedWhereAPivotVanishes(
+            string matrixRaw, string a, int aVal, string b, int bVal, int expected)
+        {
+            Matrix m = matrixRaw;
+            var det = TestExtensions.AsNotNull(m.Determinant);
+            Assert.Equal(expected, det.Substitute(a, aVal).Substitute(b, bVal).EvalNumerical());
+        }
+
+        [Theory]
+        // both of these are singular, and both defeated the divided form: the first
+        // has a zero in the pivot position, the second a vanishing 2x2 leading minor
+        [InlineData("[[0, 1, 2], [3, 4, 5], [6, 7, 8]]", 0)]
+        [InlineData("[[1, 2, 3], [2, 4, 6], [1, 1, 1]]", 0)]
+        [InlineData("[[1, 2, 3], [4, 5, 6], [7, 8, 10]]", -3)]
+        [InlineData("[[2, 1, 0], [1, 2, 1], [0, 1, 2]]", 4)]
+        public void GeneralSymbolicDeterminantAgreesWithTheNumericOne(string matrixRaw, int expected)
+        {
+            Matrix numeric = matrixRaw;
+            Matrix symbolic = "[[a, b, c], [d, e, f], [g, h, i]]";
+            var det = TestExtensions.AsNotNull(symbolic.Determinant);
+            foreach (var (name, value) in new[] { "a", "b", "c", "d", "e", "f", "g", "h", "i" }
+                .Select((name, index) => (name, numeric[index / 3, index % 3])))
+                det = det.Substitute(name, value);
+            Assert.Equal(expected, det.EvalNumerical());
+            Assert.Equal(expected, TestExtensions.AsNotNull(numeric.Determinant).EvalNumerical());
+        }
+
+        [Theory]
+        [InlineData("[[x, 1], [2, y]]", "x * y - 2")]
+        [InlineData("[[a, b], [c, d]]", "a * d - b * c")]
+        public void SymbolicDeterminantCarriesNoCondition(string matrixRaw, string expected)
+        {
+            Matrix m = matrixRaw;
+            var det = TestExtensions.AsNotNull(m.Determinant).Simplify();
+            Assert.DoesNotContain(det.Nodes, node => node is Providedf);
+            Assert.Equal(0, ((Entity)expected - det).Simplify().EvalNumerical());
+        }
     }
 }