From 31b29201ede7ce562c40abaac900fd8b71ebff69 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Thu, 22 Oct 2026 13:52:18 +0000 Subject: [PATCH] Integrate numerically by Gauss-Legendre, and give the digits it knows DefiniteIntegral was the mean of steps + 1 equally weighted samples, its bounds among them: a first-order rule, off in the third digit on x^2 over [0, 1] at the default 100 steps and printed to a hundred digits, and NaN for a convergent integral with a singular bound, sin(x)/x and ln(x) over [0, 1] among them. It is now the five-point Gauss-Legendre rule on steps panels, which never evaluates a bound and is exact on each panel for degree nine, run 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 -- so where the differences between the three do not shrink above the rounding, the integral is NaN. Where they shrink, the rest of their geometric series is added and is the error: for x^-0.9 the last difference is a fourteenth of it, and the value without it was 7 where it is 10. The rounding is measured against the integral of the modulus, so x^3 over [-1, 1] is 0 and not NaN. The value is rounded to the place of the error, at most fifteen digits, with the downcasting off so that -0.9999 does not print as -9999/10000. A singularity inside the segment is not seen, and the docs say so. The MCP server ran the old rule at 2000 and 4000 steps and reported the digits the two agreed on, and its self-test documented the first-order rule, which no longer reproduces. It now asks once and reports the digits of the value it gets, with half a unit in the last place for the error; its caveat says what the rule does at a bound, and its README lists the workarounds that remain. Fixes #1675. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 24 ++++ .../Integration/Integration.Definition.cs | 118 +++++++++++++++--- Sources/MCP/AngouriMath.Mcp/Numeric.cs | 66 +++------- Sources/MCP/AngouriMath.Mcp/SelfTest.cs | 15 --- Sources/MCP/AngouriMath.Mcp/Tools.cs | 20 ++- Sources/MCP/README.md | 10 +- .../DefiniteIntegralByGaussLegendreTest.cs | 67 ++++++++++ 7 files changed, 227 insertions(+), 93 deletions(-) create mode 100644 Sources/Tests/UnitTests/Calculus/DefiniteIntegralByGaussLegendreTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index d50c7c089..74c10497c 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -40,6 +40,30 @@ 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` | + +--- + ## 2.6.0 — since 2.5.0 ### A high half-odd power of a quadratic beside a power of another is reduced in exact rationals 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 913e7af76..f8751b54d 100644 --- a/Sources/MCP/AngouriMath.Mcp/SelfTest.cs +++ b/Sources/MCP/AngouriMath.Mcp/SelfTest.cs @@ -91,21 +91,6 @@ private static readonly (string Name, Func StillBroken, string Documented) return true; } }, "x + y = 3, x - y = 1 has no solution when written as equalities; each equation must be given in '= 0' form"), - ("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 b260d487a..b65d1498b 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 2a037db77..46435da11 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,9 +151,7 @@ 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 four workarounds here wait for it: +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)), -`a = b` rewritten to `a - b` for a system ([#1673](https://github.com/asc-community/AngouriMath/issues/1673)), -and the numeric definite integral run twice to find its good digits -([#1675](https://github.com/asc-community/AngouriMath/issues/1675)). +and `a = b` rewritten to `a - b` for a system ([#1673](https://github.com/asc-community/AngouriMath/issues/1673)). 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}"); + } + } +}