From d574b0b68b83cc648b068789745202bd8050facb Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 3 Nov 2026 13:20:42 +0000 Subject: [PATCH] Decide a symbolic pivot before choosing it The row echelon forms were GenericTensor's, whose elimination asked whether a pivot was the literal 0. An entry the elimination had made itself, like c + -c / a * a, is 0 without being written so: it was chosen and divided by, so the reduced form of [[a, b], [c, d]] simplified to NaN, and the rank of [[a, b], [2a, 2b]], whose second row is twice the first, was 2. 2.6.0 does the same. The elimination is written over the entries, with GenericTensor's steps: the first entry of a column that is not 0, from the current row down, for its pivot, rows swapped up, and the reduced form's pivot rows divided by their pivot. Each entry it makes is normalised, a rational function over the rationals to RationalFunction's canonical quotient, where 0 is the literal 0, and anything else is simplified when asked whether it is 0. A symbolic matrix's rank is its rank for general values of its symbols: [[a, b], [2a, 2b]] has rank 1, and [[a, b, c], [d, f, g], [a + d, b + f, c + g]] rank 2. Matrices of numbers come out as they did: the 177 matrix tests pass unchanged, and matrixcheck is clean on 1200 random ones. Part of #2345 Co-Authored-By: Claude Opus 5.5 --- BREAKING-CHANGES.md | 17 +++ .../Core/Entity/Omni/Entity.Matrix.cs | 117 +++++++++++++++--- ...bolicPivotIsDecidedBeforeItIsChosenTest.cs | 67 ++++++++++ 3 files changed, 185 insertions(+), 16 deletions(-) create mode 100644 Sources/Tests/UnitTests/Algebra/ASymbolicPivotIsDecidedBeforeItIsChosenTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 2072155fb..f7c7fe93d 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -445,6 +445,23 @@ value too. No issue was filed; a check of compiled expressions against evaluatio | `"binomial(10000, 1.43)".EvalNumerical()` | throws `CannotEvalException` | `414189.539...` | | `"binomial(pi, 2000)".EvalNumerical()` | `NaN` | `2.10439...E-14` | +### A symbolic matrix's rank and row echelon forms are decided pivot by pivot + +**Silent: a different answer.** The elimination asked whether a pivot was the literal 0. An entry it had made itself, +like `c + -c / a * a`, is 0 without being written so, and it was chosen and divided by: the rank of `[[a, b], [2a, 2b]]`, +whose second row is twice the first, was 2, and the reduced form of every symbolic matrix in the table simplified to +`NaN`. Each entry the elimination makes is now normalised, a rational function over the rationals to its canonical +quotient, where 0 is the literal 0, and anything else is simplified before it is taken for a pivot. A symbolic +matrix's rank is its rank for general values of its symbols. Matrices of numbers are reduced as before. +https://github.com/asc-community/AngouriMath/issues/2345 + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `[[a, b], [2a, 2b]].Rank` | `2` | `1` | +| `[[a, b, c], [d, f, g], [a + d, b + f, c + g]].Rank` | `3` | `2` | +| `[[x, x^2], [1, x]].Rank` | `2` | `1` | +| `[[a, b], [c, d]].ReducedRowEchelonForm.Simplify()` | `NaN` | `[[1, 0], [0, 1]]` | + ### The row echelon forms divide exactly **Not silent for numbers: what threw or did not return now answers. Silent for symbols: the forms are written diff --git a/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs b/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs index fcd071646..dd713cab7 100644 --- a/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs +++ b/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs @@ -583,41 +583,126 @@ public IEnumerator GetEnumerator() /// Matrix's form, transformed via Gaussian elimination. /// /// - /// Both forms divide as they go. GenericTensor's SafeDivision variants keep every entry as a numerator and a - /// denominator that are never reduced against each other, which is for element types whose division truncates, - /// and an Entity divides exactly. Unreduced, the pairs grew with every row operation and the back substitution - /// squared them: one 7x7 of small entries took two minutes, and an 8x8 did not return in ten, and the 4x4 - /// of https://github.com/asc-community/AngouriMath/issues/2270 ended with a denominator near 1e-246, which - /// https://github.com/asc-community/AngouriMath/issues/2277 made a NaN. + /// Both forms divide as they go, each entry an exact quotient: GenericTensor's SafeDivision variants kept + /// every entry as a numerator and a denominator never reduced against each other, the pairs grew with every + /// row operation, and the 4x4 of https://github.com/asc-community/AngouriMath/issues/2270 ended with a + /// denominator near 1e-246, which https://github.com/asc-community/AngouriMath/issues/2277 made a NaN. A + /// pivot is the first entry of its column, from the current row down, that is not zero, and a symbolic + /// entry is simplified to say so; see . /// public Matrix RowEchelonForm => @ref.GetValue(static @this => - OfTheShapeOf(@this, new Matrix(@this.InnerMatrix.RowEchelonFormSimple()).InnerSimplified), + OfTheShapeOf(@this, MathS.Matrix(Eliminate(@this, reduced: false).Entries).InnerSimplified), this); private LazyPropertyA @ref; /// /// Reduced row echelon form via Gaussian elimination. /// + /// + /// Each pivot row is divided by its pivot, so that where a pivot is symbolic the form holds where that + /// pivot is not 0, and the quotients say so as their simplification does. + /// public Matrix ReducedRowEchelonForm => rref.GetValue(static @this => - OfTheShapeOf(@this, new Matrix(@this.InnerMatrix.ReducedRowEchelonFormSimple()).InnerSimplified), + OfTheShapeOf(@this, MathS.Matrix(Eliminate(@this, reduced: true).Entries).InnerSimplified), this); private LazyPropertyA rref; /// - /// The number of linearly independent rows + /// The number of linearly independent rows: of the elimination's pivots. For a symbolic matrix it is the + /// rank for general values of the symbols, as over the field of rational functions in them: the rank of + /// [[a, b], [2 a, 2 b]] is 1, which is 0 where a = b = 0. /// public int Rank => - rank.GetValue(static @this => + rank.GetValue(static @this => Eliminate(@this, reduced: false).Pivots, this); + private LazyPropertyA rank; + + /// + /// Gaussian elimination over the entries, to the row echelon form, or with to the + /// reduced one, and the number of pivots it found. + /// + /// + /// GenericTensor's elimination asked whether a pivot was the literal 0, and an entry the elimination had + /// itself made, like c + -c / a * a, is 0 without being written so: it was chosen and divided by, and + /// the reduced form of [[a, b], [c, d]] simplified to NaN and the rank of [[a, b], [2 a, 2 b]] + /// was 2. Each entry the elimination makes is normalised, a rational function over the rationals to its + /// canonical quotient (), where being 0 is being the literal 0, and + /// the third row of [[a, b, c], [d, e, f], [a + d, b + e, c + f]] is 0 as it should be; anything + /// else is simplified when it is asked whether it is 0. + /// https://github.com/asc-community/AngouriMath/issues/2345 + /// + private static (Entity[,] Entries, int Pivots) Eliminate(Matrix matrix, bool reduced) + { + var (rows, columns) = (matrix.RowCount, matrix.ColumnCount); + var entries = new Entity[rows, columns]; + for (var r = 0; r < rows; r++) + for (var c = 0; c < columns; c++) + entries[r, c] = matrix[r, c]; + var row = 0; + for (var column = 0; column < columns && row < rows; column++) { - for (int i = 0; i < @this.RowCount; i++) - if (@this.ReducedRowEchelonForm.InnerMatrix.RowGetLeadingElement(i) is null) - return i; - return @this.RowCount; + var pivotRow = -1; + for (var r = row; r < rows && pivotRow < 0; r++) + if (IsZeroEntry(entries[r, column])) + entries[r, column] = Number.Integer.Zero; + else + pivotRow = r; + if (pivotRow < 0) + continue; + if (pivotRow != row) + for (var c = 0; c < columns; c++) + (entries[row, c], entries[pivotRow, c]) = (entries[pivotRow, c], entries[row, c]); + var pivot = entries[row, column]; + if (reduced) + { + for (var c = column + 1; c < columns; c++) + entries[row, c] = Normalised(entries[row, c] / pivot); + entries[row, column] = pivot = Number.Integer.One; + } + for (var r = reduced ? 0 : row + 1; r < rows; r++) + { + if (r == row || entries[r, column] is Number.Complex { IsZero: true }) + continue; + var factor = Normalised(entries[r, column] / pivot); + for (var c = column + 1; c < columns; c++) + entries[r, c] = Normalised(entries[r, c] - factor * entries[row, c]); + entries[r, column] = Number.Integer.Zero; + } + row++; } - , this); - private LazyPropertyA rank; + return (entries, row); + } + + /// + /// An entry the elimination made, normalised: a number as it is, a rational function over the rationals as + /// its canonical quotient, which keeps the entries from growing with each row operation, and anything else + /// as its inner simplification. + /// + private static Entity Normalised(Entity entry) + { + var inner = entry.InnerSimplified; + if (inner is Number) + return inner; + return Functions.RationalFunction.TryCanonicalize(inner, out var canonical, wholeCoefficients: true) ? canonical : inner; + } + + /// + /// Whether an entry is 0: a number exactly, a rational function over the rationals by its canonical quotient, + /// and anything else once simplified, a condition on the zero set aside, since an entry that is 0 wherever + /// it has a value is no pivot to divide by. + /// + private static bool IsZeroEntry(Entity entry) + { + if (entry is Number.Complex number) + return number.IsZero; + if (Functions.RationalFunction.TryCanonicalize(entry, out var canonical)) + return canonical is Number.Complex { IsZero: true }; + var simplified = entry.Simplify(); + while (simplified is Providedf(var expression, _)) + simplified = expression; + return simplified is Number.Complex { IsZero: true }; + } /// /// Adjugate form of a matrix diff --git a/Sources/Tests/UnitTests/Algebra/ASymbolicPivotIsDecidedBeforeItIsChosenTest.cs b/Sources/Tests/UnitTests/Algebra/ASymbolicPivotIsDecidedBeforeItIsChosenTest.cs new file mode 100644 index 000000000..b2b01512b --- /dev/null +++ b/Sources/Tests/UnitTests/Algebra/ASymbolicPivotIsDecidedBeforeItIsChosenTest.cs @@ -0,0 +1,67 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using AngouriMath.Extensions; +using Xunit; +using static AngouriMath.Entity; + +namespace AngouriMath.Tests.Algebra +{ + /// + /// The row echelon forms asked whether a pivot was the literal 0, and an entry the elimination had made, like + /// c + -c / a * a, is 0 without being written so: it was chosen and divided by. The reduced form of + /// [[a, b], [c, d]] simplified to NaN, and the rank of [[a, b], [2 a, 2 b]], whose second row is + /// twice the first, was 2. The rank of a symbolic matrix is the rank for general values of its symbols. + /// https://github.com/asc-community/AngouriMath/issues/2345 + /// + [Trait("Area", "Algebra")] + public sealed class ASymbolicPivotIsDecidedBeforeItIsChosenTest + { + [Theory] + [InlineData("[[a, b], [2a, 2b]]", 1)] + [InlineData("[[a, b, c], [d, f, g], [a + d, b + f, c + g]]", 2)] + [InlineData("[[a, b, c], [d, f, g], [2a - d, 2b - f, 2c - g]]", 2)] + [InlineData("[[x, x^2], [1, x]]", 1)] + [InlineData("[[sin(t)^2 + cos(t)^2, 1], [1, 1]]", 1)] + // Independent, to the end. + [InlineData("[[a, b], [c, d]]", 2)] + [InlineData("[[x, 1, 0], [1, x, 1], [0, 1, x]]", 3)] + public void HasItsRank(string matrix, int rank) + => Assert.Equal(rank, ((Matrix)matrix.ToEntity()).Rank); + + [Theory] + [InlineData("[[a, b], [c, d]]")] + [InlineData("[[1, b], [c, d]]")] + public void AnIndependentMatrixReducesToTheIdentity(string matrix) + => Assert.Equal(MathS.IdentityMatrix(2), ((Matrix)matrix.ToEntity()).ReducedRowEchelonForm.Simplify()); + + [Fact] + public void ADependentRowReducesToZeros() + { + var reduced = (Matrix)((Matrix)"[[a, b], [2a, 2b]]".ToEntity()).ReducedRowEchelonForm.Simplify(); + Assert.Equal(Number.Integer.One, reduced[0, 0]); + // b / a, which holds where a is not 0, and says so. + var difference = (reduced[0, 1] - "b / a".ToEntity()).Simplify(); + while (difference is Providedf(var inner, _)) + difference = inner; + Assert.Equal(Number.Integer.Zero, difference); + Assert.Equal(Number.Integer.Zero, reduced[1, 0]); + Assert.Equal(Number.Integer.Zero, reduced[1, 1]); + } + + [Theory] + [InlineData("[[a, b], [c, d]]")] + [InlineData("[[a, b], [2a, 2b]]")] + [InlineData("[[a, b, c], [d, f, g], [a + d, b + f, c + g]]")] + public void NoFormHasANaN(string matrix) + { + var m = (Matrix)matrix.ToEntity(); + Assert.DoesNotContain("NaN", m.RowEchelonForm.Simplify().Stringize()); + Assert.DoesNotContain("NaN", m.ReducedRowEchelonForm.Simplify().Stringize()); + } + } +}