diff --git a/Sources/AngouriMath/Functions/Continuous/Solvers/NumericalSolving/NewtonSolver.cs b/Sources/AngouriMath/Functions/Continuous/Solvers/NumericalSolving/NewtonSolver.cs index 26cf2bbf2..f83042bfe 100644 --- a/Sources/AngouriMath/Functions/Continuous/Solvers/NumericalSolving/NewtonSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Solvers/NumericalSolving/NewtonSolver.cs @@ -66,11 +66,20 @@ NumericsComplex ChooseGood() => var res = new HashSet(); var df = WithoutConditions(expr.Differentiate(v).Simplify()).Compile(v); var f = WithoutConditions(expr.Simplify()).Compile(v); + // The shares are worked out by CtxDivide, which rounds to the working precision, + // rather than by EDecimal's own operator, which has no context and answers NaN + // whenever the quotient does not terminate in base ten. Only step counts of the + // form 2^a * 5^b divide exactly, so for 3, 6, 7, 9, 12, 21 and most other values + // every share but the first came out NaN, and with it every starting point: the + // grid quietly shrank to the single corner where x and y are both zero. The + // default of 10 divides exactly, which is why this went unseen -- and why asking + // for a finer search made the answer worse rather than better. x^3 - 2x gives all + // three of its roots at 10 steps and only one at 21. for (int x = 0; x < settings.StepCount.Re; x++) for (int y = 0; y < settings.StepCount.Im; y++) { - var xShare = ((EDecimal)x) / settings.StepCount.Re; - var yShare = ((EDecimal)y) / settings.StepCount.Im; + var xShare = CtxDivide((EDecimal)x, settings.StepCount.Re); + var yShare = CtxDivide((EDecimal)y, settings.StepCount.Im); var value = Complex.Create( settings.From.Re * xShare + settings.To.Re * (1 - xShare), settings.From.Im * yShare + settings.To.Im * (1 - yShare)); diff --git a/Sources/Tests/UnitTests/Algebra/SolveTest/NewtonGridTest.cs b/Sources/Tests/UnitTests/Algebra/SolveTest/NewtonGridTest.cs new file mode 100644 index 000000000..bd676222c --- /dev/null +++ b/Sources/Tests/UnitTests/Algebra/SolveTest/NewtonGridTest.cs @@ -0,0 +1,80 @@ +// +// 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.Extensions; +using PeterO.Numbers; +using System.Collections.Generic; +using System.Linq; +using Xunit; + +namespace AngouriMath.Tests.Algebra.PolynomialSolverTests +{ + /// + /// The grid starts from was laid out with + /// EDecimal's contextless division, which answers NaN for a quotient that does not + /// terminate in base ten. Every step count outside 2^a * 5^b therefore produced one usable + /// starting point and a grid of NaN, so the search covered a single corner of the region + /// it was asked for and quietly returned fewer roots the finer it was told to look. + /// + public sealed class NewtonGridTest + { + private static HashSet Roots(string expression, int steps) + { + using var _ = MathS.Settings.NewtonSolver.Set(new MathS.Settings.NewtonSetting + { + From = (-10, -10), + To = (10, 10), + StepCount = (steps, steps), + Precision = 30, + }); + return expression.ToEntity().SolveNt("x"); + } + + private static double[] RealRoots(string expression, int steps) => + Roots(expression, steps) + .Where(root => root.ImaginaryPart.EDecimal.Abs().CompareTo(EDecimal.Create(1, -8)) < 0) + .Select(root => System.Math.Round(root.RealPart.EDecimal.ToDouble(), 6)) + .Distinct().OrderBy(value => value).ToArray(); + + /// + /// 8, 10, 16, 20 and 40 divide exactly and always worked; 3, 6, 7, 9, 12, 21 and 41 did + /// not, and used to give back at most the one root reachable from the corner. + /// + [Theory] + [InlineData(8)] + [InlineData(9)] + [InlineData(10)] + [InlineData(12)] + [InlineData(20)] + [InlineData(21)] + [InlineData(41)] + public void EveryStepCountSearchesTheWholeRegion(int steps) => + Assert.Equal(new[] { -1.414214d, 0d, 1.414214d }, RealRoots("x^3 - 2*x", steps)); + + /// + /// Asking for a finer search must not lose roots a coarser one found. This is the + /// property the defect broke, stated without naming any particular step count. + /// + [Fact] + public void AFinerSearchKeepsWhatACoarserOneFound() + { + var coarse = RealRoots("x^3 - 2*x", 10); + foreach (var steps in new[] { 11, 13, 17, 21, 33 }) + Assert.All(coarse, root => Assert.Contains(root, RealRoots("x^3 - 2*x", steps))); + } + + /// + /// A step count that does not divide exactly must lay out the same spread of starting + /// points as one that does, rather than collapsing onto a corner: at 21 steps over + /// [-10, 10] the search has to reach the roots of sin(x) - x/2 near +-1.8954942. + /// + [Fact] + public void AStepCountThatDoesNotDivideExactlyStillReachesTheOutlyingRoots() => + Assert.Equal(new[] { -1.895494d, 0d, 1.895494d }, RealRoots("sin(x) - x / 2", 21)); + } +}