diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index b5268044d..24a085356 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -552,6 +552,23 @@ power of the variable every term holds comes out the same way. | `"x^2/((a*x + a)*(1 - x^2))".Integrate("x")` | unevaluated | the partial fractions over `(x + 1)^2 (1 - x)`, over `a` | | `"cot(x)^3/(a + b*csc(x))".Integrate("x")` | unevaluated | `-1/(b sin(x)) - a ln(sin(x))/b^2 + (b^2 - a^2) ln(a sin(x) + b)/(a b^2)`, once `a u^4 + b u^3` is read as `u^3 (a u + b)` | +### Symbolic factors of a denominator that share a factor are written over it + +**Improvement, not silent.** Partial fractions over factors with symbols in them read the written +factors as coprime and squarefree, and the refactoring after them reads rational coefficients only. +A quadratic that shares a root with a linear beside it, or that is a square, was read as an +irreducible quadratic and declined. The written factors are now taken apart over their greatest +common divisors, in every variable, with one another and each with its derivative, and a factor +shared with the numerator is cancelled. Rubi's 1.2.1.2 +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"1/((d+e*x)*(a*d*e+(c*d^2+a*e^2)*x+c*d*e*x^2)^2)".Integrate("x")` | unevaluated | the partial fractions over `(d + e x)^3 (a e + c d x)^2` | +| `"(d+e*x)^6/(a*d*e+(c*d^2+a*e^2)*x+c*d*e*x^2)^4".Integrate("x")` | no answer in 90 s | the antiderivative of `(d + e x)^2/(a e + c d x)^4` | +| `"1/((d+e*x)^2*(c*d^2+2*c*d*e*x+c*e^2*x^2)^3)".Integrate("x")` | unevaluated | `-1/(7 c^3 e^8 (d/e + x)^7)` | +| `"x^2/((a+b*x)*(a^2+2*a*b*x+b^2*x^2))".Integrate("x")` | unevaluated | the partial fractions over `(a + b x)^3` | + ### Fractional powers of `a ± a sin` are integrated by the half angle at which they are squares **Improvement, not silent.** `(A + C sin(y)^2)/((c - c sin(y))^(3/2) sqrt(a + a sin(y)))` and diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index c56b465f2..6dbddf974 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -665,6 +665,19 @@ is var (multiple, leftover) ?? Integration.ComputeIndefiniteIntegral(numerator / primitive, x, integrateByParts)) is { } overPrimitives) return overPrimitives; + // Written factors with a symbol in them that share a factor, or repeat one inside, + // written over what they share: `(d + e x)(a d e + (c d^2 + a e^2) x + c d e x^2)^2` + // is `(d + e x)^3 (a e + c d x)^2`, and `a^2 + 2 a b x + b^2 x^2` beside `a + b x` + // is a third power of it. The splits below read the written factors as coprime and + // squarefree, and the refactoring after them reads rational coefficients only: + // Rubi's 1.2.1.2 `1/((d + e x)(a d e + (c d^2 + a e^2) x + c d e x^2)^2)` was declined + // after four seconds, where written so it is answered in a tenth of one. Once: what + // this writes is coprime and squarefree already. + if (OverSymbolicFactorsWrittenApart(numerator, denominator, x) is var (above, below) + && (SolveByPartialFractions(above / below, x, integrateByParts) + ?? Integration.ComputeIndefiniteIntegral(above / below, x, integrateByParts)) is { } overWhatTheyShare) + return overWhatTheyShare; + // Written linear factors with symbols in their coefficients, two or more, one of // them to a power: decomposed over the written factors, the coefficients read // off derivatives at the roots and each a line. First, before the respellings @@ -10275,6 +10288,221 @@ private static bool IsASumOfMonomials(Entity expr, Entity.Variable x) return changed ? product : null; } + /// + /// over , the denominator's + /// written factors taken apart over their greatest common divisors -- with one another, + /// and each with its own derivative in -- until they are coprime and + /// squarefree as written, and what the numerator shares with them cancelled; + /// where no factor holds a symbol, they are coprime and squarefree + /// already, or one is not a polynomial. + /// + /// + /// The divisors are taken over the rationals in every variable, so a factor of the symbols + /// alone can come with one: that goes in front, as a constant, and each factor kept is + /// primitive in . `c d^2 + 2 c d e x + c e^2 x^2` shares `c (d + e x)` + /// with its derivative, and is `c (d + e x)^2`. + /// + private static (Entity Numerator, Entity Denominator)? OverSymbolicFactorsWrittenApart(Entity numerator, Entity denominator, Entity.Variable x) + { + var written = new List<(Entity Base, int Power)>(); + Entity constant = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(denominator)) + { + if (!factor.ContainsNode(x)) + { + constant = constant == Number.Integer.One ? factor : constant * factor; + continue; + } + var (@base, power) = factor is Powf(var b, Number.Integer e) && e.EInteger.Sign > 0 && e.EInteger.CanFitInInt32() + ? (b, e.EInteger.ToInt32Unchecked()) + : (factor, 1); + // Over one bar, the constant below it in front: the integrator writes a quadratic + // monic, `a/c + (c d^2 + a e^2)/(c d e) x + x^2`, and it is a polynomial in the + // symbols over `c d e`. + var (overOneBar, belowTheBar) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(@base)); + if (belowTheBar.ContainsNode(x)) + return null; + if (belowTheBar != Number.Integer.One) + { + var front = power == 1 ? belowTheBar : MathS.Pow(belowTheBar, power); + constant = constant == Number.Integer.One ? 1 / front : constant / front; + } + written.Add((overOneBar, power)); + } + if (!written.Any(pair => pair.Base.Vars.Any(v => v != x))) + return null; + var (numeratorOverOneBar, belowTheNumerator) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(numerator)); + if (!belowTheNumerator.ContainsNode(x)) + { + numerator = numeratorOverOneBar; + if (belowTheNumerator != Number.Integer.One) + constant = constant == Number.Integer.One ? belowTheNumerator : constant * belowTheNumerator; + } + var variables = numerator.Vars.Concat(denominator.Vars).Distinct().OrderBy(v => v.Name, System.StringComparer.Ordinal).ToList(); + if (variables.Count > MultivariatePolynomial.MaxVariables) + return null; + var indices = new Dictionary(); + for (var i = 0; i < variables.Count; i++) + indices[variables[i]] = i; + var all = Enumerable.Range(0, variables.Count).ToList(); + var at = indices[x]; + var rest = all.Where(i => i != at).ToList(); + var parts = new List<(MultivariatePolynomial Poly, int Power)>(); + var degree = 0; + foreach (var (@base, power) in written) + { + if (MultivariatePolynomial.TryParse(@base, indices) is not { } poly || poly.DegreeIn(at) == 0) + return null; + degree += poly.DegreeIn(at); + parts.Add((poly, power)); + } + if (degree > 12 || parts.Count < 2 && parts[0].Poly.DegreeIn(at) < 2) + return null; + + // A power of x is what a factor shares with x^k beside it, or with its own derivative + // where its lowest term is x^k, and taking it out is the content's, before this. Split + // here, `b x^2 + c x^4` came apart as `x^2 (b/c + x^2) c`, and over that monic quadratic + // with a symbol below its bar Rubi's `x^12 (A + B x^2)/(b x^2 + c x^4)^3` had an answer + // of four million characters, where with the content taken out it has seven thousand. + bool IsAPowerOfX(MultivariatePolynomial poly) + => poly.Monomials.Take(2).Count() == 1 + && rest.All(variable => MultivariatePolynomial.PowerOfMonomial(poly.Monomials.First(), variable) == 0); + + // What is free of x in a factor goes in front, to the factor's power. + void Front(MultivariatePolynomial free, int power) + { + if (free.IsConstant && free.SameAs(MultivariatePolynomial.One(free.VariableCount))) + return; + var front = free.ToEntity(variables); + var raised = power == 1 ? front : MathS.Pow(front, power); + constant = constant == Number.Integer.One ? raised : constant * raised; + } + bool Keep(MultivariatePolynomial poly, int power) + { + // A part free of x goes in front whole, its sign with it: a divisor is found up to + // a unit, so taking -(b u^2 - 1)^2 apart leaves a part -1, whose content is 1. + if (poly.DegreeIn(at) == 0) + { + Front(poly, power); + return true; + } + if (Functions.PolynomialGcd.ContentIn(poly, at, rest, 0) is not { } content || poly.DivideExact(content) is not { } primitive) + return false; + Front(content, power); + parts.Add((primitive, power)); + return true; + } + + var changed = false; + for (var round = 0; round < 64; round++) + { + var split = false; + for (var i = 0; i < parts.Count && !split; i++) + { + var (poly, power) = parts[i]; + // A factor repeated inside: what it shares with its derivative. + if (poly.DegreeIn(at) >= 2 && Functions.PolynomialGcd.Gcd(poly, poly.DerivativeIn(at), all, 0) is { } repeated + && repeated.DegreeIn(at) >= 1 && !IsAPowerOfX(repeated) && poly.DivideExact(repeated) is { } remaining) + { + parts.RemoveAt(i); + if (!Keep(repeated, power) || !Keep(remaining, power)) + return null; + split = true; + break; + } + // Two factors that share one. + for (var j = i + 1; j < parts.Count; j++) + { + if (Functions.PolynomialGcd.Gcd(poly, parts[j].Poly, all, 0) is not { } common || common.DegreeIn(at) == 0 || IsAPowerOfX(common)) + continue; + if (poly.DivideExact(common) is not { } left || parts[j].Poly.DivideExact(common) is not { } right) + return null; + var other = parts[j].Power; + parts.RemoveAt(j); + parts.RemoveAt(i); + if (!Keep(common, power + other) || !Keep(left, power) || !Keep(right, other)) + return null; + split = true; + break; + } + } + if (!split) + break; + changed = true; + } + + // A factor that shares one with the numerator, taken apart the same way: + // `(d + e x)^6/(a d e + (c d^2 + a e^2) x + c d e x^2)^4` is `(d + e x)^2/(a e + c d x)^4`. + var top = MultivariatePolynomial.TryParse(numerator, indices); + for (var round = 0; top is not null && round < 16; round++) + { + var split = false; + for (var i = 0; i < parts.Count; i++) + { + var (poly, power) = parts[i]; + if (poly.DegreeIn(at) < 2 || Functions.PolynomialGcd.Gcd(top, poly, all, 0) is not { } common + || common.DegreeIn(at) == 0 || IsAPowerOfX(common) || common.DegreeIn(at) == poly.DegreeIn(at) + || poly.DivideExact(common) is not { } remaining) + continue; + parts.RemoveAt(i); + if (!Keep(common, power) || !Keep(remaining, power)) + return null; + split = true; + break; + } + if (!split) + break; + changed = true; + } + if (!changed) + return null; + + // What the numerator shares with the factors, cancelled. + Entity above = numerator; + if (top is not null) + { + var cancelled = false; + for (var i = 0; i < parts.Count; i++) + while (parts[i].Power > 0 && top.DivideExact(parts[i].Poly) is { } fewer) + { + top = fewer; + parts[i] = (parts[i].Poly, parts[i].Power - 1); + cancelled = true; + } + if (cancelled) + above = top.ToEntity(variables); + } + // Each factor monic in x, the leading coefficient to its power in front, which is how + // the integrator writes a factor: `d + e x` is `e (d/e + x)`, and written the other + // way the next normalisation took `e` out again and the answer carried `e^(-8) e^8`. + var monicFactors = new List(); + foreach (var (poly, power) in parts) + { + if (power == 0) + continue; + var lead = poly.LeadingCoefficientIn(at).ToEntity(variables); + Entity monic = Number.Integer.Zero; + foreach (var pair in poly.CoefficientsIn(at).OrderBy(pair => pair.Key)) + { + var coefficient = pair.Key == poly.DegreeIn(at) ? Number.Integer.One : Reduced(pair.Value.ToEntity(variables) / lead); + var term = pair.Key == 0 ? coefficient : coefficient * (pair.Key == 1 ? x : MathS.Pow(x, pair.Key)); + monic = monic == Number.Integer.Zero ? term : monic + term; + } + monicFactors.Add(power == 1 ? monic : MathS.Pow(monic, power)); + var raised = power == 1 ? lead : MathS.Pow(lead, power); + constant = constant == Number.Integer.One ? raised : constant * raised; + } + // The constant in lowest terms: the bars and the contents taken out of the factors + // and the leading coefficients put back leave `e^(-8) e^8` in it as written. + Entity below = Reduced(constant); + foreach (var factor in monicFactors) + below = below == Number.Integer.One ? factor : below * factor; + return (above, below); + + static Entity Reduced(Entity constant) + => constant.Vars.Any() ? Functions.PartialFractions.InLowestTermsOverTheSymbols(constant) : constant.InnerSimplified; + } + /// /// with every written factor that has rational coefficients /// and splits over the rationals written in its factors, where another written factor diff --git a/Sources/Tests/UnitTests/Calculus/SymbolicFactorsPartialFractionsTest.cs b/Sources/Tests/UnitTests/Calculus/SymbolicFactorsPartialFractionsTest.cs index fed7dadbc..860e2f961 100644 --- a/Sources/Tests/UnitTests/Calculus/SymbolicFactorsPartialFractionsTest.cs +++ b/Sources/Tests/UnitTests/Calculus/SymbolicFactorsPartialFractionsTest.cs @@ -278,6 +278,35 @@ public void TheCoefficientsAreInLowestTerms(string integrand, int atMost) public void ProportionalLinearFactorsAreOneFactor(string integrand) => DifferentiatesBack(integrand, "x", ("a", 1.7), ("b", 2.3), ("c", 0.4), ("d", 1.1), ("g", 0.7)); + /// + /// Written factors that share a factor, or repeat one inside, are written over what they + /// share, the greatest common divisors taken in every variable: Rubi's 1.2.1.2, + /// (d + e x)^m (a d e + (c d^2 + a e^2) x + c d e x^2)^p, whose quadratic is + /// (d + e x)(a e + c d x), and the square c d^2 + 2 c d e x + c e^2 x^2 + /// beside a power of its linear, or with the numerator sharing it. Each was declined or + /// ran past the corpus's budget. + /// + [Theory] + [InlineData("1/((d + h*x)*(a*d*h + (c*d^2 + a*h^2)*x + c*d*h*x^2)^2)")] + [InlineData("(d + h*x)^6/(a*d*h + (c*d^2 + a*h^2)*x + c*d*h*x^2)^4")] + [InlineData("1/((d + h*x)^2*(c*d^2 + 2*c*d*h*x + c*h^2*x^2)^3)")] + [InlineData("x^2/((a + b*x)*(a^2 + 2*a*b*x + b^2*x^2))")] + public void FactorsThatShareOneAreWrittenOverIt(string integrand) + => DifferentiatesBack(integrand, "x", ("a", 1.7), ("b", 2.3), ("c", 0.4), ("d", 1.1), ("h", 0.7)); + + /// + /// The negative of a square, -(b x^2 - 1)^2, is the square with its sign in front: + /// a divisor is found up to a unit, and the part the split left was -1. Rubi's + /// 1.1.2.2 x^4/(a + b x^2)^(3/2) is this under t = x/sqrt(a + b x^2), and was + /// answered with the wrong sign. + /// + [Theory] + [InlineData("(2*x^2/b - 1/b^2)/(-b^2*x^4 + 2*b*x^2 - 1)")] + [InlineData("x^4/(a + b*x^2)^(3/2)")] + [InlineData("(c + d*x^3)^2/(a + b*x^3)^(4/3)")] + public void TheNegativeOfASquareKeepsItsSign(string integrand) + => DifferentiatesBack(integrand, "x", ("a", 1.7), ("b", 2.3), ("c", 0.4), ("d", 1.1)); + /// /// A written power of a linear is integrated as the power it is: expanded, /// (f/g + x)^2 is a quadratic whose discriminant is zero in a spelling the