Solve a polynomial through its factors rather than by its degree (#272)

Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/AnalyticalEquationSolver.cs

@@ -140,10 +140,53 @@ static Entity TryDowncast(Entity equation, Variable x, Entity root)
                     return Solve(expression, x, compensateSolving).Filter(predicate, x);
                 case Piecewise p:
                     return EquationSolver.SolvePiecewise(p, x, (e, x) => Solve(e, x, compensateSolving));
+                // A product is zero exactly where one of its factors is. Without this the
+                // factors are expanded back into a polynomial and the equation is answered
+                // by the general formula for whatever degree that turns out to be, so
+                // (x - 1)(x^2 - 3) = 0 -- handed over already factored -- came back as two
+                // nested cube roots of 26 + 18i rather than as 1 and +-sqrt(3).
+                // https://github.com/asc-community/AngouriMath/issues/272
+                //
+                // A factor free of x contributes no roots. It could still be zero itself,
+                // which would make every x a solution, but that is a reading the solver
+                // does not take anywhere else either: a * (x^2 - 3) = 0 already answered
+                // +-sqrt(3) by dividing a out, so nothing changes here.
+                case Mulf(var left, var right) when left.ContainsNode(x) || right.ContainsNode(x):
+                {
+                    var factors = new List<Entity>(2);
+                    if (left.ContainsNode(x)) factors.Add(left);
+                    if (right.ContainsNode(x)) factors.Add(right);
+                    var product = expr;
+                    // Zeroing one factor makes the product zero only where the others are
+                    // defined: x * (1/x) is undefined at 0, not zero there. This has to be
+                    // checked here rather than left to the verification every solver shares,
+                    // because that deliberately keeps a root whose residual does not
+                    // evaluate to a number at all -- a root is only dropped on positive
+                    // evidence, and NaN is not evidence (see SolveStatement.IsSpurious).
+                    bool ProductIsDefinedAt(Entity root)
+                        => root.Vars.Any()
+                           || product.Substitute(x, root).Evaled is not Complex value
+                           || value.IsFinite;
+                    var fromFactors = factors.Select(factor => Solve(factor, x, compensateSolving)).Unite();
+                    // Collapsed first: uniting sets builds a Unionf, and the roots cannot be
+                    // looked at one at a time until it is a finite set again.
+                    return fromFactors.InnerSimplified is FiniteSet found
+                        ? found.Where(ProductIsDefinedAt).ToSet()
+                        : fromFactors;
+                }
                 default:
                     break;
             }
 
+            // The same equation written out as a sum. Where a rational root can be divided
+            // out, what is left is a polynomial of lower degree, and the factors are then
+            // solved one at a time by the case above -- which is how x^3 - x^2 - 3x + 3 is
+            // answered exactly, and how a quintic is answered at all.
+            // https://github.com/asc-community/AngouriMath/issues/272
+            if (Functions.PolynomialFactoring.TrySplitOffRationalRoots(expr, x, out var split)
+                && split is Mulf)
+                return Solve(split, x, compensateSolving);
+
             if (PolynomialSolver.SolveAsPolynomial(expr, x, out var isIdentity) is { } poly)
                 return poly.Select(e => TryDowncast(expr, x, e.InnerSimplified)).ToSet();
             if (isIdentity)

Sources/AngouriMath/Functions/Simplification/PolynomialFactoring.cs

@@ -105,6 +105,62 @@ internal static bool TryFactor(Entity expr, Variable x, [NotNullWhen(true)] out
             return degree > 1;
         }
 
+        /// <summary>
+        /// <paramref name="expr"/> written as the linear factors of its rational roots
+        /// times whatever will not divide, or <see langword="false"/> where it has no
+        /// rational root at all.
+        /// </summary>
+        /// <remarks>
+        /// <para>
+        /// <see cref="TryFactor"/> answers a different question and answers it more
+        /// strictly: it offers a factorization only when the polynomial splits completely
+        /// into whole linear factors, because a partly factored sum is not obviously a
+        /// better way of *writing* the same thing. Solving is not a question of how to
+        /// write it. <c>(x - 1)(x^2 - 3)</c> is no nicer to read than
+        /// <c>x^3 - x^2 - 3x + 3</c>, but it is the difference between answering with
+        /// <c>1</c> and <c>+-sqrt(3)</c> and answering with two nested cube roots of
+        /// <c>26 + 18i</c> -- and above degree four, between answering at all and handing
+        /// the equation to the numeric solver.
+        /// https://github.com/asc-community/AngouriMath/issues/272
+        /// </para>
+        /// <para>
+        /// Fractional roots are kept here where <see cref="TryFactor"/> declines them, for
+        /// the same reason: a root of <c>1/2</c> is as exact an answer as any, whatever
+        /// the factored form looks like on the page.
+        /// </para>
+        /// </remarks>
+        internal static bool TrySplitOffRationalRoots(
+            Entity expr, Variable x, [NotNullWhen(true)] out Entity? factored)
+        {
+            factored = null;
+            if (!TryGetRationalCoefficients(expr, x, out var coefficients))
+                return false;
+
+            // Two terms is a x^n + b, which is answered whole by inverting x^n = -b/a, and
+            // that gives the n roots in the form a + bi. Splitting it would divide out the
+            // rational ones and leave the rest to be dug out of a quotient: x^3 - 8 reads
+            // as 2 and (-1/2 +- i*sqrt(3)/2)*2 that way, and as 2 and (-2 -+ sqrt(-12))/2
+            // this way. Nothing here improves on that, so it is left alone.
+            if (coefficients.Count(coefficient => !coefficient.IsZero) <= 2)
+                return false;
+
+            var factors = new List<Entity>();
+            foreach (var candidate in RootCandidates(coefficients))
+                while (coefficients.Length > 1 && Evaluate(coefficients, candidate).IsZero)
+                {
+                    factors.Add(LinearFactor(candidate, x));
+                    coefficients = DivideByLinear(coefficients, candidate);
+                }
+            if (factors.Count == 0)
+                return false;
+
+            var remainder = BuildPolynomial(coefficients, x);
+            if (remainder != Integer.One)
+                factors.Add(remainder);
+            factored = factors.Aggregate((left, right) => left * right);
+            return true;
+        }
+
         /// <summary>
         /// One step of a partial fraction decomposition, at a rational root of the
         /// denominator: <c>N/D</c> becomes <c>A/(x - r) + R/Q</c>, where <c>D = (x - r)Q</c>.

Sources/Tests/UnitTests/Algebra/SolveTest/FactoredPolynomialSolveTest.cs

@@ -0,0 +1,175 @@
+//
+// 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.Linq;
+using AngouriMath;
+using AngouriMath.Extensions;
+using Xunit;
+
+namespace AngouriMath.Tests.Algebra
+{
+    /// <summary>
+    /// A polynomial equation that has a factor the solver can answer exactly must be
+    /// answered through that factor, rather than by handing the whole polynomial to the
+    /// general formula for its degree -- or, above degree four, to the numeric solver.
+    /// https://github.com/asc-community/AngouriMath/issues/272
+    /// </summary>
+    public sealed class FactoredPolynomialSolveTest
+    {
+        private static Entity.Set.FiniteSet Roots(string equation) =>
+            (Entity.Set.FiniteSet)equation.ToEntity().Solve("x");
+
+        private static void AssertRootsAre(string equation, params string[] expected)
+        {
+            var roots = Roots(equation);
+            foreach (var root in expected)
+                Assert.Contains(root.ToEntity().InnerSimplified, roots);
+            Assert.Equal(expected.Length, roots.Count);
+        }
+
+        /// <summary>
+        /// The clearest statement of the defect, and the one that does not depend on how a
+        /// radical is spelled: an equation whose roots are all real must not be answered
+        /// with expressions built out of complex numbers. Cardano's formula reaches the
+        /// real roots of x^3 - x^2 - 3x + 3 through (26 + 18i)^(1/3).
+        /// </summary>
+        private static void AssertNoImaginaryPartsInTheWorking(string equation)
+        {
+            foreach (var root in Roots(equation))
+                foreach (var node in root.Nodes)
+                    if (node is Entity.Number.Complex and not Entity.Number.Real)
+                        Assert.True(false,
+                            $"{equation} has only real roots, but {root.Stringize()} reaches one " +
+                            $"through {node.Stringize()}");
+        }
+
+        /// <summary>A root that is a decimal is a root the solver gave up on.</summary>
+        private static void AssertEveryRootIsExact(string equation)
+        {
+            foreach (var root in Roots(equation))
+                foreach (var node in root.Nodes)
+                    if (node is Entity.Number.Real and not Entity.Number.Rational)
+                        Assert.True(false,
+                            $"{equation} was answered with the decimal {node.Stringize()} " +
+                            $"inside {root.Stringize()}");
+        }
+
+        // Handed the factors outright, the solver still expanded them into a cubic and ran
+        // Cardano over it, so the answer to an equation whose roots are 1 and +-sqrt(3) came
+        // back as two nested cube roots of 26 + 18i.
+        [Fact]
+        public void AProductIsSolvedThroughItsFactors() =>
+            AssertRootsAre("(x - 1) * (x ^ 2 - 3) = 0", "1", "sqrt(3)", "-sqrt(3)");
+
+        [Fact]
+        public void AProductOfSeveralFactorsIsSolvedThroughThemToo() =>
+            AssertRootsAre("(x - 1) * (x - 2) * (x + 2) * (x ^ 2 - 3) = 0",
+                "1", "2", "-2", "sqrt(3)", "-sqrt(3)");
+
+        // The same equation written out. A rational root is there to be found, and dividing
+        // it out leaves a quadratic that is answered exactly.
+        [Fact]
+        public void AnExpandedCubicWithARationalRootIsSplitAtIt() =>
+            AssertRootsAre("x ^ 3 - x ^ 2 - 3 * x + 3 = 0", "1", "sqrt(3)", "-sqrt(3)");
+
+        // Above degree four there is no general formula, so the whole equation used to fall
+        // through to the numeric solver -- which returned decimals, one of them carrying a
+        // spurious -2.3e-15i.
+        [Fact]
+        public void AQuinticThatSplitsIntoSolvableFactorsIsSolvedExactly() =>
+            AssertRootsAre("x ^ 5 - x ^ 4 - 7 * x ^ 3 + 7 * x ^ 2 + 12 * x - 12 = 0",
+                "1", "2", "-2", "sqrt(3)", "-sqrt(3)");
+
+        [Theory]
+        [InlineData("(x - 1) * (x ^ 2 - 3) = 0")]
+        [InlineData("x ^ 3 - x ^ 2 - 3 * x + 3 = 0")]
+        [InlineData("x ^ 3 - 6 * x ^ 2 + 11 * x - 6 = 0")]
+        public void RealRootsAreNotReachedThroughComplexNumbers(string equation) =>
+            AssertNoImaginaryPartsInTheWorking(equation);
+
+        [Theory]
+        [InlineData("x ^ 3 - x ^ 2 - 3 * x + 3 = 0")]
+        [InlineData("x ^ 5 - x ^ 4 - 7 * x ^ 3 + 7 * x ^ 2 + 12 * x - 12 = 0")]
+        [InlineData("(x - 1) * (x - 2) * (x + 2) * (x ^ 2 - 3) = 0")]
+        public void RootsThatCanBeWrittenExactlyAre(string equation) =>
+            AssertEveryRootIsExact(equation);
+
+        private static bool IsNear(Entity root, double re, double im)
+        {
+            if (root.Vars.Any())
+                return false;
+            var value = root.EvalNumerical();
+            return System.Math.Abs(value.RealPart.EDecimal.ToDouble() - re) < 1e-9
+                && System.Math.Abs(value.ImaginaryPart.EDecimal.ToDouble() - im) < 1e-9;
+        }
+
+        /// <summary>
+        /// Worse than an unreadable answer: an incomplete one. Above degree four the whole
+        /// equation went to the numeric solver, which searches for real roots, so
+        /// (x - 1)(x^2 - 2)(x^2 + x + 1) = 0 came back as { 1, sqrt(2), -sqrt(2) } -- the
+        /// two roots of its x^2 + x + 1 factor were missing, and a set that is short of two
+        /// roots does not look any different from one that is not.
+        /// Splitting off the rational root leaves a quartic, which is solved whole.
+        /// </summary>
+        [Theory]
+        [InlineData("x ^ 5 - 2 * x ^ 3 - x ^ 2 + 2 = 0")]
+        [InlineData("(x - 1) * (x ^ 2 - 2) * (x ^ 2 + x + 1) = 0")]
+        public void AQuinticKeepsTheRootsOfItsIrreducibleFactor(string equation)
+        {
+            var roots = Roots(equation);
+            var sqrt2 = System.Math.Sqrt(2);
+            var halfSqrt3 = System.Math.Sqrt(3) / 2;
+            foreach (var (re, im) in new[] { (1.0, 0.0), (sqrt2, 0.0), (-sqrt2, 0.0),
+                                             (-0.5, halfSqrt3), (-0.5, -halfSqrt3) })
+                Assert.Contains(roots, root => IsNear(root, re, im));
+            Assert.Equal(5, roots.Count);
+        }
+
+        // What must not change. A polynomial with no rational root has nothing to split off,
+        // and its roots are what they always were.
+        [Fact]
+        public void APolynomialWithoutARationalRootIsUnaffected() =>
+            AssertRootsAre("x ^ 2 - 2 = 0", "sqrt(2)", "-sqrt(2)");
+
+        // A two-term polynomial is answered whole by inverting x^n = c, which gives its
+        // roots as a + bi. Splitting off the rational ones instead leaves the others to be
+        // dug out of the quotient, and x^3 - 8 comes back as (-2 -+ sqrt(-12))/2 rather
+        // than as (-1/2 +- i*sqrt(3)/2)*2. Nothing is gained, so those are left alone.
+        [Theory]
+        [InlineData("x ^ 3 - 8 = 0")]
+        [InlineData("x ^ 6 - 1 = 0")]
+        public void ABinomialKeepsTheAnswerInversionGivesIt(string equation)
+        {
+            foreach (var root in Roots(equation))
+                foreach (var node in root.Nodes)
+                    if (node is Entity.Powf(Entity.Number.Real { IsNegative: true } radicand, Entity.Number.Rational))
+                        Assert.True(false,
+                            $"{equation} was answered through a root of the negative " +
+                            $"{radicand.Stringize()}, in {root.Stringize()}");
+        }
+
+        // Repeated factors are one root, not several.
+        [Fact]
+        public void ARepeatedFactorContributesItsRootOnce() =>
+            AssertRootsAre("(x - 1) ^ 2 * (x + 3) = 0", "1", "-3");
+
+        // A factor free of the variable has no roots to contribute, and must not make the
+        // equation look unsolvable either.
+        [Fact]
+        public void AConstantFactorContributesNothing() =>
+            AssertRootsAre("5 * (x ^ 2 - 3) = 0", "sqrt(3)", "-sqrt(3)");
+
+        // Solving factor by factor must not smuggle in a value that is not a root of the
+        // whole equation. Here x = 0 zeroes the first factor but makes the second undefined.
+        [Fact]
+        public void AFactorsRootThatTheWholeEquationRejectsIsDropped()
+        {
+            var roots = Roots("x * (1 / x) = 0");
+            Assert.DoesNotContain(Entity.Number.Integer.Create(0), roots);
+        }
+    }
+}