@@ -99,6 +99,8 @@ internal static class IndefiniteIntegralSolver
9999 // corpus for nothing.
100100 if (!Integration.AnsweringTheQuestionAsked)
101101 return null;
102+ if (integrateByParts && ByPartsAgainstThePowersOfX(splitted, x) is { } byPartsOverTheSum)
103+ return byPartsOverTheSum;
102104 // What is polynomial in a term is what has no root of the variable in it: a
103105 // logarithm or an exponential of x is an indeterminate as far as gathering and
104106 // cancelling go. Bronstein's
@@ -143,6 +145,276 @@ internal static class IndefiniteIntegralSolver
143145 });
144146 }
145147
148+ /// <summary>
149+ /// A sum whose terms are not elementary apart, integrated by parts against the power
150+ /// of <paramref name="x"/> in each term that has one, with what parts leaves and the
151+ /// other terms asked together as one integrand: <c>x^k f + g</c> is
152+ /// <c>x^k F + ∫(g - k x^(k - 1) F)</c> for <c>F</c> the antiderivative of <c>f</c>.
153+ /// </summary>
154+ /// <remarks>
155+ /// <para>
156+ /// <c>sqrt(csc(x)) (x cos(x) - 4 sec(x) tan(x))</c>, Rubi's 4.7.7, is
157+ /// <c>(2 x sqrt(sin(x)) - 4 sin(x)^(3/2)/cos(x))'</c>. Split into its two terms, neither is
158+ /// elementary: parts on the first leaves <c>-2 sqrt(sin(x))</c>, and that and the second
159+ /// term are elementary only together, as <c>sqrt(sin(x)) (-2 - 4/cos(x)^2)</c>. So the
160+ /// remainder is gathered into one quotient with the factors its terms share taken out,
161+ /// which is the spelling the rules read; and where it is zero, as for
162+ /// <c>x cos(x) e^sin(x) + e^sin(x)</c>, the parts are the answer.
163+ /// </para>
164+ /// <para>
165+ /// Only where some term with a power of <paramref name="x"/> beside something else in
166+ /// <paramref name="x"/> has an elementary <c>F</c>, and only at the question asked: the
167+ /// remainder is a new question, and the powers of <paramref name="x"/> of the terms
168+ /// integrated are lower in it, so asking it again descends.
169+ /// </para>
170+ /// </remarks>
171+ private static Entity? ByPartsAgainstThePowersOfX(IReadOnlyList<Entity> terms, Entity.Variable x)
172+ {
173+ // Once: the cofactors and the remainder are asked at the top, where this runs again,
174+ // and `(c + d x) erfc(a + b x)^2` went from a 13-second decline to a hang that way.
175+ if (integratingASumByParts)
176+ {
177+ Integration.DeclinedForItsScope();
178+ return null;
179+ }
180+ integratingASumByParts = true;
181+ try
182+ {
183+ return ByPartsAgainstThePowersOfXOnce(terms, x);
184+ }
185+ finally
186+ {
187+ integratingASumByParts = false;
188+ }
189+ }
190+
191+ /// <summary>
192+ /// Whether a term of <paramref name="antiderivative"/> has a factor that is a polynomial
193+ /// in <paramref name="x"/> of positive degree.
194+ /// </summary>
195+ private static bool HasAPolynomialFactorIn(Entity antiderivative, Entity.Variable x)
196+ => Sumf.LinearChildren(antiderivative).Any(term => Mulf.LinearChildren(term).Any(factor =>
197+ factor.ContainsNode(x) && TreeAnalyzer.TryGetPolynomial(factor, x, out var powers)
198+ && powers.Keys.All(power => power.Sign >= 0) && powers.Values.All(coefficient => !coefficient.ContainsNode(x))));
199+
200+ /// <summary>
201+ /// Whether every occurrence of <paramref name="x"/> in <paramref name="expr"/> is inside
202+ /// the argument of a function that is not algebraic: <c>cos(x)/sqrt(sin(x))</c>, not
203+ /// <c>sqrt(x^2 + 3)</c> or <c>atanh(x)/(1 - x^2)</c>.
204+ /// </summary>
205+ private static bool OnlyInsideTranscendentalFunctions(Entity expr, Entity.Variable x)
206+ => !expr.ContainsNode(x)
207+ || expr is not (Sumf or Minusf or Mulf or Divf or Powf or Entity.Variable or Number)
208+ || expr is not Entity.Variable && expr.DirectChildren.All(child => OnlyInsideTranscendentalFunctions(child, x));
209+
210+ /// <summary>Set while <see cref="ByPartsAgainstThePowersOfX"/> asks its questions.</summary>
211+ [System.ThreadStatic] private static bool integratingASumByParts;
212+
213+ private static Entity? ByPartsAgainstThePowersOfXOnce(IReadOnlyList<Entity> terms, Entity.Variable x)
214+ {
215+ Entity parts = Number.Integer.Zero;
216+ var remainder = new List<Entity>();
217+ var any = false;
218+ foreach (var term in terms)
219+ {
220+ var power = EInteger.Zero;
221+ Entity rest = Number.Integer.One;
222+ foreach (var factor in Mulf.LinearChildren(term))
223+ if (factor == x)
224+ power = power.Add(EInteger.One);
225+ else if (factor is Powf(var @base, Number.Integer { EInteger.Sign: > 0 } exponent) && @base == x)
226+ power = power.Add(exponent.EInteger);
227+ else
228+ rest *= factor;
229+ // Only a cofactor with x in it nowhere but inside transcendental functions. A power
230+ // of x beside an algebraic one is what the algebraic reductions answer, and asking
231+ // each such cofactor and the remainder took `(x^2 + 3)^2/(x + sqrt(x^2 + 3))` from a
232+ // half-second decline past two minutes; and beside a rational function of x as
233+ // well, `x^3 atanh(c x^2)/(1 - c^2 x^4)` in what parts leaves of
234+ // `(a + b atanh(c x^2))^2/x^5`, the cofactor is a slow decline that doubled its time.
235+ if (power.IsZero || !rest.ContainsNode(x) || !OnlyInsideTranscendentalFunctions(rest, x))
236+ {
237+ remainder.Add(term);
238+ continue;
239+ }
240+ // A term whose cofactor is not elementary goes into the remainder as it is: in
241+ // `x^2 cos(x) e^sin(x) + 2 x e^sin(x)` the second is what parts on the first leaves.
242+ // And one with a polynomial in x beside something in its antiderivative: what parts
243+ // leaves then has the power it took away, and asking it does not descend --
244+ // `x erfc(a + b x)^2`, whose cofactor's antiderivative has `x erfc(a + b x)^2` in it.
245+ if (Integration.ComputeAsAQuestionOfItsOwn(rest, x, integrateByParts: true) is not { } antiderivative
246+ || HasAPolynomialFactorIn(antiderivative, x))
247+ {
248+ remainder.Add(term);
249+ continue;
250+ }
251+ var monomial = power.Equals(EInteger.One) ? (Entity)x : MathS.Pow(x, Number.Integer.Create(power));
252+ parts += monomial * antiderivative;
253+ // No `x^0` for a first power: it would bring its `provided not x = 0` along.
254+ var lower = power.Equals(EInteger.One) ? Number.Integer.One : MathS.Pow(x, Number.Integer.Create(power.Subtract(EInteger.One)));
255+ remainder.Add(-Number.Integer.Create(power) * lower * antiderivative);
256+ any = true;
257+ }
258+ if (!any)
259+ return null;
260+ var left = remainder.Aggregate((l, r) => l + r).InnerSimplified;
261+ if (TreeAnalyzer.IsZero(left))
262+ return parts;
263+ var gathered = OverTheLeastPowersOfTheSineAndTheCosine(left, x, out var closed);
264+ if (closed is not null)
265+ return parts + closed;
266+ if (gathered is null)
267+ {
268+ var (above, below) = Functions.SingleQuotient.Of(left);
269+ gathered = WithTheFactorsEveryTermSharesOut(above.InnerSimplified, x) / below;
270+ }
271+ return Integration.ComputeAsAQuestionOfItsOwn(gathered, x, integrateByParts: true) is { } integrated
272+ ? parts + integrated
273+ : null;
274+ }
275+
276+ /// <summary>
277+ /// <paramref name="sum"/> with each term read as a number times powers of
278+ /// <c>sin(u)</c> and <c>cos(u)</c> for one argument <c>u</c>, and written as the least
279+ /// power of each times the sum of what is left; <see langword="null"/> where a term
280+ /// is not such a product.
281+ /// </summary>
282+ /// <remarks>
283+ /// What parts leaves is spelled the way the antiderivative came back, and the other
284+ /// terms the way they were written: <c>(1/sin(x))^(-1/2)</c> beside
285+ /// <c>sqrt(csc(x)) sec(x) tan(x)</c> share nothing as written, and are
286+ /// <c>sqrt(sin(x))</c> and <c>sqrt(sin(x))/cos(x)^2</c> as powers. Reading
287+ /// <c>(1/f)^p</c> as <c>f^(-p)</c> holds wherever <c>f</c> is positive or <c>p</c> is
288+ /// whole; where neither, the integrand is not real.
289+ /// <para>
290+ /// And <paramref name="antiderivative"/> where the sum is the derivative of one such
291+ /// product, which the rules otherwise answer in the half angle's tangent at four times
292+ /// the size: <c>(K sin(u)^(a + 1) cos(u)^(b + 1))'</c> is
293+ /// <c>K u' sin(u)^a cos(u)^b ((a + b + 2) cos(u)^2 - (b + 1))</c>, and
294+ /// <c>K u' sin(u)^a cos(u)^b ((a + 1) - (a + b + 2) sin(u)^2)</c> the other way, so
295+ /// <c>sqrt(sin(x)) (-2 - 4/cos(x)^2)</c> is <c>(-4 sin(x)^(3/2)/cos(x))'</c>.
296+ /// </para>
297+ /// </remarks>
298+ private static Entity? OverTheLeastPowersOfTheSineAndTheCosine(Entity sum, Entity.Variable x, out Entity? antiderivative)
299+ {
300+ antiderivative = null;
301+ Entity? argument = null;
302+ var read = new List<(Entity Coefficient, Number.Rational Sine, Number.Rational Cosine)>();
303+ foreach (var term in TreeAnalyzer.GatherLinearChildrenOverSumAndExpand(sum, e => e.ContainsNode(x)) ?? new List<Entity> { sum })
304+ if (Read(term) is { } powers)
305+ read.Add(powers);
306+ else
307+ return null;
308+ if (argument is null || read.Count < 2)
309+ return null;
310+ var leastSine = read.Select(r => r.Sine).Aggregate((l, r) => r < l ? r : l);
311+ var leastCosine = read.Select(r => r.Cosine).Aggregate((l, r) => r < l ? r : l);
312+ Entity Power(Entity @base, Number.Rational power)
313+ => power == Number.Integer.Zero ? Number.Integer.One : power == Number.Integer.One ? @base : MathS.Pow(@base, power);
314+ var slope = argument.Differentiate(x).InnerSimplified;
315+ var two = Number.Integer.Create(2);
316+ if (!slope.ContainsNode(x))
317+ {
318+ Entity Total(System.Func<(Entity Coefficient, Number.Rational Sine, Number.Rational Cosine), bool> which)
319+ => read.Where(which).Select(r => r.Coefficient).Aggregate((Entity)Number.Integer.Zero, (l, r) => l + r).InnerSimplified;
320+ var (a, b) = (leastSine, leastCosine);
321+ Entity? k = null;
322+ if (read.All(r => r.Sine == a && (r.Cosine == b || r.Cosine == b + two)) && b + Number.Integer.One != Number.Integer.Zero)
323+ {
324+ var (squared, alone) = (Total(r => r.Cosine != b), Total(r => r.Cosine == b));
325+ var candidate = (-alone / ((b + Number.Integer.One) * slope)).InnerSimplified;
326+ if (TreeAnalyzer.IsZero((squared - candidate * slope * (a + b + two)).InnerSimplified))
327+ k = candidate;
328+ }
329+ else if (read.All(r => r.Cosine == b && (r.Sine == a || r.Sine == a + two)) && a + Number.Integer.One != Number.Integer.Zero)
330+ {
331+ var (squared, alone) = (Total(r => r.Sine != a), Total(r => r.Sine == a));
332+ var candidate = (alone / ((a + Number.Integer.One) * slope)).InnerSimplified;
333+ if (TreeAnalyzer.IsZero((squared + candidate * slope * (a + b + two)).InnerSimplified))
334+ k = candidate;
335+ }
336+ if (k is not null)
337+ {
338+ var (sinePower, cosinePower) = (a + Number.Integer.One, b + Number.Integer.One);
339+ Entity product = k;
340+ if (sinePower < Number.Integer.Zero) product /= Power(MathS.Sin(argument), -sinePower); else product *= Power(MathS.Sin(argument), sinePower);
341+ if (cosinePower < Number.Integer.Zero) product /= Power(MathS.Cos(argument), -cosinePower); else product *= Power(MathS.Cos(argument), cosinePower);
342+ antiderivative = product;
343+ }
344+ }
345+ var left = read.Select(r => r.Coefficient * Power(MathS.Sin(argument), r.Sine - leastSine) * Power(MathS.Cos(argument), r.Cosine - leastCosine))
346+ .Aggregate((l, r) => l + r);
347+ Entity above = left, below = Number.Integer.One;
348+ if (leastSine < Number.Integer.Zero) below *= Power(MathS.Sin(argument), -leastSine); else above = Power(MathS.Sin(argument), leastSine) * above;
349+ if (leastCosine < Number.Integer.Zero) below *= Power(MathS.Cos(argument), -leastCosine); else above = Power(MathS.Cos(argument), leastCosine) * above;
350+ return below == Number.Integer.One ? above : above / below;
351+
352+ (Entity Coefficient, Number.Rational Sine, Number.Rational Cosine)? Read(Entity node)
353+ {
354+ if (!node.ContainsNode(x))
355+ return (node, Number.Integer.Zero, Number.Integer.Zero);
356+ switch (node)
357+ {
358+ case Sinf(var u): return Base(u, Number.Integer.One, Number.Integer.Zero);
359+ case Cosf(var u): return Base(u, Number.Integer.Zero, Number.Integer.One);
360+ case Secantf(var u): return Base(u, Number.Integer.Zero, Number.Integer.MinusOne);
361+ case Cosecantf(var u): return Base(u, Number.Integer.MinusOne, Number.Integer.Zero);
362+ case Tanf(var u): return Base(u, Number.Integer.One, Number.Integer.MinusOne);
363+ case Cotanf(var u): return Base(u, Number.Integer.MinusOne, Number.Integer.One);
364+ case Mulf(var l, var r):
365+ return (Read(l), Read(r)) is ({ } a, { } b) ? ((a.Coefficient * b.Coefficient).InnerSimplified, a.Sine + b.Sine, a.Cosine + b.Cosine) : null;
366+ case Divf(var l, var r):
367+ return (Read(l), Read(r)) is ({ } n, { } d) ? ((n.Coefficient / d.Coefficient).InnerSimplified, n.Sine - d.Sine, n.Cosine - d.Cosine) : null;
368+ case Powf(var @base, Number.Rational power) when Read(@base) is { } inner:
369+ // A whole power of anything read; a fractional one only of a single
370+ // power of the sine or of the cosine with nothing beside it, so that
371+ // no root of a product is split.
372+ if (power is Number.Integer)
373+ return (MathS.Pow(inner.Coefficient, power).InnerSimplified, inner.Sine * power, inner.Cosine * power);
374+ if (inner.Coefficient == Number.Integer.One && (inner.Sine == Number.Integer.Zero || inner.Cosine == Number.Integer.Zero))
375+ return (Number.Integer.One, inner.Sine * power, inner.Cosine * power);
376+ return null;
377+ default:
378+ return null;
379+ }
380+ }
381+
382+ (Entity Coefficient, Number.Rational Sine, Number.Rational Cosine)? Base(Entity u, Number.Rational sine, Number.Rational cosine)
383+ {
384+ if (argument is not null && argument != u)
385+ return null;
386+ argument = u;
387+ return (Number.Integer.One, sine, cosine);
388+ }
389+ }
390+
391+ /// <summary>
392+ /// <paramref name="sum"/> as the factors with <paramref name="x"/> in them that every
393+ /// term has, as written, times the sum of what is left; as written where there are none.
394+ /// </summary>
395+ private static Entity WithTheFactorsEveryTermSharesOut(Entity sum, Entity.Variable x)
396+ {
397+ var terms = Sumf.LinearChildren(sum).Select(term => Mulf.LinearChildren(term).ToList()).ToList();
398+ if (terms.Count < 2)
399+ return sum;
400+ var shared = new List<Entity>();
401+ foreach (var factor in terms[0])
402+ if (factor.ContainsNode(x) && terms.Skip(1).All(other => other.Contains(factor)))
403+ {
404+ shared.Add(factor);
405+ foreach (var other in terms.Skip(1))
406+ other.Remove(factor);
407+ }
408+ if (shared.Count == 0)
409+ return sum;
410+ var firstLeft = new List<Entity>(terms[0]);
411+ foreach (var factor in shared)
412+ firstLeft.Remove(factor);
413+ Entity Product(List<Entity> factors) => factors.Count == 0 ? Number.Integer.One : factors.Aggregate((l, r) => l * r);
414+ var left = new[] { Product(firstLeft) }.Concat(terms.Skip(1).Select(Product)).Aggregate((l, r) => l + r);
415+ return Product(shared) * left;
416+ }
417+
146418 /// <summary>The largest term a sum is split over as written, before the expanding gather.</summary>
147419 private const int LargestTermTakenAsWritten = 60;
148420
0 commit comments