diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 2d8c7ac2b..44b28ce24 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -383,6 +383,23 @@ taken apart into pieces each closed the same way | `"(2 + x)/((2 + 4*x - 3*x^2)*(1 + 3*x + 2*x^2)^(3/2))".ToEntity().Integrate("x")` | `integral(...)` | an antiderivative | | `"1/((x^2 + 1)*sqrt(x^2 + x + 1))".ToEntity().Integrate("x")` | `integral(...)` | a logarithm and an arctangent | +### A polynomial over a power of a quadratic beside the root of another is integrated + +`P/(A^k sqrt(B))`, `A` and `B` two different quadratics and `k` at least two, was left +unevaluated. So were Rubi's `sqrt(a + a sec(x))/(c + d sec(x))^2` and its kin, which the +half-angle tangent writes as `(1 - t^2)^(3/2)/((1 + t^2)((c + d) + (d - c) t^2)^2)`. They are +reduced a power at a time: `L/(A^j sqrt(B))`, with `L` linear, is the derivative of +`(p + q x) sqrt(B)/A^(j - 1)` plus a quadratic over `A^(j - 1) sqrt(B)`, from five linear +equations. The reduction goes down to the closed form over `A` of the entry above. With symbols +in both quadratics, the cube's coefficients grow past what the reduction takes, and it declines +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(g + h*x)/((d + k*x + f*x^2)^2*sqrt(a + b*x + c*x^2))".ToEntity().Integrate("x")` | `integral(...)` | a linear times the root over the quadratic, and two arctangents | +| `"(7 + 13*x)/((5 + x + 2*x^2)^3*sqrt(2 + x + 3*x^2))".ToEntity().Integrate("x")` | `integral(...)` | an antiderivative | +| `"sqrt(a + a*sec(x))/(c + d*sec(x))^2".ToEntity().Integrate("x")` | `integral(...)` | an antiderivative in `tan(x/2)` | + ### A function comes out of a fractional power of its even power with its sign **Improvement, not silent.** The entry two above made `(sin(x)^2)^(3/2)` the modulus `|sin(x)|^3` diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 3abeafc79..8f809d452 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -2661,7 +2661,12 @@ Entity Step(ERational power) integrand = integrand.InnerSimplified; if (integrand.ContainsNode(x) || integrand.Nodes.Any(node => node == MathS.NaN)) return null; - if (Integration.ComputeAsAQuestionOfItsOwn(integrand, t, integrateByParts) is not { } result) + // A rational function beside the root, which is what this lands on, is taken apart by + // the rule for it before the chain is asked: asked as a question, the substitution + // search ahead of that rule spent the budget on `(1 - t^2)^(3/2)/((1 + t^2)(p + q t^2)^2)`, + // which `sqrt(a + a sec(y))/(c + d sec(y))^2` is. + if ((SolveARationalFunctionBesideTheRootOfAQuadratic(integrand, t) + ?? Integration.ComputeAsAQuestionOfItsOwn(integrand, t, integrateByParts)) is not { } result) return null; var back = secantKind ? MathS.Tan(argument / 2) : MathS.Tan(MathS.pi / 4 - argument / 2); var answer = (signs == Number.Integer.One ? result : signs * result).Substitute(t, back); @@ -11306,6 +11311,184 @@ Dictionary PowerOfQ(int j) return answer; } + /// + /// The largest the linear p + q x of one step of + /// is let grow, + /// its two coefficients' complexities together: past it the reduction declines rather than + /// hand the next step a numerator of tens of thousands of nodes. With six symbols the + /// cube of d + k x + f x^2 beside sqrt(a + b x + c x^2) reached 10,700 at + /// its second step and took over two minutes to answer and check; the square reaches a + /// few hundred, and Rubi's `sqrt(a + a sec(x))/(c + d sec(x))^3` 1,010. + /// + private const int LargestReducedCoefficients = 4000; + + /// + /// A polynomial over a power of a quadratic beside the square root of another quadratic, + /// P/(A^k sqrt(B)) with k at least two and P of degree below + /// 2k, one power of A at a time down to the first, which + /// closes. + /// + /// + /// + /// Over A^j, P is L + A R with L linear, and R goes a + /// power down as it is. For L = g + h x, the derivative of + /// (p + q x) sqrt(B)/A^(j - 1) is + /// (2 q A B + (p + q x) W)/(2 A^j sqrt(B)) with W = B' A - 2 (j - 1) B A', + /// so L/(A^j sqrt(B)) is that derivative and S/(A^(j - 1) sqrt(B)) wherever + /// 2 L = 2 q A B + (p + q x) W + 2 A S. That is five linear equations, one for each + /// power of x up to the fourth, in p, q and the three coefficients of a + /// quadratic S, and they have a solution where A and B have no root + /// in common: Hermite's reduction, as Rubi's 1.2.1.6 takes these. Over A itself, + /// what is left is sigma A + L, sigma over the root alone and L over + /// A beside it. + /// + /// + /// Beside a half-odd power of the secant or the cosecant, the half-angle tangent leaves + /// these: sqrt(a + a sec(x))/(c + d sec(x))^2 has (c + d) + (d - c) t^2 + /// squared below the bar beside sqrt(1 - t^2). + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + /// + internal static Entity? SolveAPolynomialOverAPowerOfAQuadraticBesideTheRootOfAnother(Entity expr, Entity.Variable x) + { + var (numerator, denominator) = Functions.SingleQuotient.Of(expr); + Entity? radicand = null; + Entity? quadratic = null; + var power = 0; + Entity constant = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(denominator)) + { + if (!factor.ContainsNode(x)) + { + constant = constant * factor; + continue; + } + if (factor is Powf(var @base, Number.Rational half) && half == Number.Rational.Create(1, 2) && radicand is null) + radicand = @base; + else if (quadratic is null && factor is Powf(var repeated, Number.Integer { EInteger.Sign: > 0 } whole) + && repeated is not Powf && whole.EInteger.CanFitInInt32() && whole.EInteger.ToInt32Unchecked() >= 2) + { + quadratic = repeated; + power = whole.EInteger.ToInt32Unchecked(); + } + else + return null; + } + if (radicand is null || quadratic is null || quadratic == radicand + || !TreeAnalyzer.TryGetPolyQuadratic(quadratic, x, out var a2, out var a1, out var a0) + || !TreeAnalyzer.TryGetPolyQuadratic(radicand, x, out var b2, out var b1, out var b0) + || !TreeAnalyzer.TryGetPolynomial(numerator, x, out var monomials) + || VanishesIdentically(a2) || VanishesIdentically(b2)) + return null; + if (monomials.Any(pair => pair.Key.Sign < 0 || pair.Value.ContainsNode(x)) + || monomials.Count > 0 && monomials.Keys.Max()!.CompareTo(EInteger.FromInt32(2 * power)) >= 0) + return null; + // A coefficient that is a number is a real one, as in every rule about a real root. + foreach (var coefficient in new[] { a0, a1, a2, b0, b1, b2, constant }.Concat(monomials.Values)) + if (coefficient.Evaled is Number.Complex and not Number.Real) + return null; + // Every coefficient without the condition a division by a symbol writes on it, + // `k/f provided not f = 0`: that is the generic case this answers in, and the linear + // solve below reads no condition. + static Entity WithoutConditions(Entity coefficient) + => coefficient.InnerSimplified.Replace(node => node is Providedf(var inner, _) ? inner : node); + (a0, a1, a2, b0, b1, b2) = (WithoutConditions(a0), WithoutConditions(a1), WithoutConditions(a2), WithoutConditions(b0), WithoutConditions(b1), WithoutConditions(b2)); + var coefficients = new Entity[2 * power]; + for (var i = 0; i < coefficients.Length; i++) + coefficients[i] = monomials.TryGetValue(EInteger.FromInt32(i), out var monomial) ? WithoutConditions(monomial) : Number.Integer.Zero; + + var root = MathS.Sqrt(radicand); + Entity answer = Number.Integer.Zero; + for (var j = power; j >= 1; j--) + { + // The polynomial over A^j as L + A R: R's coefficients and L's. + var rest = (Entity[])coefficients.Clone(); + var quotient = new Entity[System.Math.Max(rest.Length - 2, 0)]; + for (var d = rest.Length - 1; d >= 2; d--) + { + var term = WithoutConditions(rest[d] / a2); + quotient[d - 2] = term; + rest[d] = Number.Integer.Zero; + rest[d - 1] = WithoutConditions(rest[d - 1] - term * a1); + rest[d - 2] = WithoutConditions(rest[d - 2] - term * a0); + } + var (g, h) = (rest.Length > 0 ? rest[0] : Number.Integer.Zero, rest.Length > 1 ? rest[1] : Number.Integer.Zero); + if (j == 1) + { + // sigma A + L over A, beside the root: sigma over the root alone, and L over A. + if (quotient.Length > 1 && quotient.Skip(1).Any(extra => !VanishesIdentically(extra))) + return null; + if (quotient.Length > 0 && !VanishesIdentically(quotient[0])) + { + if (Integration.ComputeIndefiniteIntegral(quotient[0] / root, x, integrateByParts: false) is not { } overTheRoot) + return null; + answer = answer + overTheRoot; + } + if (!VanishesIdentically(g) || !VanishesIdentically(h)) + { + if (SolveALinearOverAQuadraticBesideTheRootOfAnother((g + h * x) / (quadratic * root), x) is not { } overTheQuadratic) + return null; + answer = answer + overTheQuadratic; + } + break; + } + // 2 L = 2 q A B + (p + q x) W + 2 A S, one equation for each power of x, in + // p, q and S = S0 + S1 x + S2 x^2. The fourth power's gives S2, and the third's + // and the second's S1 and S0, each linear in p and q; the first's and the + // constant's are then two equations in p and q alone, solved by their + // determinant, which is not zero where A and B have no root in common. + var m = j - 1; + var w3 = WithoutConditions(2 * a2 * b2 * (1 - 2 * m)); + var w2 = WithoutConditions(2 * a1 * b2 + a2 * b1 - 2 * m * (a1 * b2 + 2 * a2 * b1)); + var w1 = WithoutConditions(2 * a0 * b2 + a1 * b1 - 2 * m * (a1 * b1 + 2 * a2 * b0)); + var w0 = WithoutConditions(a0 * b1 - 2 * m * a1 * b0); + // S2 = sigma2 q, S1 = alpha1 q + beta1 p, S0 = alpha0 q + beta0 p. + var sigma2 = WithoutConditions(2 * b2 * (m - 1)); + var alpha1 = WithoutConditions(-(2 * (a2 * b1 + a1 * b2) + w2 + 2 * a1 * sigma2) / (2 * a2)); + var beta1 = WithoutConditions(-w3 / (2 * a2)); + var alpha0 = WithoutConditions(-(2 * (a2 * b0 + a1 * b1 + a0 * b2) + w1 + 2 * a0 * sigma2 + 2 * a1 * alpha1) / (2 * a2)); + var beta0 = WithoutConditions(-(w2 + 2 * a1 * beta1) / (2 * a2)); + // p and q from the constant's equation and the first power's. + var pOfConstant = WithoutConditions(w0 + 2 * a0 * beta0); + var qOfConstant = WithoutConditions(2 * a0 * b0 + 2 * a0 * alpha0); + var pOfFirst = WithoutConditions(w1 + 2 * a0 * beta1 + 2 * a1 * beta0); + var qOfFirst = WithoutConditions(2 * (a1 * b0 + a0 * b1) + w0 + 2 * a0 * alpha1 + 2 * a1 * alpha0); + var determinant = WithoutConditions(pOfConstant * qOfFirst - qOfConstant * pOfFirst); + if (VanishesIdentically(determinant)) + return null; + var pValue = WithoutConditions((2 * g * qOfFirst - 2 * h * qOfConstant) / determinant); + var qValue = WithoutConditions((2 * h * pOfConstant - 2 * g * pOfFirst) / determinant); + if (pValue.Complexity + qValue.Complexity > LargestReducedCoefficients) + return null; + var values = new[] + { + pValue, qValue, + WithoutConditions(alpha0 * qValue + beta0 * pValue), + WithoutConditions(alpha1 * qValue + beta1 * pValue), + WithoutConditions(sigma2 * qValue), + }; + answer = answer + (values[0] + values[1] * x) * root / MathS.Pow(quadratic, m); + // What is left over A^(j - 1): R, and S. + var next = new Entity[2 * m]; + for (var i = 0; i < next.Length; i++) + { + Entity sum = i < quotient.Length ? quotient[i] : Number.Integer.Zero; + if (i < 3) + sum = sum + values[2 + i]; + next[i] = WithoutConditions(sum); + } + // S is a quadratic, and over A^(j - 1) with j - 1 = 1 its square term is sigma's. + if (m == 1 && !VanishesIdentically(values[4])) + next = new[] { next[0], next[1], WithoutConditions(values[4]) }; + coefficients = next; + } + answer = (answer / constant).InnerSimplified; + if (answer.Nodes.Any(node => node is Number.Complex { IsNaN: true }) + || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x)) + return null; + return answer; + } + /// /// A polynomial over a power of a linear beside the square root of a quadratic, /// P/((x - p)^k sqrt(Q)) with P of degree below k, by the reciprocal @@ -11456,11 +11639,12 @@ Dictionary PowerOfQ(int j) // numeric, and as written otherwise. var factors = new List(); Entity constant = Number.Integer.One; - // Or with a quadratic among them, a power of a linear, or a power of the radicand: - // the partial fractions split over those as well, and each piece is closed beside - // the root -- over a quadratic by the rule before this one, over a power of a linear - // by its reciprocal, over a power of the radicand by the reduction for a polynomial - // beside a half-odd power. `tan(x)^5 sqrt(a + b tan(x) + c tan(x)^2)` is + // Or with a quadratic among them, a power of one, a power of a linear, or a power of + // the radicand: the partial fractions split over those as well, and each piece is + // closed beside the root -- over a quadratic by the rule before this one, over a power + // of a quadratic a power at a time down to that, over a power of a linear by its + // reciprocal, over a power of the radicand by the reduction for a polynomial beside a + // half-odd power. `tan(x)^5 sqrt(a + b tan(x) + c tan(x)^2)` is // `t^5 sqrt(a + b t + c t^2)/(1 + t^2)` under the tangent, a polynomial over the root // and a linear over `1 + t^2` beside it. var beyondTheLinears = false; @@ -11482,11 +11666,14 @@ Dictionary PowerOfQ(int j) && TreeAnalyzer.TryGetPolyLinear(repeated, x, out var repeatedSlope, out _) && repeatedSlope.Evaled is not Number.Complex { IsZero: true }; var isAPowerOfTheRadicand = factor is Powf(var raised, Number.Integer { EInteger.Sign: > 0 }) && raised == radicand; var isAQuadratic = factor is not Powf && TreeAnalyzer.TryGetPolyQuadratic(factor, x, out var leading, out _, out _) && !VanishesIdentically(leading); - if (!isAPowerOfALinear && !isAPowerOfTheRadicand && !isAQuadratic) + var isAPowerOfAQuadratic = factor is Powf(var repeatedQuadratic, Number.Integer { EInteger.Sign: > 0 }) && repeatedQuadratic != radicand + && repeatedQuadratic is not Powf && TreeAnalyzer.TryGetPolyQuadratic(repeatedQuadratic, x, out var repeatedLeading, out _, out _) + && !VanishesIdentically(repeatedLeading); + if (!isAPowerOfALinear && !isAPowerOfTheRadicand && !isAQuadratic && !isAPowerOfAQuadratic) return null; factors.Add(factor); beyondTheLinears = true; - aQuadratic |= isAQuadratic || isAPowerOfTheRadicand; + aQuadratic |= isAQuadratic || isAPowerOfTheRadicand || isAPowerOfAQuadratic; } // Powers of linears only beside a quadratic, which is what the split is for: over // the linears alone the reciprocal writes a sign of each, `sgn(sqrt(1 + x) - 1)` once @@ -11508,8 +11695,11 @@ Dictionary PowerOfQ(int j) // P/D as a polynomial plus a proper part, the proper part over each linear. Entity polynomialPart = Number.Integer.Zero; Entity properNumerator = above; + // A zero is one with the condition the division writes on it as well, `0 provided not + // 1 - t^2 = 0`, which read as a piece asks for its integral and declines the whole. + static bool IsZero(Entity e) => Functions.PartialFractions.Bare(e).Evaled is Number.Complex { IsZero: true }; if (TreeAnalyzer.PolynomialLongDivision(above, linears, genericCase: true, inTermsOf: x) is var (quotient, proper) - && quotient.Evaled is not Number.Complex { IsZero: true }) + && !IsZero(quotient)) { polynomialPart = quotient; var (properTop, properBottom) = Functions.SingleQuotient.Of(proper); @@ -11520,9 +11710,9 @@ Dictionary PowerOfQ(int j) return null; } var pieces = new List(); - if (polynomialPart.Evaled is not Number.Complex { IsZero: true }) + if (!IsZero(polynomialPart)) pieces.Add(polynomialPart); - if (properNumerator.Evaled is not Number.Complex { IsZero: true }) + if (!IsZero(properNumerator)) { if (factors.Count == 1) pieces.Add(properNumerator / factors[0]); @@ -11532,7 +11722,7 @@ Dictionary PowerOfQ(int j) var (terms, over) = decomposition is Divf(var splitTop, var splitBottom) && !splitBottom.ContainsNode(x) ? (splitTop, splitBottom) : (decomposition, Number.Integer.One as Entity); foreach (var term in Sumf.LinearChildren(terms)) - if (term.Evaled is not Number.Complex { IsZero: true }) + if (!IsZero(term)) pieces.Add(term / over); } else @@ -11558,6 +11748,7 @@ Dictionary PowerOfQ(int j) var integrated = SolveALinearBesideTheRootOfAQuadratic(overTheRoot, x) ?? (beyondTheLinears ? SolveALinearOverAQuadraticBesideTheRootOfAnother(overTheRoot, x) + ?? SolveAPolynomialOverAPowerOfAQuadraticBesideTheRootOfAnother(overTheRoot, x) ?? SolveAPolynomialOverAPowerOfALinearBesideTheRoot(overTheRoot, x) ?? SolveAPolynomialTimesAnOddHalfPowerOfAQuadratic( OverAPowerOfTheRadicand((piece / constant).InnerSimplified, radicand, x) ?? overTheRoot, x) diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 73efda8ca..92dc8b563 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -926,6 +926,8 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // And over a quadratic, closed as well: beside the linear's, since each of the two is // a piece the rational function below is taken apart into. if ((answer = IndefiniteIntegralSolver.SolveALinearOverAQuadraticBesideTheRootOfAnother(expr, x)) is { }) return answer; + // And over a power of the quadratic, a power at a time down to that one. + if ((answer = IndefiniteIntegralSolver.SolveAPolynomialOverAPowerOfAQuadraticBesideTheRootOfAnother(expr, x)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionBesideTheRootOfAQuadratic(expr, x)) is { }) return answer; // A linear below the bar that divides the radicand, written over it: after the two // rules above, which read a linear beside the root as a pole and find nothing to take diff --git a/Sources/Tests/UnitTests/Calculus/QuadraticBesideARootIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/QuadraticBesideARootIntegralTest.cs index 842e4a9e6..e80bfa78b 100644 --- a/Sources/Tests/UnitTests/Calculus/QuadraticBesideARootIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/QuadraticBesideARootIntegralTest.cs @@ -95,5 +95,30 @@ public void UnderTheTangent(string integrand) [InlineData("x^2/((a + b*x + c*x^2)^(1/2)*(d + k*x + f*x^2))")] public void APolynomialPartBeside(string integrand) => DifferentiatesBackPinned(integrand, new[] { 0.3, 0.7, 1.2, -0.5 }, Pins); + + /// + /// Over a power of the quadratic, a power at a time: the derivative of + /// (p + q x) sqrt(B)/A^(k - 1) takes one away, and what it leaves goes down to the + /// closed form over the first. The cube is numeric: with six symbols its coefficients + /// grow past the bound the reduction declines at. + /// + [Theory] + [InlineData("1/((1 + x^2)^2*sqrt(1 - x^2))")] + [InlineData("(g + h*x)/((d + k*x + f*x^2)^2*sqrt(a + b*x + c*x^2))")] + [InlineData("(7 + 13*x)/((5 + x + 2*x^2)^3*sqrt(2 + x + 3*x^2))")] + [InlineData("x^3/((d + f*x^2)^2*sqrt(a + b*x + c*x^2))")] + public void OverAPowerOfTheQuadratic(string integrand) + => DifferentiatesBackPinned(integrand, new[] { 0.3, 0.7, 0.9, -0.5 }, Pins); + + /// + /// Beside a half-odd power of a + a sec(x), the half-angle tangent leaves a power + /// of (c + d) + (d - c) t^2 below the bar, which ran out Rubi's budget in 4.5.2.1. + /// + [Theory] + [InlineData("sqrt(a + a*sec(x))/(c + d*sec(x))^2")] + [InlineData("(a + a*sec(x))^(3/2)/(c + d*sec(x))^2")] + [InlineData("sqrt(a + a*sec(x))/(c + d*sec(x))^3")] + public void BesideAHalfOddPowerOfOnePlusASecant(string integrand) + => DifferentiatesBackPinned(integrand, new[] { 0.3, 0.7, 1.1, -0.4 }, ("a", 1.3), ("c", 2), ("d", 0.7)); } }