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()); + } + } +}