diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 9f1a36a05..8e9e914a1 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -42,6 +42,26 @@ A not-equal written `<>` is read as `not (a = b)`, and printed that way. ## Unreleased — since 2.6.0 +### A numerical definite integral is Gauss–Legendre, with the digits it knows + +**Silent: wrong answers, and values where there was NaN.** `Entity.DefiniteIntegral(x, a, b, steps)` was the mean of +`steps + 1` equally weighted samples, its bounds among them: a first-order rule, whose value on `x^2` over `[0, 1]` was +off in the third digit at the default 100 steps and printed to a hundred digits. A convergent integral with a singular +bound had no value. It is now the five-point Gauss–Legendre rule on `steps` panels, which never evaluates a bound, run +on the panels, twice as many and four times as many. Where the differences between the three do not shrink, the +integral diverges at a bound and is `NaN`. Otherwise the rest of their geometric series is added and taken as the +error, and the value is rounded to its place, at most fifteen digits, a double's. A singularity inside the segment is +not seen, and a symmetric divergent one comes to its principal value +([#1675](https://github.com/asc-community/AngouriMath/issues/1675)). + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"x^2".ToEntity().DefiniteIntegral("x", 0, 1)` | `0.334999999999999938…`, a hundred digits | `0.333333333333333` | +| `"e^x".ToEntity().DefiniteIntegral("x", 0, 2)` | `6.409068937660210827…` | `6.38905609893066` | +| `"sin(x)/x".ToEntity().DefiniteIntegral("x", 0, 1)` | `NaN` | `0.946083070367182` | +| `"ln(x)".ToEntity().DefiniteIntegral("x", 0, 1)` | `NaN` | `-1.00000` | +| `"1/x".ToEntity().DefiniteIntegral("x", 0, 1)` | `NaN` | `NaN` | + ### An equality in an equation or a system is read as a difference equal to zero **Silent: wrong answers.** `MathS.SolveEquation`, `Entity.SolveEquation` and `MathS.Equations` read each equation as an diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 150528623..62ca25bb0 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -139,7 +139,12 @@ private static Entity AtTheBound(Entity antiderivative, Variable x, Entity bound /// The variable to integrate over. /// The lower bound. /// The upper bound. - /// How many steps to take; more is more accurate and slower. + /// + /// How many panels to cut the segment into, each integrated by the five-point Gauss–Legendre rule; more + /// is more accurate and slower. Neither bound is evaluated, so an integrable singularity at one is + /// allowed, and one that diverges is NaN. The value has the digits the rule agrees on as the panels are + /// doubled, at most fifteen, a double's. A singularity inside the segment is not seen. + /// /// The value of the integral. /// /// This is the instance method the obsolete MathS.Compute.DefiniteIntegral @@ -156,7 +161,12 @@ public Number.Complex DefiniteIntegral(Variable x, Number.Complex from, Number.C /// The variable to integrate over. /// The lower bound. /// The upper bound. - /// How many steps to take; more is more accurate and slower. + /// + /// How many panels to cut the segment into, each integrated by the five-point Gauss–Legendre rule; more + /// is more accurate and slower. Neither bound is evaluated, so an integrable singularity at one is + /// allowed, and one that diverges is NaN. The value has the digits the rule agrees on as the panels are + /// doubled, at most fifteen, a double's. A singularity inside the segment is not seen. + /// /// The value of the integral. public Number.Complex DefiniteIntegral(Variable x, EDecimal from, EDecimal to, int stepCount = 100) => Integration.IntegrateNumerically(this, x, @@ -1825,26 +1835,106 @@ internal static bool ASumOverRootsWouldAnswer return null; } /// - /// Returns the approximate numeric value of a definite integral of a function. Only works for one-variable functions. - /// Accuracy is limited to the number specified (default is 100). - /// See more at + /// The integral of a function of one variable along the segment from to + /// , by the five-point Gauss–Legendre rule on each of + /// equal panels, in double precision. See + /// . /// + /// + /// It was the mean of panels + 1 equally weighted samples, the endpoints among them: a first-order + /// rule, whose error on x^2 over [0, 1] halved with each doubling, 1.67e-4 at 1000 steps, + /// and printed to a hundred digits. And a convergent integral with a singular endpoint had no value: + /// sin(x)/x and ln(x) over [0, 1] were NaN. Gauss–Legendre never evaluates an + /// endpoint, and the five-point rule is exact on each panel for a polynomial of degree nine. It is run on + /// the panels, twice as many and four times as many: where the differences do not shrink, the integral + /// diverges at a bound and is NaN, and otherwise the rest of their geometric series is added and is the + /// error, and the value is rounded to its place, fifteen digits at most. A singularity inside the + /// segment is not seen, and a symmetric divergent one, as in 1/x over [-1, 1], comes to its principal + /// value. + /// https://github.com/asc-community/AngouriMath/issues/1675 + /// /// Expression to integrate /// Variable to integrate over /// The complex lower bound for integrating /// The complex upper bound for integrating - /// Accuracy (for now, number of iterations) - internal static Complex IntegrateNumerically(Entity expr, Entity.Variable x, Complex from, Complex to, int accuracy = 100) + /// How many panels the segment is cut into, five evaluations each + internal static Complex IntegrateNumerically(Entity expr, Entity.Variable x, Complex from, Complex to, int panels = 100) { - System.Numerics.Complex res = 0; - var cfunc = expr.Compile(x); - for (int i = 0; i <= accuracy; i++) + var function = expr.Compile(x); + var start = from.ToNumerics(); + var length = to.ToNumerics() - start; + panels = System.Math.Max(1, panels); + // The integral on so many panels, and the integral of the modulus, which the rounding is measured + // against: an integral that comes to zero, x^3 over [-1, 1], is a sum of terms that do not. + (System.Numerics.Complex Value, double OfTheModulus) OnPanels(int count) { - var share = ((EDecimal)i) / accuracy; - var tmp = Complex.Create(from.RealPart.EDecimal * share + to.RealPart.EDecimal * (1 - share), from.ImaginaryPart.EDecimal * share + to.ImaginaryPart.EDecimal * (1 - share)); - res += cfunc.Substitute(tmp.ToNumerics()); + System.Numerics.Complex sum = 0; + var ofTheModulus = 0.0; + for (var i = 0; i < count; i++) + { + var centre = start + length * ((i + 0.5) / count); + for (var k = 0; k < GaussLegendreNodes.Length; k++) + { + var term = GaussLegendreWeights[k] * function.Substitute(centre + length * (GaussLegendreNodes[k] / (2.0 * count))); + sum += term; + ofTheModulus += System.Numerics.Complex.Abs(term); + } + } + var width = System.Numerics.Complex.Abs(length) / (2.0 * count); + return (sum * length / (2.0 * count), ofTheModulus * width); } - return res.ToNumber() / (accuracy + 1) * (to - from); + // On the panels, twice as many and four times as many. Never evaluating a bound, the rule gives a + // divergent integral a finite value: 1/x over [0, 1] came to about 10. Where it converges the + // differences between the three shrink, by r = 2^-(1 - p) for an endpoint singularity like x^-p, + // p < 1; for 1/x they stay at ln(2), and they grow for anything worse. So differences that do not + // shrink, above the rounding, are NaN; and where they shrink, the rest of the geometric series, + // d r/(1 - r) for the last difference d, is both the correction and the error: for x^-0.9 the last + // difference is a fourteenth of the error, and the value without it was 7 where it is 10. + var (coarse, _) = OnPanels(panels); + var (middle, _) = OnPanels(2 * panels); + var (fine, ofTheModulus) = OnPanels(4 * panels); + if (double.IsNaN(fine.Real) || double.IsNaN(fine.Imaginary) || double.IsInfinity(fine.Real) || double.IsInfinity(fine.Imaginary)) + return fine.ToNumber(); + var rounding = 1e-14 * System.Math.Max(ofTheModulus, System.Numerics.Complex.Abs(fine)); + if (rounding == 0) + return Integer.Zero; + var first = System.Numerics.Complex.Abs(middle - coarse); + var last = System.Numerics.Complex.Abs(fine - middle); + if (last > rounding && last >= 0.99 * first) + return Real.NaN; + var value = fine; + var error = System.Math.Max(last, rounding); + if (last > rounding && first > rounding) + { + var ratio = last / first; + var tail = (fine - middle) * (ratio / (1 - ratio)); + value = fine + tail; + error = System.Math.Max(System.Numerics.Complex.Abs(tail), rounding); + } + // Then rounded to the place of that error, and to no more digits than a double holds: the value + // was the double's exact binary expansion, fifty digits of which most were not known. + var magnitude = System.Numerics.Complex.Abs(value); + var quantum = (int)System.Math.Floor(System.Math.Log10(error)); + if (magnitude > 0) + quantum = System.Math.Max(quantum, (int)System.Math.Floor(System.Math.Log10(magnitude)) - MostDigitsOfADouble + 1); + // A part that rounds to zero is zero, which prints as 0 and not as -0E-15; and a value with no + // imaginary part is real. + EDecimal Rounded(double part) + => EDecimal.FromDouble(part).Quantize(quantum, ERounding.HalfEven) is var rounded && rounded.IsZero ? EDecimal.Zero : rounded; + var (realPart, imaginaryPart) = (Rounded(value.Real), Rounded(value.Imaginary)); + using var _ = MathS.Settings.DowncastingEnabled.Set(false); + return imaginaryPart.IsZero ? Real.Create(realPart) : Complex.Create(realPart, imaginaryPart); } + + private const int MostDigitsOfADouble = 15; + + // The five-point Gauss–Legendre rule on [-1, 1]: the roots of the fifth Legendre polynomial, + // 0 and ±sqrt(5 ∓ 2 sqrt(10/7))/3, and their weights. + [ConstantField] + private static readonly double[] GaussLegendreNodes = + { 0.0, -0.5384693101056831, 0.5384693101056831, -0.9061798459386640, 0.9061798459386640 }; + [ConstantField] + private static readonly double[] GaussLegendreWeights = + { 0.5688888888888889, 0.4786286704993665, 0.4786286704993665, 0.2369268850561891, 0.2369268850561891 }; } } diff --git a/Sources/MCP/AngouriMath.Mcp/Numeric.cs b/Sources/MCP/AngouriMath.Mcp/Numeric.cs index 7964bf784..4e12eceae 100644 --- a/Sources/MCP/AngouriMath.Mcp/Numeric.cs +++ b/Sources/MCP/AngouriMath.Mcp/Numeric.cs @@ -179,64 +179,36 @@ public static Entity Strip(Entity e) => /// A quadrature result: the value, and how far it can be trusted. public readonly record struct Quadrature(double Value, double ErrorEstimate, int Digits); - // Step counts for the two runs. AngouriMath's rule is first order — measured, by - // doubling: the error of `∫[0,1] e^(x^2)` halves at every doubling of the step count - // rather than falling by four — so the gap between a run and its double IS the error of - // the coarser one, to a constant. Two runs at these counts cost about a third of a - // second together, which is inside the request budget; going further buys one digit per - // quadrupling and is not worth the wall clock. - private const int CoarseSteps = 2_000; - private const int FineSteps = 4_000; + // Panels for the library's rule: five evaluations each, on these, twice as many and four times + // as many, which is seven thousand evaluations of a compiled function. + private const int Panels = 200; /// - /// Numeric definite integration via `Entity.DefiniteIntegral`, run TWICE so that the - /// answer arrives with an error estimate rather than as a hundred digits of false - /// precision. https://github.com/asc-community/AngouriMath/issues/1675 + /// Numeric definite integration via `Entity.DefiniteIntegral`, with the digits it reports. The + /// library's rule is five-point Gauss-Legendre, run on its panels, twice as many and four times + /// as many, and its value is rounded to the place of the error that leaves, so the digits of the + /// value it returns are the digits it knows: `∫[0,1] ln(x)` is -1.000000, six of them. + /// https://github.com/asc-community/AngouriMath/issues/1675 /// - /// That precision is the whole reason this is not a one-line call. The library returns - /// an arbitrary-precision decimal whatever the step count, so `∫[0,1] e^(x^2)` prints as - /// 1.4666222366253148581140206684127892598 when the true value is 1.46265174590718 — a - /// figure that is wrong in its third decimal and reads as converged to thirty. Reporting - /// only the digits two runs agree on is exactly the job this server exists to do. - /// - /// Returns null when the integrand cannot be sampled: the rule evaluates at the - /// endpoints, so a singularity at either — `∫[0,1] sin(x)/x`, `∫[0,1] ln(x)`, both of - /// which converge — comes back NaN and must be declined rather than reported. + /// Returns null where there is nothing to report: a NaN, which is how the library answers an + /// integral that diverges at a bound, or a complex value for a real integrand over a real + /// interval, which means the rule left the domain. /// public static Quadrature? TryDefiniteIntegral(Entity integrand, Variable variable, double from, double to) - { - var coarse = Sample(integrand, variable, from, to, CoarseSteps); - var fine = Sample(integrand, variable, from, to, FineSteps); - if (coarse is not { } c || fine is not { } f) return null; - - var error = Math.Abs(f - c); - - // Digits justified by the error, relative to the magnitude of the answer: an error - // of 2e-4 on a value near 1.46 pins three decimals, and the fourth is noise. One - // digit is dropped on purpose — the estimate is itself an estimate. - var scale = Math.Max(Math.Abs(f), double.Epsilon); - var digits = error <= 0 - ? 15 - : (int)Math.Floor(-Math.Log10(error / scale)) - 1; - digits = Math.Clamp(digits, 0, 15); - - return new Quadrature(f, error, digits); - } - - private static double? Sample(Entity integrand, Variable variable, - double from, double to, int steps) { try { - var value = integrand.DefiniteIntegral(variable, from, to, steps); - var real = value.RealPart.EDecimal.ToDouble(); + var value = integrand.DefiniteIntegral(variable, from, to, Panels); + var real = value.RealPart.EDecimal; var imaginary = value.ImaginaryPart.EDecimal.ToDouble(); - if (double.IsNaN(real) || double.IsInfinity(real)) return null; - // A real integrand integrated over a real interval that comes back complex means - // the rule wandered off the domain; that is not an answer to report. + if (real.IsNaN() || real.IsInfinity()) return null; if (double.IsNaN(imaginary) || Math.Abs(imaginary) > 1e-12) return null; - return real; + // A value the library rounded to zero has no significant digit to count, and was + // rounded at the fifteenth of the integral of the modulus at most. + var digits = real.IsZero ? 15 : Math.Min(15, real.Mantissa.Abs().ToString().Length); + var error = 0.5 * Math.Pow(10, real.Exponent.ToInt32Checked()); + return new Quadrature(real.ToDouble(), error, digits); } catch { diff --git a/Sources/MCP/AngouriMath.Mcp/SelfTest.cs b/Sources/MCP/AngouriMath.Mcp/SelfTest.cs index 5dad347e9..600a25504 100644 --- a/Sources/MCP/AngouriMath.Mcp/SelfTest.cs +++ b/Sources/MCP/AngouriMath.Mcp/SelfTest.cs @@ -75,21 +75,6 @@ private static readonly (string Name, Func StillBroken, string Documented) [ ("Simplify(sqrt(x^2))", () => Value("sqrt(x^2)") != "abs(x)", "left as written rather than reduced to abs(x)"), - ("DefiniteIntegral is a first-order rule", () => - { - // Halving the step count roughly doubles the error iff the rule is first order. - // https://github.com/asc-community/AngouriMath/issues/1675 - // Numeric.TryDefiniteIntegral's step counts, its error estimate and the caveat - // am_integrate prints all assume that; a rule of higher order would make the - // reported error bound far too pessimistic and the wording wrong. - var e = "e^(x^2)".ToEntity(); - var v = MathS.Var("x"); - var coarse = e.DefiniteIntegral(v, 0, 1, 500).RealPart.EDecimal.ToDouble(); - var fine = e.DefiniteIntegral(v, 0, 1, 1000).RealPart.EDecimal.ToDouble(); - var finer = e.DefiniteIntegral(v, 0, 1, 2000).RealPart.EDecimal.ToDouble(); - var ratio = Math.Abs(fine - coarse) / Math.Abs(finer - fine); - return ratio is > 1.7 and < 2.3; - }, "error halves per doubling, so ~4 digits at 4000 steps; Simpson would give more"), ]; /// Runs every check, writes what it found, and returns the exit code: 1 when an identity fails. diff --git a/Sources/MCP/AngouriMath.Mcp/Tools.cs b/Sources/MCP/AngouriMath.Mcp/Tools.cs index 89b5359f7..40c933e54 100644 --- a/Sources/MCP/AngouriMath.Mcp/Tools.cs +++ b/Sources/MCP/AngouriMath.Mcp/Tools.cs @@ -857,10 +857,8 @@ private static JsonObject Integrate(JsonObject args) if (r.Item4 is { } q) { - // Rounded to the digits two step counts agree on. Printing the library's - // full arbitrary-precision output here would be a hundred digits of a - // number correct to four, which is the failure this server exists to - // prevent — the caller cannot see the step count and would read it as exact. + // To the digits the library's rule knows, which are the digits of the value it + // returned: the caller cannot see the rule and would read more as exact. response["numeric_definite_value"] = Math.Round(q.Value, q.Digits) .ToString("G" + Math.Max(q.Digits + 1, 2), System.Globalization.CultureInfo.InvariantCulture); @@ -868,13 +866,13 @@ private static JsonObject Integrate(JsonObject args) .ToString("G3", System.Globalization.CultureInfo.InvariantCulture); response["numeric_definite_caveat"] = "APPROXIMATE, and produced only because the exact route gave nothing " + - "usable. It is quadrature over 4000 steps, cross-checked against 2000, " + - "and reported to the digits the two agree on — every digit shown is " + - "believed, and there are no more to be had from this call. Unlike the " + - "symbolic answer it is NOT verified by differentiating back. It also " + - "assumes the integrand is finite across the closed interval: the rule " + - "samples both endpoints, so a convergent integral with a singular " + - "endpoint returns nothing here rather than a value."; + "usable. It is five-point Gauss-Legendre quadrature on 200 panels, run " + + "again on 400 and 800, and reported to the digits the runs leave known — " + + "every digit shown is believed, and there are no more to be had from this " + + "call. Unlike the symbolic answer it is NOT verified by differentiating " + + "back. The rule never samples an endpoint, so a convergent integral with " + + "a singular endpoint has a value, and one that diverges there returns " + + "nothing; a singularity inside the interval is not detected."; if ((string?)response["status"] == "declined") response["note"] = (string?)response["note"] + diff --git a/Sources/MCP/README.md b/Sources/MCP/README.md index cfd5afc61..856d6ccba 100644 --- a/Sources/MCP/README.md +++ b/Sources/MCP/README.md @@ -104,8 +104,8 @@ raised `NotSufficientlySupportedException`. `unchanged` means no progress, not " **A printed `NaN` is screened.** An answer that exists has none in it. **An approximate answer says how approximate it is.** The numeric fallback for a declined definite -integral, `Entity.DefiniteIntegral`, is a first-order rule that returns a hundred digits whatever -its accuracy. It is run at two step counts, and only the digits they agree on are reported. +integral, `Entity.DefiniteIntegral`, returns the digits its rule knows, and only those are reported, +with the error they leave. **Every call has a budget and its own stack.** Each runs under the library's cancellation token on a thread with a 64 MB stack, abandoned rather than killed on timeout, so neither a runaway search nor @@ -151,8 +151,5 @@ existing tool, then a prompt, and add a tool last. A tool's description is a rou documentation: say when to call it, and where a model would wrongly trust itself, say so. This is an adapter, not a second computer algebra system. Mathematics belongs in the library, where -the terminal and the notebooks get it too, and three workarounds here wait for it: -eigenvalues by the characteristic polynomial ([#1676](https://github.com/asc-community/AngouriMath/issues/1676)), -a numeric comparison with a tolerance ([#1674](https://github.com/asc-community/AngouriMath/issues/1674)), -and the numeric definite integral run twice to find its good digits -([#1675](https://github.com/asc-community/AngouriMath/issues/1675)). +the terminal and the notebooks get it too. Where the server has to work around the library, the +workaround links the library's issue in a comment beside it, and the fix for that issue removes it. diff --git a/Sources/Tests/UnitTests/Calculus/DefiniteIntegralByGaussLegendreTest.cs b/Sources/Tests/UnitTests/Calculus/DefiniteIntegralByGaussLegendreTest.cs new file mode 100644 index 000000000..ade5912ca --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/DefiniteIntegralByGaussLegendreTest.cs @@ -0,0 +1,67 @@ +// +// 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; +using AngouriMath.Extensions; +using PeterO.Numbers; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// is composite five-point + /// Gauss–Legendre, where it was a first-order rule printed to a hundred digits: x^2 over [0, 1] + /// was off by 1.67e-4 at 1000 steps. It never evaluates a bound, so a convergent integral with a singular + /// endpoint has a value, where sin(x)/x and ln(x) over [0, 1] were NaN; one that + /// diverges there is NaN; and the value has the digits the rule agrees on, fifteen at most. + /// #1675 + /// + [Trait("Area", "Calculus")] + public sealed class DefiniteIntegralByGaussLegendreTest + { + private static double Integral(string integrand, double from, double to, int steps = 100) + => integrand.ToEntity().DefiniteIntegral("x", EDecimal.FromDouble(from), EDecimal.FromDouble(to), steps).RealPart.EDecimal.ToDouble(); + + [Theory] + [InlineData("x^2", 0, 1, 1.0 / 3, 1e-14)] + [InlineData("e^x", 0, 2, 6.38905609893065, 1e-12)] + [InlineData("sin(x)", 0, 3.141592653589793, 2, 1e-12)] + [InlineData("x^3", -1, 1, 0, 1e-14)] + // A removable singularity at a bound: Si(1). + [InlineData("sin(x)/x", 0, 1, 0.946083070367183, 1e-12)] + // Integrable singularities at a bound, where the rule converges slowly and the rest of the + // series of differences is added. + [InlineData("ln(x)", 0, 1, -1, 1e-4)] + [InlineData("1/sqrt(x)", 0, 1, 2, 1e-2)] + public void IsTheIntegral(string integrand, double from, double to, double expected, double tolerance) + { + var value = Integral(integrand, from, to); + Assert.False(double.IsNaN(value), $"{integrand} over [{from}, {to}] is NaN"); + Assert.True(Math.Abs(value - expected) <= tolerance, $"{integrand} over [{from}, {to}] is {value}, where it is {expected}"); + } + + [Theory] + [InlineData("1/x")] + [InlineData("1/x^2")] + public void DivergesAtABoundToNaN(string integrand) + => Assert.True(double.IsNaN(Integral(integrand, 0, 1)), $"{integrand} over [0, 1] diverges, and came to {Integral(integrand, 0, 1)}"); + + [Fact] + public void IntegratesAlongAComplexSegment() + { + var value = "x".ToEntity().DefiniteIntegral("x", 0, Entity.Number.Complex.ImaginaryOne, 10); + Assert.True((value + Entity.Number.Rational.Create(1, 2)).Abs().EDecimal.ToDouble() < 1e-14, $"x over [0, i] is {value}"); + } + + [Fact] + public void HasNoMoreDigitsThanADouble() + { + var text = "x^2".ToEntity().DefiniteIntegral("x", 0, 1, 1000).Stringize(); + Assert.True(text.TrimStart('-', '0', '.').Length <= 15, $"x^2 over [0, 1] printed as {text}"); + } + } +}