@@ -19544,6 +19544,137 @@ void Add(int j, int k, Entity[] poly)
1954419544 return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
1954519545 }
1954619546
19547+ /// <summary>
19548+ /// A power of a linear that is not whole, times a linear, over a whole power of a quadratic,
19549+ /// with the power of the quadratic lowered one at a time:
19550+ /// <c>integral L^m N Q^p = (alpha + beta x) L^(m + 1) Q^(p + 1) + integral L^m W Q^(p + 1)</c>
19551+ /// for <c>p <= -2</c>, where <c>alpha</c> and <c>beta</c> make
19552+ /// <c>N - (alpha + beta x) (p + 1) L Q'</c> a multiple of <c>Q</c> and
19553+ /// <c>W = (N - (alpha + beta x) (p + 1) L Q')/Q - beta L - (alpha + beta x) (m + 1) e</c> is
19554+ /// linear again; at <c>p = -1</c> the integrator answers what is left.
19555+ /// </summary>
19556+ /// <remarks>
19557+ /// The derivative of the first term is <c>L^m Q^p</c> times
19558+ /// <c>beta L Q + (alpha + beta x) (m + 1) e Q + (alpha + beta x) L (p + 1) Q'</c>, which with
19559+ /// <c>L^m Q^(p + 1) W</c> added is <c>L^m Q^p N</c>. <c>(alpha + beta x) L Q'</c> modulo
19560+ /// <c>Q</c> is <c>(alpha + beta x) (r1 x + r0)</c> with <c>r1 = 2 c d - b e</c> and
19561+ /// <c>r0 = b d - 2 a e</c>, for <c>L = d + e x</c> and <c>Q = a + b x + c x^2</c>, and
19562+ /// matching it to <c>N</c> is two equations in two unknowns. Rubi's 1.2.1.2
19563+ /// <c>(d + e x)^(7/2)/(a + b x + c x^2)^2</c> ran past a minute and
19564+ /// <c>sqrt(d + e x)/(a + b x + c x^2)^3</c> took thirty seconds for 26,437 characters.
19565+ /// https://github.com/asc-community/AngouriMath/issues/718
19566+ /// </remarks>
19567+ internal static Entity? SolveAPowerOfALinearOverAPowerOfAQuadratic(Entity expr, Entity.Variable x, bool integrateByParts)
19568+ {
19569+ if (!Integration.AnsweringTheQuestionAsked)
19570+ return null;
19571+ var (numerator, denominator) = Functions.SingleQuotient.Of(expr);
19572+ Entity? quadratic = null;
19573+ var power = 0;
19574+ Entity? linear = null;
19575+ Entity exponent = Number.Integer.Zero;
19576+ Entity? other = null;
19577+ Entity constant = Number.Integer.One;
19578+ foreach (var (side, below) in new[] { (numerator, false), (denominator, true) })
19579+ foreach (var factor in Mulf.LinearChildren(side))
19580+ {
19581+ if (!factor.ContainsNode(x))
19582+ {
19583+ constant = below ? constant / factor : constant * factor;
19584+ continue;
19585+ }
19586+ var (@base, raised) = factor is Powf(var b0, var raisedTo) && !raisedTo.ContainsNode(x) ? (b0, raisedTo) : (factor, (Entity)Number.Integer.One);
19587+ if (below && quadratic is null && raised is Number.Integer { EInteger: var whole } && whole.CompareTo(EInteger.FromInt32(2)) >= 0
19588+ && whole.CanFitInInt32() && TreeAnalyzer.TryGetPolyQuadratic(@base, x, out var leading, out _, out _) && !VanishesIdentically(leading))
19589+ {
19590+ quadratic = @base;
19591+ power = whole.ToInt32Unchecked();
19592+ continue;
19593+ }
19594+ if (!TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out _) || VanishesIdentically(slope))
19595+ return null;
19596+ // The linear to a power that is not whole, either side; one more linear above, once.
19597+ if (raised is not Number.Integer && (linear is null || linear == @base))
19598+ {
19599+ linear = @base;
19600+ exponent = below ? exponent - raised : exponent + raised;
19601+ continue;
19602+ }
19603+ if (!below && other is null && raised == Number.Integer.One)
19604+ {
19605+ other = @base;
19606+ continue;
19607+ }
19608+ return null;
19609+ }
19610+ exponent = exponent.InnerSimplified;
19611+ if (quadratic is null || linear is null || exponent is Number.Integer
19612+ || !TreeAnalyzer.TryGetPolyQuadratic(quadratic, x, out var c, out var b, out var a)
19613+ || !TreeAnalyzer.TryGetPolyLinear(linear, x, out var e, out var d)
19614+ || VanishesIdentically(b * b - 4 * a * c))
19615+ return null;
19616+ // Not a positive power of x itself, which u = sqrt(x) answers sooner -- the
19617+ // reduction took `x^(9/2)/(a + b x + c x^2)^3` from 2.6 seconds to 8.3; not a quadratic
19618+ // over the imaginary unit; and not one that factors over linears, b^2 - 4 a c the
19619+ // square of a polynomial, which the partial fractions take apart cheaply.
19620+ if (linear == x && exponent is Number.Rational { IsPositive: true } || Functions.PartialFractions.HoldsTheImaginaryUnit(quadratic))
19621+ return null;
19622+ var opposite = Functions.PartialFractions.Bare((b * b - 4 * a * c).Expand().InnerSimplified);
19623+ var symbols = opposite.Vars.Distinct().OrderBy(v => v.Name, System.StringComparer.Ordinal).ToArray();
19624+ var indices = new Dictionary<Variable, int>();
19625+ for (var i = 0; i < symbols.Length; i++)
19626+ indices[symbols[i]] = i;
19627+ if (symbols.Length <= MultivariatePolynomial.MaxVariables
19628+ && MultivariatePolynomial.TryParse(opposite, indices) is { } parsed && parsed.TrySquareRoot() is not null)
19629+ return null;
19630+ Entity nb = Number.Integer.Zero, na = constant;
19631+ if (other is not null)
19632+ {
19633+ if (!TreeAnalyzer.TryGetPolyLinear(other, x, out var ob, out var oa))
19634+ return null;
19635+ nb = (constant * ob).InnerSimplified;
19636+ na = (constant * oa).InnerSimplified;
19637+ }
19638+ var r1 = (2 * c * d - b * e).InnerSimplified;
19639+ var r0 = (b * d - 2 * a * e).InnerSimplified;
19640+ Entity answer = Number.Integer.Zero;
19641+ for (var p = -power; p <= -2; p++)
19642+ {
19643+ // (p + 1) [r1 alpha + (r0 - r1 b/c) beta] = B, (p + 1) [r0 alpha - (r1 a/c) beta] = A.
19644+ var m11 = r1;
19645+ var m12 = r0 - r1 * b / c;
19646+ var m21 = r0;
19647+ var m22 = -r1 * a / c;
19648+ var determinant = Functions.PartialFractions.InLowestTermsOverTheSymbols((p + 1) * (m11 * m22 - m12 * m21));
19649+ if (VanishesIdentically(determinant))
19650+ return null;
19651+ var alpha = Functions.PartialFractions.InLowestTermsOverTheSymbols((nb * m22 - m12 * na) / determinant);
19652+ var beta = Functions.PartialFractions.InLowestTermsOverTheSymbols((m11 * na - m21 * nb) / determinant);
19653+ // N - (alpha + beta x) (p + 1) L Q', a cubic, over Q: its quotient is linear, read off
19654+ // the two highest coefficients, the remainder being zero by the choice of alpha, beta.
19655+ var c3 = -(p + 1) * beta * e * 2 * c;
19656+ var c2 = -(p + 1) * (beta * (d * 2 * c + e * b) + alpha * e * 2 * c);
19657+ var q1 = (c3 / c).InnerSimplified;
19658+ var q0 = ((c2 - q1 * b) / c).InnerSimplified;
19659+ // W = quotient - beta L - (alpha + beta x)(m + 1) e.
19660+ var wb = Functions.PartialFractions.InLowestTermsOverTheSymbols(q1 - beta * e - beta * (exponent + 1) * e);
19661+ var wa = Functions.PartialFractions.InLowestTermsOverTheSymbols(q0 - beta * d - alpha * (exponent + 1) * e);
19662+ answer = answer + (alpha + beta * x) * MathS.Pow(linear, (exponent + 1).InnerSimplified) * MathS.Pow(quadratic, Number.Integer.Create(p + 1));
19663+ nb = wb;
19664+ na = wa;
19665+ }
19666+ // What is left asked with its two coefficients named: written out they are long, and the
19667+ // integrator took past a minute over them where with two symbols it takes a third of a
19668+ // second.
19669+ var namedA = Variable.CreateUnique(expr, "A_l");
19670+ var namedB = Variable.CreateUnique(expr + namedA, "B_l");
19671+ var left = (namedA + namedB * x) * MathS.Pow(linear, exponent) / quadratic;
19672+ if (Integration.ComputeAsAQuestionOfItsOwn(left, x, integrateByParts) is not { } rest)
19673+ return null;
19674+ answer = (answer + rest.Substitute(namedA, na).Substitute(namedB, nb)).InnerSimplified;
19675+ return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
19676+ }
19677+
1954719678 /// <summary>
1954819679 /// A power of a linear beside the derivative of a quadratic, over a whole power of that
1954919680 /// quadratic, by parts once: <c>Q' L^m/Q^n</c> with <c>u = L^m</c> and
0 commit comments