diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 38c1dc016..59087d3a0 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -154,6 +154,18 @@ integrand is then `x^(2n + 1) (c + d x)^n (2 c + 3 d x)`, the derivative of `x^( | `"x^(2*n)*(c+d*x)^n*(2*c*x+3*d*x^2)".ToEntity().Integrate("x")` | `integral(...)` | 61 characters | | `"x^(n-1)*(c+d*x^2)^n*(c*x+3*d*x^3)".ToEntity().Integrate("x")` | `integral(...)` | 65 characters | +### A power of a quotient of linears beside whole powers of them is integrated under the quotient + +**Answers where there were none.** `e^(n acoth(a x))/(c - a^2 c x^2)^2` was declined, and so were its `^3` and `^4`, Rubi's +7.4.2. The parser folds `e^(n acoth(a x))` into `((a x + 1)/(a x - 1))^(n/2)`, and under `t = (a x + 1)/(a x - 1)` -- the +Moebius substitution, `p + q x = t (r + s x)` and `dx = (q r - p s)/(q - s t)^2 dt` -- a whole power of the quadratic +is a polynomial in `t` times a power of it ([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"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` | +| `"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` | + ### The last digits of an evaluated sine, cosine, tangent, arcsine or arctangent **Silent, and in the hundredth digit.** `EvalNumerical` at the default hundred digits gives the diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 127392b67..c372e4b2b 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -709,6 +709,109 @@ private static bool IsADerivativeOfAFactorBelowTheBar(Entity expr, Entity.Variab return (constant * (MathS.Pow(x, raised) * logarithm - 2 * slope * integral) / raised).InnerSimplified; } + /// + /// A power of a quotient of two linears beside whole powers of the two, under the Moebius + /// substitution t = (p + q x)/(r + s x): then p + q x = t (r + s x), + /// r + s x = D/(q - s t) and dx = D/(q - s t)^2 dt with D = q r - p s, so + /// t^alpha (p + q x)^j (r + s x)^k dx is D^(j + k + 1) t^(alpha + j) (q - s t)^(-(j + k + 2)) dt, + /// a polynomial in t times a power of it wherever j + k <= -2. A whole power of + /// a quadratic that is a multiple of the two linears' product counts as a power of each. + /// e^(n acoth(a x)) arrives folded as ((a x + 1)/(a x - 1))^(n/2), and over + /// (c - a^2 c x^2)^m it is Rubi's 7.4.2. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + /// Only for a power alpha that is not whole: a whole one makes a rational function, + /// which the partial fractions answer. t^gamma is differentiated to + /// gamma t^(gamma - 1) t' on the principal branch wherever t is not 0, the + /// negative half-line included, so the answer written back in the quotient holds where + /// the integrand is defined. A half-odd power of the quadratic is left alone: it is not + /// the product of the linears' powers on every side of their roots. + /// + internal static Entity? SolveByTheMoebiusSubstitutionOfAQuotientOfLinears(Entity expr, Entity.Variable x) + { + Powf? power = null; + Entity constant = Number.Integer.One; + var others = new List<(Entity Base, EInteger Power)>(); + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (!underneath && power is null && factor is Powf(Divf(var above, var below), var alpha) + && !alpha.ContainsNode(x) && alpha.Evaled is not Number.Integer + && above.ContainsNode(x) && below.ContainsNode(x)) + { + power = (Powf)factor; + continue; + } + var (@base, whole) = factor is Powf(var b, Number.Integer n) && n.EInteger.CanFitInInt32() ? (b, n.EInteger) : (factor, EInteger.One); + others.Add((@base, underneath ? -whole : whole)); + } + if (power is null || others.Count == 0 || power.Base is not Divf(var top, var bottom) + || !TreeAnalyzer.TryGetPolyLinear(top, x, out var q, out var p) || TreeAnalyzer.IsZero(q) + || !TreeAnalyzer.TryGetPolyLinear(bottom, x, out var s, out var r) || TreeAnalyzer.IsZero(s)) + return null; + var determinant = (q * r - p * s).InnerSimplified; + if (VanishesIdentically(determinant)) + return null; + // Each factor a multiple of the top, of the bottom, or of their product. + var j = EInteger.Zero; + var k = EInteger.Zero; + foreach (var (@base, n) in others) + { + Entity multiple; + if (TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out var intercept)) + { + if (VanishesIdentically(slope * p - intercept * q)) + { + multiple = (slope / q).InnerSimplified; + j += n; + } + else if (VanishesIdentically(slope * r - intercept * s)) + { + multiple = (slope / s).InnerSimplified; + k += n; + } + else + return null; + } + else if (TreeAnalyzer.TryGetPolyQuadratic(@base, x, out var a2, out var b2, out var c2) && !TreeAnalyzer.IsZero(a2)) + { + multiple = (a2 / (q * s)).InnerSimplified; + if (!VanishesIdentically(b2 - multiple * (p * s + q * r)) || !VanishesIdentically(c2 - multiple * p * r)) + return null; + j += n; + k += n; + } + else + return null; + constant *= n.Equals(EInteger.One) ? multiple : MathS.Pow(multiple, Number.Integer.Create(n)); + } + var degree = (j + k + EInteger.FromInt32(2)).Negate(); + if (degree.Sign < 0 || degree.CompareTo(EInteger.FromInt32(24)) > 0 || !j.CanFitInInt32()) + return null; + // D^(j + k + 1) t^(alpha + j) (q - s t)^degree, term by term in t. + var m = degree.ToInt32Checked(); + Entity sum = Number.Integer.Zero; + Entity binomial = Number.Integer.One; + for (var i = 0; i <= m; i++) + { + var offset = j + EInteger.FromInt32(i + 1); + var raised = (offset.IsZero ? power.Exponent : power.Exponent + Number.Integer.Create(offset)).InnerSimplified; + if (raised.Evaled is Number.Complex { IsZero: true }) + return null; + // Free of x and small, so simplified: the multiples are otherwise left as + // `(1/(a (-a)))^(-2) (a - (-a))^(-3)`. + var coefficient = (binomial * MathS.Pow(q, Number.Integer.Create(m - i)) * MathS.Pow(-s, Number.Integer.Create(i)) / raised).Simplify(); + sum += coefficient * MathS.Pow(power.Base, raised); + binomial = binomial * (m - i) / (i + 1); + } + return ((constant * MathS.Pow(determinant, Number.Integer.Create(j + k + EInteger.One))).Simplify() * sum).InnerSimplified; + } + /// /// A quotient by the square of a base, by parts against the base's reciprocal, where the /// base's derivative makes the rest exact: N/(M v^2) is g v'/v^2 with diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index f5645285e..3ceb0c6f3 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -784,6 +784,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) if ((answer = IndefiniteIntegralSolver.SolveAPowerOverAPowerOfALogarithm(expr, x)) is { }) return answer; // And x^m artanh(A x/sqrt(d + A^2 x^2)), by parts, the logarithm's derivative read as the root it is. if ((answer = IndefiniteIntegralSolver.SolveByPartsAgainstAnArtanhWhoseDerivativeIsARoot(expr, x, integrateByParts)) is { }) return answer; + // A power of a quotient of linears beside whole powers of them, which is how e^(n acoth(a x)) + // over a power of c - a^2 c x^2 arrives folded, under t = the quotient. + if ((answer = IndefiniteIntegralSolver.SolveByTheMoebiusSubstitutionOfAQuotientOfLinears(expr, x)) is { }) return answer; // A power of an exponential with a positive base is the exponential of the product, // exactly, and only that spelling is one the exponential rules read. if ((answer = IndefiniteIntegralSolver.SolveByFlatteningAPowerOfAnExponential(expr, x, integrateByParts)) is { }) return answer; diff --git a/Sources/Tests/UnitTests/Calculus/APowerOfAQuotientOfLinearsBesideThemIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerOfAQuotientOfLinearsBesideThemIntegralTest.cs new file mode 100644 index 000000000..96ce18423 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APowerOfAQuotientOfLinearsBesideThemIntegralTest.cs @@ -0,0 +1,54 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A power of a quotient of two linears beside whole powers of them, under + /// t = (p + q x)/(r + s x): e^(n acoth(a x))/(c - a^2 c x^2)^m arrives as + /// ((a x + 1)/(a x - 1))^(n/2) over the power of the quadratic, and is a polynomial in + /// t times a power of it. For m past 1 it was declined. Rubi's 7.4.2. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class APowerOfAQuotientOfLinearsBesideThemIntegralTest + { + [Theory] + [InlineData("e^(n*acoth(a*x))/(c - a^2*c*x^2)^2")] + [InlineData("e^(n*acoth(a*x))/(c - a^2*c*x^2)^3")] + [InlineData("e^(n*atanh(a*x))/(c - a^2*c*x^2)^2")] + [InlineData("((a*x + 1)/(a*x - 1))^(n/2)/((a*x + 1)^2*(a*x - 1)^2)")] + [InlineData("((x + 2)/(3*x - 1))^(1/3)/((x + 2)*(3*x - 1)^2)")] + public void UnderTheQuotient(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("c", 0.7).Substitute("n", 0.9); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + // On both sides of the roots, -0.77 and 0.77 for a = 1.3, and of 1/3 and -2. + foreach (var at in new[] { -2.9, -1.4, -0.4, 0.2, 0.6, 1.3, 2.9 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +}