From e7cbc4eb8ad7ba27d641332bcf195254fe717f71 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 5 Aug 2026 02:34:48 +0000 Subject: [PATCH 1/6] Pull the perfect powers out from under a radical, and out of a denominator sqrt(12) was left as written. The reduced form was generated as a Simplify candidate but lost to the shorter one on node count, so it only ever appeared where two like radicals had to be collected -- sqrt(12) + sqrt(27) came to 5*sqrt(3), while sqrt(12) on its own did not move. That is the wrong judge for this: a reduced radicand is what makes two surds comparable without simplifying their difference, and it is what sympy, Mathematica and Maxima all print. So the reduction now happens where the power is built. Trial division, with two bounds chosen so the extraction is complete rather than best-effort. Once every factor below 10^4 is out, what is left is a prime, a product of two primes above the bound, or the square of one -- three primes above 10^4 would already exceed the 10^12 radicand bound -- and only the last of those has anything to give up, which the perfect-power check at the end takes. So nothing within the bound is missed, and the ordinary case costs a handful of divisions. A negative exponent is the same question asked of the reciprocal, so it is turned round rather than declined, and a radical below the line is answered the same way where the division is built. Without that, 1/sqrt(2) and sqrt(1/2) printed differently -- one number, two forms, only because one of them was refused. A negative radicand is declined. The principal cube root of -8 is 1 + i*sqrt(3), not -2; which branch a negative base takes is settled elsewhere in the library, and reducing the radicand under an even root would be making that same choice by accident. sqrt(12) -> 2 * sqrt(3) sqrt(1/2) -> sqrt(2) / 2 1 / sqrt(2) -> sqrt(2) / 2 1 / sqrt(12) -> sqrt(3) / 6 1 / cbrt(2) -> 4^(1/3) / 2 (2^3 * 5^7)^(1/3) -> 50 * 5^(1/3) 12^(5/2) -> 288 * sqrt(3) A surd against a *sum* -- (5 - sqrt(3)) / (5 + sqrt(3)) -- still needs the conjugate, and is the part of 205 this does not close. Two existing tests moved with it, both because the answer got better: Issue205_StandaloneSurdsKeepTheShorterForm pinned the behaviour this replaces, and Issue662_ExactPathStaysSymbolic expected sqrt(369), where 369 is 9 * 41. https://github.com/asc-community/AngouriMath/issues/281 https://github.com/asc-community/AngouriMath/issues/205 --- ...aluation.Continuous.Arithmetics.Classes.cs | 12 +- .../Functions/NumberTheory/RootExtraction.cs | 168 ++++++++++++++++++ .../Common/SimplificationRegressionTest.cs | 110 ++++++++++-- .../UnitTests/Common/SolverRegressionTest.cs | 6 +- 4 files changed, 283 insertions(+), 13 deletions(-) create mode 100644 Sources/AngouriMath/Functions/NumberTheory/RootExtraction.cs diff --git a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs index 1ef8bca20..4e8d8c827 100644 --- a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs +++ b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs @@ -96,6 +96,13 @@ n0.Evaled is Complex c (_, Integer(0)) => Real.NaN, (Complex n1, Complex n2) when !isExact => n1 / n2, (var n1, Integer(1)) => n1, + // A radical in the denominator is the reciprocal power asked for outright, + // and answering it there is what makes 1/sqrt(2) and sqrt(1/2) the same + // thing on the page as well as in value. + // https://github.com/asc-community/AngouriMath/issues/205 + (var n1, Powf(Rational radicand, Rational exponent)) when + RootExtraction.PullOutOfDenominator(radicand, exponent) is { } rationalized + => (n1 * rationalized).InnerSimplified(isExact), _ => null }, (@this, a, b) => ((Divf)@this).New(a, b), isExact); @@ -182,11 +189,14 @@ x.Evaled is Complex c : 1 : new Providedf(1, !x.EqualTo(0)), (var n1, Integer(1)) => n1, + // https://github.com/asc-community/AngouriMath/issues/281 + (Rational radicand, Rational exponent) when + RootExtraction.PullOutOfRadical(radicand, exponent) is { } pulled => pulled.InnerSimplified(isExact), _ => null }, (@this, a, b) => ((Powf)@this).New(a, b), isExact); } - + public partial record Logf { // Logarithm is undefined when: diff --git a/Sources/AngouriMath/Functions/NumberTheory/RootExtraction.cs b/Sources/AngouriMath/Functions/NumberTheory/RootExtraction.cs new file mode 100644 index 000000000..987cc9b82 --- /dev/null +++ b/Sources/AngouriMath/Functions/NumberTheory/RootExtraction.cs @@ -0,0 +1,168 @@ +// +// Copyright (c) 2019-2022 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using PeterO.Numbers; + +namespace AngouriMath.Functions +{ + using static Entity; + using static Entity.Number; + + /// + /// Writes a rational raised to a rational power as a whole part times a radical that has + /// nothing left to give up, so that sqrt(12) is 2 * sqrt(3) and sqrt(1/2) is sqrt(2) / 2. + /// https://github.com/asc-community/AngouriMath/issues/281 + /// + internal static class RootExtraction + { + // Trial division is the whole method, and these two bounds are what keep it from costing + // more than the tidier answer is worth: a power is simplified over and over, and most of + // them have nothing to give up, so the ordinary case has to be a handful of divisions. + // + // They are also chosen so that the extraction is *complete* rather than best-effort. Once + // every factor below the divisor bound has been taken out, what is left is either one, a + // prime, or a product of two primes above the bound -- three primes above 10^4 would + // already exceed 10^12. A product of two distinct primes has nothing to extract for any + // root, and the one remaining case, a square of a single prime, is what the perfect-power + // check at the end is for. So nothing within the bound is missed. + [ConstantField] private static readonly EInteger maxRadicand = EInteger.FromInt64(1_000_000_000_000L); + [ConstantField] private static readonly EInteger maxDivisor = EInteger.FromInt32(10_000); + + /// + /// as a bit count, rounded up, so that a radicand can be + /// rejected before it is raised to the exponent's numerator. + /// + private const int MaxRadicandBits = 40; + + /// + /// A root of an order beyond this is not one anybody wrote down, and the radicand is + /// raised to the numerator before anything else happens, so both have to stay small. + /// + private const int MaxExponentPart = 64; + + /// + /// The reciprocal of ^ , written + /// without a radical below the line, or null when the exponent is not a root at all. + /// + internal static Entity? PullOutOfDenominator(Rational radicand, Rational exponent) + => PullOutOfRadical(radicand, + Rational.Create(exponent.Numerator.EInteger.Negate(), exponent.Denominator.EInteger)); + + /// + /// Rewrites ^ as an exact factor + /// times a radical, or returns null when there is nothing to pull out. + /// + internal static Entity? PullOutOfRadical(Rational radicand, Rational exponent) + { + // A negative radicand is declined for correctness and not out of caution: the + // principal cube root of -8 is 1 + i * sqrt(3), not -2, so pulling the sign out + // would quietly pick the real branch over the one the rest of the library uses. + if (radicand.Numerator.EInteger.Sign <= 0) + return null; + if (!exponent.Numerator.EInteger.CanFitInInt32() || !exponent.Denominator.EInteger.CanFitInInt32()) + return null; + var power = exponent.Numerator.EInteger.ToInt32Checked(); + var rootIndex = exponent.Denominator.EInteger.ToInt32Checked(); + // A whole exponent is somebody else's case. + if (rootIndex < 2 || power is 0 || System.Math.Abs(power) > MaxExponentPart || rootIndex > MaxExponentPart) + return null; + + // A negative exponent is the same question asked of the reciprocal: 2^(-1/2) is + // (1/2)^(1/2). Turning it round rather than declining it is what keeps 1/sqrt(2) + // and sqrt(1/2) from printing differently, which they did. + var (above, below) = power < 0 + ? (radicand.Denominator.EInteger, radicand.Numerator.EInteger) + : (radicand.Numerator.EInteger, radicand.Denominator.EInteger); + power = System.Math.Abs(power); + + // Checked before the power is taken rather than after. A thousand-digit rational + // raised to the 64th is a sixty-thousand-digit integer, and building one only to + // find it over the bound and throw it away is the whole cost of the call. + if (TooLargeRaised(above, power) || TooLargeRaised(below, power)) + return null; + + var numerator = above.Pow(power); + var denominator = below.Pow(power); + if (numerator.CompareTo(maxRadicand) > 0 || denominator.CompareTo(maxRadicand) > 0) + return null; + + Split(numerator, rootIndex, out var wholeAbove, out var underAbove); + Split(denominator, rootIndex, out var wholeBelow, out var underBelow); + + // A radical left in the denominator is worse than the form we started from, so it is + // moved upstairs the usual way: 1 / b^(1/n) is b^((n-1)/n) / b. What arrives up there + // may itself have a whole part -- 1/8 under a cube root leaves 8^2, which is 4 cubed + // -- so it goes back through the same split rather than being multiplied on as is. + if (!underBelow.Equals(EInteger.One)) + { + var raised = underBelow.Pow(rootIndex - 1); + if (raised.CompareTo(maxRadicand) > 0) + return null; + var moved = underAbove.Multiply(raised); + if (moved.CompareTo(maxRadicand) > 0) + return null; + Split(moved, rootIndex, out var wholeMoved, out underAbove); + wholeAbove = wholeAbove.Multiply(wholeMoved); + wholeBelow = wholeBelow.Multiply(underBelow); + underBelow = EInteger.One; + } + + if (wholeAbove.Equals(EInteger.One) && wholeBelow.Equals(EInteger.One) + && underAbove.Equals(numerator) && underBelow.Equals(denominator)) + return null; + + var whole = Rational.Create(wholeAbove, wholeBelow); + if (underAbove.Equals(EInteger.One)) + return whole; + return whole * new Powf(Integer.Create(underAbove), Rational.Create(EInteger.One, EInteger.FromInt32(rootIndex))); + } + + /// + /// Whether raised to would exceed the + /// radicand bound, decided from the bit length so that the power is never taken. + /// + private static bool TooLargeRaised(EInteger value, int power) + => value.GetUnsignedBitLengthAsEInteger() + .Multiply(power) + .CompareTo(EInteger.FromInt32(MaxRadicandBits)) > 0; + + /// + /// Splits into the part that comes out from under a root of + /// order and the part that stays under it. + /// + private static void Split(EInteger value, int rootIndex, out EInteger whole, out EInteger under) + { + whole = EInteger.One; + under = EInteger.One; + var remaining = value; + for (var divisor = EInteger.FromInt32(2); + divisor.CompareTo(maxDivisor) <= 0 && divisor.Multiply(divisor).CompareTo(remaining) <= 0; + divisor = divisor.Add(EInteger.One)) + { + var multiplicity = 0; + while (remaining.Remainder(divisor).IsZero) + { + remaining = remaining.Divide(divisor); + multiplicity++; + } + if (multiplicity is 0) + continue; + whole = whole.Multiply(divisor.Pow(multiplicity / rootIndex)); + under = under.Multiply(divisor.Pow(multiplicity % rootIndex)); + } + if (remaining.Equals(EInteger.One)) + return; + // What survives trial division is prime or a product of two primes, so the only thing + // it can still hide is being an exact power itself. + var root = remaining.Root(rootIndex); + if (root.Pow(rootIndex).Equals(remaining)) + whole = whole.Multiply(root); + else + under = under.Multiply(remaining); + } + } +} diff --git a/Sources/Tests/UnitTests/Common/SimplificationRegressionTest.cs b/Sources/Tests/UnitTests/Common/SimplificationRegressionTest.cs index f9a9120ad..e69719a8d 100644 --- a/Sources/Tests/UnitTests/Common/SimplificationRegressionTest.cs +++ b/Sources/Tests/UnitTests/Common/SimplificationRegressionTest.cs @@ -44,16 +44,41 @@ public void Issue205_LikeRadicalsCollect(string input, string expected) Assert.DoesNotContain("-", simplified.Stringize()); } - // A standalone surd is a separate matter: Simplify picks between candidate - // forms by a complexity metric, and by that metric `sqrt(12)` beats - // `2 * sqrt(3)`, so the reduced form is generated but not selected. Changing - // that is a change to the metric, not to this rule. Pinning the current - // behaviour so the distinction is deliberate rather than accidental. - [Theory] - [InlineData("sqrt(12)")] - [InlineData("sqrt(18)")] - public void Issue205_StandaloneSurdsKeepTheShorterForm(string input) => - Assert.Equal(input.ToEntity(), input.ToEntity().Simplify()); + // https://github.com/asc-community/AngouriMath/issues/281 + // A standalone surd used to keep the form it was written in, because Simplify + // picks between candidates by a complexity metric and by that metric `sqrt(12)` + // beats `2 * sqrt(3)`. The metric is the wrong judge here: the reduced radicand + // is the canonical form every other system prints -- sympy, Mathematica and + // Maxima all answer 2*sqrt(3) -- and it is what makes two surds comparable + // without simplifying their difference. So the reduction now happens where the + // power itself is built rather than being offered as a candidate. + [Theory] + [InlineData("sqrt(12)", "2 * sqrt(3)")] + [InlineData("sqrt(18)", "3 * sqrt(2)")] + [InlineData("sqrt(1/2)", "sqrt(2) / 2")] + [InlineData("(2 ^ 3 * 5 ^ 7) ^ (1/3)", "50 * 5 ^ (1/3)")] + [InlineData("12 ^ (5/2)", "288 * sqrt(3)")] + // A negative exponent is the same question asked of the reciprocal, and it has to give + // the same answer: 1/sqrt(2) and sqrt(1/2) are one number and printed differently only + // because one of them was declined. + [InlineData("1 / sqrt(2)", "sqrt(2) / 2")] + [InlineData("2 ^ (-1/2)", "sqrt(2) / 2")] + [InlineData("1 / sqrt(12)", "sqrt(3) / 6")] + [InlineData("1 / sqrt(8)", "sqrt(2) / 4")] + [InlineData("1 / cbrt(2)", "4 ^ (1/3) / 2")] + [InlineData("x / sqrt(2)", "x * sqrt(2) / 2")] + [InlineData("8 ^ (-1/3)", "1/2")] + public void Issue281_PerfectPowersComeOutFromUnderTheRoot(string input, string expected) => + Assert.Equal(Entity.Number.Integer.Create(0), + (input.ToEntity().Simplify() - expected.ToEntity()).Simplify()); + + // A negative radicand is not reduced: the extraction would have to choose a branch, + // and which branch a negative base takes is already settled elsewhere -- the odd + // roots take the real one, so (-8)^(1/3) is -2 and not 1 + i*sqrt(3). Reducing the + // radicand under an even root would be picking that same choice by accident. + [Fact] + public void Issue281_NegativeRadicandsAreNotReduced() => + Assert.Equal("sqrt(-12)", "(-12) ^ (1/2)".ToEntity().Simplify().Stringize()); // Whatever form comes out, it must still be the same number. [Theory] @@ -215,5 +240,70 @@ public void Issue178_FactorizationPreservesValue(string input) [InlineData("x * ln(x)")] public void AntiderivativesNeverContainNaN(string integrand) => Assert.DoesNotContain("NaN", integrand.ToEntity().Integrate("x").Stringize()); + + // https://github.com/asc-community/AngouriMath/issues/569 + // https://github.com/asc-community/AngouriMath/issues/179 + // The trigonometric tables were only ever read forwards, so sin(pi/4) was sqrt(2)/2 + // but arcsin(sqrt(2)/2) stayed as written. arctan was the exception, and only at + // +-1, which is why arctan(1) already answered pi/4. + [Theory] + [InlineData("arcsin(1/2)", "pi / 6")] + [InlineData("arcsin(sqrt(2) / 2)", "pi / 4")] + [InlineData("arcsin(sqrt(3) / 2)", "pi / 3")] + [InlineData("arcsin(1)", "pi / 2")] + [InlineData("arcsin(0)", "0")] + [InlineData("arccos(1/2)", "pi / 3")] + [InlineData("arccos(0)", "pi / 2")] + [InlineData("arccos(1)", "0")] + [InlineData("arctan(sqrt(3))", "pi / 3")] + [InlineData("arctan(sqrt(3) / 3)", "pi / 6")] + [InlineData("arccotan(1)", "pi / 4")] + [InlineData("arccotan(sqrt(3))", "pi / 6")] + // The 5ths, 8ths, 10ths and 12ths of the circle, which are in the forward table too + [InlineData("arcsin((sqrt(5) - 1) / 4)", "pi / 10")] + [InlineData("arcsin(sqrt(2 - sqrt(2)) / 2)", "pi / 8")] + [InlineData("arcsin((sqrt(6) - sqrt(2)) / 4)", "pi / 12")] + [InlineData("arctan(2 - sqrt(3))", "pi / 12")] + [InlineData("arctan(sqrt(2) - 1)", "pi / 8")] + [InlineData("arctan(sqrt(5 - 2 * sqrt(5)))", "pi / 5")] + // Both functions are odd, and arccos/arccotan take the complement + [InlineData("arcsin(-1/2)", "-pi / 6")] + [InlineData("arctan(-sqrt(3))", "-pi / 3")] + [InlineData("arccos(-1/2)", "2 * pi / 3")] + [InlineData("arccos(-1)", "pi")] + public void Issue569_InverseTrigonometryReadsTheTableBackwards(string input, string expected) => + Assert.Equal(Entity.Number.Integer.Create(0), + (input.ToEntity().Simplify() - expected.ToEntity()).Simplify()); + + // A value that is merely near a table entry is not that entry, and one that is not in + // the table at all has no closed form to give. `arcsin(0.4999999)` in particular must + // not come back pi/6: the double comparison is a sieve, not the answer. + [Theory] + [InlineData("arcsin(4999999/10000000)", "arcsin")] + [InlineData("arcsin(1/3)", "arcsin")] + [InlineData("arctan(2)", "arctan")] + [InlineData("arcsin(x)", "arcsin")] + [InlineData("arcsin(2)", "arcsin")] // outside the real branch + public void Issue569_NonTableValuesAreLeftAlone(string input, string stillThere) + { + var simplified = input.ToEntity().Simplify().Stringize(); + Assert.Contains(stillThere, simplified); + Assert.DoesNotContain("pi", simplified); + } + + // And the angle that comes back must be the angle that was asked for. + [Theory] + [InlineData("arcsin(1/2)")] + [InlineData("arccos(1/2)")] + [InlineData("arctan(sqrt(3))")] + [InlineData("arccotan(1)")] + [InlineData("arcsin(-sqrt(2) / 2)")] + [InlineData("arccos(-sqrt(3) / 2)")] + public void Issue569_RecognisedAnglesKeepTheirValue(string input) + { + var before = input.ToEntity().EvalNumerical().RealPart.EDecimal.ToDouble(); + var after = input.ToEntity().Simplify().EvalNumerical().RealPart.EDecimal.ToDouble(); + Assert.Equal(before, after, 12); + } } } diff --git a/Sources/Tests/UnitTests/Common/SolverRegressionTest.cs b/Sources/Tests/UnitTests/Common/SolverRegressionTest.cs index 51fd59251..8583ab220 100644 --- a/Sources/Tests/UnitTests/Common/SolverRegressionTest.cs +++ b/Sources/Tests/UnitTests/Common/SolverRegressionTest.cs @@ -102,10 +102,12 @@ public void Issue662_EvaledIsANumber() => Assert.IsAssignableFrom("(|[12,15]|)".ToEntity().Evaled); // The exact path must keep its symbolic answer -- this fix is only about the - // numeric one. + // numeric one. The radicand comes out reduced (369 = 9 * 41) since perfect powers + // are pulled out from under a root: + // https://github.com/asc-community/AngouriMath/issues/281 [Fact] public void Issue662_ExactPathStaysSymbolic() => - Assert.Equal("sqrt(369)".ToEntity(), "(|[12,15]|)".ToEntity().InnerSimplified); + Assert.Equal("3 * sqrt(41)".ToEntity(), "(|[12,15]|)".ToEntity().InnerSimplified); // Not filed upstream; found while measuring solver coverage. // `lim x->0 (1+x)^(1/x)` answered 1 instead of e. The second remarkable limit From a5c34c56cd568dfb5ef1c2be3f7d3f88a0e13f63 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 5 Aug 2026 02:34:58 +0000 Subject: [PATCH 2/6] Read the trigonometric tables backwards sin(pi/4) was sqrt(2)/2 but arcsin(sqrt(2)/2) stayed as written: the tables were only ever read forwards. arctan was the exception and only at +-1, which is why arctan(1) already answered pi/4. The inverse table is written out rather than derived from the forward one at run time, because inverting that one needs the principal branch and it holds every angle around the circle: sin(2pi/3) and sin(2pi/6) are the same number and only the second is an answer arcsin may give. arccos and arccotan take the complement, which is the same reading of arccotan the existing arctan(x) + arccotan(x) rule already takes. A double comparison picks the candidate and the exact value confirms it, so arcsin(0.4999999) comes back unanswered rather than as pi/6. arcsin(1/2) -> pi / 6 arcsin(sqrt(2) / 2) -> pi / 4 arccos(1/2) -> pi / 3 arccos(-1/2) -> 2 * pi / 3 arctan(sqrt(3)) -> pi / 3 arctan(2 - sqrt(3)) -> pi / 12 arccotan(1) -> pi / 4 arcsin(sqrt(2 - sqrt(2)) / 2) -> pi / 8 https://github.com/asc-community/AngouriMath/issues/569 https://github.com/asc-community/AngouriMath/issues/179 --- ...luation.Continuous.Trigonometry.Classes.cs | 4 + .../InverseTrigonometricTableValues.cs | 134 ++++++++++++++++++ 2 files changed, 138 insertions(+) create mode 100644 Sources/AngouriMath/Functions/Simplification/InverseTrigonometricTableValues.cs diff --git a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Trigonometry.Classes.cs b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Trigonometry.Classes.cs index 0cf0fc667..b85a34c66 100644 --- a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Trigonometry.Classes.cs +++ b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Trigonometry.Classes.cs @@ -192,6 +192,7 @@ protected override Entity InnerSimplify(bool isExact) => a => a switch { Complex n when !isExact => Number.Arcsin(n), + { Evaled: Complex n } when isExact && InverseTrigonometricTableValues.PullArcsin(n, out var res) => res, _ => null }, (@this, a) => ((Arcsinf)@this).New(a), isExact); @@ -209,6 +210,7 @@ protected override Entity InnerSimplify(bool isExact) => a => a switch { Complex n when !isExact => Number.Arccos(n), + { Evaled: Complex n } when isExact && InverseTrigonometricTableValues.PullArccos(n, out var res) => res, _ => null }, (@this, a) => ((Arccosf)@this).New(a), isExact); @@ -232,6 +234,7 @@ Real r when r.EDecimal.IsNegativeInfinity() => -MathS.pi / 2, // pi/4 stopped one step short of it. Integer(1) => MathS.pi / 4, Integer(-1) => -MathS.pi / 4, + { Evaled: Complex n } when isExact && InverseTrigonometricTableValues.PullArctan(n, out var res) => res, _ => null }, (@this, a) => ((Arctanf)@this).New(a), isExact); @@ -247,6 +250,7 @@ protected override Entity InnerSimplify(bool isExact) => a => a switch { Complex n when !isExact => Number.Arccotan(n), + { Evaled: Complex n } when isExact && InverseTrigonometricTableValues.PullArccotan(n, out var res) => res, _ => null }, (@this, a) => ((Arccotanf)@this).New(a), isExact); diff --git a/Sources/AngouriMath/Functions/Simplification/InverseTrigonometricTableValues.cs b/Sources/AngouriMath/Functions/Simplification/InverseTrigonometricTableValues.cs new file mode 100644 index 000000000..726319e53 --- /dev/null +++ b/Sources/AngouriMath/Functions/Simplification/InverseTrigonometricTableValues.cs @@ -0,0 +1,134 @@ +// +// Copyright (c) 2019-2022 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System.Diagnostics.CodeAnalysis; + +namespace AngouriMath.Functions +{ + using static Entity; + using static Entity.Number; + + /// + /// The trigonometric tables read backwards: given a value, the angle in the principal + /// range whose sine or tangent it is. + /// https://github.com/asc-community/AngouriMath/issues/569 + /// https://github.com/asc-community/AngouriMath/issues/179 + /// + internal static class InverseTrigonometricTableValues + { + /// + /// The double is a sieve and nothing more -- every hit is confirmed against the exact + /// value -- so this only has to be tighter than the gap between two table entries, + /// which is about a hundredth. + /// + private const double SieveWidth = 1e-9; + + private static Entity Sqrt(Entity a) => MathS.Sqrt(a); + [ConstantField] private static readonly Entity pi = MathS.pi; + + /// + /// Written out rather than inverted from at run + /// time, because inverting it needs the principal branch and that table holds every + /// angle around the circle: sin(2pi/3) and sin(2pi/6) are the same number, and only the + /// second of them is an answer arcsin may give. + /// + [ConstantField] private static readonly List<(double approximation, Entity value, Entity angle)> arcsinTable = new() + { + (0d, 0, (Entity)0), + (0.25881904510252074d, (Sqrt(6) - Sqrt(2)) / 4, pi / 12), + (0.3090169943749474d, (Sqrt(5) - 1) / 4, pi / 10), + (0.3826834323650898d, Sqrt(2 - Sqrt(2)) / 2, pi / 8), + (0.5d, Rational.Create(1, 2), pi / 6), + (0.5877852522924731d, Sqrt(10 - 2 * Sqrt(5)) / 4, pi / 5), + (0.7071067811865476d, Sqrt(2) / 2, pi / 4), + (0.8090169943749475d, (Sqrt(5) + 1) / 4, 3 * pi / 10), + (0.8660254037844386d, Sqrt(3) / 2, pi / 3), + (0.9238795325112867d, Sqrt(2 + Sqrt(2)) / 2, 3 * pi / 8), + (0.9510565162951535d, Sqrt(10 + 2 * Sqrt(5)) / 4, 2 * pi / 5), + (0.9659258262890683d, (Sqrt(6) + Sqrt(2)) / 4, 5 * pi / 12), + (1d, 1, pi / 2) + }; + + [ConstantField] private static readonly List<(double approximation, Entity value, Entity angle)> arctanTable = new() + { + (0d, 0, (Entity)0), + (0.2679491924311227d, 2 - Sqrt(3), pi / 12), + (0.3249196962329063d, Sqrt(25 - 10 * Sqrt(5)) / 5, pi / 10), + (0.41421356237309503d, Sqrt(2) - 1, pi / 8), + (0.5773502691896257d, Sqrt(3) / 3, pi / 6), + (0.7265425280053609d, Sqrt(5 - 2 * Sqrt(5)), pi / 5), + (1d, 1, pi / 4), + (1.3763819204711736d, Sqrt(25 + 10 * Sqrt(5)) / 5, 3 * pi / 10), + (1.7320508075688772d, Sqrt(3), pi / 3), + (2.414213562373095d, Sqrt(2) + 1, 3 * pi / 8), + (3.0776835371752527d, Sqrt(5 + 2 * Sqrt(5)), 2 * pi / 5), + (3.7320508075688776d, 2 + Sqrt(3), 5 * pi / 12) + }; + + internal static bool PullArcsin(Complex arg, [NotNullWhen(true)] out Entity? res) + => TryPull(arcsinTable, arg, out res); + + internal static bool PullArctan(Complex arg, [NotNullWhen(true)] out Entity? res) + => TryPull(arctanTable, arg, out res); + + /// + /// arccos(x) = pi/2 - arcsin(x), which is also the range arccos is stated over. + /// + internal static bool PullArccos(Complex arg, [NotNullWhen(true)] out Entity? res) + => Complement(PullArcsin(arg, out var arcsin), arcsin, out res); + + /// + /// arccotan(x) = pi/2 - arctan(x), the same reading of it the simplification rule for + /// arctan(x) + arccotan(x) already takes. + /// + internal static bool PullArccotan(Complex arg, [NotNullWhen(true)] out Entity? res) + => Complement(PullArctan(arg, out var arctan), arctan, out res); + + private static bool Complement(bool found, Entity? angle, [NotNullWhen(true)] out Entity? res) + { + if (!found || angle is null) + { + res = null; + return false; + } + res = (pi / 2 - angle).InnerSimplified; + return true; + } + + private static bool TryPull(List<(double approximation, Entity value, Entity angle)> table, + Complex arg, [NotNullWhen(true)] out Entity? res) + { + res = null; + // Only the real branch. The table says nothing about arcsin(2), and answering the + // complex cases from it would mean choosing a branch cut here rather than where the + // numeric functions already choose one. + if (arg is not Real { EDecimal: var given } || !given.IsFinite) + return false; + var approximation = given.ToDouble(); + var negative = approximation < 0; + if (negative) + approximation = -approximation; + foreach (var (candidate, value, angle) in table) + { + if (System.Math.Abs(candidate - approximation) > SieveWidth) + continue; + // Near a table value is not the same as being one: 0.4999999 must come back + // unanswered rather than as pi/6, so the sieve hit is checked against the exact + // value at the precision the caller asked for. + if (value.Evaled is not Real { EDecimal: var exact }) + return false; + if (!Number.IsZero(negative ? exact + given : exact - given)) + return false; + // Both functions in the table are odd, so a negative argument is the same angle + // the other way round. + res = negative ? (-angle).InnerSimplified : angle; + return true; + } + return false; + } + } +} From d3f4c9eac21a19753f597fe597951e1b4364cfd0 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 5 Aug 2026 02:34:58 +0000 Subject: [PATCH 3/6] Give a square matrix its characteristic polynomial det(x*I - A), whose roots are the eigenvalues. The determinant was already there, so this is the shift and the sign convention. Monic -- det(x*I - A) rather than det(A - x*I). The two are the same polynomial up to a sign for an odd-sized matrix, and the systems disagree: Mathematica prints det(A - x*I), sympy's charpoly is monic. Monic is what most linear algebra texts state and it is the form that composes with the polynomial solver without a stray -1 in front, so it is the one documented and tested here. [[1, 2], [3, 4]] -> x^2 - 5x - 2 [[a, b], [c, d]] -> x^2 - (a + d)x + ad - bc [[1, 0, 0], [0, 1, 0], [0, 0, 1]] -> x^3 - 3x^2 + 3x - 1 https://github.com/asc-community/AngouriMath/issues/381 --- Sources/AngouriMath/Convenience/MathS.cs | 24 +++++++ .../Core/Entity/Omni/Entity.Matrix.cs | 33 +++++++++ .../Core/CharacteristicPolynomialTest.cs | 72 +++++++++++++++++++ 3 files changed, 129 insertions(+) create mode 100644 Sources/Tests/UnitTests/Core/CharacteristicPolynomialTest.cs diff --git a/Sources/AngouriMath/Convenience/MathS.cs b/Sources/AngouriMath/Convenience/MathS.cs index ea15c8657..d44d97fa8 100644 --- a/Sources/AngouriMath/Convenience/MathS.cs +++ b/Sources/AngouriMath/Convenience/MathS.cs @@ -4419,6 +4419,30 @@ public static Rational CreateRational(EInteger numerator, EInteger denominator) public static Entity? Det(Matrix m) => m.Determinant; + /// + /// Returns the characteristic polynomial det(x*I - m) of a square matrix, whose roots + /// are its eigenvalues. https://github.com/asc-community/AngouriMath/issues/381 + /// + /// A square matrix + /// The variable the polynomial is written in + /// + /// + /// Matrix a = "[[1, 2], [3, 4]]"; + /// Console.WriteLine(CharacteristicPolynomial(a, "x")); + /// Console.WriteLine(CharacteristicPolynomial(a, "x").SolveEquation("x")); + /// + /// Prints + /// + /// (1 - x) * (4 - x) - 6 + /// { 5/2 - sqrt(33) / 2, 5/2 + sqrt(33) / 2 } + /// + /// + /// + /// Thrown for a non-square matrix + /// + public static Entity CharacteristicPolynomial(Matrix m, Variable variable) + => m.CharacteristicPolynomial(variable); + /// Creates an instance of . /// /// A two-dimensional array of values. diff --git a/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs b/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs index ff8786a1d..24da6f814 100644 --- a/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs +++ b/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs @@ -250,6 +250,39 @@ public Entity AsScalar() }, this); private LazyPropertyA inverse; + /// + /// The characteristic polynomial det(x*I - A), whose roots are the eigenvalues + /// of the matrix. https://github.com/asc-community/AngouriMath/issues/381 + /// + /// + /// The variable the polynomial is written in + /// + /// + /// A polynomial of degree , inner simplified. It is monic: + /// det(x*I - A) rather than det(A - x*I), which are the same polynomial up to a + /// sign for an odd-sized matrix. The eigenvalues are the roots either way, and the + /// monic form is the one sympy and most linear algebra texts state, so it is the + /// one that composes with the polynomial solver without a stray -1 in front. + /// + /// + /// Thrown for a non-square matrix, which has no characteristic polynomial + /// + public Entity CharacteristicPolynomial(Variable variable) + { + if (!IsSquare) + throw new InvalidMatrixOperationException( + "Only a square matrix has a characteristic polynomial"); + var shifted = InnerMatrix.Copy(true); + for (var i = 0; i < RowCount; i++) + shifted[i, i] = (variable - shifted[i, i]).InnerSimplified; + for (var row = 0; row < RowCount; row++) + for (var column = 0; column < ColumnCount; column++) + if (row != column) + shifted[row, column] = (-shifted[row, column]).InnerSimplified; + return new Matrix(shifted).Determinant + ?? throw new AngouriBugException("A square matrix always has a determinant"); + } + /// /// The Add operator. Performs an active operation diff --git a/Sources/Tests/UnitTests/Core/CharacteristicPolynomialTest.cs b/Sources/Tests/UnitTests/Core/CharacteristicPolynomialTest.cs new file mode 100644 index 000000000..1d6e79f96 --- /dev/null +++ b/Sources/Tests/UnitTests/Core/CharacteristicPolynomialTest.cs @@ -0,0 +1,72 @@ +// +// Copyright (c) 2019-2022 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using AngouriMath; +using AngouriMath.Core.Exceptions; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Core +{ + /// + /// https://github.com/asc-community/AngouriMath/issues/381 + /// + public sealed class CharacteristicPolynomialTest + { + private static Entity Expanded(Entity.Matrix matrix, string variable) + => MathS.CharacteristicPolynomial(matrix, variable).Expand().Simplify(); + + // det(x*I - A) written out by hand for each of these. + [Theory] + [InlineData("[[1, 2], [3, 4]]", "x ^ 2 - 5 * x - 2")] + [InlineData("[[2, 0], [0, 3]]", "x ^ 2 - 5 * x + 6")] + [InlineData("[[0, 1], [-1, 0]]", "x ^ 2 + 1")] + [InlineData("[[1, 0, 0], [0, 1, 0], [0, 0, 1]]", "x ^ 3 - 3 * x ^ 2 + 3 * x - 1")] + [InlineData("[[2, 1], [1, 2]]", "x ^ 2 - 4 * x + 3")] + public void MatchesTheHandComputedPolynomial(string matrix, string expected) => + Assert.Equal(Entity.Number.Integer.Create(0), + (Expanded((Entity.Matrix)matrix.ToEntity(), "x") - expected.ToEntity()).Expand().Simplify()); + + // The point of the polynomial: its roots are the eigenvalues. A diagonal matrix wears + // them on its face, so it is the case that can be checked without solving anything. + [Theory] + [InlineData("[[5, 0], [0, 7]]", 5, 7)] + [InlineData("[[-1, 0], [0, 4]]", -1, 4)] + public void DiagonalEntriesAreRoots(string matrix, int first, int second) + { + var polynomial = MathS.CharacteristicPolynomial((Entity.Matrix)matrix.ToEntity(), "x"); + Assert.Equal(Entity.Number.Integer.Create(0), polynomial.Substitute("x", first).Simplify()); + Assert.Equal(Entity.Number.Integer.Create(0), polynomial.Substitute("x", second).Simplify()); + } + + // Monic, which is what distinguishes det(x*I - A) from det(A - x*I): for an odd-sized + // matrix the other convention would put a -1 in front of the leading term. + [Theory] + [InlineData("[[1, 2], [3, 4]]", 2)] + [InlineData("[[1, 2, 3], [4, 5, 6], [7, 8, 10]]", 3)] + public void IsMonic(string matrix, int degree) + { + var polynomial = Expanded((Entity.Matrix)matrix.ToEntity(), "x"); + // The leading coefficient is what survives when the polynomial is divided by x^n + // and x is sent to infinity. + Assert.Equal(Entity.Number.Integer.Create(1), + (polynomial / MathS.Pow("x", degree)).Limit("x", Entity.Number.Real.PositiveInfinity).Simplify()); + } + + // Symbolic entries are as good as numeric ones: this is det(x*I - A) either way. + [Fact] + public void WorksOnSymbolicEntries() => + Assert.Equal(Entity.Number.Integer.Create(0), + (Expanded((Entity.Matrix)"[[a, b], [c, d]]".ToEntity(), "x") + - "x ^ 2 - (a + d) * x + a * d - b * c".ToEntity()).Expand().Simplify()); + + [Fact] + public void NonSquareMatrixIsRefused() => + Assert.Throws( + () => MathS.CharacteristicPolynomial((Entity.Matrix)"[[1, 2, 3], [4, 5, 6]]".ToEntity(), "x")); + } +} From 9bb8fe17423a52d83bfe404a8b838887c042fdda Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 5 Aug 2026 02:34:58 +0000 Subject: [PATCH 4/6] Collect coverage from the F# and Interactive suites too Only the C# suite reported to codecov, so the wrappers were invisible in it. Both now collect on the ubuntu leg of their matrix and upload with a flag, which keeps the three distinguishable rather than merging into one total; the other two operating systems keep running the tests without paying for the collector or racing each other's upload, which is the shape the C# workflow already used. C++ is deliberately not included. AngouriMath.CPP is a native wrapper over the AOT-exported library, so its coverage would have to come from gcov or llvm-cov and would measure the wrapper rather than the library. That is a separate decision about whether a coverage build is worth having in CI at all, so the issue stays open for it. https://github.com/asc-community/AngouriMath/issues/397 --- .github/workflows/CSharpTest.yml | 1 + .github/workflows/FSharpTest.yml | 13 +++++++++++++ .github/workflows/InteractiveTest.yml | 17 ++++++++++++++++- 3 files changed, 30 insertions(+), 1 deletion(-) diff --git a/.github/workflows/CSharpTest.yml b/.github/workflows/CSharpTest.yml index 325ae3526..bd97620b2 100644 --- a/.github/workflows/CSharpTest.yml +++ b/.github/workflows/CSharpTest.yml @@ -62,3 +62,4 @@ jobs: uses: codecov/codecov-action@v5 with: file: ./Sources/Tests/UnitTests/coverage.opencover.xml + flags: csharp diff --git a/.github/workflows/FSharpTest.yml b/.github/workflows/FSharpTest.yml index c09fd4fdb..384d5a38e 100644 --- a/.github/workflows/FSharpTest.yml +++ b/.github/workflows/FSharpTest.yml @@ -34,5 +34,18 @@ jobs: dotnet build Sources/Tests/FSharpWrapperUnitTests - name: Test + if: ${{ matrix.os != 'ubuntu-latest' }} run: | dotnet test "Sources/Tests/FSharpWrapperUnitTests" -c Release -l GitHubActions + + - name: Test with coverage + if: ${{ matrix.os == 'ubuntu-latest' }} + run: | + dotnet test "Sources/Tests/FSharpWrapperUnitTests" -c Release -l GitHubActions --collect:"XPlat Code Coverage" --results-directory ./coverage -- DataCollectionRunSettings.DataCollectors.DataCollector.Configuration.Format=opencover + + - name: 'Send to codecov' + if: ${{ matrix.os == 'ubuntu-latest' }} + uses: codecov/codecov-action@v5 + with: + directory: ./coverage + flags: fsharp diff --git a/.github/workflows/InteractiveTest.yml b/.github/workflows/InteractiveTest.yml index d301aab3a..37b8b4555 100644 --- a/.github/workflows/InteractiveTest.yml +++ b/.github/workflows/InteractiveTest.yml @@ -33,9 +33,24 @@ jobs: dotnet build Sources/Tests/TerminalUnitTests -c Release - name: Test Interactive Wrapper + if: ${{ matrix.os != 'ubuntu-latest' }} run: | dotnet test "Sources/Tests/InteractiveWrapperUnitTests" -c Release - + - name: Test Terminal Lib + if: ${{ matrix.os != 'ubuntu-latest' }} run: | dotnet test "Sources/Tests/TerminalUnitTests" -c Release + + - name: Test with coverage + if: ${{ matrix.os == 'ubuntu-latest' }} + run: | + dotnet test "Sources/Tests/InteractiveWrapperUnitTests" -c Release --collect:"XPlat Code Coverage" --results-directory ./coverage -- DataCollectionRunSettings.DataCollectors.DataCollector.Configuration.Format=opencover + dotnet test "Sources/Tests/TerminalUnitTests" -c Release --collect:"XPlat Code Coverage" --results-directory ./coverage -- DataCollectionRunSettings.DataCollectors.DataCollector.Configuration.Format=opencover + + - name: 'Send to codecov' + if: ${{ matrix.os == 'ubuntu-latest' }} + uses: codecov/codecov-action@v5 + with: + directory: ./coverage + flags: interactive From e38103a087d46e6b28a71c3c768a1b777b027dea Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 5 Aug 2026 02:34:58 +0000 Subject: [PATCH 5/6] Note in AGENTS.md that the tracker is older than the code Eleven open issues turned out to be already fixed when re-measured on 2026-08-05, and re-measuring cost minutes where planning fixes for them would have cost days. Worth saying next to the rule about reproducing before claiming something is broken, since it is the same measurement taken earlier. --- AGENTS.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/AGENTS.md b/AGENTS.md index b66a3aba4..d450b40b3 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -83,6 +83,12 @@ go, named for the issue. **Reproduce on a stock `master` build before claiming anything is broken.** More than one confident report here has turned out to describe behaviour that was never broken. +The same measurement is worth taking before you plan a fix at all. An open issue names the version +its reporter was on, and the tracker is older than the code: on 2026-08-05, eleven of the open +issues turned out to be already fixed, and re-measuring them cost minutes where planning fixes for +them would have cost days. Run the reporter's expression first. Close it with the measurement if it +answers, and say which build you measured. + ## Aim high The ambition is to answer the hard problems: olympiad and competition algebra, calculus and number From 045ad780c1acf876a1644a716d9d840acba4f628 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 5 Aug 2026 03:33:59 +0000 Subject: [PATCH 6/6] Take the characteristic polynomial back out det(x*I - A) came back NaN at particular values of x rather than the number the polynomial takes there. The determinant is computed by Gaussian elimination, which divides, so it carries the conditions under which those divisions were legal -- and on the shifted matrix those pivots are polynomials in x that vanish at isolated points. Every one of them is a removable singularity of something defined everywhere, so the conditions are right of the determinant and wrong of the characteristic polynomial, and clearing them is not a matter of stripping them off afterwards. Getting that right means either a Laplace expansion for this case or cancelling the conditions against the factors that produced them. Neither is a line of code, and a polynomial that is NaN at some of its own points is worse than not having one, so it comes out rather than shipping. https://github.com/asc-community/AngouriMath/issues/381 --- Sources/AngouriMath/Convenience/MathS.cs | 24 ------- .../Core/Entity/Omni/Entity.Matrix.cs | 33 --------- .../Core/CharacteristicPolynomialTest.cs | 72 ------------------- 3 files changed, 129 deletions(-) delete mode 100644 Sources/Tests/UnitTests/Core/CharacteristicPolynomialTest.cs diff --git a/Sources/AngouriMath/Convenience/MathS.cs b/Sources/AngouriMath/Convenience/MathS.cs index d44d97fa8..ea15c8657 100644 --- a/Sources/AngouriMath/Convenience/MathS.cs +++ b/Sources/AngouriMath/Convenience/MathS.cs @@ -4419,30 +4419,6 @@ public static Rational CreateRational(EInteger numerator, EInteger denominator) public static Entity? Det(Matrix m) => m.Determinant; - /// - /// Returns the characteristic polynomial det(x*I - m) of a square matrix, whose roots - /// are its eigenvalues. https://github.com/asc-community/AngouriMath/issues/381 - /// - /// A square matrix - /// The variable the polynomial is written in - /// - /// - /// Matrix a = "[[1, 2], [3, 4]]"; - /// Console.WriteLine(CharacteristicPolynomial(a, "x")); - /// Console.WriteLine(CharacteristicPolynomial(a, "x").SolveEquation("x")); - /// - /// Prints - /// - /// (1 - x) * (4 - x) - 6 - /// { 5/2 - sqrt(33) / 2, 5/2 + sqrt(33) / 2 } - /// - /// - /// - /// Thrown for a non-square matrix - /// - public static Entity CharacteristicPolynomial(Matrix m, Variable variable) - => m.CharacteristicPolynomial(variable); - /// Creates an instance of . /// /// A two-dimensional array of values. diff --git a/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs b/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs index 24da6f814..ff8786a1d 100644 --- a/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs +++ b/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs @@ -250,39 +250,6 @@ public Entity AsScalar() }, this); private LazyPropertyA inverse; - /// - /// The characteristic polynomial det(x*I - A), whose roots are the eigenvalues - /// of the matrix. https://github.com/asc-community/AngouriMath/issues/381 - /// - /// - /// The variable the polynomial is written in - /// - /// - /// A polynomial of degree , inner simplified. It is monic: - /// det(x*I - A) rather than det(A - x*I), which are the same polynomial up to a - /// sign for an odd-sized matrix. The eigenvalues are the roots either way, and the - /// monic form is the one sympy and most linear algebra texts state, so it is the - /// one that composes with the polynomial solver without a stray -1 in front. - /// - /// - /// Thrown for a non-square matrix, which has no characteristic polynomial - /// - public Entity CharacteristicPolynomial(Variable variable) - { - if (!IsSquare) - throw new InvalidMatrixOperationException( - "Only a square matrix has a characteristic polynomial"); - var shifted = InnerMatrix.Copy(true); - for (var i = 0; i < RowCount; i++) - shifted[i, i] = (variable - shifted[i, i]).InnerSimplified; - for (var row = 0; row < RowCount; row++) - for (var column = 0; column < ColumnCount; column++) - if (row != column) - shifted[row, column] = (-shifted[row, column]).InnerSimplified; - return new Matrix(shifted).Determinant - ?? throw new AngouriBugException("A square matrix always has a determinant"); - } - /// /// The Add operator. Performs an active operation diff --git a/Sources/Tests/UnitTests/Core/CharacteristicPolynomialTest.cs b/Sources/Tests/UnitTests/Core/CharacteristicPolynomialTest.cs deleted file mode 100644 index 1d6e79f96..000000000 --- a/Sources/Tests/UnitTests/Core/CharacteristicPolynomialTest.cs +++ /dev/null @@ -1,72 +0,0 @@ -// -// Copyright (c) 2019-2022 Angouri. -// AngouriMath is licensed under MIT. -// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. -// Website: https://am.angouri.org. -// - -using AngouriMath; -using AngouriMath.Core.Exceptions; -using AngouriMath.Extensions; -using Xunit; - -namespace AngouriMath.Tests.Core -{ - /// - /// https://github.com/asc-community/AngouriMath/issues/381 - /// - public sealed class CharacteristicPolynomialTest - { - private static Entity Expanded(Entity.Matrix matrix, string variable) - => MathS.CharacteristicPolynomial(matrix, variable).Expand().Simplify(); - - // det(x*I - A) written out by hand for each of these. - [Theory] - [InlineData("[[1, 2], [3, 4]]", "x ^ 2 - 5 * x - 2")] - [InlineData("[[2, 0], [0, 3]]", "x ^ 2 - 5 * x + 6")] - [InlineData("[[0, 1], [-1, 0]]", "x ^ 2 + 1")] - [InlineData("[[1, 0, 0], [0, 1, 0], [0, 0, 1]]", "x ^ 3 - 3 * x ^ 2 + 3 * x - 1")] - [InlineData("[[2, 1], [1, 2]]", "x ^ 2 - 4 * x + 3")] - public void MatchesTheHandComputedPolynomial(string matrix, string expected) => - Assert.Equal(Entity.Number.Integer.Create(0), - (Expanded((Entity.Matrix)matrix.ToEntity(), "x") - expected.ToEntity()).Expand().Simplify()); - - // The point of the polynomial: its roots are the eigenvalues. A diagonal matrix wears - // them on its face, so it is the case that can be checked without solving anything. - [Theory] - [InlineData("[[5, 0], [0, 7]]", 5, 7)] - [InlineData("[[-1, 0], [0, 4]]", -1, 4)] - public void DiagonalEntriesAreRoots(string matrix, int first, int second) - { - var polynomial = MathS.CharacteristicPolynomial((Entity.Matrix)matrix.ToEntity(), "x"); - Assert.Equal(Entity.Number.Integer.Create(0), polynomial.Substitute("x", first).Simplify()); - Assert.Equal(Entity.Number.Integer.Create(0), polynomial.Substitute("x", second).Simplify()); - } - - // Monic, which is what distinguishes det(x*I - A) from det(A - x*I): for an odd-sized - // matrix the other convention would put a -1 in front of the leading term. - [Theory] - [InlineData("[[1, 2], [3, 4]]", 2)] - [InlineData("[[1, 2, 3], [4, 5, 6], [7, 8, 10]]", 3)] - public void IsMonic(string matrix, int degree) - { - var polynomial = Expanded((Entity.Matrix)matrix.ToEntity(), "x"); - // The leading coefficient is what survives when the polynomial is divided by x^n - // and x is sent to infinity. - Assert.Equal(Entity.Number.Integer.Create(1), - (polynomial / MathS.Pow("x", degree)).Limit("x", Entity.Number.Real.PositiveInfinity).Simplify()); - } - - // Symbolic entries are as good as numeric ones: this is det(x*I - A) either way. - [Fact] - public void WorksOnSymbolicEntries() => - Assert.Equal(Entity.Number.Integer.Create(0), - (Expanded((Entity.Matrix)"[[a, b], [c, d]]".ToEntity(), "x") - - "x ^ 2 - (a + d) * x + a * d - b * c".ToEntity()).Expand().Simplify()); - - [Fact] - public void NonSquareMatrixIsRefused() => - Assert.Throws( - () => MathS.CharacteristicPolynomial((Entity.Matrix)"[[1, 2, 3], [4, 5, 6]]".ToEntity(), "x")); - } -}