@@ -19588,6 +19588,95 @@ void Add(int j, int k, Entity[] poly)
1958819588 return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
1958919589 }
1959019590
19591+ /// <summary>
19592+ /// A power of <c>x</c> over a sextic that is a sum of two cubes,
19593+ /// <c>d6 ((x^2 + s)^3 + q x^3)</c>, split over the cube roots of unity:
19594+ /// <c>1/(U^3 + V^3) = sum_j w_j/(3 V^2 (U + w_j V))</c> for <c>U = x^2 + s</c>,
19595+ /// <c>V = m x</c>, <c>m^3 = q</c>, so that <c>x^k/(U^3 + V^3)</c> is
19596+ /// <c>sum_j w_j x^(k - 2)/(3 m^2 (x^2 + w_j m x + s))</c>, three quadratics.
19597+ /// </summary>
19598+ /// <remarks>
19599+ /// The sextic is read as <c>d6 x^6 + d4 x^4 + d3 x^3 + d2 x^2 + d0</c> with
19600+ /// <c>d4^2 = 3 d6 d2</c> and <c>d2^2 = 3 d4 d0</c>, which is what makes
19601+ /// <c>d6 x^6 + d4 x^4 + d2 x^2 + d0</c> the cube <c>d6 (x^2 + s)^3</c> with
19602+ /// <c>s = d4/(3 d6)</c>. The identity is the partial fractions of <c>1/(U^3 + V^3)</c> in
19603+ /// <c>U</c>, whose poles <c>-w_j V</c> have residues <c>1/(3 w_j^2 V^2) = w_j/(3 V^2)</c>; any
19604+ /// cube root <c>m</c> of <c>q</c> gives it. The three quadratics are one integral
19605+ /// <c>x^n/(x^2 + B x + S)</c> asked with <c>B</c> and <c>S</c> named. Rubi's 1.3.1
19606+ /// <c>x^k/(27 a^3 + 27 a^2 b x^2 + 27 a^2 c x^3 + 9 a b^2 x^4 + b^3 x^6)</c> was declined,
19607+ /// the integrator's resultant needing numbers for coefficients.
19608+ /// https://github.com/asc-community/AngouriMath/issues/718
19609+ /// </remarks>
19610+ internal static Entity? SolveAPowerOfXOverASumOfTwoCubes(Entity expr, Entity.Variable x, bool integrateByParts)
19611+ {
19612+ if (!Integration.AnsweringTheQuestionAsked || !expr.Vars.Any(symbol => symbol != x))
19613+ return null;
19614+ var (numerator, denominator) = Functions.SingleQuotient.Of(expr);
19615+ Entity? sextic = null;
19616+ var power = 0;
19617+ Entity constant = Number.Integer.One;
19618+ foreach (var (side, below) in new[] { (numerator, false), (denominator, true) })
19619+ foreach (var factor in Mulf.LinearChildren(side))
19620+ {
19621+ if (!factor.ContainsNode(x))
19622+ {
19623+ constant = below ? constant / factor : constant * factor;
19624+ continue;
19625+ }
19626+ if (factor == x || factor is Powf(var bx, Number.Integer { EInteger: var wx }) && bx == x && wx.CanFitInInt32())
19627+ {
19628+ var whole = factor == x ? 1 : ((Number.Integer)((Powf)factor).Exponent).EInteger.ToInt32Unchecked();
19629+ power += below ? -whole : whole;
19630+ continue;
19631+ }
19632+ if (below && sextic is null)
19633+ {
19634+ sextic = factor;
19635+ continue;
19636+ }
19637+ return null;
19638+ }
19639+ if (sextic is null || Functions.PartialFractions.HoldsTheImaginaryUnit(sextic)
19640+ || !TreeAnalyzer.TryGetPolynomial(sextic, x, out var read) || read.Keys.Any(key => key.Sign < 0 || key.CompareTo(EInteger.FromInt32(6)) > 0))
19641+ return null;
19642+ Entity At(int degree) => read.TryGetValue(EInteger.FromInt32(degree), out var coefficient) ? coefficient.InnerSimplified : Number.Integer.Zero;
19643+ var (d6, d5, d4, d3, d2, d1, d0) = (At(6), At(5), At(4), At(3), At(2), At(1), At(0));
19644+ if (!VanishesIdentically(d5) || !VanishesIdentically(d1)
19645+ || VanishesIdentically(d6) || VanishesIdentically(d4) || VanishesIdentically(d3) || VanishesIdentically(d0)
19646+ || !VanishesIdentically((d4 * d4 - 3 * d6 * d2).Expand()) || !VanishesIdentically((d2 * d2 - 3 * d4 * d0).Expand()))
19647+ return null;
19648+ var shift = Functions.PartialFractions.InLowestTermsOverTheSymbols(d4 / (3 * d6));
19649+ var cube = Functions.PartialFractions.InLowestTermsOverTheSymbols(d3 / d6);
19650+ var root = MathS.Pow(cube, Number.Rational.Create(1, 3));
19651+ var namedB = Variable.CreateUnique(expr, "B_c");
19652+ var namedS = Variable.CreateUnique(expr + namedB, "S_c");
19653+ var term = MathS.Pow(x, Number.Integer.Create(power - 2)) / (MathS.Pow(x, Number.Integer.Create(2)) + namedB * x + namedS);
19654+ if (Integration.ComputeAsAQuestionOfItsOwn(term, x, integrateByParts) is not { } each)
19655+ return null;
19656+ // The integral over a real quadratic splits on the sign of 4 S - B^2, which means nothing
19657+ // for the complex B of two of the three: its first arm, the arctangent over
19658+ // sqrt(4 S - B^2), differentiates back to the integrand for any B with 4 S != B^2, and
19659+ // the answer is checked below.
19660+ each = each.Replace(node => node switch
19661+ {
19662+ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expression,
19663+ Providedf provided => provided.Expression,
19664+ _ => node,
19665+ });
19666+ var half = Number.Rational.Create(1, 2);
19667+ var unities = new Entity[]
19668+ {
19669+ Number.Integer.One,
19670+ -half + half * MathS.Sqrt(3) * MathS.i,
19671+ -half - half * MathS.Sqrt(3) * MathS.i,
19672+ };
19673+ Entity answer = Number.Integer.Zero;
19674+ foreach (var unity in unities)
19675+ answer = answer + unity * each.Substitute(namedB, unity * root).Substitute(namedS, shift);
19676+ answer = (constant / (3 * d6 * MathS.Pow(root, Number.Integer.Create(2))) * answer).InnerSimplified;
19677+ return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
19678+ }
19679+
1959119680 /// <summary>
1959219681 /// A power of a linear that is not whole, times a linear, over a whole power of a quadratic,
1959319682 /// with the power of the quadratic lowered one at a time:
0 commit comments