diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index fcbd8c337..be5134275 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -4370,6 +4370,141 @@ Entity Integral(ERational power) return roots < 2 ? answer : answer.Provided((secantKind ? MathS.Cos(argument) : MathS.Sin(argument)) >= Number.Integer.Zero); } + /// + /// A half-odd power of a ± a sec(y) beside one of c + d sec(y), integrated in + /// the half angle's sine, or its cosine for the minus: a + a sec(y) is + /// 2a cos(y/2)^2/cos(y) and c + d sec(y) is (c cos(y) + d)/cos(y), the two + /// powers of cos(y) make a whole one, and in w = sin(y/2), where + /// cos(y) = 1 - 2w^2 and dy = 2 dw/cos(y/2), what is left is rational beside + /// the one root of c + d - 2c w^2. Either sum may be in the cosine instead, and the + /// cosecant's, with the sine, are the same by the complement. + /// + /// + /// Rubi's sec(e + f x) sqrt(a + a sec(e + f x))/sqrt(c + d sec(e + f x)) and eight more + /// of 4.5.2.1 and 4.5.2.3 were declined: the half-angle tangent reads a root of the one sum + /// and not of the other. The constants are not written: the answer is the integrand times the + /// antiderivative in w over what that differentiates back to, a quotient whose square + /// is one; and it is kept only where that square is one at the sampled points. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveTwoHalfOddPowersOfSecantSumsByTheHalfAngleSine(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAskedOrOneBelow) + return null; + Entity? argument = null; + bool? ofTheSecant = null; + foreach (var node in expr.Nodes) + { + if (TrigonometricArgument(node) is not { } thisArgument || !thisArgument.ContainsNode(x)) + continue; + if (argument is null) + argument = thisArgument; + else if (argument != thisArgument) + return null; + var isSecant = node switch { Secantf or Cosf => true, Cosecantf or Sinf => false, _ => (bool?)null }; + if (isSecant is not { } kind || ofTheSecant is { } seen && seen != kind) + return null; + ofTheSecant = kind; + } + if (argument is null || ofTheSecant is not { } secantKind + || !TreeAnalyzer.TryGetPolyLinear(argument, x, out var rate, out _) || rate.ContainsNode(x) || TreeAnalyzer.IsZero(rate)) + return null; + if (!expr.Nodes.All(node => !node.ContainsNode(x) + || node is Variable or Sumf or Minusf or Mulf or Divf or Sinf or Cosf or Secantf or Cosecantf + || node is Powf(_, Number.Rational))) + return null; + // In y, the argument for the secant and its complement for the cosecant. + var y = Variable.CreateUnique(expr, "y_half_sine"); + var secant = MathS.Sec(y); + var cosine = MathS.Cos(y); + var inY = expr.Replace(node => node switch + { + Secantf(var a) when a == argument => secant, + Cosf(var a) when a == argument => cosine, + Cosecantf(var a) when a == argument => secant, + Sinf(var a) when a == argument => cosine, + _ => node, + }); + if (inY.ContainsNode(x)) + return null; + (Entity Radicand, Number.Rational Exponent)? square = null, other = null; + Entity rest = Number.Integer.One; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(inY)) + { + if (factor is Powf(var radicand, Number.Rational exponent) && exponent is not Number.Integer && radicand.ContainsNode(y)) + { + if (!exponent.ERational.Denominator.Equals(EInteger.FromInt32(2)) || !exponent.ERational.Numerator.CanFitInInt32()) + return null; + var power = underneath ? Number.Rational.Create(exponent.ERational.Negate()) : exponent; + if (square is null && ReadAsOnePlusMinusAFunction(radicand, secant, cosine) is not null) + square = (radicand, power); + else if (other is null) + other = (radicand, power); + else + return null; + continue; + } + rest = underneath ? rest / factor : rest * factor; + } + if (square is not var (squareRadicand, p) || other is not var (otherRadicand, q) + || ReadAsOnePlusMinusAFunction(squareRadicand, secant, cosine) is not var (a, plus, squareInTheSecant) + || rest.Nodes.Any(node => node is Powf(var @base, Number.Rational exponent) && exponent is not Number.Integer && @base.ContainsNode(y))) + return null; + // The other as c + d f, for f the secant or the cosine. + var f = Variable.CreateUnique(expr, "f_half_sine"); + (Entity C, Entity D, bool InTheSecant)? linear = null; + foreach (var (function, inTheSecant) in new[] { (secant, true), (cosine, false) }) + { + var read = otherRadicand.Replace(node => node == function ? f : node); + if (!read.ContainsNode(y) && TreeAnalyzer.TryGetPolyLinear(read, f, out var d, out var c) && !d.ContainsNode(f) && !c.ContainsNode(f)) + { + linear = (c, d, inTheSecant); + break; + } + } + if (linear is not var (cOther, dOther, otherInTheSecant)) + return null; + // The powers of cos(y) the two sums bring, which have to make a whole one. + var cosinePower = ERational.Zero; + if (squareInTheSecant) + cosinePower = cosinePower.Subtract(p.ERational); + if (otherInTheSecant) + cosinePower = cosinePower.Subtract(q.ERational); + if (!cosinePower.IsInteger()) + return null; + // Plus: a (1 + sec(y)) is 2a cos(y/2)^2/cos(y), in w = sin(y/2), cos(y) = 1 - 2w^2 and + // dy = 2 dw/cos(y/2). Minus: -2a sin(y/2)^2/cos(y), in w = cos(y/2), cos(y) = 2w^2 - 1 + // and dy = -2 dw/sin(y/2). Either way the half angle's power over it is a whole power of + // 1 - w^2. + var w = Variable.CreateUnique(expr, "w_half_sine"); + var cosineInW = plus ? 1 - 2 * MathS.Sqr(w) : 2 * MathS.Sqr(w) - 1; + var halfPower = p.ERational.Multiply(ERational.FromInt32(2)).Subtract(ERational.One).Divide(ERational.FromInt32(2)); + if (!halfPower.IsInteger()) + return null; + var restInW = rest.Replace(node => node == secant ? 1 / cosineInW : node == cosine ? cosineInW : node); + var otherInW = otherInTheSecant ? cOther * cosineInW + dOther : cOther + dOther * cosineInW; + var inW = (restInW + * MathS.Pow(1 - MathS.Sqr(w), Number.Integer.Create(halfPower.ToEInteger())) + * MathS.Pow(otherInW, q) + * MathS.Pow(cosineInW, Number.Integer.Create(cosinePower.ToEInteger()))).InnerSimplified; + // Bare: a zeroth power of 1 - w^2 comes back as one provided it is not zero. + inW = Functions.PartialFractions.Bare(inW); + if (inW.ContainsNode(y) || inW.ContainsNode(x) || inW.Nodes.Any(node => node == MathS.NaN)) + return null; + var angle = secantKind ? argument : MathS.pi / 2 - argument; + var wOfX = plus ? MathS.Sin(angle / 2) : MathS.Cos(angle / 2); + var inX = inW.Substitute(w, wOfX) * wOfX.Differentiate(x); + // The integrand is that times a constant and a quotient of roots whose square is one; the + // square of the constant is 4 (±2a)^(2p) over the rate's. + var constantSquared = MathS.Pow(plus ? 2 * a : -2 * a, Number.Integer.Create(p.ERational.Numerator)) * 4 / MathS.Sqr(rate); + if (!Functions.PartialFractions.HoldsAtSampledPoints(constantSquared * MathS.Sqr(inX), MathS.Sqr(expr), x)) + return null; + if (Integration.ComputeAsAQuestionOfItsOwn(inW, w, integrateByParts) is not { } inTermsOfW + || inTermsOfW.Nodes.Any(node => node == MathS.NaN)) + return null; + return expr * inTermsOfW.Substitute(w, wOfX) / inX; + } + /// /// tan(y) tan(2y) written sec(2y) - 1, and the integrand asked again in the /// one argument 2y. diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index f5f23bba3..c408d4161 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -905,6 +905,9 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // And `tan(y) tan(2y)` written `sec(2y) - 1` first, so that the rule reads one argument. if ((answer = IndefiniteIntegralSolver.SolveByWritingATangentTimesThatOfItsDoubleThroughTheSecant(expr, x, integrateByParts)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveByTheHalfAngleTangentBesideAHalfOddPowerOfOnePlusASecant(expr, x, integrateByParts)) is { }) return answer; + // And beside a half-odd power of `c + d sec(y)`, which the half-angle tangent leaves as a + // second root: in sin(y/2), where the two powers of cos(y) they bring make a whole one. + if ((answer = IndefiniteIntegralSolver.SolveTwoHalfOddPowersOfSecantSumsByTheHalfAngleSine(expr, x, integrateByParts)) is { }) return answer; // And `a ± a cosh(y)` under a fractional power: `2a cosh(y/2)^2`, `-2a sinh(y/2)^2`. if ((answer = IndefiniteIntegralSolver.SolveByTheHalfAngleWhereOnePlusAHyperbolicCosineIsASquare(expr, x, integrateByParts)) is { }) return answer; // The hyperbolic sine's, off the real line: 1 + i sinh(y) is (cosh(y/2) + i sinh(y/2))^2. diff --git a/Sources/Tests/UnitTests/Calculus/TwoRootsOfSecantSumsIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/TwoRootsOfSecantSumsIntegralTest.cs new file mode 100644 index 000000000..61d8ada47 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/TwoRootsOfSecantSumsIntegralTest.cs @@ -0,0 +1,52 @@ +// +// 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 half-odd power of a ± a sec(y) beside one of c + d sec(y), in the half + /// angle's sine, for a symbolic slope and shift too. Rubi's 4.5.2.1 and 4.5.2.3. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class TwoRootsOfSecantSumsIntegralTest + { + [Theory] + [InlineData("sec(g + f*x)*sqrt(a + a*sec(g + f*x))/sqrt(c + d*sec(g + f*x))")] + [InlineData("sqrt(c + d*sec(g + f*x))/sqrt(a + a*sec(g + f*x))")] + [InlineData("sqrt(a + a*sec(g + f*x))/(c + d*sec(g + f*x))^(3/2)")] + [InlineData("sqrt(a - a*sec(x))/sqrt(c + d*sec(x))")] + [InlineData("csc(x)*sqrt(a + a*csc(x))/sqrt(c + d*csc(x))")] + public void InTheHalfAngleSine(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 5000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("c", 0.9).Substitute("d", 0.4).Substitute("g", 0.2).Substitute("f", 0.7); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -1.2, -0.7, 0.3, 0.8, 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}"); + } + } +}