From 064f6ec642366586673964226b1c88fb158dc523 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 30 Sep 2026 16:35:28 +0000 Subject: [PATCH] A rational integral over the roots of its denominator The Rothstein-Trager resultant declined a rational function whose residues lie in a field of degree above two: 1/(x^3 + x + 1) and #1285's sqrt(x)/(1 + x + x^4), 2u^2/(1 + u^2 + u^8) under u = sqrt(x), were left unevaluated. Their logarithms are written now as a sum over the roots of the factor of the denominator they belong to, sum(A(r)/D'(r) ln(x - r), r in { r : E(r) = 0 }), the form #1285 asked for, with the node #1624 added. The poles whose residues are roots of R are the common roots of D and D'^n R(A/D'), so E is their gcd over the rationals, and no arithmetic in the residues' field is needed. Only as a last resort. Asked beside the other rules, the sum answered each of Jeffrey's three terms of (-1 + 4 cos(x) + 5 cos(x)^2)/(-1 - 4 cos(x) - 3 cos(x)^2 + 4 cos(x)^3) with a sum over the roots of a sextic, where the whole is 2 arctan(tan(x/2)^3 - 2 tan(x/2)). So the rule declines as before and notes that a sum would answer, and only where the question comes back unanswered and that was noted is it asked again with sums allowed, under a memo of its own; the answer found then is remembered for the question. The integrator checks an answer in double intervals and then in decimals, and neither read the sum node, so every such answer would have been declined by its own check. PreciseEvaluation encloses each root now: Durand-Kerner's approximations, each widened by its Weierstrass correction in intervals, n |w_k| about z_k, which holds exactly one root where the rectangles are apart (Braess and Hadeler), and the summand worked out on each rectangle. Measured on the Rubi corpus, master at ad05a79b and this change on it, run side by side: the independent suites 1753 -> 1756 of 1814 (Bondarenko 22 and 23, and Hearn 66, 1/(1 - x^4 + x^8)); family 1 at ten problems a file, 295 -> 301 of 377 (1/(x^2 (1 - x^3 + x^6)), (1 + x^4)/(1 - x^4 + x^8), (1 + x^4)/(1 - 5x^4 + x^8), (-1 + 2x^4 + sqrt(3))/(1 - x^4 + x^8), (1 - x^4)/(1 - x^4 + x^8) and x^3 (5 + x + 3x^2 + 2x^3)/(2 + x + 5x^2 + x^3 + 2x^4)); families 2 to 7 at the usual sample, 663 of 722 on both. 0 wrong everywhere. Run alone, those nine are 0 of 9 on master and 9 of 9 here, and the tenth problem that moved, 1.1.4.3 #231, is 17.9 s unevaluated on both: its master timeout was the load of the side-by-side run. A question still declined costs what it did: the problems unevaluated in both arms took 14.2 s against 13.4 s in the independent suites, 65.6 s against 64.7 s in families 2 to 7 and 89.9 s against 92.3 s in the family 1 sample. The performance gate passes on 34c548a8, whose library code is this commit's: allocation is what the baseline says on all 19 gated benchmarks. RothsteinTragerTest: the cubic, quartic and quintic rows that were pinned as declined are sums now, and #1285's integral in PowerSubstitutionIntegralTest; that test's declined witnesses are the elliptic x^2/sqrt(x^4 + x + 1) and x^2/sqrt(x^4 + x^3 + 1), which nothing answers. BinomialDenominatorIntegralTest, PartialFractionsTest and RationalIntegralsTest pinned 1/(x^3 + x + 1), 1/(x^3 + x^2 + x + 2), 1/(x^4 + x + 1) and 1/(x^4 + x^3 + 1) as declined: each is the sum now, checked by differentiating back, and the boundary pinned in their place is the same kind of denominator with a symbol among its coefficients, which the sum is not written for. DecliningStaysCheap keeps its ten-second bound as FindingThatNothingSplitsStaysCheap: the second pass starts only once the first has declined, so the decline it guards is still inside the bound. Part of #1285. Co-Authored-By: Claude Opus 5.5 --- BREAKING-CHANGES.md | 25 +++++ .../Integration/Integration.Definition.cs | 48 ++++++++- .../Continuous/Integration/RothsteinTrager.cs | 99 ++++++++++++++++--- .../AngouriMath/Numerics/PreciseEvaluation.cs | 97 ++++++++++++++++++ .../BinomialDenominatorIntegralTest.cs | 17 +++- .../Calculus/PartialFractionsTest.cs | 37 ++++--- .../Calculus/PowerSubstitutionIntegralTest.cs | 34 +++---- .../Calculus/RationalIntegralsTest.cs | 20 +++- .../UnitTests/Calculus/RothsteinTragerTest.cs | 18 +++- 9 files changed, 329 insertions(+), 66 deletions(-) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 841ea252a..7cf40b66e 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1699,6 +1699,31 @@ where it is worked out, is not affected. | `"sum(1/(k - x), k, 1, n)".ToEntity().DomainCondition` | `not k - x = 0` | `forall k in { k : k in ZZ and 1 <= k and k <= n } : not k - x = 0` | | `"integral(1/(t - x), t, 0, 1)".ToEntity().DomainCondition` | `not t - x = 0` | `forall t in [0; 1] : not t - x = 0` | +### A rational function whose residues are of degree three or more is integrated over its poles + +**Wider.** The logarithmic part of a rational integral is a logarithm at each pole, weighted by the +residue there. Where the residues are rational or roots of a quadratic they were written out, as +logarithms and arctangents. Where they lie in a field of degree three or more, the rational +integrator declined, and `1/(x^3 + x + 1)` was left unevaluated. It is written now as a sum over the +poles, `sum(A(r)/D'(r) ln(x - r), r in { r : E(r) = 0 })`, with the node above: `E` is the factor +of the denominator whose poles have those residues, the gcd of `D` and `D'^n R(A/D')` over the +rationals. The integrand of [#1285](https://github.com/asc-community/AngouriMath/issues/1285), +`sqrt(x)/(1 + x + x^4)`, is `2u^2/(1 + u^2 + u^8)` under `u = sqrt(x)`, and is answered through it. + +**Only where nothing else answers.** The integral is asked again with sums over roots allowed once +the whole of it has come back unevaluated, so an integrand another rule writes in closed form keeps +that form. Asked with the other rules, the sum answered each of Jeffrey's three terms of +`(-1 + 4 cos(x) + 5 cos(x)^2)/(-1 - 4 cos(x) - 3 cos(x)^2 + 4 cos(x)^3)` with a sum over the roots of +a sextic, where the whole is one arctangent. A denominator with a symbol among its coefficients is +still declined. Both columns measured on a build, `v2.5.0` against this change. + +| | Was (2.5.0) | Is | +|---|---|---| +| `"1/(x^3 + x + 1)".Integrate("x")` | `integral(1 / (x ^ 3 + x + 1), x)` | `sum(1 / (3 * r ^ 2 + 1) * ln(x - r), r in { r : r ^ 3 + r + 1 = 0 }) + C` | +| `"sqrt(x)/(1 + x + x^4)".Integrate("x")` | `integral(sqrt(x) / (1 + x + x ^ 4), x)` | `sum(2 * r ^ 2 / (8 * r ^ 7 + 2 * r) * ln(x ^ (1/2) - r), r in { r : r ^ 8 + r ^ 2 + 1 = 0 }) + C` | +| `"1/(1 - x^4 + x^8)".Integrate("x")` | `integral(1 / (1 - x ^ 4 + x ^ 8), x)` | `sum(1 / (8 * r ^ 7 + (-4) * r ^ 3) * ln(x - r), r in { r : r ^ 8 + -r ^ 4 + 1 = 0 }) + C` | +| `"1/(x^4 + a*x + 1)".Integrate("x")` | `integral(1 / (x ^ 4 + a * x + 1), x)` | the same | + ### `(a + b asech(c x))/(d + e x)^2` is integrated, and a root written apart with `|x|` no longer needs a parity Rubi's 7.5.1 with a symbolic linear below the bar ran for ten minutes without an answer. By diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 92dc8b563..a2770fba7 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -446,14 +446,54 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => descentTruncated = true; return null; } - if (descentDepth == 0) + if (descentDepth != 0) + return ComputeIndefiniteIntegralGuarded(expr, x, integrateByParts); + descentTruncated = false; + inProgress?.Clear(); + ASumOverRootsWouldAnswer = false; + var answer = ComputeIndefiniteIntegralGuarded(expr, x, integrateByParts); + if (answer is not null || !ASumOverRootsWouldAnswer || SumsOverRootsAllowed) + return answer; + // Once more with sums over roots, and with a memo of its own, since the first pass + // remembered every one of those rational functions as declined. Last, and only where + // nothing else answered: asked with the rest, a sum over roots answered each of + // Jeffrey's three terms, where the whole is one arctangent the rational integrator + // writes when it is asked the whole. + // https://github.com/asc-community/AngouriMath/issues/1285 + var (memo, stamp) = (answered, answeredUnder); + SumsOverRootsAllowed = true; + answered = null; + descentTruncated = false; + inProgress?.Clear(); + try { - descentTruncated = false; - inProgress?.Clear(); + answer = ComputeIndefiniteIntegralGuarded(expr, x, integrateByParts); } - return ComputeIndefiniteIntegralGuarded(expr, x, integrateByParts); + finally + { + SumsOverRootsAllowed = false; + (answered, answeredUnder) = (memo, stamp); + } + // In place of the decline the first pass remembered, which a question asked again + // would otherwise be answered with, without the rule ever being reached to say that a + // second pass could answer it. + if (answer is not null && memo is not null && stamp is not null && SettingsState.StillHolds(stamp)) + memo[(Normalized(expr, x), x, integrateByParts, true)] = answer; + return answer; } + /// + /// Whether the rational integrator may answer with a sum over the roots of a factor of the + /// denominator: only in the second pass of a question nothing else answered. + /// + [System.ThreadStatic] internal static bool SumsOverRootsAllowed; + + /// + /// Set where the rational integrator declined only because the residues lie in a field of + /// degree above two, which a sum over roots would have written. + /// + [System.ThreadStatic] internal static bool ASumOverRootsWouldAnswer; + /// /// past its entry check: the cycle guard, the /// depth, and the memo. diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/RothsteinTrager.cs b/Sources/AngouriMath/Functions/Continuous/Integration/RothsteinTrager.cs index cfefb6a75..6cd535791 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/RothsteinTrager.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/RothsteinTrager.cs @@ -39,8 +39,10 @@ namespace AngouriMath.Functions.Algebra /// a ± w with w^2 rational, and the gcd is computed in Q(w): for /// w real the pair is two logarithms with a root in them, and for w = i b it /// is a ln(P^2 + b^2 Q^2) + b LogToAtan(P, b Q), Rioboo's continuous arctangent - /// form. A factor of higher degree would put the residues in a field this does not do - /// arithmetic in, and is declined. + /// form. A factor of higher degree puts the residues in a field this does no arithmetic in, + /// and those logarithms are written as a sum over the roots of the factor of the denominator + /// they belong to, sum(A(r)/D'(r) ln(x - r), r in { r : E(r) = 0 }) + /// (https://github.com/asc-community/AngouriMath/issues/1285). /// /// /// Bronstein, Symbolic Integration I, §2.2 (HermiteReduce), §2.5 @@ -64,8 +66,8 @@ internal static class RothsteinTrager /// /// The integral of over , /// both polynomials in with rational coefficients and the fraction - /// proper; where they are not, or where a residue lies in a field - /// of degree above two. + /// proper; where they are not, or where the answer does not + /// differentiate back to the integrand at the sampled points. /// internal static Entity? Integrate(Entity numerator, Entity denominator, Variable x) { @@ -148,23 +150,54 @@ private static (Entity Rational, RationalPolynomial Above, RationalPolynomial Be if (SquareFree(resultantPolynomial) is not { } residueParts) return null; - // Every residue field first, and the residues only then: a residue in a field of - // degree above two declines the whole, and the conjugate pair of a quadratic factor - // beside it was twelve seconds of gcds over the extension before the quartic factor - // was reached and declined it. - var factored = new List<(IReadOnlyList Factors, int Multiplicity)>(); + // Every residue field first, and the residues only then. A residue in a field of degree + // above two is written as a sum over roots, below, and where there is one, every + // residue that is not rational goes into that sum with it: the conjugate pair of a + // quadratic factor beside a quartic one was twelve seconds of gcds over the extension, + // and a sum over roots needs no arithmetic in any extension at all. + var factors = new List<(IReadOnlyList Factors, int Multiplicity)>(); foreach (var part in residueParts) { - var factors = PolynomialFactorization.FactorPrimitive(ToInteger(part.Factor).PrimitivePart()); - if (factors is null || factors.Any(irreducible => irreducible.Factor.Degree > 2)) + if (PolynomialFactorization.FactorPrimitive(ToInteger(part.Factor).PrimitivePart()) is not { } irreducibles) return null; - factored.Add((factors, part.Multiplicity)); + factors.Add((irreducibles, part.Multiplicity)); + } + var overRoots = factors.Any(part => part.Factors.Any(irreducible => irreducible.Factor.Degree > 2)); + // Only once nothing else has answered the question: see + // Integration.SumsOverRootsAllowed. Until then this declines, as it always did. + if (overRoots && !Integration.SumsOverRootsAllowed) + { + Integration.ASumOverRootsWouldAnswer = true; + return null; + } + var factored = new List<(IReadOnlyList Factors, int Multiplicity)>(); + var summed = RationalPolynomial.One; + var summedRoots = 0; + foreach (var (irreducibles, multiplicity) in factors) + { + if (!overRoots) + { + factored.Add((irreducibles, multiplicity)); + continue; + } + factored.Add((irreducibles.Where(irreducible => irreducible.Factor.Degree == 1).ToList(), multiplicity)); + foreach (var irreducible in irreducibles.Where(irreducible => irreducible.Factor.Degree > 1)) + { + summed = summed.Multiply(RationalPolynomial.FromInteger(irreducible.Factor)); + summedRoots += irreducible.Factor.Degree * multiplicity; + } } Entity total = Integer.Create(0); - foreach (var (factors, multiplicity) in factored) + if (overRoots) { - var part = (Multiplicity: multiplicity, Factors: factors); - foreach (var irreducible in factors) + if (OverTheRoots(summed, summedRoots, above, below, derivative, x) is not { } sum) + return null; + total = sum; + } + foreach (var (kept, multiplicity) in factored) + { + var part = (Multiplicity: multiplicity, Factors: kept); + foreach (var irreducible in kept) { var f = irreducible.Factor; Entity? term = f.Degree switch @@ -181,6 +214,42 @@ private static (Entity Rational, RationalPolynomial Above, RationalPolynomial Be return total; } + /// + /// The logarithms at the poles whose residues are roots of , + /// as a sum over those poles: sum(A(r)/D'(r) ln(x - r), r in { r : E(r) = 0 }), + /// where E is the factor of D they are the roots of. + /// where E does not have the the residues' multiplicities + /// count. + /// + /// + /// The residue at a simple pole b of A/D is A(b)/D'(b), a root of the + /// resultant, and D'(b) is not zero there. So the poles whose residues are roots of + /// R are the common roots of D and the polynomial D'^n R(A/D'), and + /// E is their gcd over the rationals: no arithmetic in the field the residues lie in + /// is needed. The sum is the one Rothstein–Trager's theorem writes over the roots of the + /// resultant, with its terms taken pole by pole. E has a factor of degree above two, + /// since a pole in a field of degree at most two has its residue there too, so the sum is + /// left standing rather than written out in radicals, unless each such factor is a + /// binomial, whose roots are written. + /// https://github.com/asc-community/AngouriMath/issues/1285 + /// + private static Entity? OverTheRoots(RationalPolynomial residues, int roots, + RationalPolynomial above, RationalPolynomial below, RationalPolynomial derivative, Variable x) + { + // D'^n R(A/D') = sum over k of r_k A^k D'^(n - k). + var n = residues.Degree; + var cleared = RationalPolynomial.Zero; + for (var k = 0; k <= n; k++) + cleared = cleared.Add(above.Pow(k).Multiply(derivative.Pow(n - k)).ScaleBy(residues[k])); + var poles = RationalPolynomial.FromInteger(IntegerPolynomial.Gcd(ToInteger(below), ToInteger(cleared)).PrimitivePart()); + if (poles.Degree != roots) + return null; + // The polynomials are in x alone, so any other name is free. + var r = MathS.Var(x.Name == "r" ? "w" : "r"); + var summand = above.ToEntity(r) / derivative.ToEntity(r) * MathS.Ln(x - r); + return MathS.Sum(summand, r, new Set.ConditionalSet(r, poles.ToEntity(r).Equalizes(Integer.Zero))); + } + /// c ln(gcd(D, A - c D')) for a rational residue . private static Entity? RationalResidue( ERational c, int degree, RationalPolynomial above, RationalPolynomial below, RationalPolynomial derivative, Variable x) diff --git a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs index 6261c61c7..262a4e4e4 100644 --- a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs +++ b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs @@ -7,6 +7,7 @@ using System; using System.Collections.Generic; +using System.Linq; using PeterO.Numbers; using static AngouriMath.Entity; @@ -76,6 +77,10 @@ internal sealed class PreciseEvaluation // every subtree free of the point, and each of those is worked out once. private readonly Dictionary known = new(ByReference.Instance); + // The roots a sum over roots is working through, each under a name of its own and as the + // rectangle that encloses it. + private Dictionary? bound; + // One set of contexts per precision, for the life of the process: the library's constant // cache and its guard-digit contexts are keyed by the context instance, so contexts // built afresh for every evaluation had pi and the logarithm tables recomputed each time. @@ -395,6 +400,10 @@ private PreciseInterval E() private PreciseComplexInterval Evaluate(Entity expr) { + // Ahead of the values remembered by node: a name bound to one root now is bound to + // another the next time round. + if (bound is not null && expr is Variable name && bound.TryGetValue(name, out var root)) + return root; if (known.TryGetValue(expr, out var value)) return value; value = EvaluateNode(expr); @@ -521,11 +530,99 @@ private PreciseComplexInterval EvaluateNode(Entity expr) return BySlope(Evaluate(argument), Number.Shi, z => Divide(Half(Subtract(Exp(z), Exp(Negate(z)))), z), cutEnd: null); case Chif(var argument): return BySlope(Evaluate(argument), Number.Chi, z => Divide(Half(Add(Exp(z), Exp(Negate(z)))), z), cutEnd: EDecimal.Zero); + // A sum over the roots of a polynomial, the form Rothstein–Trager answers in: each + // root enclosed, the summand worked out on its rectangle, and the terms added. + // https://github.com/asc-community/AngouriMath/issues/1285 + case SumOverSetf(var summand, Variable name, Set.ConditionalSet { Var: Variable w, Predicate: Equalsf(var left, var right) }): + return OverTheRoots(summand, name, left - right, w); default: return undefined; } } + /// + /// added up over the roots of in + /// , taken by in turn; undefined where a root + /// cannot be enclosed apart from the others. + /// + private PreciseComplexInterval OverTheRoots(Entity summand, Variable name, Entity polynomial, Variable w) + { + if (RootsEnclosed(polynomial, w) is not { } roots) + return undefined; + var total = Real(PreciseInterval.Exactly(EDecimal.Zero)); + foreach (var root in roots) + { + bound ??= new(); + var own = Variable.CreateTemp(summand.Vars.Concat(bound.Keys).Append(name)); + bound[own] = root; + total = Add(total, Evaluate(summand.Substitute(name, own))); + } + return total; + } + + /// + /// A rectangle about each root of in , + /// holding that root and no other; null where the polynomial does not have rational + /// coefficients, or the rectangles cannot be kept apart. + /// + /// + /// Durand–Kerner's approximations z_k at this precision, each enclosed by its + /// Weierstrass correction w_k = p(z_k)/(a prod_{j != k} (z_k - z_j)), worked out in + /// intervals. Every root lies in a disk about some z_k - w_k of radius + /// (n - 1)|w_k|, and a disk apart from the others holds exactly one (Braess and + /// Hadeler). That disk is inside the one of radius n|w_k| about z_k, whose + /// square is the rectangle, so rectangles apart from each other hold one root each. + /// + private List? RootsEnclosed(Entity polynomial, Variable w) + { + if (AngouriMath.Functions.SumOverSet.SquareFreeParts(polynomial, w) is not { } parts) + return null; + var enclosed = new List(); + foreach (var part in parts) + { + var p = part.Factor; + var n = p.Degree; + if (n < 1) + continue; + if (AngouriMath.Functions.Algebra.NumericalSolving.DurandKerner.Roots(p, near) is not { } approximations + || approximations.Count != n) + return null; + var points = new PreciseComplexInterval[n]; + for (var k = 0; k < n; k++) + points[k] = new(PreciseInterval.Exactly(approximations[k].RealPart.EDecimal), + PreciseInterval.Exactly(approximations[k].ImaginaryPart.EDecimal)); + var lead = Real(PreciseInterval.Exactly(EDecimal.FromEInteger(p[n]))); + var rectangles = new PreciseComplexInterval[n]; + for (var k = 0; k < n; k++) + { + var value = lead; + for (var i = n - 1; i >= 0; i--) + value = Add(Multiply(value, points[k]), Real(PreciseInterval.Exactly(EDecimal.FromEInteger(p[i])))); + var product = lead; + for (var j = 0; j < n; j++) + if (j != k) + product = Multiply(product, Subtract(points[k], points[j])); + var correction = Divide(value, product); + if (!correction.IsFinite) + return null; + var radius = EDecimal.FromInt32(n).Multiply(correction.Re.Magnitude.Add(correction.Im.Magnitude, up), up); + var re = points[k].Re.Low; + var im = points[k].Im.Low; + rectangles[k] = new(new(re.Subtract(radius, down), re.Add(radius, up)), new(im.Subtract(radius, down), im.Add(radius, up))); + } + for (var j = 0; j < n; j++) + for (var k = j + 1; k < n; k++) + if (Overlap(rectangles[j], rectangles[k])) + return null; + enclosed.AddRange(rectangles); + } + return enclosed; + } + + private static bool Overlap(PreciseComplexInterval a, PreciseComplexInterval b) + => a.Re.Low.CompareTo(b.Re.High) <= 0 && b.Re.Low.CompareTo(a.Re.High) <= 0 + && a.Im.Low.CompareTo(b.Im.High) <= 0 && b.Im.Low.CompareTo(a.Im.High) <= 0; + /// /// A special function over the rectangle: its value at the middle, which the library /// works out with ten digits to spare, widened by the most the function can move within diff --git a/Sources/Tests/UnitTests/Calculus/BinomialDenominatorIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/BinomialDenominatorIntegralTest.cs index 79186ef82..4d33c2ed7 100644 --- a/Sources/Tests/UnitTests/Calculus/BinomialDenominatorIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/BinomialDenominatorIntegralTest.cs @@ -128,14 +128,21 @@ public void AnExactFactorisationIsStillTakenFirst(string integrand, string expec } /// - /// What is not this shape and must stay declined by it rather than mis-answered: a - /// trinomial with no rational root, which no rule reads and - /// #1285 is about. + /// What is not this shape is not this rule's to claim: a trinomial with no rational root. + /// It is answered by the sum over its roots that + /// #1285 is about, + /// so the answer is that sum rather than this rule's cosines of roots of unity, and it + /// differentiates back. /// [Theory] [InlineData("1/(x^3 + x + 1)")] - public void OutsideTheShapeNothingIsClaimed(string integrand) - => Assert.Contains("integral(", integrand.ToEntity().Integrate("x").Stringize()); + public void OutsideTheShapeTheSumOverTheRootsAnswers(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.Contains(integral.Nodes, node => node is Entity.SumOverSetf); + Assert.DoesNotContain("cos(", integral.Stringize()); + DifferentiatesBack(integrand); + } /// /// A repeated binomial is the Hermite reduction's first, and what that leaves over the diff --git a/Sources/Tests/UnitTests/Calculus/PartialFractionsTest.cs b/Sources/Tests/UnitTests/Calculus/PartialFractionsTest.cs index 845d55173..cf3a323b0 100644 --- a/Sources/Tests/UnitTests/Calculus/PartialFractionsTest.cs +++ b/Sources/Tests/UnitTests/Calculus/PartialFractionsTest.cs @@ -140,28 +140,34 @@ public void ABiquadraticWithTwoRealRootsInTheSquare(string integrand, double[] p AssertIsAntiderivative(integrand, points); /// - /// What is still left unevaluated rather than answered wrongly: a denominator that is - /// irreducible and not biquadratic. x^2/(x^4 + 1) used to be on this list and is - /// now answered above; the step over the reals reaches a biquadratic only, so a - /// quartic with an odd power in it stays here. 1/(x^4 + 2x^2 + 1) was here too, - /// as a power of a single irreducible with no coprime pair to split into -- it is - /// (x^2 + 1)^2, which the Hermite reduction answers once the repeated factor is - /// written, and the denominator is now written that way first; see - /// RationalIntegralsTest.ARepeatedFactorTheSpellingHides. + /// What has nothing to split into: a denominator that is irreducible and not biquadratic. + /// The step over the reals reaches a biquadratic only, so a quartic with an odd power in + /// it is answered whole, its logarithms in a sum over its roots + /// (#1285). + /// x^2/(x^4 + 1) used to be on this list and is answered above by that step. + /// 1/(x^4 + 2x^2 + 1) was here too, as a power of a single irreducible with no + /// coprime pair to split into -- it is (x^2 + 1)^2, which the Hermite reduction + /// answers once the repeated factor is written, and the denominator is now written that + /// way first; see RationalIntegralsTest.ARepeatedFactorTheSpellingHides. /// [Theory] [InlineData("1 / (x ^ 3 + x ^ 2 + x + 2)")] [InlineData("1 / (x ^ 4 + x ^ 3 + 1)")] [InlineData("1 / (x ^ 4 + x + 1)")] - public void WhatCannotBeSplitIsLeftAlone(string integrand) => - Assert.Contains("integral(", integrand.ToEntity().Integrate("x").Stringize()); + public void WhatCannotBeSplitIsASumOverItsRoots(string integrand) + { + Assert.Contains(integrand.ToEntity().Integrate("x").Nodes, node => node is Entity.SumOverSetf); + AssertIsAntiderivative(integrand, 0.3, 1.7, 3.2, -2.4); + } /// - /// The guard that keeps declining cheap, which nothing widening the split must undo: + /// The guard that keeps the search cheap, which nothing widening the split must undo: /// this factorises into x^4 + x + 1, an irreducible quartic with an odd power in - /// it that no rule reads, and the whole point of reading the factorisation is that + /// it that no split reads, and the whole point of reading the factorisation is that /// finding that out costs one factorisation rather than a search of every half of every - /// split. + /// split. The quartic's roots do answer it, but in a second pass that starts only once + /// the first has declined, so the first pass's decline is still inside the time asserted + /// here. /// /// /// The integrand here used to be (1 - x^4)/(1 + x^4 + x^8), whose quartic is @@ -171,12 +177,13 @@ public void WhatCannotBeSplitIsLeftAlone(string integrand) => /// that actually reaches it. /// [Fact] - public void DecliningStaysCheap() + public void FindingThatNothingSplitsStaysCheap() { var clock = System.Diagnostics.Stopwatch.StartNew(); var answer = "(1 - x ^ 4) / ((1 + x ^ 2) * (x ^ 4 + x + 1))".ToEntity().Integrate("x"); - Assert.Contains("integral(", answer.Stringize()); Assert.True(clock.Elapsed < System.TimeSpan.FromSeconds(10), $"took {clock.Elapsed}"); + Assert.Contains(answer.Nodes, node => node is Entity.SumOverSetf); + AssertIsAntiderivative("(1 - x ^ 4) / ((1 + x ^ 2) * (x ^ 4 + x + 1))", 0.3, 1.7, 3.2, -2.4); } /// diff --git a/Sources/Tests/UnitTests/Calculus/PowerSubstitutionIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/PowerSubstitutionIntegralTest.cs index 01ffe674a..b412d93c3 100644 --- a/Sources/Tests/UnitTests/Calculus/PowerSubstitutionIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/PowerSubstitutionIntegralTest.cs @@ -99,15 +99,12 @@ public void AFractionalPowerIsAlsoASubstitution(string integrand) => DifferentiatesBack(integrand); /// - /// Where a fractional substitution reaches and the rest of the chain does not, as a - /// note rather than a pin: sqrt(x)/(1 + x + x^4) becomes - /// 2u^2/(1 + u^2 + u^8), whose denominator is irreducible over the rationals, - /// not a biquadratic and not a binomial. Its antiderivative is a sum over the eight - /// roots of w^8 + w^2 + 1 of w ln(sqrt(x) - w)/(4 w^6 + 1), which this - /// library has no node to write, and that is - /// #1285 rather - /// than a verdict to record here — a test that pins the decline would have to be - /// falsified to close the issue. + /// Where a fractional substitution reaches and the rational integrator then answers with + /// a sum over roots: sqrt(x)/(1 + x + x^4) becomes 2u^2/(1 + u^2 + u^8), + /// whose denominator is irreducible over the rationals, not a biquadratic and not a + /// binomial. Its antiderivative is a sum over the eight roots of w^8 + w^2 + 1 of + /// w ln(sqrt(x) - w)/(4 w^6 + 1), the form + /// #1285 asked for. /// /// /// sqrt(x)/(x + 1) and sqrt(x)/(1 + x^4) were pinned here in turn, for an @@ -115,13 +112,10 @@ public void AFractionalPowerIsAlsoASubstitution(string integrand) /// moved up into the theory above when its rule arrived. /// [Fact] - public void WhereTheChainStopsIsAnIssueNotAPin() + public void TheSubstitutionEndsInASumOverRoots() { - var integral = "sqrt(x)/(1 + x + x^4)".ToEntity().Integrate("x"); - // Either answer is acceptable here: unevaluated today, and a correct antiderivative - // once #1285 gives it a form. What is not acceptable is a wrong one. - if (!integral.Stringize().Contains("integral(")) - DifferentiatesBack("sqrt(x)/(1 + x + x^4)"); + DifferentiatesBack("sqrt(x)/(1 + x + x^4)"); + Assert.Contains("sqrt(x)/(1 + x + x^4)".ToEntity().Integrate("x").Nodes, node => node is Entity.SumOverSetf); } /// @@ -182,12 +176,16 @@ public void WhatIntegratedBeforeStillDoes(string integrand) /// step learning to factor a biquadratic denominator over the reals. A test that the /// substitution declines something needs an integrand nothing else answers either, or it /// stops testing the substitution the moment a neighbouring capability arrives. These - /// carry an odd power, which puts them out of reach of that step as well. + /// carry an odd power, which puts them out of reach of that step as well. Then + /// x^2/(x^4 + x + 1) and x^2/(x^4 + x^3 + 1) were answered too, as sums over + /// the roots of their denominators + /// (), and the + /// witnesses are their square roots now, elliptic integrals nothing answers. /// . /// [Theory] - [InlineData("x^2/(x^4 + x + 1)")] - [InlineData("x^2/(x^4 + x^3 + 1)")] + [InlineData("x^2/sqrt(x^4 + x + 1)")] + [InlineData("x^2/sqrt(x^4 + x^3 + 1)")] public void AnIntegralThisDoesNotReachIsStillDeclined(string integrand) => Assert.Contains("integral(", integrand.ToEntity().Integrate("x").Stringize()); diff --git a/Sources/Tests/UnitTests/Calculus/RationalIntegralsTest.cs b/Sources/Tests/UnitTests/Calculus/RationalIntegralsTest.cs index 159bd5853..78930012f 100644 --- a/Sources/Tests/UnitTests/Calculus/RationalIntegralsTest.cs +++ b/Sources/Tests/UnitTests/Calculus/RationalIntegralsTest.cs @@ -209,9 +209,12 @@ public void AHighPowerOfAWrittenBase(string integrand, double[] points) } } - // What is out of reach is a denominator that does not factor over Q and is not a - // biquadratic either -- an odd power puts it past the step that factors over the reals. - // Recorded so the boundary is visible rather than inferred from an absence. + // A denominator that does not factor over Q and is not a biquadratic either -- an odd + // power puts it past the step that factors over the reals -- has its logarithms in a sum + // over its roots (https://github.com/asc-community/AngouriMath/issues/1285). What is out + // of reach is the same denominator with a symbol among its coefficients, which the sum + // is not written for. Recorded so the boundary is visible rather than inferred from an + // absence. // // x^2/(x^4 + 1) was the first entry here, on the grounds that x^4 + 1 is irreducible // over Q and only factors once real coefficients are allowed. Allowing them is what the @@ -221,7 +224,16 @@ public void AHighPowerOfAWrittenBase(string integrand, double[] points) [Theory] [InlineData("1 / (x ^ 4 + x + 1)")] [InlineData("1 / (x ^ 4 + x ^ 3 + 1)")] - public void ADenominatorThatDoesNotFactorIsStillDeclined(string integrand) => + public void ADenominatorThatDoesNotFactorIsASumOverItsRoots(string integrand) + { + Assert.Contains(integrand.ToEntity().Integrate("x").Nodes, node => node is Entity.SumOverSetf); + AssertIsAntiderivative(integrand, 0.3, 1.7, 3.2, -2.4); + } + + [Theory] + [InlineData("1 / (x ^ 4 + a * x + 1)")] + [InlineData("1 / (x ^ 3 + x + a)")] + public void ADenominatorWithASymbolThatDoesNotFactorIsStillDeclined(string integrand) => Assert.Contains("integral(", integrand.ToEntity().Integrate("x").Stringize()); } } diff --git a/Sources/Tests/UnitTests/Calculus/RothsteinTragerTest.cs b/Sources/Tests/UnitTests/Calculus/RothsteinTragerTest.cs index 13575b242..2ae06ec6b 100644 --- a/Sources/Tests/UnitTests/Calculus/RothsteinTragerTest.cs +++ b/Sources/Tests/UnitTests/Calculus/RothsteinTragerTest.cs @@ -34,7 +34,7 @@ namespace AngouriMath.Tests.Calculus [Trait("Area", "Calculus")] public sealed class RothsteinTragerTest { - private static void DifferentiatesBack(string integrand, double[] points) + private static Entity DifferentiatesBack(string integrand, double[] points) { var integral = integrand.ToEntity().Integrate("x"); Assert.DoesNotContain("integral(", integral.Stringize()); @@ -59,6 +59,7 @@ private static void DifferentiatesBack(string integrand, double[] points) } Assert.True(compared >= 5, $"only {compared} of {points.Length} points were comparable for {integrand}"); + return integral; } private static readonly double[] TheLine = { -2.3, -1.1, -0.4, 0.3, 0.9, 1.7, 2.6 }; @@ -103,13 +104,20 @@ public void TheArctangentsAreContinuous() } /// - /// A residue in a field of degree three is declined rather than answered in a field - /// this does no arithmetic in, and quickly. + /// A residue in a field of degree three or more puts its logarithms in a sum over the poles + /// it belongs to, sum(A(r)/D'(r) ln(x - r), r in { r : E(r) = 0 }), where these + /// were declined. x^5 + x + 1 is (x^2 + x + 1)(x^3 - x^2 + 1), and the + /// quadratic factor's residues go into the sum beside the cubic's. + /// #1285 /// [Theory] [InlineData("1/(x^3 + x + 1)")] + [InlineData("x/(x^3 - x + 1)")] [InlineData("(x^2 + 3)/(x^5 + x + 1)")] - public void DeclinedWhereTheResiduesAreCubic(string integrand) - => Assert.Contains("integral(", integrand.ToEntity().Integrate("x").Stringize()); + [InlineData("1/(x^5 - x + 1)")] + [InlineData("x^2/(x^4 + x + 1)")] + [InlineData("x^2/(x^4 + x^3 + 1)")] + public void AResidueOfDegreeThreeIsASumOverRoots(string integrand) + => Assert.Contains(DifferentiatesBack(integrand, TheLine).Nodes, node => node is Entity.SumOverSetf); } }