diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 7d8b0e0ac..812dbbc3f 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -3826,6 +3826,26 @@ simplified as before | `"sec(x)^2/(a + b*cos(x))^4".ToEntity().Integrate("x")` | `integral(...)` | the same | | `"sec(c + d*x)^2/(a + b*cos(c + d*x))^4".ToEntity().Integrate("x")` | `integral(...)` | the same in `tan((c + d x)/2)` | +### A polynomial times powers of polynomials is integrated by raising each power by one + +**Answers where there were none.** `(a g + e x - b g x^4)/(a + b x^4)^(3/2)`, +`(d + f x^2)^p (2 c d f + 2 b f^2 (3 + 2 p) x + 2 c f^2 (3 + 2 p) x^2)` and +`x^2 (1 - a x)^(-1 - n (1 + n)/2) (1 + a x)^(-1 - (n - 1) n/2)`, Rubi's 1.1 and 1.2, were declined. A polynomial `P` +times `R_1^(s_1) ... R_m^(s_m)`, an exponent among them not a whole number, is now `Q R_1^(s_1 + 1) ... R_m^(s_m + 1)` +for a polynomial `Q` where there is one: that derivative equal to the integrand is a linear system in the +coefficients of `Q`, solved with the symbols pinned first, and the answer is differentiated back at sampled points. +It holds wherever no base is zero. And a factor `(c x)^t` is read as `x^t` with `(c x)^t x^(-t)` in front, constant +on each interval for either sign of `c`, so `(c x)^(-1 - 2n/3)/(a + b x^n)` is answered too +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(a*g + k*x - b*g*x^4)/(a + b*x^4)^(3/2)".ToEntity().Integrate("x")` | `integral(...)` | `(g x + k x^2/(2 a))/sqrt(a + b x^4)`, 60 characters | +| `"(d + f*x^2)^p*(2*c*d*f + 2*b*f^2*(3 + 2*p)*x + 2*c*f^2*(3 + 2*p)*x^2)".ToEntity().Integrate("x")` | `integral(...)` | `(b f (3 + 2 p)/(p + 1) + 2 c f x) (d + f x^2)^(p + 1)`, 87 characters | +| `"(3*a + b*x^2)^2/(a - b*x^2)^(7/3)".ToEntity().Integrate("x")` | `integral(...)` | `(9 a x - 3 b x^3) (a - b x^2)^(-4/3)`, 59 characters | +| `"x^2*(1 - a*x)^(-1 - 1/2*n*(1 + n))*(1 + a*x)^(-1 - 1/2*(-1 + n)*n)".ToEntity().Integrate("x")` | `integral(...)` | a linear times both powers raised, 167 characters | +| `"(c*x)^(-1 - 2/3*n)/(a + b*x^n)".ToEntity().Integrate("x")` | `integral(...)` | 581 characters | + ### A polynomial times an odd half power of a quadratic with compound coefficients is answered shortly **Answers where there were none, and short ones where the unreleased master's were not.** diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index f289ef6fc..9f8a5af4e 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -15317,6 +15317,192 @@ Entity AtUnit(int which) return answer; } + /// + /// A polynomial times powers of polynomials, P R_1^(s_1) ... R_m^(s_m) with an + /// exponent that is not a whole number -- a fraction or a symbol -- answered as + /// Q R_1^(s_1 + 1) ... R_m^(s_m + 1) for a polynomial Q where there is one: + /// with R = R_1 ... R_m, the derivative of that is the integrand's powers times + /// Q' R + Q sum_i (s_i + 1) R_i' R/R_i, and that equal to P is a linear + /// system in the coefficients of Q. + /// + /// + /// + /// Rubi's (a g + e x - b g x^4)/(a + b x^4)^(3/2) is + /// (2 a g x + e x^2)/(2 a sqrt(a + b x^4)); its + /// (d + f x^2)^p (2 c d f + 2 b f^2 (3 + 2 p) x + 2 c f^2 (3 + 2 p) x^2) is + /// (b f (3 + 2 p)/(1 + p) + 2 c f x) (d + f x^2)^(p + 1); and + /// x^2 (1 - a x)^(-1 - n (1 + n)/2) (1 + a x)^(-1 - (n - 1) n/2) is a linear times both + /// powers raised. Each was declined: + /// reads one base, holding a function of x, and solves for rational coefficients; + /// wants a constant for Q; and the + /// rules for the root of a quadratic read no quartic and no symbolic exponent. + /// + /// + /// The degree of Q. The leading coefficient of the bracket is that of Q + /// times lc(R) (deg Q + sum_i (s_i + 1) deg R_i), so deg Q is + /// deg P - deg R + 1, or sigma = -sum_i (s_i + 1) deg R_i where that is a whole + /// number above it, at which the top equation leaves the top coefficient free: the quartic + /// above has sigma = 2 against 1. + /// + /// + /// Where it holds. Wherever no base is zero: (z^w)' = w z^(w - 1) z' and + /// z^(w + 1) = z z^w hold for the principal power at every z but zero, on the + /// negative real axis too. The system is solved with the symbols pinned first, which + /// declines one with no solution at once, and the answer is differentiated back at sampled + /// points. A whole exponent of -1 is a logarithm's, and a rational function, with no + /// exponent that is not whole, is the rational integrator's; both are declined. Asked late, + /// beside the base's rule, so that the substitutions and the rules for the roots of + /// quadratics answer what they read in their own forms first. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveAPolynomialTimesPowersOfPolynomialsByRaisingThePowers(Entity expr, Entity.Variable x) + { + // The exponent of each base, a factor with x in it that is no power counted as its + // base to the first. + Entity constant = Number.Integer.One; + var powers = new List<(Entity Base, Entity Exponent)>(); + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + var (@base, exponent) = factor is Powf(var b, var e) && !e.ContainsNode(x) ? (b, e) : (factor, Number.Integer.One); + if (underneath) + exponent = -exponent; + var index = powers.FindIndex(pair => pair.Base == @base); + if (index >= 0) + powers[index] = (@base, powers[index].Exponent + exponent); + else if (powers.Count == MaximumBasesRaised + 4) + return null; + else + powers.Add((@base, exponent)); + } + + // A whole positive power is a factor of P; the rest are raised. + Entity polynomial = Number.Integer.One; + var raised = new List<(Entity Base, Entity Exponent, Dictionary Coefficients, int Degree)>(); + foreach (var (@base, written) in powers) + { + var exponent = Functions.PartialFractions.Bare(written.InnerSimplified); + if (exponent.Evaled is Number.Integer whole && !exponent.Vars.Any()) + { + if (whole.EInteger.Sign > 0) + { + polynomial *= MathS.Pow(@base, whole); + continue; + } + if (whole.EInteger.Sign == 0) + continue; + if (whole == Number.Integer.MinusOne) + return null; + exponent = whole; + } + if (ReadAPolynomial(@base) is not var (coefficients, degree) || degree == 0 || raised.Count == MaximumBasesRaised) + return null; + raised.Add((@base, exponent, coefficients, degree)); + } + var degreeOfR = raised.Sum(power => power.Degree); + if (!raised.Any(power => power.Exponent.Evaled is not Number.Integer || power.Exponent.Vars.Any()) + || degreeOfR > MaximumDegreeRaised || ReadAPolynomial(polynomial) is not var (p, degreeOfP)) + return null; + + // R and the sum of (s_i + 1) R_i' R/R_i, the bracket being Q' R + Q times that. + var product = new Dictionary { [EInteger.Zero] = Number.Integer.One }; + foreach (var power in raised) + product = PolynomialProduct(product, power.Coefficients); + var beside = new Dictionary(); + Entity sigma = Number.Integer.Zero; + for (var i = 0; i < raised.Count; i++) + { + var raisedBy = (raised[i].Exponent + 1).InnerSimplified; + sigma -= raisedBy * raised[i].Degree; + var term = new Dictionary(); + foreach (var pair in PolynomialDerivative(raised[i].Coefficients)) + term[pair.Key] = (raisedBy * pair.Value).InnerSimplified; + for (var j = 0; j < raised.Count; j++) + if (j != i) + term = PolynomialProduct(term, raised[j].Coefficients); + beside = PolynomialSum(beside, term); + } + var degreeOfQ = degreeOfP - degreeOfR + 1; + if (LowestOverTheSymbols(Functions.PartialFractions.Bare(sigma.InnerSimplified)) is var sigmaRead && !sigmaRead.Vars.Any() + && sigmaRead.Evaled is Number.Integer free && free.EInteger.CanFitInInt32() && free.EInteger.ToInt32Unchecked() > degreeOfQ) + degreeOfQ = free.EInteger.ToInt32Unchecked(); + if (degreeOfQ < 0 || degreeOfQ > MaximumDegreeRaised) + return null; + + // The bracket of each monomial of Q is a column, one row a power of x. + var columns = new List>(); + for (var k = 0; k <= degreeOfQ; k++) + { + var monomial = new Dictionary { [EInteger.FromInt32(k)] = Number.Integer.One }; + columns.Add(PolynomialSum(PolynomialProduct(PolynomialDerivative(monomial), product), PolynomialProduct(monomial, beside))); + } + var rows = System.Math.Max(degreeOfP, degreeOfQ + degreeOfR - 1) + 1; + var matrix = new Entity[rows][]; + var rhs = new Entity[rows]; + for (var row = 0; row < rows; row++) + { + var power = EInteger.FromInt32(row); + matrix[row] = new Entity[columns.Count]; + for (var column = 0; column < columns.Count; column++) + matrix[row][column] = columns[column].TryGetValue(power, out var entry) ? entry : Number.Integer.Zero; + rhs[row] = p.TryGetValue(power, out var wanted) ? wanted : Number.Integer.Zero; + } + if (!Functions.PartialFractions.TrySolveLinearWithSymbols(matrix, rhs, out var values)) + return null; + Entity q = Number.Integer.Zero; + for (var k = 0; k <= degreeOfQ; k++) + { + var value = LowestOverTheSymbols(Functions.PartialFractions.Bare(values[k].InnerSimplified)); + if (value.Evaled is Number.Complex { IsZero: true }) + continue; + q += value * (k == 0 ? Number.Integer.One : k == 1 ? x : MathS.Pow(x, k)); + } + if (q == Number.Integer.Zero || q.Nodes.Any(node => node == MathS.NaN)) + return null; + Entity answer = constant * q; + foreach (var power in raised) + { + var raisedTo = (power.Exponent + 1).InnerSimplified; + if (raisedTo.Evaled is not Number.Complex { IsZero: true }) + answer *= raisedTo == Number.Integer.One ? power.Base : MathS.Pow(power.Base, raisedTo); + } + answer = answer.InnerSimplified; + return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; + + // The coefficients from x^0 up and the degree; null where it is not a polynomial in x + // or is of a degree past what is solved for. + (Dictionary Coefficients, int Degree)? ReadAPolynomial(Entity e) + { + if (!TreeAnalyzer.TryGetPolynomial(e, x, out var read)) + return null; + var coefficients = new Dictionary(); + var degree = 0; + foreach (var pair in read) + { + if (pair.Key.Sign < 0 || pair.Key.CompareTo(EInteger.FromInt32(2 * MaximumDegreeRaised)) > 0 || pair.Value.ContainsNode(x)) + return null; + var simplified = pair.Value.InnerSimplified; + if (simplified.Evaled is Number.Complex { IsZero: true }) + continue; + coefficients[pair.Key] = simplified; + degree = System.Math.Max(degree, pair.Key.ToInt32Unchecked()); + } + return (coefficients, degree); + } + } + + /// + /// The most bases + /// raises, and the largest degree of their product and of the polynomial it solves for. + /// + private const int MaximumBasesRaised = 4; + private const int MaximumDegreeRaised = 12; + /// /// An integrand e^h R, with R rational in x and h rational /// in x or absent, integrated by the ansatz F = e^h N/D: D read @@ -34987,6 +35173,13 @@ private static bool IsARationalFunctionOfExponentials(Entity expr, Entity.Variab /// x^(4n), and the unit is each exponent left, divided by up to six. /// https://github.com/asc-community/AngouriMath/issues/718 /// + /// + /// A factor (c x)^t is read as x^t, with (c x)^t x^(-t) in front of the + /// answer: its derivative is zero wherever it is defined, so it is a constant on each + /// interval of the domain, for either sign of c. Read as it was written, Rubi's + /// (c x)^(-1 - 2n/3)/(a + b x^n) left an x outside every power of x and was declined. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// /// internal static Entity? SolveByAPowerOfXWithASymbolicExponent(Entity expr, Entity.Variable x, bool integrateByParts) { @@ -34999,6 +35192,28 @@ private static bool IsARationalFunctionOfExponentials(Entity expr, Entity.Variab if (expr.Nodes.Any(node => node.ContainsNode(x) && node is not (Sumf or Minusf or Mulf or Divf or Variable) && !(node is Powf(_, var exponent) && !exponent.ContainsNode(x)))) return null; + // A factor that is a power of a multiple of x, `(c x)^t` with t not whole, is `x^t` times + // `(c x)^t x^(-t)`, whose derivative is zero wherever it is defined, since + // `t c (c x)^(t - 1)` is `t (c x)^t/x` for the principal power: that goes in front, + // constant on each interval of the domain, with no sign assumed for c. Rubi's + // `(c x)^(-1 - 2n/3)/(a + b x^n)` was declined, where `x^(-1 - 2n/3)/(a + b x^n)` is answered. + Entity inFront = Number.Integer.One; + Entity withPowersOfX = Number.Integer.One; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + var written = factor; + if (factor is Powf(var multiple, var t) && multiple != x && !t.ContainsNode(x) && t.Evaled is not Number.Integer + && TreeAnalyzer.TryGetPolyLinear(multiple, x, out var slope, out var intercept) + && intercept.Evaled is Number.Complex { IsZero: true } && !slope.ContainsNode(x)) + { + written = MathS.Pow(x, t); + var ratio = MathS.Pow(multiple, t) * MathS.Pow(x, -t); + inFront = underneath ? inFront / ratio : inFront * ratio; + } + withPowersOfX = underneath ? withPowersOfX / written : withPowersOfX * written; + } + if (inFront != Number.Integer.One) + expr = withPowersOfX; // Over one bar first, so that an x inside `(x (a + b x^n))^(-1)` is a factor to gather. var g = Patterns.GatherPowersOfOneBase(Functions.SingleQuotient.Combine(x * expr)); var exponents = g.Nodes.OfType().Where(power => power.Base == x && power.Exponent is not Number && !power.Exponent.ContainsNode(x)) @@ -35041,7 +35256,7 @@ private static bool IsARationalFunctionOfExponentials(Entity expr, Entity.Variab return null; foreach (var (name, value) in named) inTermsOfU = inTermsOfU.Substitute(name, value); - var answer = (d / unit * inTermsOfU).Substitute(u, MathS.Pow(x, step)).InnerSimplified; + var answer = (inFront * d / unit * inTermsOfU).Substitute(u, MathS.Pow(x, step)).InnerSimplified; return answer.Nodes.Any(node => node is Number.Complex { IsNaN: true }) ? null : answer; // Arithmetic and whole powers taken apart; a root of something with u in it keeps diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 779ce6448..d55feb5a3 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -1396,6 +1396,11 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // answered as a polynomial times the next power of the base: a linear system in // the polynomial's coefficients, exact, and volunteered at any depth. if ((answer = IndefiniteIntegralSolver.SolveByAPolynomialTimesAPowerOfTheBase(expr, x)) is { }) return answer; + // And a polynomial times powers of polynomials, one at least not whole, as another + // polynomial times each power raised by one: Rubi's + // `(a g + e x - b g x^4)/(a + b x^4)^(3/2)` and `(d + f x^2)^p (...)`, which no rule for + // a root of a quadratic reads. After the substitutions, which answer their shapes first. + if ((answer = IndefiniteIntegralSolver.SolveAPolynomialTimesPowersOfPolynomialsByRaisingThePowers(expr, x)) is { }) return answer; // And an exponential times a rational function of the sine and cosine, closed by // the same kind of ansatz in the half-angle tangent: `e^x (1 - sin(x))/(1 - cos(x))` // is `-e^x cot(x/2)`. diff --git a/Sources/Tests/UnitTests/Calculus/APolynomialTimesPowersOfPolynomialsByRaisingThePowersIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APolynomialTimesPowersOfPolynomialsByRaisingThePowersIntegralTest.cs new file mode 100644 index 000000000..6c42d6335 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APolynomialTimesPowersOfPolynomialsByRaisingThePowersIntegralTest.cs @@ -0,0 +1,88 @@ +// +// 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 Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A polynomial times powers of polynomials, one at least not whole, whose antiderivative is + /// another polynomial times each power raised by one: Rubi's + /// (a g + e x - b g x^4)/(a + b x^4)^(3/2) and + /// (d + f x^2)^p (2 c d f + 2 b f^2 (3 + 2 p) x + 2 c f^2 (3 + 2 p) x^2), which were declined. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class APolynomialTimesPowersOfPolynomialsByRaisingThePowersIntegralTest + { + [Theory] + // Rubi's 1.1.3.8 and 1.2.2.5: a quartic under the root. + [InlineData("(a*g + k*x - b*g*x^4)/(a + b*x^4)^(3/2)")] + [InlineData("(a*g + f*x^3 - b*g*x^4)/(a + b*x^4)^(3/2)")] + [InlineData("(a*g + k*x + f*x^3 - b*g*x^4)/(a + b*x^4)^(3/2)")] + [InlineData("(a*g + k*x - c*g*x^4)/(a + b*x^2 + c*x^4)^(3/2)")] + [InlineData("(a*g + f*x^3 - c*g*x^4)/(a + b*x^2 + c*x^4)^(3/2)")] + [InlineData("(a*g + k*x + f*x^3 - c*g*x^4)/(a + b*x^2 + c*x^4)^(3/2)")] + // Rubi's 1.2.1.9: a symbolic exponent. + [InlineData("(d + f*x^2)^p*(2*c*d*f + 2*b*f^2*(3 + 2*p)*x + 2*c*f^2*(3 + 2*p)*x^2)")] + [InlineData("(d + k*x + f*x^2)^p*(-2*c*k^2 + 2*c*d*f - c*k^2*p + 2*c*f^2*(3 + 2*p)*x^2)")] + [InlineData("(d + k*x + f*x^2)^p*(-2*c*k^2 + 2*c*d*f + 3*b*k*f - c*k^2*p + 2*b*k*f*p + 2*b*f^2*(3 + 2*p)*x + 2*c*f^2*(3 + 2*p)*x^2)")] + // Rubi's 1.1.2.3: a cube root. + [InlineData("(3*a + b*x^2)^2/(a - b*x^2)^(7/3)")] + // Rubi's 1.1.1.3: two linears, each raised. + [InlineData("x^2*(1 - a*x)^(-1 - 1/2*n*(1 + n))*(1 + a*x)^(-1 - 1/2*(-1 + n)*n)")] + public void AsAPolynomialTimesThePowersRaised(string integrand) => AssertAnAntiderivative(integrand); + + [Theory] + // Rubi's 1.1.3.8 and 1.2.3.5: the same under u = x^(n/4) or u = x^(n/2), which the + // substitution for a symbolic power of x asks. + [InlineData("(-a*h*x^(-1 + 1/4*n) + b*f*x^(-1 + 1/2*n) + b*g*x^(-1 + n) + b*h*x^(-1 + 5/4*n))/(a + b*x^n)^(3/2)")] + [InlineData("(-a*h*x^(-1 + 1/2*n) + c*f*x^(-1 + n) + c*g*x^(-1 + 2*n) + c*h*x^(-1 + 5/2*n))/(a + b*x^n + c*x^(2*n))^(3/2)")] + [InlineData("x^(-1 + 1/4*n)*(-a*h + c*f*x^(1/4*n) + c*g*x^(3/4*n) + c*h*x^n)/(a + c*x^n)^(3/2)")] + [InlineData("x^(-1 + 1/2*n)*(-a*h + c*f*x^(1/2*n) + c*g*x^(3/2*n) + c*h*x^(2*n))/(a + b*x^n + c*x^(2*n))^(3/2)")] + public void UnderAPowerOfX(string integrand) => AssertAnAntiderivative(integrand); + + [Theory] + // Rubi's 1.2.3.5 and 1.1.3.2: a power of a multiple of x beside the powers of x, read as + // the power of x with the ratio, constant on each interval, in front -- for a negative + // multiple too, where both are complex and their ratio is still a constant. + [InlineData("(d*x)^(-1 + 1/4*n)*(-a*h + c*f*x^(1/4*n) + c*g*x^(3/4*n) + c*h*x^n)/(a + c*x^n)^(3/2)")] + [InlineData("(d*x)^(-1 + 1/2*n)*(-a*h + c*f*x^(1/2*n) + c*g*x^(3/2*n) + c*h*x^(2*n))/(a + b*x^n + c*x^(2*n))^(3/2)")] + [InlineData("(c*x)^(-1 - n)/(a + b*x^n)")] + [InlineData("(c*x)^(-1 - 2/3*n)/(a + b*x^n)")] + [InlineData("(c*x)^(-1 - 5/4*n)/(a + b*x^n)")] + [InlineData("(-c*x)^(-1 - n)/(a + b*x^n)")] + public void BesideAPowerOfAMultipleOfX(string integrand) => AssertAnAntiderivative(integrand); + + private static void AssertAnAntiderivative(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 2000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.4).Substitute("d", 0.3) + .Substitute("f", 0.9).Substitute("g", 1.1).Substitute("h", 0.6).Substitute("k", 0.5).Substitute("p", 0.35).Substitute("n", 2.3); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { 0.3, 0.55, 0.8, 1.3, 1.7, 2.9 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +}