Answer a polynomial with its roots, not with one per starting point

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

@@ -79,12 +79,67 @@ NumericsComplex ChooseGood() =>
                         MathS.Settings.PrecisionErrorCommon.Value)
                         res.Add(root);
                 }
-            // Rounded first, then verified: rounding collapses the duplicates the grid
-            // search produces, so there is less to verify, and what gets verified is the
-            // values actually handed back.
-            var rounded = WithoutIterationNoise(res);
-            rounded.RemoveWhere(root => !Satisfies(expr, v, root));
-            return rounded;
+            // Rounded and collapsed first, then verified: both cut down the duplicates the
+            // grid search produces, so there is less to verify, and what gets verified is
+            // the values actually handed back.
+            var distinct = OnePerRoot(WithoutIterationNoise(res), f);
+            distinct.RemoveWhere(root => !Satisfies(expr, v, root));
+            return distinct;
+        }
+
+        /// <summary>
+        /// How far apart two candidates may be, relative to their own size, and still be
+        /// the same root.
+        /// </summary>
+        /// <remarks>
+        /// The search starts from a grid and iterates in double precision, so the same root
+        /// reached from different starting points comes back agreeing to about sixteen
+        /// significant digits and differing after that. Left as they were, those were
+        /// counted as separate answers: x^5 + 3x + 1 came back with 28 roots and x^6 + x + 1
+        /// with 23, four of the 28 being -0.83907243306660750, -0.83907243306660761,
+        /// -0.83907243306660773 and -0.83907243306660784.
+        ///
+        /// The measured room is wide. Candidates for one root lie within 1e-15 of each other
+        /// relative to their size, while the nearest two distinct roots over the twelve
+        /// polynomials tried are 0.42 apart. This is 1e-13: a hundred times above the noise
+        /// and still twelve orders below the closest two roots that were ever told apart.
+        /// </remarks>
+        [ConstantField] private static readonly EDecimal SameRootTolerance = EDecimal.Create(1, -13);
+
+        private static bool AreTheSameRoot(Complex a, Complex b)
+        {
+            var apart = EDecimal.Max(
+                a.RealPart.EDecimal.Subtract(b.RealPart.EDecimal).Abs(),
+                a.ImaginaryPart.EDecimal.Subtract(b.ImaginaryPart.EDecimal).Abs());
+            var size = EDecimal.Max(EDecimal.One, EDecimal.Max(
+                a.RealPart.EDecimal.Abs(), a.ImaginaryPart.EDecimal.Abs()));
+            return !apart.GreaterThan(size.Multiply(SameRootTolerance));
+        }
+
+        /// <summary>
+        /// One candidate per root, keeping whichever of each group leaves the equation
+        /// closest to zero. Ordered first, so that which candidate is compared against
+        /// which does not depend on the order the grid happened to produce them in.
+        /// </summary>
+        private static HashSet<Complex> OnePerRoot(HashSet<Complex> roots, FastExpression f)
+        {
+            EDecimal Residual(Complex root)
+            {
+                try { return f.Call(root.ToNumerics()).ToNumber().Abs().EDecimal; }
+                catch (System.Exception) { return EDecimal.PositiveInfinity; }
+            }
+            var kept = new List<Complex>();
+            foreach (var root in roots
+                .OrderBy(root => root.RealPart.EDecimal)
+                .ThenBy(root => root.ImaginaryPart.EDecimal))
+            {
+                var same = kept.FindIndex(other => AreTheSameRoot(other, root));
+                if (same < 0)
+                    kept.Add(root);
+                else if (Residual(root).CompareTo(Residual(kept[same])) < 0)
+                    kept[same] = root;
+            }
+            return new HashSet<Complex>(kept);
         }
 
         /// <summary>

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

@@ -0,0 +1,92 @@
+//
+// 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;
+using System.Linq;
+using AngouriMath;
+using AngouriMath.Extensions;
+using Xunit;
+
+namespace AngouriMath.Tests.Algebra.SolveTest
+{
+    /// <summary>
+    /// A polynomial whose roots are found numerically must be answered with its roots and
+    /// not with one entry per starting point the search happened to use. No issue is filed
+    /// for this; it was found while checking the answers of another fix.
+    /// </summary>
+    public sealed class DuplicateRootTest
+    {
+        /// <summary>
+        /// The search starts from a grid, so the same root reached from different starting
+        /// points came back agreeing to sixteen significant digits and differing after
+        /// that, and each of those counted as an answer of its own: x^5 + 3x + 1 was
+        /// answered with 28 roots, four of them -0.83907243306660750, -0.83907243306660761,
+        /// -0.83907243306660773 and -0.83907243306660784.
+        /// </summary>
+        [Theory]
+        [InlineData("x ^ 5 + 3 * x + 1", 5)]
+        [InlineData("x ^ 6 + x + 1", 6)]
+        [InlineData("x ^ 5 + x + 1", 5)]
+        [InlineData("x ^ 5 - 5 * x + 2", 5)]
+        [InlineData("x ^ 6 - 3 * x ^ 2 + 1", 6)]
+        [InlineData("x ^ 9 - x - 1", 9)]
+        public void APolynomialHasAsManyRootsAsItsDegree(string expression, int degree) =>
+            Assert.Equal(degree, RootsOf(expression).Count);
+
+        // The polynomials answered exactly, by radicals rather than by iteration, are not
+        // touched by any of this and must keep the roots they had.
+        [Theory]
+        [InlineData("x ^ 2 - 4", 2)]
+        [InlineData("x ^ 3 - 1", 3)]
+        [InlineData("x ^ 5 - 1", 5)]
+        [InlineData("x ^ 8 - 1", 8)]
+        [InlineData("x ^ 7 - 2", 7)]
+        [InlineData("x ^ 4 - 10 * x ^ 3 + 35 * x ^ 2 - 50 * x + 24", 4)]
+        public void TheOnesSolvedExactlyAreUnaffected(string expression, int degree) =>
+            Assert.Equal(degree, RootsOf(expression).Count);
+
+        // Whichever candidate is kept from a group has to be one, so collapsing them may
+        // not hand back something that fails the equation.
+        [Theory]
+        [InlineData("x ^ 5 + 3 * x + 1")]
+        [InlineData("x ^ 6 + x + 1")]
+        [InlineData("x ^ 9 - x - 1")]
+        public void EveryRootKeptIsStillARoot(string expression)
+        {
+            var expr = expression.ToEntity();
+            foreach (var root in RootsOf(expression))
+            {
+                var residual = expr.Substitute("x", root).EvalNumerical();
+                Assert.True(
+                    residual.Abs().EDecimal.ToDouble() < 1e-8,
+                    $"{root.Stringize()} leaves {residual.Stringize()}");
+            }
+        }
+
+        // No two of the answers may be the same root written twice.
+        [Theory]
+        [InlineData("x ^ 5 + 3 * x + 1")]
+        [InlineData("x ^ 6 + x + 1")]
+        [InlineData("x ^ 9 - x - 1")]
+        [InlineData("x ^ 5 - 5 * x + 2")]
+        public void NoTwoAnswersAreTheSameNumber(string expression)
+        {
+            var values = RootsOf(expression)
+                .Select(root => (Entity.Number.Complex)root.EvalNumerical())
+                .Select(value => (Re: value.RealPart.EDecimal.ToDouble(), Im: value.ImaginaryPart.EDecimal.ToDouble()))
+                .ToList();
+            for (var i = 0; i < values.Count; i++)
+                for (var j = i + 1; j < values.Count; j++)
+                    Assert.True(
+                        Math.Max(Math.Abs(values[i].Re - values[j].Re), Math.Abs(values[i].Im - values[j].Im)) > 1e-9,
+                        $"{values[i]} and {values[j]} are the same root of {expression}");
+        }
+
+        private static Entity.Set.FiniteSet RootsOf(string expression) =>
+            Assert.IsType<Entity.Set.FiniteSet>(expression.ToEntity().SolveEquation("x"));
+    }
+}

Sources/Tests/UnitTests/Convenience/FromStringTest.cs

@@ -71,7 +71,13 @@ public void TestFormula8()
         {
             // Only needed for Mac
             using var _ = Settings.PrecisionErrorZeroRange.Set(2e-16m);
-            Assert.Equal(i, FromString("x^2+1").SolveNt(x).First());
+            var roots = FromString("x^2+1").SolveNt(x);
+            // Both roots, rather than whichever of them a set happens to hand back first.
+            // That order was never the solver's to decide and it changed once the search
+            // began collapsing its duplicates in a fixed order.
+            Assert.Equal(2, roots.Count);
+            Assert.Contains(Entity.Number.Complex.Create(0, 1), roots);
+            Assert.Contains(Entity.Number.Complex.Create(0, -1), roots);
         }
         [Fact] public void TestFormula9() => Assert.Equal(1, FromString("cos(sin(0))").EvalNumerical());
         [Fact] public void TestFormula10() => Assert.Equal(Entity.Number.Complex.Create(4, 1), FromString("2i + 2 * 2 - 1i").EvalNumerical());