diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index e5e1db970..245032122 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1281,6 +1281,54 @@ which the rules for a root of an even power already answer with its sign | `"cos(a+b*ln(c*x^n))".Integrate("x")` | left unevaluated | `x(cos(L) + b n sin(L))/(1 + b^2 n^2)` | | `"(c*x^n)^b".Integrate("x")` | left unevaluated | `c^b x^(n b + 1)/(n b + 1) provided c > 0` | +### A constant comes out of a fractional power beside another power of the same function + +`sqrt(b sec(x))/sec(x)^(7/2)` was left as written. The two powers of the secant are the same +function to a total power of `-3`, which the rules for a power of a secant answer -- but they +cannot see it while one of them is written over `b sec(x)` rather than over `sec(x)`. + +`(c f)^b` is `c^b f^b` for a **positive** real `c` and any `f`: both sides pick up the same phase +where `f` is negative, so the rewrite needs nothing assumed about `f`, and where `c` carries +symbols the answer says `provided c > 0`. That condition is not decoration. 2.5.0 answered +`sec(x)^(3/2)/(b sec(x))^(5/2)` without it, and that answer is **wrong wherever `b` and the +secant are both negative** -- at `b = -3` its derivative is `+0.0267i` against the integrand's +`-0.0267i` at `x = 2`, while the two agree at `x = 0.5`. +[#1388](https://github.com/asc-community/AngouriMath/pull/1388) withdrew the unconditional +reading for that reason, and this brings the shape back with the condition it owes. + +Where it is asked, and what it takes, is what keeps it from costing anything elsewhere: + +- **After the rules that answer the same shapes for any real constant.** A function comes out + of its even power with its sign -- `(a sin(x)^2)^(5/2)` is `a^(5/2) sgn(sin x) sin(x)^5` for + every real `a`, an even power being never negative -- and `sqrt(a + a sin(x))` is integrated + by the half angle at which `1 + sin(x)` is a square, again for every real `a`. Asked before + them, this answered both `provided a > 0` and lost the negative half of the line. It is asked + after them, and before the substitution search, which spends the budget on the roots. +- **A single factor in which `x` enters only through trigonometric functions.** The constant + comes out only where the rest of the base is one factor like `sec(x)` or `1 - sin(x)^2`, with + the constant a sum writes in every term read as well: `a - a sin(x)^2` is `a (1 - sin(x)^2)`. + The rewritten integrand is asked again at the same depth, so a rewrite that does not bring the + answer closer multiplies the search at every level below it. Over a product it would rewrite + the question into one no easier -- `(cos(x)^11 sin(x)^13)^(-1/4)` measured at ten seconds + against forty-seven -- and over a hyperbolic function, which is written in exponentials, it + finds a `2` in `csch(x) = 1/((e^x - e^(-x))/2)` under every substitution below it: + `x/csch(x)^(3/2)`, declined in seven seconds, was not declined in ninety. + +The power of the variable is unchanged and still splits only where its exponent is not whole +(the entry above), because `(2 u^3)^(3/2)` is `2^(3/2) sgn(u) u^(9/2)` and the sign belongs to +the rules for a root of an even power. + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sqrt(b*sec(c+d*x))/sec(c+d*x)^(7/2)".Integrate("x")` | left unevaluated | `sqrt(b)(sin(y) - sin(y)^3/3)/d provided b > 0` | +| `"sec(c+d*x)^(3/2)/(b*sec(c+d*x))^(5/2)".Integrate("x")` | `sin(y)/(b^(5/2) d)`, wrong for a negative `b` where the secant is negative | `sin(y)/(b^(5/2) d) provided b > 0` | +| `"sqrt(a-a*sin(x)^2)*tan(x)^6".Integrate("x")` | left unevaluated | `sqrt(a)` times the antiderivative of `sqrt(1 - sin(x)^2) tan(x)^6`, `provided a > 0` | +| `"(a+a*sin(x))^(3/2)/(c+d*sin(x))^(5/2)".Integrate("x")` | left unevaluated | `a^(3/2)` times a closed form, `provided a > 0` | + +`y` is `c + d x`. Rubi's 4.1.0, 4.1.2, 4.1.7, 4.2.0, 4.3.0 and 4.5.0: family 4 goes from 290 +to 309 of 422 with five fewer timeouts, and family 6 from 377 to 378; no row is lost in families +1, 4, 5, 6 or 7, and the 1774-problem suite stays at 1707 with no wrong answer. + ### `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 8f08b74cd..d1a5049e0 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -14655,10 +14655,98 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( } /// - /// A power of a monomial with the power distributed: (c x^n)^b is - /// c^b x^(n b), which the power rule answers where the monomial as written is read - /// by nothing -- x^2 sin(a + b ln(c x^n)) folds to a power of c x^n and - /// stops there. + /// 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 + /// a (1 + sin(x)). A factor counts only where every term writes it, compared as + /// written, so 2 + 4 sin(x) keeps its 2. Null where no factor is common. + /// + private static (Entity Content, Entity Remainder)? TakeTheConstantWrittenInEveryTermOut(Entity sum, Entity.Variable x) + { + var terms = Sumf.LinearChildren(sum); + if (terms.Count < 2) + return null; + List? common = null; + foreach (var term in terms) + { + var constants = new List(); + foreach (var factor in Mulf.LinearChildren(term)) + if (!factor.ContainsNode(x) && factor != Number.Integer.One) + constants.Add(factor); + if (common is null) + common = constants; + else + { + var kept = new List(); + foreach (var candidate in common) + if (constants.Remove(candidate)) + kept.Add(candidate); + common = kept; + } + if (common.Count == 0) + return null; + } + Entity content = Number.Integer.One; + foreach (var factor in common!) + content = content == Number.Integer.One ? factor : content * factor; + Entity rest = Number.Integer.Zero; + foreach (var term in terms) + { + var factors = new List(Mulf.LinearChildren(term)); + foreach (var factor in common) + factors.Remove(factor); + Entity reduced = Number.Integer.One; + foreach (var factor in factors) + reduced = reduced == Number.Integer.One ? factor : reduced * factor; + rest = rest == Number.Integer.Zero ? reduced : rest + reduced; + } + return (content, rest); + } + + /// + /// Whether occurs in at all, and then only + /// inside a sine, cosine, tangent, cotangent, secant or cosecant: 1 + sin(x) and + /// sec(2x)^3 are such, e^x, x sin(x) and a hyperbolic function -- + /// written as exponentials -- are not. + /// + private static bool EntersOnlyThroughTrigonometricFunctions(Entity expr, Entity.Variable x) + { + return expr.ContainsNode(x) && Through(expr, x); + + static bool Through(Entity node, Entity.Variable x) + { + if (node is Sinf or Cosf or Tanf or Cotanf or Secantf or Cosecantf || !node.ContainsNode(x)) + return true; + if (node is Variable) + return false; + foreach (var child in node.DirectChildren) + if (!Through(child, x)) + return false; + return true; + } + } + + /// + /// A fractional or symbolic power of a monomial, distributed: (c x^n)^b is + /// c^b x^(n b) for a positive c and a n that is not whole. + /// + internal static Entity? SolveByDistributingAPowerOfAMonomial(Entity expr, Entity.Variable x, bool integrateByParts) + => DistributeAFractionalPower(expr, x, integrateByParts, constantAlone: false); + + /// + /// A constant taken out of a fractional power of a factor in which x enters only + /// through trigonometric functions: sqrt(b sec(x)) is sqrt(b) sqrt(sec(x)) + /// and sqrt(a - a sin(x)^2) is sqrt(a) sqrt(1 - sin(x)^2), for a positive + /// constant, which the answer says where the constant is a symbol. + /// + internal static Entity? SolveByTakingAConstantOutOfAFractionalPower(Entity expr, Entity.Variable x, bool integrateByParts) + => DistributeAFractionalPower(expr, x, integrateByParts, constantAlone: true); + + /// + /// A constant taken out of a fractional power, and a power of the variable split off + /// with it: (c g)^b is c^b g^b for a positive c, and (c x^n)^b + /// is c^b x^(n b) where n is not whole. sqrt(b sec(y))/sec(y)^(7/2) + /// is sqrt(b) cos(y)^3 and was read by nothing; x^2 sin(a + b ln(c x^n)) + /// folds to a power of c x^n and stopped there. /// /// /// Exact on the domain the integrand has. With a symbolic n, x^n is real @@ -14668,9 +14756,12 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( /// provided c > 0, unless c is a positive number already; the /// (e x)^(n - 1) of Rubi's 6.5.2 carries the same one. A whole exponent is /// 's, and needs no condition. + /// The two are asked from two places: the monomial where it always was, and the constant + /// alone after the rules that answer its shapes for any real constant, since this one + /// owes provided c > 0 and, asked first, answered them weaker. /// https://github.com/asc-community/AngouriMath/issues/718 /// - internal static Entity? SolveByDistributingAPowerOfAMonomial(Entity expr, Entity.Variable x, bool integrateByParts) + private static Entity? DistributeAFractionalPower(Entity expr, Entity.Variable x, bool integrateByParts, bool constantAlone) { Entity assumed = Entity.Boolean.True; var changed = false; @@ -14679,19 +14770,36 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( if (node is not Powf(var @base, var exponent) || exponent.ContainsNode(x) || exponent is Number.Integer || !@base.ContainsNode(x)) return node; - // The base as `c x^n`, with `c` and `n` free of the variable and `n` not whole -- - // a whole one is the distributing rule's, which owes no condition. + // The base as a constant times the rest: `(c g)^b` is `c^b g^b` for a positive + // real `c` and any `g`, on the principal branch and with nothing assumed about + // `g` -- both sides pick up the same phase where `g` is negative. The rest is + // read as `x^n` where it is one, since `(x^n)^b` is `x^(n b)` under the + // condition below and that is what makes the power rule answer it. Entity constant = Number.Integer.One; + Entity rest = Number.Integer.One; Entity? degree = null; + Entity? single = null; + var restFactors = 0; + // Each product starts at its first factor rather than at `1 * factor`: a + // quotient `1/s` arrives as the factors `1` and `s^(-1)`, and a constant built as + // `1 * 1` is not `1` as a tree, so the rule would take out nothing, believe it had + // taken something, and ask the same question again. foreach (var factor in Mulf.LinearChildren(@base)) { if (!factor.ContainsNode(x)) { - constant *= factor; + if (factor != Number.Integer.One) + constant = constant == Number.Integer.One ? factor : constant * factor; continue; } + rest = rest == Number.Integer.One ? factor : rest * factor; + restFactors++; + single = factor; if (degree is not null) - return node; + { + degree = null; + continue; + } switch (factor) { case Variable bare when bare == x: @@ -14700,18 +14808,65 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( case Powf(Variable bare, var power) when bare == x && !power.ContainsNode(x): degree = power; break; - default: - return node; } } - // A whole `n` only where the identity needs nothing: for a negative `x`, - // `x^n` is a real number when `n` is whole, and `(c x^n)^b` is then `|x|`'s - // power, not `x`'s -- `(2 u^3)^(3/2)` is `2^(3/2) sgn(u) u^(9/2)`, and - // distributing it without the sign is a wrong answer where `u < 0`, which the - // rules for a root of an even power answer correctly and this must leave to - // them. With a fractional or symbolic `n` the integrand is real only for a - // positive `x`, and there the distribution is exact. - if (degree is null || degree is Number.Integer || degree.Evaled is Number.Integer) + if (rest == Number.Integer.One) + return node; + // A sum is read for the constant written in every one of its terms: `a - a sin(x)^2` + // is `a (1 - sin(x)^2)`. The substitution for a linear argument happens to write it + // so, and without this `sqrt(a - a sin(c + d x)^2)` came out where + // `sqrt(a - a sin(x)^2)` did not. + if (constantAlone && restFactors == 1 && single is Sumf or Minusf + && TakeTheConstantWrittenInEveryTermOut(single, x) is var (content, withoutIt)) + { + constant = constant == Number.Integer.One ? content : constant * content; + single = withoutIt; + rest = withoutIt; + degree = null; + } + // The power of `x` is split only where `n` is not whole: for a negative `x`, + // `x^n` is a real number when `n` is whole, and `(x^n)^b` is then `|x|`'s power, + // not `x`'s -- `(2 u^3)^(3/2)` is `2^(3/2) sgn(u) u^(9/2)`, and splitting it + // without the sign is a wrong answer where `u < 0`, which the rules for a root + // of an even power answer correctly and this must leave to them. With a + // fractional or symbolic `n` the integrand is real only for a positive `x`, and + // there the split is exact. The constant comes out either way. + // An even whole power of a function is not this rule's: `(a sin(x)^2)^(5/2)` is + // `a^(5/2) sgn(sin x) sin(x)^5` for *any* real `a`, because an even power is + // never negative, and the rule that takes the function out of it with its sign + // says so unconditionally -- at the top, where it is asked first; below the top + // it is not asked, and taking the constant out here would answer the same shape + // `provided a > 0` and carry that up. Read through the nesting and past the + // sign: `1/sqrt(a cot(x)^2)` arrives at this rule, below the tangent + // substitution, as `(a/u^2)^(-1/2)`, + // whose factor is `(u^2)^(-1)` -- an outer exponent of `-1` over an even one, + // and no less non-negative for either. + var evenness = EInteger.One; + for (var peeled = single; peeled is Powf(var peeledBase, Number.Integer step); peeled = peeledBase) + evenness = evenness.Multiply(step.EInteger); + if (evenness.IsEven && !evenness.IsZero) + return node; + var split = degree is not null && degree is not Number.Integer && degree.Evaled is not Number.Integer; + // Each entry point does one of the two: the power of `x` is split where the + // monomial rule is asked, and a constant alone taken out where this one is -- after + // the rules that answer its shapes for any real constant. + if (split == constantAlone) + return node; + // Taking the constant out rewrites the question, and the rewritten question is + // asked at the same depth, so the rewrite has to be one that brings the answer + // closer or it multiplies the search at every level. It does where the rest is a + // single factor in which `x` enters only through trigonometric functions: `(c f)^b` + // then becomes a power of `f` that can meet another power of `f` beside it -- + // `sqrt(b sec x)/sec(x)^(7/2)` -- or a root the half-angle rules read, as + // `sqrt(a + a sin x)` is `sqrt(a) sqrt(1 + sin x)`. Over a product it is no easier: + // `(cos(x)^11 sin(x)^13)^(-1/4)` went from 10 to 47 seconds. And the hyperbolic + // functions are written as exponentials, where the same rewrite finds a `2` in + // `csch(x) = 1/((e^x - e^(-x))/2)` under every substitution below it: + // `x/csch(x)^(3/2)`, declined in seven seconds, was not declined in ninety. Below a + // substitution the trigonometric functions are gone, so this also keeps the + // rewrite to the levels where it was measured. + if (constantAlone && (constant == Number.Integer.One || restFactors > 1 + || !EntersOnlyThroughTrigonometricFunctions(single!, x))) return node; if (constant != Number.Integer.One && constant.Evaled is not Number.Real { IsPositive: true }) { @@ -14721,12 +14876,20 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( assumed = assumed == Entity.Boolean.True ? positive : assumed & positive; } changed = true; - var distributed = MathS.Pow(x, (degree * exponent).InnerSimplified); + var distributed = split ? MathS.Pow(x, (degree! * exponent).InnerSimplified) : MathS.Pow(rest, exponent); return constant == Number.Integer.One ? distributed : MathS.Pow(constant, exponent) * distributed; }); - if (!changed || written == expr) + if (!changed) + return null; + // A rewrite of the same question has to change the question. The question is asked + // again without consuming depth, so an integrand rewritten into itself -- which is + // what `(1/s)^b` was, read as `1^b (1/s)^b` -- is asked again at every level of the + // search below it, without bound: `x/csch(x)^(3/2)` declined in seven seconds before + // this rule and did not come back at all with it. + var rewritten = written.InnerSimplified; + if (rewritten == expr || rewritten == expr.InnerSimplified) return null; - if (Integration.ComputeAsTheSameQuestion(written.InnerSimplified, x, integrateByParts) is not { } answer) + if (Integration.ComputeAsTheSameQuestion(rewritten, x, integrateByParts) is not { } answer) return null; return assumed == Entity.Boolean.True ? answer : answer.Provided(assumed); } diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 6df2872e3..831916973 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -629,6 +629,14 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => if ((answer = IndefiniteIntegralSolver.SolveByTheHalfAngleWhereOnePlusASineIsASquare(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 + // `sqrt(b) sqrt(sec(x))` for a positive `b`, which meets the other powers of the + // secant beside it. After the rules that answer the same shapes for any real + // constant -- a function out of its even power with its sign, the half angle of + // `a +- a sin(y)` -- since this one owes `provided b > 0` and, asked first, answered + // `sqrt(a + a sin(x))` for a positive `a` only. Before the substitution search, which + // spends the budget on the roots. + if ((answer = IndefiniteIntegralSolver.SolveByTakingAConstantOutOfAFractionalPower(expr, x, integrateByParts)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveBySubstitution(expr, x, integrateByParts)) is { }) return answer; // A logarithmic derivative the substitution above could not read for a symbol in // an exponent: `(x^(n-1) - 1)/(x^n - n x)`. diff --git a/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs index 804a4638a..001b6f5db 100644 --- a/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs @@ -1,4 +1,4 @@ -// +// // Copyright (c) 2019-2026 Angouri. // AngouriMath is licensed under MIT. // Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. @@ -155,5 +155,55 @@ public void APowerOfAMonomialIsDistributed() Assert.True(Math.Abs((double)(got - want).RealPart) < 1e-9, $"at x = {at}: {got} for {want}"); } } + + /// + /// (c f)^b is c^b f^b for a positive real c and any f: both + /// sides take the same branch where f is negative, so the rewrite is exact and + /// not conditional on the sign of f. It is worth doing where the rest of the base + /// is a single factor, because the power it leaves can then meet another power of the + /// same function beside it -- sqrt(b sec x)/sec(x)^(7/2) is the secant to the + /// power -3, which the rules for those answer and cannot see while one of the two + /// powers is written over b sec(x); and sqrt(a + a sin x) is + /// sqrt(a) sqrt(1 + sin x), a root the half-angle rules read. Over a product of + /// several factors it would only rewrite the question into one no easier, at the price + /// of a whole descent, and over anything in which x enters other than through a + /// trigonometric function -- a hyperbolic function, written in exponentials -- it + /// multiplied the search at every level below it. + /// Rubi's 4.1.0, 4.1.2, 4.1.7, 4.2.0, 4.3.0 and 4.5.0. + /// + [Theory] + [InlineData("sqrt(b*sec(x))/sec(x)^(7/2)", "b=1.7")] + [InlineData("sec(x)^(3/2)/(b*sec(x))^(5/2)", "b=1.7")] + [InlineData("(b*sec(x))^(3/2)*sec(x)^(1/2)", "b=1.7")] + [InlineData("(3*sin(x))^(5/2)/sin(x)^(3/2)", "")] + [InlineData("(a+a*sin(x))^(3/2)/(c+d*sin(x))^(5/2)", "a=1.7,c=2.3,d=0.4")] + [InlineData("sqrt(a-a*sin(x)^2)*tan(x)^6", "a=1.7")] + public void AConstantComesOutOfAPowerOfASingleFactor(string integrand, string pins) + { + var integral = integrand.ToEntity().Integrate("x").Substitute("C", 0); + Assert.DoesNotContain("integral(", integral.Stringize()); + Assert.DoesNotContain("NaN", integral.Stringize()); + Entity Pin(Entity e) + { + foreach (var pin in pins.Split(',', StringSplitOptions.RemoveEmptyEntries)) + { + var parts = pin.Split('='); + e = e.Substitute(parts[0], double.Parse(parts[1], System.Globalization.CultureInfo.InvariantCulture)); + } + return e; + } + var derivative = Pin(integral).Differentiate("x"); + var original = Pin(integrand.ToEntity()); + 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(); + Assert.False(got.IsNaN, $"the antiderivative of {integrand} differentiates to NaN at x = {at}"); + 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}"); + } + } } }