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); } }