A radical of a constant is a coefficient, and a polynomial with real coefficients has all its roots

BREAKING-CHANGES.md

@@ -169,6 +169,27 @@ answered 5 numbers that are not roots. Now 2 roots are missing and none is wrong
 | `"sqrt(x) - 1/x - 1 = 0".Solve("x")` | `{ 2.1478… }` | that root in radicals |
 | `"(3 - x)*sqrt(x) - x = 0".Solve("x")` | `{ 0, (7 - sqrt(13)) / 2 }` | `{ 0, (-7 + sqrt(13)) / (-2) }` |
 
+### A radical of a constant is a coefficient
+
+**Not silent: an equation that ended the process with a stack overflow is answered, and a polynomial with a radical
+or another real constant among its coefficients has all its roots.** The solver raises an equation to a power to remove
+a radical, and it took a radical of a constant, such as `sqrt(2)`, for one. Raising `x^5 + x - sqrt(2) = 0` gives
+`(x^5 + x)^2 = 2`, and the substitution that solves that gives back `x^5 + x = sqrt(2)`, the equation it started from,
+so the solver recursed until the stack ran out. A radical of a constant is now a coefficient. No issue was filed; a
+check of systems built from known solutions found it, eliminating `y` from `sqrt(y) - x = sqrt(2) - 1, |x| + y = 3`.
+
+A polynomial with such a coefficient and no roots in closed form went to Newton's method, which starts inside a
+bounded region and misses a root outside it. It is now solved as one with rational coefficients is, by the
+Durand–Kerner iteration, which answers all its roots or none.
+
+| Input | Was (2.6.0) | Now |
+|---|---|---|
+| `"x^5 + x - sqrt(2) = 0".Solve("x")` | stack overflow | five roots |
+| `"x^5 + x^2 - 3 * sqrt(2) = 0".Solve("x")` | three roots of five | five roots |
+| `"x^5 - 30 * x^4 + sqrt(2) = 0".Solve("x")` | one root of five | five roots |
+| `"x^5 - 30 * x^4 + pi = 0".Solve("x")` | eight numbers, none of them the root near 30 | five roots |
+| `"2 * x^5 + 37 * x^4 - x^3 - 4 * x^2 - 2 * sqrt(2) * x + 4 = 0".Solve("x")` | four roots of five | five roots |
+
 ### A cube shifted is solved
 
 **Silent: a cubic `a(x + h)^3 = k` written out, with `b^2 = 3ac`, is solved where it was answered `{ }`.** Cardano's

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

@@ -411,7 +411,12 @@ bool ProductIsDefinedAt(Entity root)
             // answered the roots it found there as though they were all: x^5 - 30x^4 + 1 = 0 has a root
             // near 30, and was answered with one root of its five. It is numerical, so it is asked
             // only where Newton's method would have been: turning that off asks for analytical answers.
-            if (MathS.Settings.AllowNewton && expr.Vars.Count == 1 && SumOverSet.NumericRoots(expr, x) is { Count: > 0 } numeric)
+            // So does one whose coefficients are real constants, sqrt(2) or pi among them, which
+            // Newton's method left the same way: x^5 + 37x^4 - 4x^2 - 2 sqrt(2) x + 4 = 0 has a root
+            // near -18.5, and was answered with four roots of its five.
+            if (MathS.Settings.AllowNewton && expr.Vars.Count == 1
+                && (SumOverSet.NumericRoots(expr, x)
+                    ?? Algebra.NumericalSolving.DurandKerner.RootsOf(expr, x, MathS.Settings.DecimalPrecisionContext.Value)) is { Count: > 0 } numeric)
                 return numeric.Select(ent => TryDowncast(expr, x, ent)).ToSet();
 
             // An identity is satisfied by every x, and not by the points a numerical search

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

@@ -39,7 +39,11 @@ internal static bool TrySolve(Entity expr, Variable x, out Set dst)
                     // whose exponent is a sum, and read as it stands that was no fractional power, so
                     // sqrt(x) - 1/x - 1 - i = 0, cleared of its denominator, went to Newton's method and
                     // was answered without its root -1.
-                    if (!(mpChild is Powf(var @base, var exponent) && exponent.InnerSimplified is Number num and not Integer)) // (x + 3) ^ 3
+                    // And only a power of something in x is a radical to raise. sqrt(2) is a coefficient:
+                    // taken for one, -2 sqrt(2) + |x| + (2 sqrt(2) - 2) x + x^2 = 0 was squared into
+                    // (|x| + ...)^2 = 8, whose substitution gave back the equation it started from, and
+                    // the solver recursed until the stack ran out.
+                    if (!(mpChild is Powf(var @base, var exponent) && @base.ContainsNode(x) && exponent.InnerSimplified is Number num and not Integer)) // (x + 3) ^ 3
                     {
                         potentialFraction.multiplier *= mpChild;
                         continue;

Sources/AngouriMath/Functions/Continuous/Solvers/NumericalSolving/DurandKerner.cs

@@ -49,14 +49,79 @@ internal static class DurandKerner
         internal static IReadOnlyList<Complex>? Roots(IntegerPolynomial squareFree, EContext context)
         {
             var degree = squareFree.Degree;
+            if (degree < 1)
+                return null;
+            var coefficients = new EDecimal[degree + 1];
+            for (var i = 0; i <= degree; i++)
+                coefficients[i] = EDecimal.FromEInteger(squareFree[i]);
+            // Distinct to half the working precision, which a square-free polynomial's roots are
+            // unless the iteration has not separated two of them.
+            return Roots(coefficients, context, context.Precision.ToInt32Checked() / 2);
+        }
+
+        /// <summary>
+        /// The roots of <paramref name="polynomial"/> in <paramref name="w"/> at the precision of
+        /// <paramref name="context"/>, where it is a polynomial whose coefficients are real
+        /// constants and not all of them rational, <c>x^5 - sqrt(2) x + 1</c> or <c>x^5 + pi</c>;
+        /// <see langword="null"/> otherwise, and where they are not settled.
+        /// </summary>
+        /// <remarks>
+        /// The coefficients are known to the working precision and not exactly, so a repeated root
+        /// cannot be told from a cluster: perturbed by that error, a root of multiplicity <c>m</c>
+        /// splits into <c>m</c> roots about its <c>m</c>-th root apart, and <c>m</c> is at most
+        /// the degree <c>n</c>. Approximations closer than <c>10^(-digits / (2 n))</c> are taken
+        /// not to be separated then, where a polynomial known to be square-free is held to half
+        /// the digits. And a coefficient smaller than the largest by more than half the digits may
+        /// be a zero that was not simplified to one, which as the leading coefficient would put
+        /// roots beyond any bound, so it answers nothing then.
+        /// </remarks>
+        internal static IReadOnlyList<Complex>? RootsOf(Entity polynomial, Entity.Variable w, EContext context)
+        {
+            if (!TreeAnalyzer.TryGetPolynomial(polynomial, w, out var monomials))
+                return null;
+            var values = new Dictionary<int, EDecimal>();
+            var rational = true;
+            foreach (var monomial in monomials)
+            {
+                if (monomial.Key.Sign < 0 || monomial.Key.CompareTo(IntegerPolynomial.MaxDegree) > 0
+                    || monomial.Value.Vars.Count > 0
+                    || monomial.Value.EvalNumerical() is not Real { EDecimal: { IsFinite: true } value })
+                    return null;
+                rational &= monomial.Value.InnerSimplified is Rational;
+                if (!value.IsZero)
+                    values[monomial.Key.ToInt32Checked()] = value;
+            }
+            if (rational || values.Count == 0)
+                return null;
+            var degree = values.Keys.Max();
             if (degree < 1)
                 return null;
             var digits = context.Precision.ToInt32Checked();
+            var largest = values.Values.Aggregate(EDecimal.Zero, static (max, c) => EDecimal.Max(max, c.Abs()));
+            var smallest = largest.Multiply(EDecimal.Create(EInteger.One, EInteger.FromInt32(-(digits / 2))), context);
+            if (values.Values.Any(c => c.Abs().CompareTo(smallest) < 0))
+                return null;
+            var coefficients = new EDecimal[degree + 1];
+            for (var i = 0; i <= degree; i++)
+                coefficients[i] = values.TryGetValue(i, out var c) ? c : EDecimal.Zero;
+            return Roots(coefficients, context, digits / (2 * degree));
+        }
+
+        /// <summary>
+        /// The roots of the polynomial with real <paramref name="coefficients"/>, lowest power first,
+        /// at the precision of <paramref name="context"/>, each approximation at least
+        /// <c>10^-separationDigits</c> from the others; <see langword="null"/> where they are not
+        /// settled or not that far apart.
+        /// </summary>
+        private static IReadOnlyList<Complex>? Roots(EDecimal[] coefficients, EContext context, int separationDigits)
+        {
+            var degree = coefficients.Length - 1;
+            var digits = context.Precision.ToInt32Checked();
             var work = context.WithPrecision(digits + GuardDigits);
-            var lead = EDecimal.FromEInteger(squareFree[degree]);
+            var lead = coefficients[degree];
             var monic = new EDecimal[degree];
             for (var i = 0; i < degree; i++)
-                monic[i] = EDecimal.FromEInteger(squareFree[i]).Divide(lead, work);
+                monic[i] = coefficients[i].Divide(lead, work);
 
             // Cauchy's bound: every root is within 1 + max |a_i / a_n| of zero. The approximations
             // start spread round that circle, turned off the axes so that no two start as a
@@ -104,9 +169,7 @@ internal static class DurandKerner
             if (!settled)
                 return null;
 
-            // Distinct to half the working precision, which a square-free polynomial's roots are
-            // unless the iteration has not separated two of them.
-            var separation = EDecimal.Create(EInteger.One, EInteger.FromInt32(-(digits / 2)));
+            var separation = EDecimal.Create(EInteger.One, EInteger.FromInt32(-separationDigits));
             for (var i = 0; i < degree; i++)
                 for (var j = i + 1; j < degree; j++)
                     if (re[i].Subtract(re[j], work).Abs().Add(im[i].Subtract(im[j], work).Abs(), work).CompareTo(separation) <= 0)

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

@@ -51,5 +51,65 @@ public void EveryRootIsExactAndSatisfiesTheEquation(string equation, string resi
             if (among is { })
                 Assert.Contains(roots, root => (root - among.ToEntity()).Simplify() == 0);
         }
+
+        /// <summary>
+        /// A radical of a constant is a coefficient, not a radical to raise. Taken for one,
+        /// <c>-2 sqrt(2) + |x| + (2 sqrt(2) - 2) x + x^2 = 0</c> was squared into
+        /// <c>(|x| + ...)^2 = 8</c>, whose substitution gave back the equation it started from, and
+        /// the solver recursed until the stack ran out: the process ended. It is met eliminating
+        /// <c>y</c> from <c>sqrt(y) - x = sqrt(2) - 1, |x| + y = 3</c>, which a check of systems built
+        /// from known solutions found; no issue was filed.
+        /// </summary>
+        [Theory]
+        [InlineData("sqrt(3 - abs(x)) - x - sqrt(2) + 1 = 0")]
+        [InlineData("abs(x) + x^2 - 2 * sqrt(2) = 0")]
+        [InlineData("x^5 + x - sqrt(2) = 0")]
+        public void ARadicalOfAConstantIsNotRaised(string equation)
+        {
+            var solving = System.Threading.Tasks.Task.Run(() => equation.ToEntity().Solve("x"));
+            Assert.True(solving.Wait(System.TimeSpan.FromMinutes(1)), $"{equation} did not finish");
+        }
+
+        /// <summary>A polynomial with sqrt(2) among its coefficients met the same loop.</summary>
+        [Fact]
+        public void APolynomialWithARadicalCoefficientIsAnswered()
+        {
+            var roots = Assert.IsType<Set.FiniteSet>("x^5 + x - sqrt(2) = 0".ToEntity().Solve("x"));
+            Assert.Equal(5, roots.Count);
+            foreach (var root in roots)
+            {
+                var value = "x^5 + x - sqrt(2)".ToEntity().Substitute("x", root).EvalNumerical();
+                Assert.True(Complex.Abs(value).EDecimal.CompareTo(EDecimal.FromString("1E-10")) < 0, $"{root} leaves {value}");
+            }
+        }
+
+        /// <summary>
+        /// And one whose coefficients are real constants has all its roots, as one with rational
+        /// coefficients does. Newton's method answered it, from starts inside a bounded region, and
+        /// missed the root outside it: 30 for the first two and -18.5 for the third.
+        /// </summary>
+        [Theory]
+        [InlineData("x^5 - 30 * x^4 + sqrt(2)")]
+        [InlineData("x^5 - 30 * x^4 + pi")]
+        [InlineData("2 * x^5 + 37 * x^4 - x^3 - 4 * x^2 - 2 * sqrt(2) * x + 4")]
+        public void APolynomialWithRealCoefficientsHasAllItsRoots(string polynomial)
+        {
+            var solving = System.Threading.Tasks.Task.Run(() => $"{polynomial} = 0".ToEntity().Solve("x"));
+            Assert.True(solving.Wait(System.TimeSpan.FromMinutes(1)), $"{polynomial} = 0 did not finish");
+            var roots = Assert.IsType<Set.FiniteSet>(solving.Result);
+            Assert.Equal(5, roots.Count);
+            foreach (var root in roots)
+            {
+                var value = polynomial.ToEntity().Substitute("x", root).EvalNumerical();
+                Assert.True(Complex.Abs(value).EDecimal.CompareTo(EDecimal.FromString("1E-10")) < 0, $"{root} leaves {value}");
+            }
+        }
+
+        [Fact]
+        public void ASystemThatEliminatesToOneFinishes()
+        {
+            var solving = System.Threading.Tasks.Task.Run(() => MathS.Equations("sqrt(y) - x - sqrt(2) + 1", "abs(x) + y - 3").Solve("x", "y"));
+            Assert.True(solving.Wait(System.TimeSpan.FromMinutes(1)), "the system did not finish");
+        }
     }
 }