Skip to content

Commit 99fca7a

Browse files
Integrate a power of x beside powers of a quotient of linears and of their product as one term, and 1/(a + b asinh(+-i + d x^2))^n (#2032)
- x^j T^alpha Q^beta, for T = (1 + a x)/(1 - a x) and Q = c - a^2 c x^2, as T^alpha Q^(beta + 1) times a polynomial of degree j - 1, found by one linear solve whose last equation is checked; where it fails the rule declines. Rubi's 7.3.6, e^(n artanh(a x)) x^2 (c - a^2 c x^2)^p, the parser having folded the exponential into T^(n/2). - The six of Rubi's 7.1.5 1/(a + b asinh(+-i + d x^2))^n by the rule for the inverse hyperbolic sine of a quadratic, with g = cosh(u/2 + c pi/4): cosh(u/2) + c sinh(u/2) is sqrt(2) times it, and its square is 1 + c sinh(u) = c d x^2. They are complex at every x, and checked there. Chain against master f7f2a7a: pocket of 2,374 rows, 2,362 -> 2,372 solved, 0 wrong; family 7's sample 269 -> 270, the others unchanged; the 10 rows that moved, alone, 0 -> 10 solved; 7,894 rows compared for allocation, none grew past 109 MB. Part of #233. Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura Co-authored-by: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
1 parent fc1d521 commit 99fca7a

5 files changed

Lines changed: 227 additions & 4 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 23 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -228,6 +228,23 @@ is a polynomial in `t` times a power of it ([#233](https://github.com/asc-commun
228228
| `"e^(n*acoth(a*x))/(c-a^2*c*x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | three powers of `(a x + 1)/(a x - 1)`, `n/2 - 1` to `n/2 + 1` |
229229
| `"e^(n*acoth(a*x))/(c-a^2*c*x^2)^4".ToEntity().Integrate("x")` | `integral(...)` | seven powers of it, `n/2 - 3` to `n/2 + 3` |
230230

231+
### A power of x beside powers of a quotient of linears and of their product is integrated as one term
232+
233+
**Answers where there were none.** `e^(n artanh(a x)) x^2 (c - a^2 c x^2)^(-1 - n^2/2)` was declined, and so were
234+
`e^(6 artanh(a x)) x^2/(c - a^2 c x^2)^19`, `e^(5 artanh(a x)) x^2/(c - a^2 c x^2)^(27/2)` and its reciprocal twin, and
235+
`e^(artanh(a x)/2) x^2/(c - a^2 c x^2)^(9/8)`, Rubi's 7.3.6. With `T = (1 + a x)/(1 - a x)` and `Q = c - a^2 c x^2`,
236+
each is `x^j T^alpha Q^beta`, and its antiderivative is `T^alpha Q^(beta + 1)` times a polynomial of degree `j - 1`,
237+
found by one linear solve whose last equation is checked; where it fails the rule declines
238+
([#233](https://github.com/asc-community/AngouriMath/issues/233)).
239+
240+
| Input | Was (2.5.0) | Now |
241+
|---|---|---|
242+
| `"e^(n*atanh(a*x))*x^2*(c-a^2*c*x^2)^(-1-1/2*n^2)".ToEntity().Integrate("x")` | `integral(...)` | 167 characters |
243+
| `"e^(6*atanh(a*x))*x^2/(c-a^2*c*x^2)^19".ToEntity().Integrate("x")` | `integral(...)` | `T^3 Q^(-18) (6 a x - 1)/(210 a^3 c)`, 119 characters |
244+
| `"e^(5*atanh(a*x))*x^2/(c-a^2*c*x^2)^(27/2)".ToEntity().Integrate("x")` | `integral(...)` | 125 characters |
245+
| `"e^(1/2*atanh(a*x))*x^2/(c-a^2*c*x^2)^(9/8)".ToEntity().Integrate("x")` | `integral(...)` | 121 characters |
246+
| `"e^(2*atanh(a*x))*x^2/(c-a^2*c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 109 characters |
247+
231248
### The last digits of an evaluated sine, cosine, tangent, arcsine or arctangent
232249

233250
**Silent, and in the hundredth digit.** `EvalNumerical` at the default hundred digits gives the
@@ -1458,7 +1475,10 @@ values, and too few points could be compared. The answer is now refused only whe
14581475
a cosine `g`, so `x/g` is constant on each side of 0; the powers `-1` and `-2` are then sine and cosine integrals of
14591476
`(a + b u)/(2 b)`, hyperbolic ones for the inverse hyperbolic cosine, and the recurrence run backwards gives the
14601477
rest. The answers hold on the real domain of the integrand
1461-
([#718](https://github.com/asc-community/AngouriMath/issues/718)).
1478+
([#718](https://github.com/asc-community/AngouriMath/issues/718)). The six of Rubi's 7.1.5
1479+
`1/(a + b asinh(+-i + d x^2))^n` are answered the same way, with `g = cosh(u/2 + c pi/4)`: `cosh(u/2) + c sinh(u/2)` is
1480+
`sqrt(2)` times it, and its square is `1 + c sinh(u) = c d x^2`. They are complex at every `x`, and are checked there
1481+
([#233](https://github.com/asc-community/AngouriMath/issues/233)).
14621482

14631483
| Input | Was (2.5.0) | Now |
14641484
|---|---|---|
@@ -1468,6 +1488,8 @@ rest. The answers hold on the real domain of the integrand
14681488
| `"1/(a+b*asin(-1+d*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 224 characters |
14691489
| `"1/(a+b*acosh(1+d*x^2))^2".ToEntity().Integrate("x")` | `integral(...)` | 608 characters |
14701490
| `"1/(a+b*acosh(-1+d*x^2))^3".ToEntity().Integrate("x")` | `integral(...)` | 633 characters |
1491+
| `"1/(a+b*asinh(i+d*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 482 characters |
1492+
| `"1/(a+b*asinh(-i+d*x^2))^3".ToEntity().Integrate("x")` | `integral(...)` | 735 characters |
14711493

14721494
### The reciprocal of a binomial in a high power of a sine, cosine or hyperbolic one is split over its roots
14731495

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

Lines changed: 131 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1057,6 +1057,124 @@ private static bool IsADerivativeOfAFactorBelowTheBar(Entity expr, Entity.Variab
10571057
return ((constant * MathS.Pow(determinant, Number.Integer.Create(j + k + EInteger.One))).Simplify() * sum).InnerSimplified;
10581058
}
10591059

1060+
/// <summary>
1061+
/// <c>x^j T^alpha Q^beta</c>, for a power of a quotient of linears <c>T = (p + q x)/(r + s x)</c>
1062+
/// and one of a multiple of their product, <c>Q = kappa (p + q x)(r + s x)</c>, where the
1063+
/// antiderivative is algebraic: <c>T^alpha Q^(beta + 1) P(x)</c> with <c>P</c> of degree
1064+
/// <c>j - 1</c>. With <c>D = q r - p s</c>, <c>T'/T = kappa D/Q</c>, so the derivative is
1065+
/// <c>T^alpha Q^beta (alpha kappa D P + (beta + 1) Q' P + Q P')</c>, and asking the bracket
1066+
/// to be <c>x^j</c> is, with <c>Q = A x^2 + B x + C</c>, the three-term system
1067+
/// <c>A (m + 2 beta + 1) p_(m-1) + (alpha kappa D + B (m + beta + 1)) p_m + C (m + 1) p_(m+1) = [m = j]</c>
1068+
/// for <c>m</c> from 0 to <c>j</c>: solved from the top, <c>p_(j-1) = 1/(A (j + 2 beta + 1))</c>,
1069+
/// down to <c>p_0</c>, and the equation at <c>m = 0</c> checked. Rubi's 7.3.6
1070+
/// <c>e^(n artanh(a x)) x^2 (c - a^2 c x^2)^(-1 - n^2/2)</c> is
1071+
/// <c>e^(n artanh(a x)) (1 - a n x)/(a^3 c n (1 - n^2) (c - a^2 c x^2)^(n^2/2))</c>.
1072+
/// https://github.com/asc-community/AngouriMath/issues/233
1073+
/// </summary>
1074+
/// <remarks>
1075+
/// <para>
1076+
/// One equation more than unknowns: where the check fails the antiderivative is not this
1077+
/// one term -- it has a logarithm, or is hypergeometric -- and the rule declines. Where it
1078+
/// holds, the rules after it ran past the budget: <c>e^(6 artanh(a x)) x^2/(c - a^2 c x^2)^19</c>
1079+
/// was written as a rational function over a denominator of degree 38 and handed to the
1080+
/// partial fractions, and is <c>(6 a x - 1)/(210 a^3 c^19 (1 - a x)^21 (1 + a x)^15)</c>;
1081+
/// <c>e^(5 artanh(a x)) x^2/(c - a^2 c x^2)^(27/2)</c> the same, and
1082+
/// <c>e^(artanh(a x)/2) x^2/(c - a^2 c x^2)^(9/8)</c> was declined. Where the denominator
1083+
/// is small the partial fractions answer it too, as a sum over its two roots:
1084+
/// <c>e^(2 artanh(a x)) x^2/(c - a^2 c x^2)^3</c> is here <c>(2 a x - 1) T Q^(-2)/(6 a^3 c)</c>.
1085+
/// The Moebius substitution after it takes a power of <c>Q</c> that fixes the degree of
1086+
/// <c>P</c>; this takes the degree from <c>x^j</c>.
1087+
/// </para>
1088+
/// <para>
1089+
/// The powers are differentiated as principal powers wherever <c>T</c> and <c>Q</c> are not 0,
1090+
/// on a cut as off it, so the answer holds where the integrand is defined, and is checked
1091+
/// at sampled points as well. Only for <c>x</c> above the bar: without it <c>P</c> has no
1092+
/// degree to take.
1093+
/// </para>
1094+
/// </remarks>
1095+
internal static Entity? SolveAPowerOfXBesidePowersOfAQuotientOfLinearsAndTheirProductAlgebraically(Entity expr, Entity.Variable x)
1096+
{
1097+
Entity? quotient = null, alpha = null, quadratic = null, beta = null;
1098+
Entity constant = Number.Integer.One;
1099+
var j = 0;
1100+
foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr))
1101+
{
1102+
if (!factor.ContainsNode(x))
1103+
{
1104+
constant = underneath ? constant / factor : constant * factor;
1105+
continue;
1106+
}
1107+
if (!underneath && APowerOfTheVariable(factor, x) is { } k && k > 0)
1108+
{
1109+
j += k;
1110+
continue;
1111+
}
1112+
if (quotient is null && AsAPowerOfAQuotient(factor, x) is var (written, exponent))
1113+
{
1114+
(quotient, alpha) = (written, underneath ? -exponent : exponent);
1115+
continue;
1116+
}
1117+
var (@base, power) = factor is Powf(var b, var e) ? (b, e) : (factor, (Entity)Number.Integer.One);
1118+
if (quadratic is null && !power.ContainsNode(x) && @base is not Divf)
1119+
{
1120+
(quadratic, beta) = (@base, underneath ? -power : power);
1121+
continue;
1122+
}
1123+
return null;
1124+
}
1125+
if (quotient is not Divf(var top, var bottom) || alpha is null || quadratic is null || beta is null || j < 1 || j > 12
1126+
|| !TreeAnalyzer.TryGetPolyLinear(top, x, out var q, out var p) || TreeAnalyzer.IsZero(q)
1127+
|| !TreeAnalyzer.TryGetPolyLinear(bottom, x, out var s, out var r) || TreeAnalyzer.IsZero(s)
1128+
|| !TreeAnalyzer.TryGetPolyQuadratic(quadratic, x, out var a2, out var b2, out var c2) || TreeAnalyzer.IsZero(a2))
1129+
return null;
1130+
var determinant = (q * r - p * s).InnerSimplified;
1131+
var kappa = Functions.PartialFractions.Bare((a2 / (q * s)).Simplify());
1132+
if (VanishesIdentically(determinant)
1133+
|| !VanishesIdentically(b2 - kappa * (p * s + q * r)) || !VanishesIdentically(c2 - kappa * p * r))
1134+
return null;
1135+
(alpha, beta) = (alpha.InnerSimplified, beta.InnerSimplified);
1136+
var alphaKappaD = (alpha * kappa * determinant).InnerSimplified;
1137+
// p_0 .. p_(j - 1), and p_j = 0, from the top down.
1138+
var coefficients = new Entity[j + 1];
1139+
coefficients[j] = Number.Integer.Zero;
1140+
for (var m = j; m >= 1; m--)
1141+
{
1142+
var divisor = (a2 * (m + 2 * beta + 1).Expand()).InnerSimplified;
1143+
if (VanishesIdentically(divisor))
1144+
return null;
1145+
Entity rest = m == j ? Number.Integer.One
1146+
: -((alphaKappaD + b2 * (m + beta + 1).Expand()) * coefficients[m] + c2 * (m + 1) * coefficients[m + 1]);
1147+
coefficients[m - 1] = Functions.PartialFractions.Bare((rest / divisor).Simplify());
1148+
}
1149+
var second = j >= 2 ? coefficients[1] : Number.Integer.Zero;
1150+
if (!VanishesIdentically((alphaKappaD + b2 * (beta + 1)) * coefficients[0] + c2 * second))
1151+
return null;
1152+
Entity polynomial = Number.Integer.Zero;
1153+
for (var m = 0; m < j; m++)
1154+
if (!VanishesIdentically(coefficients[m]))
1155+
polynomial += coefficients[m] * (m == 0 ? Number.Integer.One : m == 1 ? x : MathS.Pow(x, m));
1156+
var answer = (constant * MathS.Pow(quotient, alpha) * MathS.Pow(quadratic, (beta + 1).Expand().InnerSimplified) * polynomial).InnerSimplified;
1157+
return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
1158+
1159+
// T^alpha as it is written, or as e^(alpha ln(T)) where it is not folded yet.
1160+
static (Entity Quotient, Entity Exponent)? AsAPowerOfAQuotient(Entity factor, Entity.Variable x)
1161+
{
1162+
if (factor is Powf(Divf(var above, var below) quotient, var exponent) && !exponent.ContainsNode(x)
1163+
&& above.ContainsNode(x) && below.ContainsNode(x))
1164+
return (quotient, exponent);
1165+
if (factor is not Powf(var e, var power) || e != MathS.e || !power.ContainsNode(x))
1166+
return null;
1167+
var logarithms = power.Nodes.OfType<Logf>().Where(node => node.Base == MathS.e && node.ContainsNode(x)).Distinct().ToList();
1168+
if (logarithms.Count != 1 || logarithms[0].Antilogarithm is not Divf(var top, var bottom) || !top.ContainsNode(x) || !bottom.ContainsNode(x))
1169+
return null;
1170+
var placeholder = Variable.CreateUnique(power, "l_quot");
1171+
if (!TreeAnalyzer.TryGetPolyLinear(power.Replace(node => node == logarithms[0] ? placeholder : node), placeholder, out var multiple, out var free)
1172+
|| multiple.ContainsNode(x) || !TreeAnalyzer.IsZero(free))
1173+
return null;
1174+
return (logarithms[0].Antilogarithm, multiple);
1175+
}
1176+
}
1177+
10601178
/// <summary>
10611179
/// A function of the hyperbolic tangent of a linear, under <c>u = tanh(a + b x)</c>: the parser
10621180
/// writes <c>tanh(L)</c> as <c>(W - 1)/(W + 1)</c> and <c>coth(L)</c> as <c>(W + 1)/(W - 1)</c>
@@ -38724,15 +38842,18 @@ static bool HoldsAPythagoreanPair(Entity sum)
3872438842
/// negative power goes through <c>u</c>: <c>w - c</c> is a constant times <c>g(u/2 + beta)^2</c>
3872538843
/// for a sine or a cosine <c>g</c>, so <c>x/g</c> is constant on each side of 0, and the
3872638844
/// powers <c>-1</c> and <c>-2</c> are sine and cosine integrals of <c>(a + b u)/(2 b)</c>
38727-
/// (hyperbolic ones for the hyperbolic cosine), the rest the recurrence run backwards. Exact on
38728-
/// the real domain; off it <c>x/g</c> can change sign across a cut.
38845+
/// (hyperbolic ones for the hyperbolic cosine and sine), the rest the recurrence run backwards.
38846+
/// Exact on the real domain; off it <c>x/g</c> can change sign across a cut. For the inverse
38847+
/// hyperbolic sine <c>g</c> is <c>cosh(u/2 + c pi/4)</c>: <c>cosh(u/2) + c sinh(u/2)</c> is
38848+
/// <c>sqrt(2)</c> times it, and its square is <c>1 + c sinh(u) = c d x^2</c>.
3872938849
/// </summary>
3873038850
/// <remarks>
3873138851
/// Each <c>u</c> has <c>u'' = e w u'^3</c>, and with <c>c^2 = 1</c> (or <c>-1</c>) the radicand of
3873238852
/// <c>u'</c> is <c>d x^2 (2 c + d x^2)</c> up to sign, so <c>H u' = 2</c> and
3873338853
/// <c>(x (a + b u)^n + 2 e n b H (a + b u)^(n - 1))'</c> is <c>(a + b u)^n</c> plus a multiple of
3873438854
/// <c>(a + b u)^(n - 2)</c>. Rubi's 5.1.5, 5.2.5 and 7.2.5 <c>(a + b acos(1 + d x^2))^4</c>
38735-
/// and the like were declined.
38855+
/// and the like were declined, and 7.1.5 <c>1/(a + b asinh(i + d x^2))^n</c>, which had no
38856+
/// row in the table below.
3873638857
/// https://github.com/asc-community/AngouriMath/issues/718
3873738858
/// </remarks>
3873838859
internal static Entity? SolveAPowerOfAnInverseFunctionOfAQuadraticByTheRecurrence(Entity expr, Entity.Variable x)
@@ -38868,6 +38989,9 @@ static bool HoldsAPythagoreanPair(Entity sum)
3886838989
(Arcsinf, _, _, false) => (false, -MathS.pi / 4, false, false, -1),
3886938990
(Logf, -1, 1, true) => (true, Number.Integer.Zero, true, true, 1),
3887038991
(Logf, -1, 1, false) => (false, Number.Integer.Zero, true, false, 1),
38992+
// asinh with c = +-i: cosh(u/2) + c sinh(u/2) is sqrt(2) cosh(u/2 + c pi/4), whose
38993+
// square is 1 + c sinh(u) = c d x^2.
38994+
(Logf, -1, -1, _) => (false, c * MathS.pi / 4, true, false, 1),
3887138995
_ => null,
3887238996
};
3887338997
if (table is not var (gIsSine, beta, hyperbolic, cosine, hSign))
@@ -38898,6 +39022,10 @@ static bool HoldsAPythagoreanPair(Entity sum)
3889839022
integrals[k - 2] = (x * MathS.Pow(a_bu, k) + 2 * e * k * b * h * MathS.Pow(a_bu, k - 1) - integrals[k])
3889939023
/ (4 * e * MathS.Sqr(b) * k * (k - 1));
3890039024
var answer = integrals[n].InnerSimplified;
39025+
// With c = +-i the integrand is real nowhere, so it is compared where it is complex:
39026+
// the argument of the logarithm stays off its cut there, its imaginary part the sign of c.
39027+
if (squareAtZero == -1)
39028+
return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
3890139029
return Functions.PartialFractions.DerivativeDiffersWhereReal(answer, expr, x) ? null : answer;
3890239030
}
3890339031
}

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

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -698,6 +698,10 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
698698
// A power of a sine, cosine or hyperbolic one over a binomial in a power of it, over the
699699
// binomial's roots, before the rules that spent the budget on `sinh(y)^4/(a + b sinh(y)^3)`.
700700
if ((answer = IndefiniteIntegralSolver.SolveAPowerOverABinomialInAPowerOfASineOverItsRoots(expr, x)) is { }) return answer;
701+
// x^j T^alpha Q^beta for a quotient of linears T and a multiple Q of their product, where the
702+
// antiderivative is T^alpha Q^(beta + 1) times a polynomial of degree j - 1, before the rewrite
703+
// below hands `e^(6 artanh(a x)) x^2/(c - a^2 c x^2)^19` to partial fractions of degree 38.
704+
if ((answer = IndefiniteIntegralSolver.SolveAPowerOfXBesidePowersOfAQuotientOfLinearsAndTheirProductAlgebraically(expr, x)) is { }) return answer;
701705
// An exponential of a whole multiple of an inverse hyperbolic tangent, written algebraically,
702706
// before the rules that timed out on `e^(artanh(a x)) x^2/(c - a^2 c x^2)^2`.
703707
if ((answer = IndefiniteIntegralSolver.SolveAnExponentialOfAnInverseHyperbolicTangentAlgebraically(expr, x)) is { }) return answer;

‎Sources/Tests/UnitTests/Calculus/APowerOfAnInverseFunctionOfAQuadraticIntegralTest.cs‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -76,6 +76,19 @@ public void ANegativePowerIsAnsweredAndDifferentiatesBack(string integrand, int
7676
DifferentiatesBack(integrand, ("a", -0.75), ("b", -0.7), ("d", 0.15 * sign));
7777
}
7878

79+
// A negative power of the inverse hyperbolic sine of c + d x^2 with c = +-i, complex at every
80+
// x, in hyperbolic sine and cosine integrals: cosh(u/2) + c sinh(u/2) is sqrt(2) cosh(u/2 + c pi/4).
81+
[Theory]
82+
[InlineData("1/(a+b*asinh(i+d*x^2))", 1)]
83+
[InlineData("1/(a+b*asinh(i+d*x^2))^2", -1)]
84+
[InlineData("1/(a+b*asinh(-i+d*x^2))^3", 1)]
85+
[InlineData("1/(a+b*asinh(-i+d*x^2))^4", -1)]
86+
public void ANegativePowerOfAnInverseHyperbolicSineIsAnsweredAndDifferentiatesBack(string integrand, int sign)
87+
{
88+
DifferentiatesBack(integrand, ("a", 0.6), ("b", 1.25), ("d", 0.3 * sign));
89+
DifferentiatesBack(integrand, ("a", -0.75), ("b", -0.7), ("d", 0.15 * sign));
90+
}
91+
7992
// Half an odd power of the inverse hyperbolic cosine, in the error functions of
8093
// sqrt(a + b u)/sqrt(2 b).
8194
[Theory]

0 commit comments

Comments
 (0)