diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index c2319c90e..0f1494e17 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1089,6 +1089,26 @@ Rubi's 6.1.5 and 6.6.3 ([#718](https://github.com/asc-community/AngouriMath/issu | `"tanh(x)^4/(a+b*csch(x))".Integrate("x")` | left unevaluated | the antiderivative | | `"csch(x)^2/(a+b*csch(x))".Integrate("x")` | left unevaluated | the antiderivative | +### Hyperbolic powers two apart whose coefficients kill the reduction's residual are integrated + +`cosh(x)^(5/2) - 3 sqrt(cosh(x))/5` was left as written, and neither of its terms has an +elementary antiderivative -- both are elliptic. Their combination has one. From the reduction +`int cosh^p = sinh cosh^(p - 1)/p + (p - 1)/p int cosh^(p - 2)`, a sum of powers two apart is +elementary exactly when the walk from the highest exponent down carries nothing past the lowest, +and the answer is what the walk accumulated; the hyperbolic sine's reduction subtracts where the +cosine's adds. The chain may be several steps long -- Rubi's +`x/sech(x)^(7/2) - 5 x sqrt(sech(x))/21` closes over two, its coefficient being `(5/7)(1/3)` -- +and one linear factor in front is taken by parts against the same antiderivative, `x F - int F` +with `int F` the next power over `p^2`. Rubi's 6.1.1, 6.2.1, 6.5.1 and 6.6.1 +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"cosh(x)^(5/2)-3/5*cosh(x)^(1/2)".Integrate("x")` | left unevaluated | `2 sinh(x) cosh(x)^(3/2)/5`, in exponentials | +| `"x/sech(x)^(3/2)-1/3*x*sqrt(sech(x))".Integrate("x")` | left unevaluated | `2 x sinh(x) sqrt(cosh(x))/3 - 4 cosh(x)^(3/2)/9` up to the form | +| `"x/sech(x)^(7/2)-5/21*x*sqrt(sech(x))".Integrate("x")` | left unevaluated | the antiderivative, over two reduction steps | +| `"x/csch(x)^(3/2)+1/3*x*sqrt(csch(x))".Integrate("x")` | left unevaluated | the antiderivative | + ### `binomial(n, k)` is a function **Addition, not silent.** The binomial coefficient is a node, `Entity.Binomialf`, spelled diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index a803e88b2..0499ba134 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -14051,10 +14051,14 @@ private static (Entity OverTheTwo, Entity Argument)? ReadTheHyperbolicFunctions( // gathering names that y: then the tangent's spelling carries e^y. Entity? argument = null; Entity overTheTwo = expr; - foreach (var halved in new[] { false, true }) - { - var one = halved ? null : (EInteger?)EInteger.One; - var two = halved ? EInteger.One : EInteger.FromInt32(2); + foreach (var scale in new[] { 1, 2, -1 }) + { + // y itself, its half (`tanh(c + d x)` alone is written with e^(2(c + d x)), and + // the gathering names that y), and twice it (`cosh(2x + 1)` with the slope kept + // whole is `cosh(2y)` for `y = x + 1/2`). + var halved = scale == 2; + var one = scale == 1 ? (EInteger?)EInteger.One : scale == 2 ? null : EInteger.FromInt32(2); + var two = scale == 1 ? EInteger.FromInt32(2) : scale == 2 ? EInteger.One : EInteger.FromInt32(4); // With y named from an exponential of the other sign, e^-y is the k = 1 one: // sinh and tanh are odd, cosh even. Entity? ReadTheTwo(Entity node) @@ -14084,7 +14088,7 @@ private static (Entity OverTheTwo, Entity Argument)? ReadTheHyperbolicFunctions( var read = expr.Replace(node => node.ContainsNode(x) && ReadTheTwo(node) is { } inTheTwo ? inTheTwo : node); if (read.ContainsNode(s) || read.ContainsNode(c)) { - argument = halved ? (y / 2).InnerSimplified : y; + argument = scale == 1 ? y : scale == 2 ? (y / 2).InnerSimplified : (2 * y).InnerSimplified; overTheTwo = read; break; } @@ -14165,6 +14169,265 @@ private static (Entity OverTheTwo, Entity Argument)? ReadTheHyperbolicFunctions( return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; } + /// + /// A sum every term of which carries the same power of one linear form in + /// , written as that power times the sum of the rest: + /// x cosh(x)^(3/2) - x sqrt(cosh(x))/3 is x (cosh(x)^(3/2) - sqrt(cosh(x))/3), + /// which parts answers and the split does not -- neither term of that pair has an + /// elementary antiderivative, and only their combination has one. + /// + /// + /// After , which answers every sum whose terms are + /// each integrable and gives the shorter answer for them; this is for the sums it + /// declines. The same question, not one of its own: what is handed on is the integrand + /// rewritten, and the rules below it -- parts among them -- see the product. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveByGatheringACommonPolynomialFactorOfASum(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (TryTakeACommonLinearFactorOfASum(expr, x) is not var (factored, inside)) + return null; + var written = (factored * inside).InnerSimplified; + return written == expr ? null : Integration.ComputeAsTheSameQuestion(written, x, integrateByParts); + } + + /// + /// The sum with the highest power of one linear form in that every + /// term carries taken out of it, as that power and what is left; + /// where the terms carry no such factor in common. Each term is rebuilt from the + /// factors that are left rather than divided, so that taking it out is exact. + /// + private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum(Entity expr, Entity.Variable x) + { + Entity? common = null; + var power = int.MaxValue; + var terms = new List<(Entity Above, Entity Below, int Power)>(); + foreach (var term in Sumf.LinearChildren(expr)) + { + Entity? termCommon = null; + var termPower = 0; + // Above the bar: `x/sech(x)^(3/2)` is one quotient node, and the factor to + // gather is in its numerator. The rest is kept as it stands, so that taking the + // factor out is exact rather than a division nothing cancels. + var (above, below) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(term)); + Entity others = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(above)) + { + var (factorBase, factorPower) = factor is Powf(var inner, Number.Integer whole) && whole.EInteger.Sign > 0 && whole.EInteger.CanFitInInt32() + ? (inner, whole.EInteger.ToInt32Unchecked()) : (factor, 1); + if (factorBase.ContainsNode(x) && TreeAnalyzer.TryGetPolyLinear(factorBase, x, out var slope, out _) && !TreeAnalyzer.IsZero(slope) + && (termCommon is null || termCommon == factorBase)) + { + termCommon = factorBase; + termPower += factorPower; + continue; + } + others *= factor; + } + if (termCommon is null || termPower == 0) + return null; + if (common is null) + common = termCommon; + else if (common != termCommon) + return null; + power = System.Math.Min(power, termPower); + terms.Add((others, below, termPower)); + } + if (common is null || terms.Count < 2 || power == int.MaxValue) + return null; + var factored = power == 1 ? common : MathS.Pow(common, Number.Integer.Create(power)); + Entity inside = Number.Integer.Zero; + foreach (var (above, below, termPower) in terms) + { + var kept = termPower == power ? above + : above * (termPower - power == 1 ? common : MathS.Pow(common, Number.Integer.Create(termPower - power))); + inside += below == Number.Integer.One ? kept : kept / below; + } + inside = inside.InnerSimplified; + return (factored, inside); + } + + /// + /// A pair of hyperbolic powers two apart whose coefficients kill the reduction's + /// residual: cosh(y)^p - (p - 1)/p cosh(y)^(p - 2) is + /// sinh(y) cosh(y)^(p - 1)/p, and sinh(y)^p + (p - 1)/p sinh(y)^(p - 2) is + /// cosh(y) sinh(y)^(p - 1)/p. Neither power has an elementary antiderivative on + /// its own for a fractional p, and the combination does, which is why the terms + /// must not be split apart before this is asked. + /// + /// + /// From the reduction int cosh^p = sinh cosh^(p - 1)/p + (p - 1)/p int cosh^(p - 2), + /// with the coefficient of the lower power chosen so that what is left to integrate is + /// nothing; differentiating the answer through sinh^2 = cosh^2 - 1 is the whole + /// proof. Rubi's 6.5.1 and 6.6.1 -- x/sech(x)^(3/2) - x sqrt(sech(x))/3 is x + /// times such a pair -- meet it beside a polynomial, which the common factor of a sum + /// gathers and parts then separates. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveAPairOfHyperbolicPowersTwoApart(Entity expr, Entity.Variable x, bool integrateByParts) + { + var s = Variable.CreateUnique(expr, "s_hyp"); + var c = Variable.CreateUnique(expr, "c_hyp"); + if (ReadTheHyperbolicFunctions(expr, x, s, c) is not var (readBack, argument)) + { + return null; + } + // One linear factor in front, at most: `(c + d x)(A f^p + B f^(p - 2))` is + // `(c + d x) F - d int F` with `F` the pair's antiderivative, and `int F` is + // `A f^p/(p^2 a^2)` -- both closed, where a second power of `x` would need + // `int f^p` itself, which is elliptic. + Entity? inFront = null; + var overTheTwo = readBack; + // A sum every term of which carries the factor -- `x A f^p + x B f^(p - 2)` -- is + // that factor times the pair, and the gathering has not happened yet: this rule + // runs before the split, which would hand each elliptic term to the whole chain. + if (readBack.ContainsNode(x) && readBack is Sumf or Minusf + && TryTakeACommonLinearFactorOfASum(readBack, x) is var (gathered, insideOfSum) + && !insideOfSum.ContainsNode(x)) + { + inFront = gathered; + overTheTwo = insideOfSum; + } + else if (readBack.ContainsNode(x)) + { + var (above, below) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(readBack)); + Entity others = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(above)) + { + if (factor.ContainsNode(x) && inFront is null + && TreeAnalyzer.TryGetPolyLinear(factor, x, out var factorSlope, out _) && !TreeAnalyzer.IsZero(factorSlope)) + inFront = factor; + else + others *= factor; + } + if (inFront is null) + return null; + overTheTwo = (below == Number.Integer.One ? others : others / below).InnerSimplified; + if (overTheTwo.ContainsNode(x)) + return null; + } + // Each term over its own bar, and each factor read with the exponent it stands at: + // `1/sech(y)^(5/2)` arrives as `(c^(-5/2))^(-1)` or as `1/(1/c)^(5/2)` depending on + // what rewrote it, and both are `c^(5/2)`. A power of a power multiplies. + var terms = new List<(Entity Coefficient, Variable Function, Number.Rational Exponent)>(); + foreach (var term in Sumf.LinearChildren(overTheTwo)) + { + Entity coefficient = Number.Integer.One; + Variable? termFunction = null; + ERational termExponent = ERational.Zero; + var (above, below) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(term)); + foreach (var (side, sign) in new[] { (above, 1), (below, -1) }) + foreach (var factor in Mulf.LinearChildren(side)) + { + if (!factor.ContainsNode(s) && !factor.ContainsNode(c)) + { + coefficient = sign == 1 ? coefficient * factor : coefficient / factor; + continue; + } + // f, f^p, or (f^p)^q, to any depth of powers. + var node = factor; + var exponent = ERational.One; + while (true) + { + if (node is Powf(var inner, Number.Rational power)) + { + exponent = exponent.Multiply(power.ERational); + node = inner; + continue; + } + // A reciprocal inside the power: `sech(y)^(5/2)` is `(1/c)^(5/2)`, + // which is `c^(-5/2)`. + if (node is Divf(Number.Integer { EInteger.IsZero: false } one, var reciprocal) && one.EInteger.Equals(EInteger.One)) + { + exponent = exponent.Negate(); + node = reciprocal; + continue; + } + break; + } + if (node is not Variable read || read != s && read != c) + return null; + if (termFunction is not null && termFunction != read) + return null; + termFunction = read; + termExponent = termExponent.Add(sign == 1 ? exponent : exponent.Negate()); + } + if (termFunction is null) + return null; + terms.Add((coefficient.InnerSimplified, termFunction, Number.Rational.Create(termExponent))); + } + if (terms.Count < 2 || terms.Any(term => term.Function != terms[0].Function)) + return null; + var function = terms[0].Function; + // The exponents are one class two apart, so the reduction walks from the highest + // down; a gap is a term with coefficient zero. + var highest = terms[0].Exponent.ERational; + var lowest = terms[0].Exponent.ERational; + foreach (var term in terms) + { + if (term.Exponent.ERational.CompareTo(highest) > 0) + highest = term.Exponent.ERational; + if (term.Exponent.ERational.CompareTo(lowest) < 0) + lowest = term.Exponent.ERational; + } + static bool IsAWholeNumberOfSteps(ERational difference) + { + // ToLowestTerms: a difference of two halves comes back as `8/2` unreduced. + var reduced = difference.ToLowestTerms(); + return reduced.Denominator.Equals(EInteger.One) && reduced.Numerator.Remainder(EInteger.FromInt32(2)).IsZero; + } + if (terms.Any(term => !IsAWholeNumberOfSteps(highest.Subtract(term.Exponent.ERational)))) + return null; + var steps = highest.Subtract(lowest).ToLowestTerms().Numerator.Divide(EInteger.FromInt32(2)); + if (steps.CompareTo(EInteger.FromInt32(8)) > 0) + return null; + if (terms.All(term => term.Exponent is Number.Integer)) + return null; + if (!TreeAnalyzer.TryGetPolyLinear(argument, x, out var slope, out _) || TreeAnalyzer.IsZero(slope)) + return null; + + // int f^e = g f^(e - 1)/(e a) + s (e - 1)/e int f^(e - 2), with g the other function, + // s = +1 for the hyperbolic cosine (sinh^2 = cosh^2 - 1) and -1 for the sine + // (cosh^2 = sinh^2 + 1), and a the argument's slope. Walking down from the highest + // exponent, each term's coefficient is its own plus what the step above carried; the + // sum is elementary exactly when nothing is carried past the lowest. + var other = function == c ? MathS.Hyperbolic.Sinh(argument) : MathS.Hyperbolic.Cosh(argument); + var self = function == c ? MathS.Hyperbolic.Cosh(argument) : MathS.Hyperbolic.Sinh(argument); + var reductionSign = function == c ? Number.Integer.One : Number.Integer.Create(-1); + Entity antiderivative = Number.Integer.Zero; + Entity integralOfAntiderivative = Number.Integer.Zero; + Entity carried = Number.Integer.Zero; + for (var at = highest; at.CompareTo(lowest) >= 0; at = at.Subtract(ERational.FromInt32(2))) + { + Entity written = Number.Rational.Create(at); + Entity coefficient = carried; + foreach (var term in terms) + if (term.Exponent.ERational.Equals(at)) + coefficient += term.Coefficient; + coefficient = Functions.PartialFractions.Bare(coefficient.InnerSimplified); + if (at.IsZero || coefficient == Number.Integer.Zero) + { + carried = Number.Integer.Zero; + if (!at.IsZero) + continue; + return null; + } + // g f^(e - 1)/(e a), and its own antiderivative f^e/(e^2 a^2), which is what a + // linear factor in front needs. + var below = (written - Number.Integer.One).InnerSimplified; + antiderivative += coefficient * other * MathS.Pow(self, below) / (written * slope); + integralOfAntiderivative += coefficient * MathS.Pow(self, written) / (written * written * slope * slope); + carried = Functions.PartialFractions.Bare((reductionSign * coefficient * (written - Number.Integer.One) / written).InnerSimplified); + } + if (carried != Number.Integer.Zero && Functions.PartialFractions.Bare(carried.Simplify()) != Number.Integer.Zero) + return null; + if (inFront is null) + return antiderivative.InnerSimplified; + if (!TreeAnalyzer.TryGetPolyLinear(inFront, x, out var frontSlope, out _)) + return null; + return (inFront * antiderivative - frontSlope * integralOfAntiderivative).InnerSimplified; + } + /// /// Bioche's first two rules for the hyperbolic functions: a rational function of /// sinh(y) and cosh(y) that is odd in the hyperbolic sine is a rational diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index a9a63adbc..886a371e0 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -714,7 +714,19 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // answered Timofeev's hyperbolic 560 in a minute through Euler's substitution on // that quartic. if ((answer = IndefiniteIntegralSolver.SolveByHyperbolicTangentSubstitution(expr, x, integrateByParts)) is { }) return answer; + // Hyperbolic powers two apart whose coefficients kill the reduction's residual: + // `cosh(y)^p - (p - 1)/p cosh(y)^(p - 2)` is `sinh(y) cosh(y)^(p - 1)/p`, where + // neither power alone is elementary. Before the split, which would hand each + // elliptic term to the whole chain and spend a minute finding that out. + if ((answer = IndefiniteIntegralSolver.SolveAPairOfHyperbolicPowersTwoApart(expr, x, integrateByParts)) is { }) + return answer; if ((answer = IndefiniteIntegralSolver.SolveBySplittingSum(expr, x, integrateByParts)) is { }) return answer; + // A sum whose terms all carry the same power of one linear form, as that power times + // the sum of the rest: `x cosh(x)^(3/2) - x sqrt(cosh(x))/3` is `x` times a pair + // neither half of which is elementary, which the split cannot answer and parts can. + if (expr is Entity.Sumf or Entity.Minusf + && (answer = IndefiniteIntegralSolver.SolveByGatheringACommonPolynomialFactorOfASum(expr, x, integrateByParts)) is { }) + return answer; // The half-angle substitution goes *after* linearity, and that is not a preference. // It fires on anything built from sines and cosines, and it answers `cos(x) + 1` with // a correct expression in tan(x/2) some forty characters long where splitting the sum diff --git a/Sources/Tests/UnitTests/Calculus/HalfAngleSquareTest.cs b/Sources/Tests/UnitTests/Calculus/HalfAngleSquareTest.cs index debcfc4e1..f77dfe4fc 100644 --- a/Sources/Tests/UnitTests/Calculus/HalfAngleSquareTest.cs +++ b/Sources/Tests/UnitTests/Calculus/HalfAngleSquareTest.cs @@ -111,5 +111,49 @@ Entity Pin(Entity e) } } } + /// + /// Hyperbolic powers two apart whose coefficients kill the reduction's residual: + /// int cosh^p = sinh cosh^(p - 1)/p + (p - 1)/p int cosh^(p - 2), so + /// cosh^p - (p - 1)/p cosh^(p - 2) is sinh cosh^(p - 1)/p although neither + /// power alone is elementary; the chain may be several steps long, as Rubi's + /// x/sech(x)^(7/2) - 5 x sqrt(sech(x))/21 is, whose coefficient is + /// (5/7)(1/3). One linear factor in front is taken by parts against the same + /// antiderivative. Rubi's 6.1.1, 6.2.1, 6.5.1 and 6.6.1. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + [Theory] + [InlineData("x/sech(x)^(3/2) - x*sqrt(sech(x))/3")] + [InlineData("x/sech(x)^(5/2) - 3*x/sqrt(sech(x))/5")] + [InlineData("x/sech(x)^(7/2) - 5*x*sqrt(sech(x))/21")] + [InlineData("x/cosh(x)^(7/2) + 3*x*sqrt(cosh(x))/5")] + [InlineData("x/csch(x)^(3/2) + x*sqrt(csch(x))/3")] + [InlineData("cosh(x)^(5/2) - 3*cosh(x)^(1/2)/5")] + [InlineData("sinh(2*x + 1)^(3/2) + sinh(2*x + 1)^(-1/2)/3")] + public void ThePairOfPowersTwoApartIsElementary(string integrand) + { + var integral = integrand.ToEntity().Integrate("x").Substitute("C", 0); + Assert.DoesNotContain("integral(", integral.Stringize()); + Assert.DoesNotContain("NaN", integral.Stringize()); + var derivative = integral.Differentiate("x"); + var original = integrand.ToEntity(); + var compared = 0; + foreach (var at in new[] { 0.3, 0.7, 1.1, 1.9, 2.6 }) + { + 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-7, $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 4, $"only {compared} of five points were comparable for {integrand}"); + } + + /// A combination that does not kill the residual is not this rule's, and is left alone. + [Fact] + public void AnotherCoefficientIsLeftAsWritten() + => Assert.Contains("integral(", "cosh(x)^(5/2) - cosh(x)^(1/2)/2".ToEntity().Integrate("x").Stringize()); } }