Skip to content

Commit f05a5b1

Browse files
Split a sum of two powers of even quadratics over the roots of its binomial (#1939)
* Split a sum of two powers of even quadratics over the roots of its binomial Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Record the sum of two powers of even quadratics in BREAKING-CHANGES Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --------- Co-authored-by: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
1 parent e04410a commit f05a5b1

3 files changed

Lines changed: 167 additions & 29 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -596,6 +596,20 @@ unevaluated, is answered by this route as well.
596596
| `"(a*x+2*b*n*ln(c*x^n))/(a*x^2+b*x*ln(c*x^n)^2)".ToEntity().Integrate("x")` | `integral(...)` | 37 characters |
597597
| `"tanh(a+2*ln(x))^2/x^3".ToEntity().Integrate("x")` | `integral(...)` | 7,796 characters |
598598

599+
### A quotient over a sum of two powers of even quadratics is split over the roots of its binomial
600+
601+
**Answers where there were none.** `(1 + x^2)^2/(a x^6 + b (1 + x^2)^3)` was left unevaluated. The rule that splits
602+
`w^m/(a + b w^n)` over the `n` roots of the binomial, for an even quadratic `w`, took a constant beside the power
603+
only. A second power of an even quadratic `v` is read as well. `w^m v^(n - 1 - m)/(a v^n + b w^n)` is the same
604+
sum in `t = w/v`, and each root's `1/(w - r v)` is again an even quadratic
605+
([#718](https://github.com/asc-community/AngouriMath/issues/718)).
606+
607+
| Input | Was (2.5.0) | Now |
608+
|---|---|---|
609+
| `"(1+x^2)^2/(a*x^6+b*(1+x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 878 characters |
610+
| `"x^2*(1+x^2)/(a*x^6+b*(1+x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 896 characters |
611+
| `"(1+x^2)/(a*x^4+b*(1+x^2)^2)".ToEntity().Integrate("x")` | `integral(...)` | 289 characters |
612+
599613
### A symbolic cubic with a repeated root is written over its roots
600614

601615
**Answers where there were none.** `1/(x^3 - 3b x + 2b^(3/2))` was left unevaluated: its coefficients are symbolic, so

‎Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs‎

Lines changed: 94 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -20066,6 +20066,14 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres
2006620066
/// <c>P</c> named, and splits on the sign of <c>P q</c>, which means nothing for a complex
2006720067
/// <c>P</c>: its first arm, the arctangent, differentiates back for any <c>P q != 0</c>, and
2006820068
/// the answer is checked below. Rubi's 1.3.1 <c>(1 - x^2)/(a + b (1 - x^2)^4)</c> was declined.
20069+
/// <para>
20070+
/// Two powers of even quadratics, <c>w^m v^(n - 1 - m)/(a v^n + b w^n)</c>, are the same sum in
20071+
/// <c>t = w/v</c>: the integrand is <c>t^m/(v (a + b t^n))</c>, and <c>1/(v (t - r_j))</c> is
20072+
/// <c>1/(w - r_j v)</c>, again an even quadratic <c>P + Q x^2</c>, whose integral
20073+
/// <c>arctan(sqrt(Q) x/sqrt(P))/(sqrt(P) sqrt(Q))</c> differentiates back for any <c>P</c> and
20074+
/// <c>Q</c> off zero. Rubi's 1.3.1 <c>(1 + x^2)^2/(a x^6 + b (1 + x^2)^3)</c>, a sum of two cubes,
20075+
/// was declined.
20076+
/// </para>
2006920077
/// https://github.com/asc-community/AngouriMath/issues/718
2007020078
/// </remarks>
2007120079
internal static Entity? SolveAPowerOfAnEvenQuadraticOverABinomialInIt(Entity expr, Entity.Variable x, bool integrateByParts)
@@ -20081,58 +20089,98 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres
2008120089
{
2008220090
if (!factor.ContainsNode(x))
2008320091
constant = under ? constant / factor : constant * factor;
20084-
else if (!under && above is null)
20085-
above = factor;
20086-
else if (under && below is null)
20092+
else if (!under)
20093+
above = above is null ? factor : above * factor;
20094+
else if (below is null)
2008720095
below = factor;
2008820096
else
2008920097
return null;
2009020098
}
2009120099
if (below is null)
2009220100
return null;
20093-
// a + b w^n: one term in x, a whole power of an even quadratic, beside terms free of it.
20094-
Entity? inner = null;
20095-
Entity scale = Number.Integer.One;
20101+
// a + b w^n: one term in x, a whole power of an even quadratic, beside terms free of it; or
20102+
// a v^n + b w^n, two such powers, a power of x^2 among them written as one of x.
20103+
var powers = new List<(Entity Scale, Entity Base, int Degree)>();
2009620104
Entity free = Number.Integer.Zero;
20097-
var degree = 0;
2009820105
foreach (var term in Sumf.LinearChildren(below))
2009920106
{
2010020107
if (!term.ContainsNode(x))
2010120108
{
2010220109
free = free + term;
2010320110
continue;
2010420111
}
20105-
if (inner is not null)
20112+
if (powers.Count == 2)
2010620113
return null;
20114+
Entity termScale = Number.Integer.One;
20115+
Entity? termBase = null;
20116+
var termDegree = 0;
2010720117
foreach (var factor in Mulf.LinearChildren(term))
2010820118
{
2010920119
if (!factor.ContainsNode(x))
20110-
scale = scale * factor;
20111-
else if (inner is null && factor is Powf(var w, Number.Integer { EInteger: var n }) && n.CompareTo(EInteger.FromInt32(2)) >= 0
20112-
&& n.CompareTo(EInteger.FromInt32(8)) <= 0)
20120+
termScale = termScale * factor;
20121+
else if (termBase is null && factor is Powf(var w, Number.Integer { EInteger: var n }) && n.CompareTo(EInteger.FromInt32(2)) >= 0
20122+
&& n.CompareTo(EInteger.FromInt32(16)) <= 0)
2011320123
{
20114-
inner = w;
20115-
degree = n.ToInt32Unchecked();
20124+
termBase = w;
20125+
termDegree = n.ToInt32Unchecked();
2011620126
}
2011720127
else
2011820128
return null;
2011920129
}
20130+
if (termBase is null)
20131+
return null;
20132+
powers.Add((termScale, termBase, termDegree));
20133+
}
20134+
if (powers.Count == 2)
20135+
{
20136+
if (!VanishesIdentically(free))
20137+
return null;
20138+
// x^(2 n) beside v^n is (x^2)^n.
20139+
for (var k = 0; k < 2; k++)
20140+
if (powers[k].Base == x && powers[k].Degree % 2 == 0 && powers[1 - k].Base != x)
20141+
powers[k] = (powers[k].Scale, MathS.Pow(x, Number.Integer.Create(2)), powers[k].Degree / 2);
2012020142
}
20121-
if (inner is null || VanishesIdentically(free) || Functions.PartialFractions.HoldsTheImaginaryUnit(expr)
20143+
if (powers.Count == 0 || powers.Any(power => power.Degree > 8) || powers.Count == 2 && powers[0].Degree != powers[1].Degree)
20144+
return null;
20145+
var (scale, inner, degree) = powers[0];
20146+
Entity? other = powers.Count == 2 ? powers[1].Base : null;
20147+
if (other is not null)
20148+
free = powers[1].Scale;
20149+
Entity otherP = Number.Integer.One, otherQ = Number.Integer.Zero;
20150+
if (VanishesIdentically(free) || Functions.PartialFractions.HoldsTheImaginaryUnit(expr)
2012220151
|| !TreeAnalyzer.TryGetPolyQuadratic(inner, x, out var q, out var linearTerm, out var p)
2012320152
|| !VanishesIdentically(linearTerm) || VanishesIdentically(q))
2012420153
return null;
20125-
// The numerator a constant multiple of w^m, m below n.
20154+
if (other is not null)
20155+
{
20156+
if (!TreeAnalyzer.TryGetPolyQuadratic(other, x, out var readQ, out var otherLinear, out var readP)
20157+
|| !VanishesIdentically(otherLinear) || VanishesIdentically(readQ))
20158+
return null;
20159+
(otherP, otherQ) = (readP, readQ);
20160+
}
20161+
// The numerator a constant multiple of w^m v^(n - 1 - m), m below n, which is w^m alone
20162+
// for a binomial: then t^m/(V (a + b t^n)) for t = w/v, split over the roots in t, and
20163+
// 1/(V (t - r)) is 1/(w - r v).
2012620164
var power = -1;
2012720165
Entity multiple = Number.Integer.One;
20128-
if (above is null)
20166+
if (above is null && other is null)
2012920167
power = 0;
20130-
else
20131-
for (var m = 1; m < degree && power < 0; m++)
20168+
else if (above is not null)
20169+
for (var m = other is null ? 1 : 0; m < degree && power < 0; m++)
2013220170
{
2013320171
var wm = MathS.Pow(inner, Number.Integer.Create(m));
20134-
if (!TreeAnalyzer.TryGetPolynomial(above, x, out var aboveRead) || !TreeAnalyzer.TryGetPolynomial(wm, x, out var wmRead)
20135-
|| !aboveRead.TryGetValue(EInteger.FromInt32(2 * m), out var lead) || !wmRead.TryGetValue(EInteger.FromInt32(2 * m), out var wLead))
20172+
if (other is not null)
20173+
{
20174+
// No zeroth powers: InnerSimplified writes w^0 as 1 provided w is not 0.
20175+
var rest = degree - 1 - m;
20176+
wm = m == 0 ? MathS.Pow(other, Number.Integer.Create(rest))
20177+
: rest == 0 ? wm : wm * MathS.Pow(other, Number.Integer.Create(rest));
20178+
}
20179+
var top = other is null ? 2 * m : 2 * (degree - 1);
20180+
// Expanded for two powers: the reader takes a sum of monomials, and v^k is not one.
20181+
var (aboveAsRead, wmAsRead) = other is null ? (above, wm) : (above.Expand(), wm.Expand());
20182+
if (!TreeAnalyzer.TryGetPolynomial(aboveAsRead, x, out var aboveRead) || !TreeAnalyzer.TryGetPolynomial(wmAsRead, x, out var wmRead)
20183+
|| !aboveRead.TryGetValue(EInteger.FromInt32(top), out var lead) || !wmRead.TryGetValue(EInteger.FromInt32(top), out var wLead))
2013620184
continue;
2013720185
var ratio = (lead / wLead).InnerSimplified;
2013820186
if (!VanishesIdentically((above - ratio * wm).Expand()))
@@ -20144,16 +20192,20 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres
2014420192
return null;
2014520193
var a = free.InnerSimplified;
2014620194
var b = scale.InnerSimplified;
20195+
Entity? each = null;
2014720196
var namedP = Variable.CreateUnique(expr, "P_b");
20148-
if (Integration.ComputeAsAQuestionOfItsOwn(1 / (namedP + q * MathS.Pow(x, Number.Integer.Create(2))), x, integrateByParts) is not { } each)
20149-
return null;
20150-
// See the remarks: the arctangent arm for every root.
20151-
each = each.Replace(node => node switch
20197+
if (other is null)
2015220198
{
20153-
Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expression,
20154-
Providedf provided => provided.Expression,
20155-
_ => node,
20156-
});
20199+
if (Integration.ComputeAsAQuestionOfItsOwn(1 / (namedP + q * MathS.Pow(x, Number.Integer.Create(2))), x, integrateByParts) is not { } asked)
20200+
return null;
20201+
// See the remarks: the arctangent arm for every root.
20202+
each = asked.Replace(node => node switch
20203+
{
20204+
Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expression,
20205+
Providedf provided => provided.Expression,
20206+
_ => node,
20207+
});
20208+
}
2015720209
var radius = MathS.Pow(Functions.PartialFractions.InLowestTermsOverTheSymbols(-a / b), Number.Rational.Create(1, degree));
2015820210
Entity sum = Number.Integer.Zero;
2015920211
for (var j = 0; j < degree; j++)
@@ -20162,7 +20214,20 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres
2016220214
// the real root, and (-1)^(2/3) is 1 (#204).
2016320215
var turn = 2 * MathS.pi * Number.Rational.Create(j, degree);
2016420216
Entity root = j == 0 ? radius : radius * (MathS.Cos(turn) + MathS.i * MathS.Sin(turn)).InnerSimplified;
20165-
sum = sum + MathS.Pow(root, Number.Integer.Create(power + 1)) * each.Substitute(namedP, p - root);
20217+
Entity over;
20218+
if (each is not null)
20219+
over = each.Substitute(namedP, p - root);
20220+
else
20221+
{
20222+
// 1/(P + Q x^2) for w - r v: arctan(sqrt(Q) x/sqrt(P))/(sqrt(P) sqrt(Q)), whose
20223+
// derivative is 1/(P + Q x^2) for any P, Q off zero, the two roots squaring back.
20224+
var atP = Functions.PartialFractions.InLowestTermsOverTheSymbols(p - root * otherP);
20225+
var atQ = Functions.PartialFractions.InLowestTermsOverTheSymbols(q - root * otherQ);
20226+
if (VanishesIdentically(atP) || VanishesIdentically(atQ))
20227+
return null;
20228+
over = MathS.Arctan(MathS.Sqrt(atQ) * x / MathS.Sqrt(atP)) / (MathS.Sqrt(atP) * MathS.Sqrt(atQ));
20229+
}
20230+
sum = sum + MathS.Pow(root, Number.Integer.Create(power + 1)) * over;
2016620231
}
2016720232
var answer = (-constant * multiple / (degree * a) * sum).InnerSimplified;
2016820233
return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
//
2+
// Copyright (c) 2019-2026 Angouri.
3+
// AngouriMath is licensed under MIT.
4+
// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
5+
// Website: https://am.angouri.org.
6+
//
7+
8+
using System;
9+
using AngouriMath.Extensions;
10+
using Xunit;
11+
12+
namespace AngouriMath.Tests.Calculus
13+
{
14+
/// <summary>
15+
/// A quotient over a sum of two whole powers of even quadratics, <c>a v^n + b w^n</c>, split over the
16+
/// roots of <c>a + b t^n</c> in <c>t = w/v</c>: Rubi's 1.3.1 <c>(1 + x^2)^2/(a x^6 + b (1 + x^2)^3)</c>,
17+
/// which was left unevaluated.
18+
/// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
19+
/// </summary>
20+
[Trait("Area", "Calculus")]
21+
public sealed class ASumOfTwoPowersOfEvenQuadraticsIntegralTest
22+
{
23+
private static void DifferentiatesBack(string integrand, params (string, double)[] pins)
24+
{
25+
var integral = integrand.ToEntity().Integrate("x");
26+
Assert.DoesNotContain("integral(", integral.Stringize());
27+
28+
var pinned = integral.Substitute("C", 0);
29+
Entity original = integrand.ToEntity();
30+
foreach (var (name, value) in pins)
31+
{
32+
pinned = pinned.Substitute(name, value);
33+
original = original.Substitute(name, value);
34+
}
35+
var derivative = pinned.Differentiate("x");
36+
foreach (var at in new[] { -2.3, -0.7, 0.4, 1.9 })
37+
{
38+
var got = derivative.Substitute("x", at).EvalNumerical();
39+
var want = original.Substitute("x", at).EvalNumerical();
40+
var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart);
41+
var scale = Math.Max(1.0, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart));
42+
Assert.True(difference / scale < 1e-8,
43+
$"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
44+
}
45+
}
46+
47+
[Theory]
48+
[InlineData("(1+x^2)^2/(a*x^6+b*(1+x^2)^3)")]
49+
[InlineData("x^2*(1+x^2)/(a*x^6+b*(1+x^2)^3)")]
50+
[InlineData("x^4/(a*x^6+b*(1+x^2)^3)")]
51+
[InlineData("(2+3*x^2)^2/(a*x^6+b*(2+3*x^2)^3)")]
52+
[InlineData("(1+x^2)/(a*x^4+b*(1+x^2)^2)")]
53+
public void IsAnsweredAndDifferentiatesBack(string integrand)
54+
{
55+
DifferentiatesBack(integrand, ("a", 2), ("b", 3));
56+
DifferentiatesBack(integrand, ("a", -5), ("b", 0.5));
57+
}
58+
}
59+
}

0 commit comments

Comments
 (0)