From bce902f53f5f4e41f93d7244596e9fa4b679c8aa Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 27 Oct 2026 15:08:02 +0000 Subject: [PATCH 1/2] A quadratic whose discriminant is a square is answered without the root SolveQuadratic wrote sqrt(b^2 - 4ac) whatever the discriminant was. With polynomial coefficients and a square discriminant that left the square under the root: a z^2 - (a^2 - bc) z - abc has discriminant (a^2 + bc)^2 and roots a and -bc/a, and came back as (a^2 - bc -/+ sqrt((a^2 - bc)^2 + 4 a^2 bc)) / (2a). So did the most ordinary one, z^2 - (p + q) z + pq, whose roots are p and q. Either square root of the discriminant gives the same pair of roots, so any polynomial r with r^2 = b^2 - 4ac serves, with no absolute value. Where a, b and c parse as polynomials and the discriminant is a rational number times the square of one (MultivariatePolynomial.TrySquareRoot, the exact square root the integrator already uses), each numerator -b -/+ r is reduced by its greatest common divisor with 2a. Where the number's root is not rational, it is the one factor left outside: a +/- sqrt(2), (p + q)/2 +/- i (p - q)/2. Anything else takes the formula as before. On 2.6.0 and on this branch, Solve: a*z^2 - (a^2 - b*c)*z - a*b*c = 0 sqrt form -> { -b * c / a, a } x^2 - (p + q)*x + p*q = 0 sqrt form -> { q, p } x^2 - 2*a*x + a^2 - 2 = 0 sqrt form -> { a - sqrt(2), a + sqrt(2) } Fourteen such equations had every root substituted back at sampled values of the other symbols, with residual 0. Unit suite 17,138 passed, 0 failed; RootCheck 596/596 clean; CasBench 120/121, its one unsolved row unchanged. Fixes https://github.com/asc-community/AngouriMath/issues/2140 Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 15 +++ .../EquationSolver/PolynomialSolver.cs | 103 ++++++++++++++++++ .../Algebra/SquareDiscriminantTest.cs | 74 +++++++++++++ 3 files changed, 192 insertions(+) create mode 100644 Sources/Tests/UnitTests/Algebra/SquareDiscriminantTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 1a1ab7b56..d6b1b44a0 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -57,6 +57,21 @@ one inside another was answered with the empty set, a claim that it has no solut | `"1/(2 + p^(-3)) = 1/3".Solve("p")` | `{ }` | `{ 1, -1/2 + i * 1/2 * sqrt(3), -1/2 + i * -1/2 * sqrt(3) }` | | `"1/(1 + x^(-2)) = y".Solve("x")` | `{ }` | `{ ((-1) / (-y) - 1) ^ (-1/2), -((-1) / (-y) - 1) ^ (-1/2) }` | +### A quadratic whose discriminant is a square is answered without the root + +**Silent: the same roots, written as what they are.** The quadratic formula wrote `sqrt(b^2 - 4ac)` whatever the +discriminant was, so a symbolic quadratic with a square discriminant came back with the square under the root. +Either square root gives the same pair of roots, so where the coefficients are polynomials and the discriminant is +a rational number times the square of one, the roots are written with that polynomial and reduced over `2a` +([#2140](https://github.com/asc-community/AngouriMath/issues/2140)). + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"a*z^2 - (a^2 - b*c)*z - a*b*c = 0".Solve("z")` | `{ (a ^ 2 + -b * c - sqrt((-(a ^ 2 + -b * c)) ^ 2 - 4 * a * -a * b * c)) / (2 * a), … }` | `{ -b * c / a, a }` | +| `"x^2 - (p + q)*x + p*q = 0".Solve("x")` | `{ (p + q - sqrt((-(p + q)) ^ 2 - 4 * p * q)) / 2, … }` | `{ q, p }` | +| `"x^2 - 2*a*x + a^2 - 2 = 0".Solve("x")` | `{ (-(-2) * a - sqrt(((-2) * a) ^ 2 - 4 * (a ^ 2 + -2))) / 2, … }` | `{ a - sqrt(2), a + sqrt(2) }` | +| `"2*x^2 - 2*(p + q)*x + p^2 + q^2 = 0".Solve("x")` | `{ (-(-2) * (p + q) - sqrt(((-2) * (p + q)) ^ 2 - 8 * (p ^ 2 + q ^ 2))) / 4, … }` | `{ 1/2 * p + 1/2 * q - i * (1/2 * p - 1/2 * q), … }` | + ### A system with a product equation is solved branch by branch **Silent: solutions that were missing are answered.** `EquationSystem.Solve` eliminated `x * y = 0` for `y` as `0`, diff --git a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs index c2150776f..d569c83f8 100644 --- a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs @@ -69,11 +69,114 @@ internal static IEnumerable SolveQuadratic(Entity a, Entity b, Entity c) { if (TreeAnalyzer.IsZero(c)) return SolveLinear(a, b).Append(0); if (TreeAnalyzer.IsZero(a)) return SolveLinear(b, c); + if (RootsOverASquareDiscriminant(a, b, c) is { } exact) + return exact; var D = MathS.Sqr(b) - 4 * a * c; return new[] { ((-b - MathS.Sqrt(D)) / (2 * a)).InnerSimplified, ((-b + MathS.Sqrt(D)) / (2 * a)).InnerSimplified }; } + /// + /// The two roots, written without a radical, where the coefficients are polynomials and + /// the discriminant is a rational number times the square of one; otherwise + /// . + /// + /// + /// + /// Either square root of the discriminant will do. The roots are (-b - r)/(2a) and + /// (-b + r)/(2a), the same pair for r and -r, so any r with + /// r^2 = b^2 - 4ac gives the same set, and a polynomial one needs no absolute value. + /// a z^2 - (a^2 - bc) z - abc has discriminant (a^2 + bc)^2 and roots + /// a and -bc/a; through sqrt they were written with the square under + /// it, and z^2 - (p + q) z + pq with sqrt((p + q)^2 - 4pq) for p - q. + /// #2140 + /// + /// + /// The numerators are polynomials then too, so each root is reduced by the greatest common + /// divisor of its numerator and 2a rather than left for a later pass to cancel. + /// The discriminant's rational factor comes out first, and where its square root is not + /// rational -- 2 (p - q)^2, or -(p - q)^2 -- the roots are written with it, + /// sqrt(2) (p - q) and i (p - q), and left as quotients. A discriminant that + /// is a number is left to the formula, whose evaluation settles it already. + /// + /// + private static Entity[]? RootsOverASquareDiscriminant(Entity a, Entity b, Entity c) + { + var names = new SortedSet(Comparer.Create(static (x, y) => string.CompareOrdinal(x.Name, y.Name))); + foreach (var coefficient in new[] { a, b, c }) + foreach (var variable in coefficient.Vars) + names.Add(variable); + if (names.Count == 0 || names.Count > MultivariatePolynomial.MaxVariables) + return null; + var variables = names.ToList(); + var indices = new Dictionary(variables.Count); + for (var i = 0; i < variables.Count; i++) + indices[variables[i]] = i; + if (MultivariatePolynomial.TryParse(a, indices) is not { } pa + || MultivariatePolynomial.TryParse(b, indices) is not { } pb + || MultivariatePolynomial.TryParse(c, indices) is not { } pc + || pb.Multiply(pb) is not { } bSquared || pa.Multiply(pc) is not { } ac) + return null; + var discriminant = bSquared.Subtract(ac.ScaleBy(ERational.FromInt32(4))); + if (discriminant.Normalized(out var scale).TrySquareRoot() is not { } root) + return null; + var denominator = pa.ScaleBy(ERational.FromInt32(2)); + var order = Enumerable.Range(0, variables.Count).ToArray(); + // A quotient over 2a with what it shares with 2a cancelled, rational factors included. + Entity Reduced(MultivariatePolynomial numerator) + { + if (numerator.IsZero) + return Integer.Zero; + var below = denominator; + if (PolynomialGcd.Gcd(numerator, denominator, order, 0) is { } common + && numerator.DivideExact(common) is { } reducedAbove + && denominator.DivideExact(common) is { } reducedBelow) + (numerator, below) = (reducedAbove, reducedBelow); + below = below.Normalized(out var unit); + numerator = numerator.ScaleBy(unit); + return below.IsConstant + ? numerator.ScaleBy(ERational.One.Divide(below.Terms.Single().Value)).ToEntity(variables) + : numerator.ToEntity(variables) / below.ToEntity(variables); + } + var minusB = MultivariatePolynomial.Zero(variables.Count).Subtract(pb); + var overScale = ERational.Create(scale.Denominator, scale.Numerator); + if (MultivariatePolynomial.Constant(variables.Count, overScale).TrySquareRoot() is { } rational) + { + var r = root.ScaleBy(rational.IsZero ? ERational.Zero : rational.Terms.Single().Value); + return new[] { Reduced(minusB.Subtract(r)), Reduced(minusB.Add(r)) }; + } + // sqrt(2) (p - q), i (p - q): sqrt(n/d) is (s/d) sqrt(k) for n d = s^2 k, the s/d goes + // into the polynomial, and sqrt(k) -- or i sqrt(k) for a negative scale -- is the one + // factor left outside it. + var (square, rest) = SquarePart(overScale.Numerator.Multiply(overScale.Denominator).Abs()); + var outside = MathS.Sqrt(overScale.Sign < 0 ? -(Entity)Integer.Create(rest) : Integer.Create(rest)).InnerSimplified; + var middle = Reduced(minusB); + var spread = (outside * Reduced(root.ScaleBy(ERational.Create(square, overScale.Denominator)))).InnerSimplified; + return new[] { (middle - spread).InnerSimplified, (middle + spread).InnerSimplified }; + } + + /// + /// as s^2 k, with the square taken out by trial division + /// by the primes up to 1000 and by whatever of k is left a square; a larger square + /// factor stays in k, which leaves sqrt(k) correct and only less reduced. + /// + private static (EInteger Root, EInteger Remaining) SquarePart(EInteger value) + { + var root = EInteger.One; + var rest = value; + for (var prime = 2; prime <= 1000 && EInteger.FromInt32(prime * prime).CompareTo(rest) <= 0; prime++) + { + var square = EInteger.FromInt32(prime * prime); + while (rest.Remainder(square).IsZero) + { + rest = rest.Divide(square); + root = root.Multiply(prime); + } + } + var whole = rest.Sqrt(); + return whole.Multiply(whole).Equals(rest) ? (root.Multiply(whole), EInteger.One) : (root, rest); + } + /// Solves ax^3 + bx^2 + cx + d /// Coefficient of x^3 /// Coefficient of x^2 diff --git a/Sources/Tests/UnitTests/Algebra/SquareDiscriminantTest.cs b/Sources/Tests/UnitTests/Algebra/SquareDiscriminantTest.cs new file mode 100644 index 000000000..01c190100 --- /dev/null +++ b/Sources/Tests/UnitTests/Algebra/SquareDiscriminantTest.cs @@ -0,0 +1,74 @@ +// +// 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 System.Linq; +using AngouriMath; +using AngouriMath.Extensions; +using Xunit; +using static AngouriMath.Entity; +using static AngouriMath.Entity.Number; +using static AngouriMath.Entity.Set; + +namespace AngouriMath.Tests.Algebra +{ + /// + /// A quadratic whose discriminant is the square of a polynomial was answered with the square + /// left under the root: a z^2 - (a^2 - bc) z - abc has discriminant + /// (a^2 + bc)^2 and roots a and -bc/a, and came back as + /// (a^2 - bc -/+ sqrt((a^2 - bc)^2 + 4a^2 bc))/(2a). + /// #2140 + /// + [Trait("Area", "Algebra")] + public sealed class SquareDiscriminantTest + { + private static Entity[] Roots(string equation, string variable) + => Assert.IsType(equation.Solve(variable)).ToArray(); + + // -bc/a less itself is 0 provided a is not 0, and the condition is not what is compared. + private static bool Equal(Entity left, Entity right) + { + var difference = (left - right).Simplify(); + while (difference is Providedf(var inner, _)) + difference = inner; + return difference == Integer.Zero; + } + + [Theory] + [InlineData("a*z^2 - (a^2 - b*c)*z - a*b*c = 0", "z", "a", "-b*c/a")] + [InlineData("x^2 - (p + q)*x + p*q = 0", "x", "p", "q")] + [InlineData("-x^2 + (p - q)*x + p*q = 0", "x", "p", "-q")] + [InlineData("x^2 + x*y - 2*y^2 = 0", "x", "y", "-2*y")] + [InlineData("x^2 - 2*a*x + a^2 - 1 = 0", "x", "a - 1", "a + 1")] + [InlineData("x^2 - 2*a*x + a^2 - b^2/9 = 0", "x", "a - b/3", "a + b/3")] + public void TheRootsAreThePolynomialsTheyAre(string equation, string variable, string first, string second) + { + var roots = Roots(equation, variable); + Assert.Equal(2, roots.Length); + Assert.Contains(roots, root => Equal(root, first.ToEntity())); + Assert.Contains(roots, root => Equal(root, second.ToEntity())); + Assert.DoesNotContain(roots, root => root.Nodes.Any(node => node is Powf(_, Rational { Denominator: var d }) && d == 2)); + } + + /// + /// Where the discriminant is a rational number times a square that is not a rational + /// square, that number's root is the only radical left, and it is over a number. + /// + [Theory] + [InlineData("x^2 - 2*a*x + a^2 - 2 = 0", "a - sqrt(2)", "a + sqrt(2)")] + [InlineData("x^2 - 2*a*x + a^2 - 2*b^2 = 0", "a - sqrt(2) * b", "a + sqrt(2) * b")] + [InlineData("x^2 - 2*a*x + a^2 + 4 = 0", "a - 2i", "a + 2i")] + [InlineData("2*x^2 - 2*(p + q)*x + p^2 + q^2 = 0", "(p + q)/2 - i*(p - q)/2", "(p + q)/2 + i*(p - q)/2")] + public void OnlyTheRootOfTheNumberIsLeft(string equation, string first, string second) + { + var roots = Roots(equation, "x"); + Assert.Equal(2, roots.Length); + Assert.Contains(roots, root => Equal(root, first.ToEntity())); + Assert.Contains(roots, root => Equal(root, second.ToEntity())); + Assert.DoesNotContain(roots, root => root.Nodes.Any(node => node is Powf(var radicand, Rational { Denominator: var d }) && d == 2 && radicand is not Number)); + } + } +} From 11b0b538fca820890a73e57e7bbd9c6a046adba3 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 27 Oct 2026 15:44:54 +0000 Subject: [PATCH 2/2] A discriminant of odd degree in a symbol is refused before any polynomial is built The performance gate failed on #2170: SolveEasyMedium, x^2 + x + a, allocated 69,690 B against a baseline of 66,661 (+4.5%, band 3%). Its discriminant 1 - 4a is of degree one in a and cannot be a square, but the three coefficients were parsed, multiplied and normalized before TrySquareRoot said so. A square has an even degree in every variable, and where 2 deg b and deg a + deg c differ the larger is the discriminant's degree, so a quadratic for which that is odd is refused first. Degrees as written bound the true ones; a coefficient that cancels inside itself can be refused when it need not be, which leaves it to the formula as before. On one machine, B per SolveEquation: master this x2 + x + a 66,976 67,064 x2 + b*x + a 70,488 70,616 x2 + x + 1 450,114 450,282 The fourteen equations still check by substitution, and the unit suite is 17,138 passed, 0 failed. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../EquationSolver/PolynomialSolver.cs | 42 +++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs index d569c83f8..021b13b52 100644 --- a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs @@ -102,6 +102,14 @@ internal static IEnumerable SolveQuadratic(Entity a, Entity b, Entity c) /// private static Entity[]? RootsOverASquareDiscriminant(Entity a, Entity b, Entity c) { + // A square has an even degree in every variable, and where 2 deg b and deg a + deg c + // differ, that larger one is the discriminant's degree. Most quadratics with a symbol + // in them are refused here, before any polynomial is built: x^2 + x + a, whose + // discriminant 1 - 4a is of degree one in a, allocated 4.6% more to be solved without + // it. Degrees as written bound the true ones, so a coefficient that cancels inside + // itself can be refused when it need not be, which only leaves it to the formula. + if (!HasEvenDegreesWherever(a, b, c)) + return null; var names = new SortedSet(Comparer.Create(static (x, y) => string.CompareOrdinal(x.Name, y.Name))); foreach (var coefficient in new[] { a, b, c }) foreach (var variable in coefficient.Vars) @@ -155,6 +163,40 @@ Entity Reduced(MultivariatePolynomial numerator) return new[] { (middle - spread).InnerSimplified, (middle + spread).InnerSimplified }; } + /// + /// Whether, for every variable in , or + /// , the degree of b^2 - 4ac is not known to be odd. + /// + private static bool HasEvenDegreesWherever(Entity a, Entity b, Entity c) + { + static bool Even(Entity a, Entity b, Entity c, Variable v) + { + var (twiceB, ac) = (2 * DegreeAsWritten(b, v), DegreeAsWritten(a, v) + DegreeAsWritten(c, v)); + return twiceB == ac || System.Math.Max(twiceB, ac) % 2 == 0; + } + foreach (var coefficient in new[] { a, b, c }) + foreach (var variable in coefficient.Vars) + if (!Even(a, b, c, variable)) + return false; + return true; + } + + /// + /// The degree of in as it is written, which + /// bounds the degree it has: of a power its base's times the exponent, of a product the + /// sum, of a sum the most. + /// + private static int DegreeAsWritten(Entity expr, Variable x) => expr switch + { + _ when !expr.ContainsNode(x) => 0, + Variable => 1, + Powf(var @base, Integer power) when power.EInteger.CanFitInInt32() => DegreeAsWritten(@base, x) * power.EInteger.ToInt32Checked(), + Mulf(var left, var right) => DegreeAsWritten(left, x) + DegreeAsWritten(right, x), + Sumf(var left, var right) => System.Math.Max(DegreeAsWritten(left, x), DegreeAsWritten(right, x)), + Minusf(var left, var right) => System.Math.Max(DegreeAsWritten(left, x), DegreeAsWritten(right, x)), + _ => 0, + }; + /// /// as s^2 k, with the square taken out by trial division /// by the primes up to 1000 and by whatever of k is left a square; a larger square