Integrate numerically by Gauss-Legendre, and give the digits it knows

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

Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs

@@ -139,7 +139,12 @@ private static Entity AtTheBound(Entity antiderivative, Variable x, Entity bound
         /// <param name="x">The variable to integrate over.</param>
         /// <param name="from">The lower bound.</param>
         /// <param name="to">The upper bound.</param>
-        /// <param name="stepCount">How many steps to take; more is more accurate and slower.</param>
+        /// <param name="stepCount">
+        /// 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.
+        /// </param>
         /// <returns>The value of the integral.</returns>
         /// <remarks>
         /// This is the instance method the obsolete <c>MathS.Compute.DefiniteIntegral</c>
@@ -156,7 +161,12 @@ public Number.Complex DefiniteIntegral(Variable x, Number.Complex from, Number.C
         /// <param name="x">The variable to integrate over.</param>
         /// <param name="from">The lower bound.</param>
         /// <param name="to">The upper bound.</param>
-        /// <param name="stepCount">How many steps to take; more is more accurate and slower.</param>
+        /// <param name="stepCount">
+        /// 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.
+        /// </param>
         /// <returns>The value of the integral.</returns>
         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;
         }
         /// <summary>
-        /// 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 <see cref="Entity.DefiniteIntegral(Entity.Variable, EDecimal, EDecimal, int)"/>
+        /// The integral of a function of one variable along the segment from <paramref name="from"/> to
+        /// <paramref name="to"/>, by the five-point Gauss–Legendre rule on each of
+        /// <paramref name="panels"/> equal panels, in double precision. See
+        /// <see cref="Entity.DefiniteIntegral(Entity.Variable, EDecimal, EDecimal, int)"/>.
         /// </summary>
+        /// <remarks>
+        /// It was the mean of <c>panels + 1</c> equally weighted samples, the endpoints among them: a first-order
+        /// rule, whose error on <c>x^2</c> over <c>[0, 1]</c> 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:
+        /// <c>sin(x)/x</c> and <c>ln(x)</c> over <c>[0, 1]</c> 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
+        /// </remarks>
         /// <param name="expr">Expression to integrate</param>
         /// <param name="x">Variable to integrate over</param>
         /// <param name="from">The complex lower bound for integrating</param>
         /// <param name="to">The complex upper bound for integrating</param>
-        /// <param name="accuracy">Accuracy (for now, number of iterations)</param>
-        internal static Complex IntegrateNumerically(Entity expr, Entity.Variable x, Complex from, Complex to, int accuracy = 100)
+        /// <param name="panels">How many panels the segment is cut into, five evaluations each</param>
+        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 };
     }
 }

Sources/MCP/AngouriMath.Mcp/Numeric.cs

@@ -179,64 +179,36 @@ public static Entity Strip(Entity e) =>
     /// <summary>A quadrature result: the value, and how far it can be trusted.</summary>
     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;
 
     /// <summary>
-    /// 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.
     /// </summary>
     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
         {

Sources/MCP/AngouriMath.Mcp/SelfTest.cs

@@ -75,21 +75,6 @@ private static readonly (string Name, Func<bool> 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"),
     ];
 
     /// <summary>Runs every check, writes what it found, and returns the exit code: 1 when an identity fails.</summary>