diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 1877c36da..a2d0e6fc2 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -405,6 +405,29 @@ one form across the zeros. | `"sqrt(a + a*sin(x))".Integrate("x")` | unevaluated (`a^2 = a^2` is not decided by evaluation) | `-2 sqrt(2a) sgn(sin(u)) cos(u)` | | `"sqrt(1 + sin(x))".Integrate("x")` | `-2 cos(x)/sqrt(1 + sin(x))` | unchanged, the closed rule's | +### Half-odd powers of `a ± a sec` and `a ± a csc` are integrated by the half-angle tangent + +Rubi's `(a + b sec)^m (d sec)^n` files with `a^2 = b^2` hold about a thousand problems with a +half-odd `m`, and none was answered, `sqrt(1 + sec(x))` included: the half angle at which +`a + a cos(y)` is a square leaves a root of the cosine below the bar here. Under +`t = tan(y/2)`, `1 + sec(y)` is `2/(1 - t^2)`, `1 - sec(y)` is `-2 t^2/(1 - t^2)` and `sec(y)` is +`(1 + t^2)/(1 - t^2)`, so beside a power of the secant and anything rational in the sine and +cosine the whole is a rational function of `t` beside one root, of `1 - t^2` or of `1 + t^2`. +Two roots, which a half-odd power of the secant beside a whole one of `1 ± sec` makes, are +elliptic and still declined. The cosecant's the same way through the complement, and a root of +`1 - sec(y)` carries `sgn(tan(y/2))`, constant between its zeros. Each root written apart is exact +where `cos(y)` is positive; beyond it two roots can make the integrand real while each has turned +its sign on its own, so an answer through two says `provided cos(y) >= 0` +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sqrt(1 + sec(x))".ToEntity().Integrate("x")` | `integral(...)` | an arctangent in `tan(x/2)` | +| `"sec(x)/sqrt(a + a*sec(x))".ToEntity().Integrate("x")` | `integral(...)` | `2 arcsin(tan(x/2))/sqrt(2a)` | +| `"1/(sec(x)^(3/2)*sqrt(1 + sec(x)))".ToEntity().Integrate("x")` | `integral(...)` | an antiderivative through a root of `1 + tan(x/2)^2` | +| `"sqrt(a - a*sec(x))".ToEntity().Integrate("x")` | `integral(...)` | a logarithm in `tan(x/2)`, times `sgn(tan(x/2))` | +| `"sqrt(1 + csc(x))".ToEntity().Integrate("x")` | `integral(...)` | an arctangent in `tan(pi/4 - x/2)` | + ### A partial-fraction coefficient with symbols in it is in lowest terms, its rational content included **Improvement, not silent.** The decomposition over written factors with symbols among their diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index cfbedaece..e3abe1f6a 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -2507,6 +2507,215 @@ Entity Step(ERational power) return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; } + /// + /// Half-odd powers of a ± a sec(y), beside powers of d sec(y) or + /// d cos(y) and anything rational in the sine and cosine of y, by the + /// half-angle tangent t = tan(y/2): 1 + sec(y) is 2/(1 - t^2), + /// 1 - sec(y) is -2 t^2/(1 - t^2) and sec(y) is + /// (1 + t^2)/(1 - t^2), so the whole is a rational function of t beside a root of + /// 1 - t^2 or of 1 + t^2 -- or of both, which is elliptic and declined. The + /// cosecant's the same way through the complement, csc(y) = sec(pi/2 - y). + /// + /// + /// + /// Rubi's (a + b sec)^m (d sec)^n files with a^2 = b^2 hold about a thousand + /// problems with a half-odd m, and none was answered, sqrt(1 + sec(x)) + /// included: the half angle at which a + a cos(y) is a square, the rule before this + /// one, leaves a root of the cosine below the bar here, since 1 + sec(y) is that + /// square over cos(y). + /// + /// + /// Exact where the integrand is real. a (1 + sec(y)) is not negative with + /// t inside (-1, 1) for a positive a and outside it for a negative + /// one, and either way (2a/(1 - t^2))^p is (2a)^p (1 - t^2)^(-p) for the + /// principal powers; so for the secant's power beside it. For 1 - sec(y) the square + /// t^2 comes out of the root as |t|, a sign constant between the zeros of + /// tan(y/2) in front. At the question asked or one below it, since it lands on the + /// chain in t. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + /// + internal static Entity? SolveByTheHalfAngleTangentBesideAHalfOddPowerOfOnePlusASecant(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAskedOrOneBelow) + return null; + Entity? argument = 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; + } + if (argument is null || !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 or Tanf or Cotanf + || node is Powf(_, Number.Rational))) + return null; + var secant = MathS.Sec(argument); + var cosecant = new Cosecantf(argument); + // The secant's, or the cosecant's by the complement: csc(y) is sec(y') for + // y' = pi/2 - y, so each function of y is the complementary one of y', and the + // half-angle tangent is t = tan(y'/2). + bool? ofTheSecant = null; + foreach (var node in expr.Nodes) + if (node is Powf(var radicand, Number.Rational exponent) && exponent is not Number.Integer + && ReadAsOnePlusMinusAFunction(radicand, secant, cosecant) is (_, _, var isSecant)) + { + if (ofTheSecant is { } kind && kind != isSecant) + return null; + ofTheSecant = isSecant; + } + if (ofTheSecant is not { } secantKind) + return null; + var primary = secantKind ? secant : cosecant; + var companion = secantKind ? MathS.Cos(argument) : MathS.Sin(argument); + var t = Variable.CreateUnique(expr, "t_half_secant"); + var oneMinus = 1 - MathS.Sqr(t); + var onePlus = 1 + MathS.Sqr(t); + Entity signs = Number.Integer.One; + var found = false; + var declined = false; + var roots = 0; + // Each half-odd power, assembled in t: of a (1 ± sec(y)), or of d sec(y) or d cos(y); + // for the cosecant, of a (1 ± csc(y)), d csc(y) or d sin(y). + var rewritten = expr.Replace(node => + { + if (declined || node is not Powf(var radicand, Number.Rational exponent) || exponent is Number.Integer || !radicand.ContainsNode(x)) + return node; + if (!exponent.ERational.Denominator.Equals(EInteger.FromInt32(2)) || !exponent.ERational.Numerator.CanFitInInt32()) + { + declined = true; + return node; + } + roots++; + if (ReadAsOnePlusMinusAFunction(radicand, secant, cosecant) is (var a, var plus, var isSecant) && isSecant == secantKind) + { + found = true; + if (plus) + return MathS.Pow(2 * a, exponent) * MathS.Pow(oneMinus, (-exponent).InnerSimplified); + // (-2a t^2/(1 - t^2))^p is (-2a)^p |t|^(2p) (1 - t^2)^(-p), and 2p is odd. + signs = signs * MathS.Signum(t); + return MathS.Pow(-2 * a, exponent) * MathS.Pow(t, Number.Integer.Create(exponent.ERational.Numerator)) * MathS.Pow(oneMinus, (-exponent).InnerSimplified); + } + // d sec(y) or d cos(y): a constant times the one function. + Entity coefficient = Number.Integer.One; + Entity? function = null; + foreach (var factor in Mulf.LinearChildren(radicand)) + { + if (!factor.ContainsNode(x)) + coefficient = coefficient * factor; + else if (function is null && (factor == primary || factor == companion)) + function = factor; + else + { + declined = true; + return node; + } + } + if (function is null) + { + declined = true; + return node; + } + var (above, below) = function == primary ? (onePlus, oneMinus) : (oneMinus, onePlus); + return MathS.Pow(coefficient, exponent) * MathS.Pow(above, exponent) * MathS.Pow(below, (-exponent).InnerSimplified); + }); + if (declined || !found) + return null; + // What is left is rational in the functions of y. + var inT = rewritten.Replace(node => node switch + { + Sinf(var a) when a == argument => secantKind ? 2 * t / onePlus : oneMinus / onePlus, + Cosf(var a) when a == argument => secantKind ? oneMinus / onePlus : 2 * t / onePlus, + Tanf(var a) when a == argument => secantKind ? 2 * t / oneMinus : oneMinus / (2 * t), + Cotanf(var a) when a == argument => secantKind ? oneMinus / (2 * t) : 2 * t / oneMinus, + Secantf(var a) when a == argument => secantKind ? onePlus / oneMinus : onePlus / (2 * t), + Cosecantf(var a) when a == argument => secantKind ? onePlus / (2 * t) : onePlus / oneMinus, + _ => node, + }); + // dy = 2 dt/(1 + t^2), and dx = dy/rate. Each power of the two quadratics and of t + // gathered to one, so that they cancel and combine: left as written, `sec(y)` beside + // `dy` reached the rules in t as `(1 + t^2)^0`, which none reads. + var bases = new[] { oneMinus, onePlus, (Entity)t }; + // dy' = -dy for the cosecant. + var (rest, exponents) = GatheredOverTheBases(inT * (secantKind ? 2 : -2) / (rate * onePlus), bases); + Entity integrand = rest.InnerSimplified; + var halfOdd = 0; + for (var i = 0; i < bases.Length; i++) + { + if (exponents[i].IsZero) + continue; + // Not normalised as it is gathered: 8/4 is a whole power. + if (!exponents[i].Numerator.Remainder(exponents[i].Denominator).IsZero) + halfOdd++; + integrand = integrand * MathS.Pow(bases[i], Number.Rational.Create(exponents[i])); + } + // A root of 1 - t^2 beside one of 1 + t^2 is elliptic: declined before it is asked. + if (halfOdd > 1) + return null; + integrand = integrand.InnerSimplified; + if (integrand.ContainsNode(x) || integrand.Nodes.Any(node => node == MathS.NaN)) + return null; + if (Integration.ComputeAsAQuestionOfItsOwn(integrand, t, integrateByParts) is not { } result) + return null; + var back = secantKind ? MathS.Tan(argument / 2) : MathS.Tan(MathS.pi / 4 - argument / 2); + var answer = (signs == Number.Integer.One ? result : signs * result).Substitute(t, back); + if (answer.Nodes.Any(node => node == MathS.NaN)) + return null; + // Each root written apart is exact inside t^2 < 1, where 1 - t^2 is positive, for any + // sign of the constants. Outside it a root whose radicand is negative there is + // imaginary, and one alone leaves the integrand complex, where the answer is nothing + // to be wrong about; but two can make it real again -- `(1 + sec(x))^(5/2) sqrt(cos(x))` + // is real where the cosine is negative -- and there each written apart may have + // turned its sign where the other did not. The answer says where it holds. + return roots < 2 ? answer : answer.Provided((secantKind ? MathS.Cos(argument) : MathS.Sin(argument)) >= Number.Integer.Zero); + } + + /// + /// as what it is a product of beside , and + /// the sum of the exponents of each base in it, read through products, quotients and + /// whole powers of them. + /// + private static (Entity Others, ERational[] Exponents) GatheredOverTheBases(Entity expr, Entity[] bases) + { + var exponents = bases.Select(_ => ERational.Zero).ToArray(); + Entity rest = Number.Integer.One; + Gather(expr, ERational.One); + return (rest, exponents); + + void Gather(Entity node, ERational power) + { + if (System.Array.IndexOf(bases, node) is var at and >= 0) + { + exponents[at] = exponents[at].Add(power); + return; + } + switch (node) + { + case Mulf(var left, var right): + Gather(left, power); + Gather(right, power); + return; + case Divf(var above, var below): + Gather(above, power); + Gather(below, power.Negate()); + return; + case Powf(var @base, Number.Rational exponent) when System.Array.IndexOf(bases, @base) >= 0 + || exponent is Number.Integer && @base is Mulf or Divf or Powf: + Gather(@base, power.Multiply(exponent.ERational)); + return; + default: + rest = power.Equals(ERational.One) ? rest * node : rest * MathS.Pow(node, Number.Rational.Create(power)); + return; + } + } + } + /// /// read as a (1 ± f) for f the given sine or /// cosine: the constant a, whether the sign is plus, and whether f is the diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 87701269e..b663405ba 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -728,6 +728,9 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // cosine, by the half angle at which they are squares: `1 + sin(y)` is `2 sin(u)^2`. Before // the substitution search, which spent twenty seconds on the radicals of the sine. if ((answer = IndefiniteIntegralSolver.SolveByTheHalfAngleWhereOnePlusASineIsASquare(expr, x, integrateByParts)) is { }) return answer; + // And a half-odd power of a +- a sec(y), which is that square over cos(y): by the half-angle + // tangent, in which the whole is rational beside one root. + if ((answer = IndefiniteIntegralSolver.SolveByTheHalfAngleTangentBesideAHalfOddPowerOfOnePlusASecant(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; // A constant out of a fractional power of a trigonometric factor: `sqrt(b sec(x))` is diff --git a/Sources/Tests/UnitTests/Calculus/OnePlusASecantHalfPowerIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/OnePlusASecantHalfPowerIntegralTest.cs new file mode 100644 index 000000000..491476cf4 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/OnePlusASecantHalfPowerIntegralTest.cs @@ -0,0 +1,110 @@ +// +// 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 +{ + /// + /// Half-odd powers of a ± a sec(y) and a ± a csc(y), beside powers of the + /// secant or cosecant and anything rational in the sine and cosine, by the half-angle + /// tangent, in which the whole is rational beside one root: Rubi's + /// (a + b sec)^m (d sec)^n files with a^2 = b^2, of which none was answered. + /// #718 + /// + /// + /// Checked by differentiating back with the symbols pinned, at points where the integrand is + /// real: inside (0, pi/2) for 1 + sec(y) and 1 + csc(y), and in + /// (pi/2, pi) for 1 - sec(y), whose root is real where the secant is negative. + /// + [Trait("Area", "Calculus")] + public sealed class OnePlusASecantHalfPowerIntegralTest + { + private static readonly (string, double)[] Pins = { ("a", 1.7), ("c", 0.9), ("A", 0.6), ("B", 1.3) }; + + private static void DifferentiatesBackPinned(string integrand, double[] points) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + + var derivative = integral.Substitute("C", 0).Differentiate("x"); + Entity original = integrand.ToEntity(); + foreach (var (name, value) in Pins) + { + derivative = derivative.Substitute(name, value); + original = original.Substitute(name, value); + } + var compared = 0; + foreach (var at in points) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + if (got.IsNaN || want.IsNaN) + continue; + compared++; + var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart); + var scale = Math.Max(1.0, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)); + Assert.True(difference / scale < 1e-9, + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, " + + $"where the integrand is {want}"); + } + Assert.True(compared >= 3, $"only {compared} points could be compared for {integrand}"); + } + + /// + /// 1 + sec(y) is 2/(1 - t^2) for t = tan(y/2): its half-odd powers + /// beside a secant, a cotangent, or a half-odd power of the secant, whose root is then of + /// 1 + t^2 alone. + /// + [Theory] + [InlineData("sqrt(1 + sec(x))")] + [InlineData("sqrt(a + a*sec(x))")] + [InlineData("sec(x)/sqrt(a + a*sec(x))")] + [InlineData("sec(x)/(a + a*sec(x))^(3/2)")] + [InlineData("(a + a*sec(x))^(3/2)")] + [InlineData("cot(x)^2/sqrt(a + a*sec(x))")] + [InlineData("tan(x)^2*sqrt(1 + sec(x))")] + [InlineData("1/(sec(x)^(3/2)*sqrt(1 + sec(x)))")] + [InlineData("sqrt(sec(x))/sqrt(1 + sec(x))")] + [InlineData("sec(x)^(3/2)*(A + B*sec(x))/(a + a*sec(x))^(5/2)")] + public void OnePlusASecant(string integrand) + => DifferentiatesBackPinned(integrand, new[] { 0.3, 0.7, 1.1, 1.4 }); + + /// + /// 1 - sec(y) is -2 t^2/(1 - t^2), and the square comes out of the root as + /// a sign of tan(y/2), constant between its zeros. + /// + [Theory] + [InlineData("sqrt(a - a*sec(x))")] + public void OneMinusASecant(string integrand) + => DifferentiatesBackPinned(integrand, new[] { 1.8, 2.2, 2.6, 3.0 }); + + /// + /// Two roots written apart are exact where cos(y) is positive, for any sign of the + /// constants, and the answer says so: beyond it each may turn its sign where the other + /// does not, while the integrand is real. Checked inside, where both roots of + /// 1 - sec(x) here are imaginary and their product real. + /// + [Fact] + public void TwoRootsSayWhereTheyHold() + { + const string integrand = "(c - c*sec(x))^(7/2)/sqrt(a - a*sec(x))"; + Assert.Contains("provided cos(x) >= 0", integrand.ToEntity().Integrate("x").Stringize()); + DifferentiatesBackPinned(integrand, new[] { 0.3, 0.7, 1.1, 1.4 }); + } + + /// The cosecant's, through the complement: csc(y) = sec(pi/2 - y). + [Theory] + [InlineData("sqrt(1 + csc(x))")] + [InlineData("csc(x)/sqrt(a + a*csc(x))")] + [InlineData("cot(x)^2/sqrt(a + a*csc(x))")] + public void OnePlusACosecant(string integrand) + => DifferentiatesBackPinned(integrand, new[] { 0.3, 0.7, 1.1, 1.4 }); + } +}