From ce381ad619445ce65ee3cebef566fdc5fc60e19b Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 20 Oct 2026 21:09:46 +0000 Subject: [PATCH] Integrate secant powers of a logarithm two apart as the derivative of a product With u = a + b ln(c x^n), B = b n and k = m + 1, the derivative of x^k sec(u)^(p - 2) (B (p - 2) tan(u) - k) is x^m (D sec(u)^p - E sec(u)^(p - 2)), where D = B^2 (p - 1)(p - 2) and E = B^2 (p - 2)^2 + k^2. So x^m (alpha sec(u)^p + beta sec(u)^(p - 2)) is a multiple of that product wherever beta D + alpha E = 0, and a lone power is one where E is zero. The cosecant's twin is the same, and k = 0 is left to the substitution. Rubi pocket (1,525 rows of family 4): 1480 -> 1487 solved, none lost. The moves re-run alone answer nine rows of 4.7.5: seven verified, and 4.7.5:337 and 395, whose integrands are complex along the whole real line, agree with Simpson's rule to 1e-13. Part of #233. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 18 +++ .../Integration/IndefiniteIntegralSolver.cs | 139 ++++++++++++++++++ .../Integration/Integration.Definition.cs | 3 + ...rtAsTheDerivativeOfAProductIntegralTest.cs | 74 ++++++++++ 4 files changed, 234 insertions(+) create mode 100644 Sources/Tests/UnitTests/Calculus/SecantPowersOfALogarithmTwoApartAsTheDerivativeOfAProductIntegralTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 977369a15..178cc5fe8 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -961,6 +961,24 @@ in Rubi's 1.1.3.4 come out shorter as well, and those that took two seconds take | `"1/(x^8*(a + b*x^3)*(c + d*x^3))".ToEntity().Integrate("x")` | `integral(...)` | 1,217 characters | | `"1/(x^6*(a + b*x^4)*(c + d*x^4))".ToEntity().Integrate("x")` | `integral(...)` | 795 characters | +### A secant or cosecant power of a logarithm beside the power two below it is integrated as the derivative of a product + +**Answers where there were none.** `-(1 + b^2 n^2) sec(a + b ln(c x^n)) + 2 b^2 n^2 sec(a + b ln(c x^n))^3` was +declined, and so were eight more of Rubi's 4.7.5, `sec(a + i ln(c x^n)/(n (p - 2)))^p` and the cosecant's twins among +them. With `u = a + b ln(c x^n)`, `B = b n` and `k = m + 1`, the derivative of `x^k sec(u)^(p - 2) (B (p - 2) tan(u) - k)` +is `x^m (D sec(u)^p - E sec(u)^(p - 2))`, where `D = B^2 (p - 1)(p - 2)` and `E = B^2 (p - 2)^2 + k^2`. So +`x^m (alpha sec(u)^p + beta sec(u)^(p - 2))` is a multiple of that product wherever `beta D + alpha E = 0`, and a lone +power is where `E` is zero; the cosecant's product is `-x^k csc(u)^(p - 2) (B (p - 2) cot(u) + k)`. Not for `k = 0`, +where the logarithm's derivative is beside it and the substitution answers +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"-(1+b^2*n^2)*sec(a+b*ln(c*x^n))+2*b^2*n^2*sec(a+b*ln(c*x^n))^3".ToEntity().Integrate("x")` | `integral(...)` | 77 characters | +| `"x^m*sec(a+2*ln(c*x^(1/2*sqrt(-(1+m)^2))))^3".ToEntity().Integrate("x")` | `integral(...)` | 189 characters | +| `"sec(a+i*ln(c*x^n)/(n*(-2+p)))^p".ToEntity().Integrate("x")` | `integral(...)` | 138 characters | +| `"csc(a+2*ln(c/x^(1/2*i)))^3".ToEntity().Integrate("x")` | `integral(...)` | 102 characters | + ### A symbolic power of a sine, cosine, secant or cosecant is integrated where it is the derivative of a product **Answers where there were none.** `sec(c + d x)^m (m - (1 + m) sec(c + d x)^2)` was declined, and so were five more diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index ccf2a3a27..ad0b910d6 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -31937,6 +31937,145 @@ private static bool IsAnImaginaryNumber(Entity expr) return answer.InnerSimplified; } + /// + /// A power of the variable times a secant or cosecant power of a logarithm, beside the + /// power two below it, where the pair is the derivative of one product. With + /// u = a + b ln(c x^n), so that u' = B/x for B = b n, and k = m + 1, + /// d/dx [x^k sec(u)^(p - 2) (B (p - 2) tan(u) - k)] = x^m (D sec(u)^p - E sec(u)^(p - 2)) + /// with D = B^2 (p - 1)(p - 2) and E = B^2 (p - 2)^2 + k^2, so + /// x^m (alpha sec(u)^p + beta sec(u)^(p - 2)) integrates to alpha/D times the + /// product wherever beta D + alpha E = 0. The cosecant's product is + /// -x^k csc(u)^(p - 2) (B (p - 2) cot(u) + k), with the same D and E. + /// + /// + /// Rubi's reduction for x^m sec(a + b ln(c x^n))^p lowers p by two and leaves + /// E/D times the lower power, which alone is a hypergeometric function; where + /// E is zero -- b^2 n^2 (p - 2)^2 + (m + 1)^2 = 0 -- or the lower power is the + /// integrand's other term, the product is the whole answer. Declined before: the + /// exponentials are written for a numerical imaginary b only, and + /// sec(a + i ln(c x^n)/(n (p - 2)))^p, x^m sec(a + 2 ln(c x^(sqrt(-(1 + m)^2)/2)))^3 + /// and -(1 + b^2 n^2) sec(u) + 2 b^2 n^2 sec(u)^3 were unanswered or ran out of time. + /// The derivative uses only z^(p - 2) z^2 = z^p, x^m x = x^k and + /// ln(c x^n)' = n/x, which hold for principal powers wherever the base is not zero, + /// so the answer holds wherever the integrand is defined. Not for k = 0, where the + /// logarithm's own derivative is beside it and the substitution answers. Rubi's 4.7.5. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveSecantPowersOfALogarithmTwoApartAsTheDerivativeOfAProduct(Entity expr, Entity.Variable x) + { + static Entity? ArgumentOf(Entity node) => node switch { Secantf(var y) => y, Cosecantf(var y) => y, _ => null }; + if (!expr.Nodes.Any(node => ArgumentOf(node) is { } y && y.Nodes.Any(inner => inner is Logf && inner.ContainsNode(x)))) + return null; + Entity? argument = null; + var isSecant = false; + var monomials = Read(expr); + if (monomials is null || monomials.Count == 0 || argument is null) + return null; + // Gathered by their powers, and those whose coefficients cancel dropped. + var gathered = new List<(Entity Coefficient, Entity OfX, Entity OfF)>(); + foreach (var (coefficient, ofX, ofF) in monomials) + { + var index = gathered.FindIndex(g => VanishesIdentically(g.OfX - ofX) && VanishesIdentically(g.OfF - ofF)); + if (index < 0) + gathered.Add((coefficient, ofX.InnerSimplified, ofF.InnerSimplified)); + else + gathered[index] = (gathered[index].Coefficient + coefficient, gathered[index].OfX, gathered[index].OfF); + } + gathered.RemoveAll(g => VanishesIdentically(g.Coefficient)); + if (gathered.Count is 0 or > 2 || gathered.Count == 2 && !VanishesIdentically(gathered[0].OfX - gathered[1].OfX)) + return null; + var (alpha, m, p) = gathered[0]; + Entity beta = Number.Integer.Zero; + if (gathered.Count == 2) + { + var (otherCoefficient, _, otherPower) = gathered[1]; + if (VanishesIdentically(otherPower - p + 2)) + beta = otherCoefficient; + else if (VanishesIdentically(p - otherPower + 2)) + (alpha, beta, p) = (otherCoefficient, alpha, otherPower); + else + return null; + } + // `x u'`, the constant `B`. + var rate = Functions.PartialFractions.Bare((argument.Differentiate(x) * x).InnerSimplified.Simplify()); + var k = (m + 1).InnerSimplified; + if (rate.ContainsNode(x) || VanishesIdentically(rate) || VanishesIdentically(k)) + return null; + var d = (rate * rate * (p - 1) * (p - 2)).InnerSimplified; + var e = (rate * rate * (p - 2) * (p - 2) + k * k).InnerSimplified; + if (VanishesIdentically(d) || !VanishesIdentically(beta * d + alpha * e)) + return null; + var lowered = (p - 2).InnerSimplified; + // Folded, as `2 b^2 n^2/((b n)^2 2)` and `i/(p - 2) (p - 2)` are not by the simplification + // the answer gets. + var scale = Functions.PartialFractions.Bare((alpha / d).Simplify()); + var slope = Functions.PartialFractions.Bare((rate * lowered).Simplify()); + var product = isSecant + ? MathS.Pow(x, k) * MathS.Pow(MathS.Sec(argument), lowered) * (slope * MathS.Tan(argument) - k) + : -MathS.Pow(x, k) * MathS.Pow(MathS.Cosec(argument), lowered) * (slope * MathS.Cotan(argument) + k); + return (scale * product).InnerSimplified; + + // The integrand as monomials `coefficient x^m f(u)^p`, one `f` of one `u` throughout. + List<(Entity Coefficient, Entity OfX, Entity OfF)>? Read(Entity node) + { + if (!node.ContainsNode(x)) + return new() { (node, Number.Integer.Zero, Number.Integer.Zero) }; + if (node == x) + return new() { (Number.Integer.One, Number.Integer.One, Number.Integer.Zero) }; + switch (node) + { + case Powf(var @base, var exponent) when @base == x && !exponent.ContainsNode(x): + return new() { (Number.Integer.One, exponent, Number.Integer.Zero) }; + case Secantf or Cosecantf: + return Takes(node) ? new() { (Number.Integer.One, Number.Integer.Zero, Number.Integer.One) } : null; + case Powf(var function, var exponent) when function is Secantf or Cosecantf && !exponent.ContainsNode(x): + return Takes(function) ? new() { (Number.Integer.One, Number.Integer.Zero, exponent) } : null; + case Mulf(var left, var right): + return Product(Read(left), Read(right)); + case Divf(var numerator, var denominator): + if (Read(denominator) is not { Count: 1 } below) + return null; + var (c, ofX, ofF) = below[0]; + return Product(Read(numerator), new() { (Number.Integer.One / c, -ofX, -ofF) }); + case Sumf(var augend, var addend): + return Sum(Read(augend), Read(addend), Number.Integer.One); + case Minusf(var minuend, var subtrahend): + return Sum(Read(minuend), Read(subtrahend), Number.Integer.MinusOne); + default: + return null; + } + } + + bool Takes(Entity function) + { + var (inner, secant) = (ArgumentOf(function)!, function is Secantf); + if (argument is null) + (argument, isSecant) = (inner, secant); + return argument == inner && isSecant == secant; + } + + static List<(Entity, Entity, Entity)>? Product(List<(Entity Coefficient, Entity OfX, Entity OfF)>? left, List<(Entity Coefficient, Entity OfX, Entity OfF)>? right) + { + if (left is null || right is null || left.Count * right.Count > 8) + return null; + var product = new List<(Entity, Entity, Entity)>(); + foreach (var l in left) + foreach (var r in right) + product.Add((l.Coefficient * r.Coefficient, l.OfX + r.OfX, l.OfF + r.OfF)); + return product; + } + + static List<(Entity, Entity, Entity)>? Sum(List<(Entity Coefficient, Entity OfX, Entity OfF)>? left, List<(Entity Coefficient, Entity OfX, Entity OfF)>? right, Entity sign) + { + if (left is null || right is null || left.Count + right.Count > 8) + return null; + var sum = new List<(Entity, Entity, Entity)>(left); + foreach (var r in right) + sum.Add((sign * r.Coefficient, r.OfX, r.OfF)); + return sum; + } + } + /// /// A sum as the constant written in every one of its terms times the sum without it: /// a - a sin(x)^2 is a (1 - sin(x)^2), and a + a sin(x) is diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index b1725c89f..9d52758e6 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -920,6 +920,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // The same with an imaginary coefficient on the logarithm, where the closed form divides // by zero: in exponentials, powers of the variable. if ((answer = IndefiniteIntegralSolver.SolveByWritingATrigonometricFunctionOfAnImaginaryLogarithmInExponentials(expr, x, integrateByParts)) is { }) return answer; + // A secant or cosecant power of a logarithm beside the power two below it, where the + // reduction closes: `sec(a + i ln(c x^n)/(n (p - 2)))^p` is the derivative of a product. + if ((answer = IndefiniteIntegralSolver.SolveSecantPowersOfALogarithmTwoApartAsTheDerivativeOfAProduct(expr, x)) is { }) return answer; // A logarithm of a quotient that cancels with the functions in it as // indeterminates, which is how an inverse hyperbolic function of a hyperbolic // one arrives: atanh(tanh(u)) is 1/2 ln(e^(2u)) once its quotient is cancelled. diff --git a/Sources/Tests/UnitTests/Calculus/SecantPowersOfALogarithmTwoApartAsTheDerivativeOfAProductIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/SecantPowersOfALogarithmTwoApartAsTheDerivativeOfAProductIntegralTest.cs new file mode 100644 index 000000000..861a1fe4b --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/SecantPowersOfALogarithmTwoApartAsTheDerivativeOfAProductIntegralTest.cs @@ -0,0 +1,74 @@ +// +// 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 +{ + /// + /// x^m (alpha sec(u)^p + beta sec(u)^(p - 2)) for u = a + b ln(c x^n), and the same + /// with cosecants, where the pair is the derivative of x^(m + 1) sec(u)^(p - 2) (b n (p - 2) tan(u) - (m + 1)): + /// the reduction's remainder vanishing, b^2 n^2 (p - 2)^2 + (m + 1)^2 = 0, or being the + /// integrand's other term. Rubi's 4.7.5. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class SecantPowersOfALogarithmTwoApartAsTheDerivativeOfAProductIntegralTest + { + private static Entity Pinned(Entity e) => e + .Substitute("a", 0.6).Substitute("b", 0.5).Substitute("c", 1.2) + .Substitute("n", 1.0).Substitute("m", 0.6).Substitute("p", 3.5); + + // Where the integral is answered, its derivative is the integrand at every point compared. + private static bool DifferentiatesBack(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + if (text.Contains("integral(")) + return false; + Assert.True(text.Length < 2000, $"{text.Length} characters of answer for {integrand}"); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { 0.4, 0.9, 1.3, 1.7, 2.2, 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-8 * 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}"); + return true; + } + + [Theory] + [InlineData("-(1+b^2*n^2)*sec(a+b*ln(c*x^n))+2*b^2*n^2*sec(a+b*ln(c*x^n))^3")] + [InlineData("-(1+b^2*n^2)*csc(a+b*ln(c*x^n))+2*b^2*n^2*csc(a+b*ln(c*x^n))^3")] + [InlineData("x^m*sec(a+2*ln(c*x^(1/2*sqrt(-(1+m)^2))))^3")] + [InlineData("x^m*csc(a+2*ln(c*x^(1/2*sqrt(-(1+m)^2))))^3")] + [InlineData("sec(a+i*ln(c*x^n)/(n*(-2+p)))^p")] + [InlineData("sec(a-i*ln(c*x^n)/(n*(-2+p)))^p")] + [InlineData("csc(a+i*ln(c*x^n)/(n*(-2+p)))^p")] + [InlineData("csc(a-i*ln(c*x^n)/(n*(-2+p)))^p")] + [InlineData("csc(a+2*ln(c/x^(1/2*i)))^3")] + public void IsAnsweredAndDifferentiatesBack(string integrand) + => Assert.True(DifferentiatesBack(integrand), $"{integrand} was left unevaluated"); + + // The remainder does not vanish: alone, each power is a hypergeometric function of x. + [Theory] + [InlineData("sec(a+b*ln(c*x^n))^p")] + [InlineData("x*sec(a+b*ln(c*x^n))^2")] + public void IsNotClaimedWhereTheReductionDoesNotClose(string integrand) + => DifferentiatesBack(integrand); + } +}