From 5b03160537fb517501e614d1790ecc18ed63a395 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 23 Sep 2026 13:43:53 +0000 Subject: [PATCH 1/2] 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, and they cannot see it while one of the two 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 nothing need be assumed about `f`, and a symbolic `c` gives the answer `provided c > 0`. That condition is not decoration: 2.5.0 answered `sec(x)^(3/2)/(b sec(x))^(5/2)` without it, and its 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`, the two agreeing at `x = 0.5`. PR #1388 withdrew the unconditional reading for that reason; this brings the shape back with the condition it owes. Two restrictions keep the rewrite from costing anything elsewhere. The constant comes out only where the rest of the base is a *single* factor, which is where it pays, the power it leaves being free to meet another power of the same function beside it; over a product it would rewrite the question into one no easier and pay a whole descent to find out, which `(cos(x)^11 sin(x)^13)^(-1/4)` measured at ten seconds against forty-seven. And an even whole power of a function is left to the rule that takes the function out of it with its sign, which answers `(a sin(x)^2)^(5/2)` for any real `a` -- read through the nesting and past the sign, since below the tangent substitution the factor arrives as `(u^2)^(-1)`. Rubi's family 4 goes from 290 to 309 of 422, no row lost and two fewer timeouts; the corpus stays at 1707 of 1774 with no wrong answer, no error and no timeout. Part of #718. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 39 ++++++++++ .../Integration/IndefiniteIntegralSolver.cs | 71 ++++++++++++++----- .../ExponentialOfALogarithmIntegralTest.cs | 46 +++++++++++- 3 files changed, 137 insertions(+), 19 deletions(-) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index e5e1db970..16047ee38 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1281,6 +1281,45 @@ 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. + +Two restrictions keep the rewrite from costing anything elsewhere: + +- **A single factor.** The constant comes out only where the rest of the base is one factor, + which is where it pays: the power it leaves can then meet another power of the same function + beside it. Over a product it would rewrite the question into one no easier and pay a whole + descent to find out -- `(cos(x)^11 sin(x)^13)^(-1/4)` measured at ten seconds against + forty-seven. +- **Not an even power.** `(a sin(x)^2)^(5/2)` is `a^(5/2) sgn(sin x) sin(x)^5` for *any* real + `a`, since an even power is never negative, and the rule that takes a function out of its even + power with its sign says so unconditionally. Taking the constant out there would answer the + same shape `provided a > 0` and lose the negative half of the line by getting there first. + +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` | + +`y` is `c + d x`. Rubi's 4.1.0, 4.2.0, 4.3.0 and 4.5.0: family 4 goes from 290 to 309 of 422, +with no row lost and two fewer timeouts. + ### `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..f018ef04f 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -14655,10 +14655,11 @@ 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 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 @@ -14679,10 +14680,16 @@ 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; foreach (var factor in Mulf.LinearChildren(@base)) { if (!factor.ContainsNode(x)) @@ -14690,8 +14697,14 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( constant *= factor; continue; } + 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 +14713,40 @@ 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; + // 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. Taking the constant out here would answer the same + // shape `provided a > 0` and lose the negative half of the line, by getting + // there first. 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; + // Taking the constant out of a base of several factors rewrites the question + // without bringing its answer closer: `(c f g)^b` becomes `c^b (f g)^b` and the + // whole search runs again on a product that is no easier. It pays where the rest + // is a *single* factor, because `(c f)^b` then becomes a power of `f` that can + // meet another power of `f` beside it -- which is what answers + // `sqrt(b sec x)/sec(x)^(7/2)`. Over a product it only costs a descent: + // `(cos(x)^11 sin(x)^13)^(-1/4)` went from 10 to 47 seconds. + if (!split && (constant == Number.Integer.One || restFactors > 1)) return node; if (constant != Number.Integer.One && constant.Evaled is not Number.Real { IsPositive: true }) { @@ -14721,7 +14756,7 @@ 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) diff --git a/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs index 804a4638a..a6521ecf4 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,49 @@ 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). Over a product of several factors it would + /// only rewrite the question into one no easier, at the price of a whole descent. + /// Rubi's 4.1.0, 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)", "")] + 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}"); + } + } } } From 2de06deb5d1d765c01c53d3bf5d5b9fd248c85da Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 23 Sep 2026 20:33:25 +0000 Subject: [PATCH 2/2] Never re-ask the question asked, and take a constant out only after the rules that need none The first head of this change hung the test run: HalfAngleSquareTest's `x/csch(x)^(3/2) + x*sqrt(csch(x))/3` started and never finished, and CI's three C# Test legs sat in the Test step for three and a half hours before they were cancelled. The library writes `csch(x)` as `1/((e^x - e^(-x))/2)`, and the rule re-asks through ComputeAsTheSameQuestion, which keeps the rewritten question at the same depth. `1/s` splits into the factors `1` and `s^(-1)`, and a constant accumulated as `1 * 1` is not `One` as a tree, so the rule believed it had taken something out and asked the identical integrand again at every level -- 38 of the first 40 traced firings rewrote an integrand into itself. With that fixed, the genuine `2` in `csch(x)` was still taken out again under each exponential substitution below it: `x/csch(x)^(3/2)`, declined in 7.4 s on master, did not return in ninety. - The rule declines where its rewrite, once simplified, is the question it was asked, and its products are built without a leading `1 *`. - A constant is taken out only of a single factor in which `x` enters only through trigonometric functions, which is the family it was written for; below a substitution those are gone, which keeps it off the levels where it compounded. A constant written in every term of such a sum counts -- `a - a sin(x)^2` is `a (1 - sin(x)^2)` -- which was the difference between `sqrt(a - a sin(c + d x)^2)` answering and `sqrt(a - a sin(x)^2)` not. - The monomial split keeps its place; the constant pull is asked after the rules that answer its shapes for any real constant, since asked first it answered `sqrt(a + a sin(x))` for a positive `a` only. The shapes that hung now take what they take on master (16.96 s against 17.12 s, 7.39 s against 7.42 s). Family 4 goes from 290 to 309 of 422 and family 6 from 377 to 378, with no row lost in families 1, 4, 5, 6 or 7; the corpus stays at 1707 of 1774 with no wrong answer, error or timeout, and the suite runs to completion, 12,677 tests. Part of #718. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 35 ++-- .../Integration/IndefiniteIntegralSolver.cs | 162 ++++++++++++++++-- .../Integration/Integration.Definition.cs | 8 + .../ExponentialOfALogarithmIntegralTest.cs | 12 +- 4 files changed, 184 insertions(+), 33 deletions(-) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 16047ee38..245032122 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1296,17 +1296,23 @@ secant are both negative** -- at `b = -3` its derivative is `+0.0267i` against t [#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. -Two restrictions keep the rewrite from costing anything elsewhere: - -- **A single factor.** The constant comes out only where the rest of the base is one factor, - which is where it pays: the power it leaves can then meet another power of the same function - beside it. Over a product it would rewrite the question into one no easier and pay a whole - descent to find out -- `(cos(x)^11 sin(x)^13)^(-1/4)` measured at ten seconds against - forty-seven. -- **Not an even power.** `(a sin(x)^2)^(5/2)` is `a^(5/2) sgn(sin x) sin(x)^5` for *any* real - `a`, since an even power is never negative, and the rule that takes a function out of its even - power with its sign says so unconditionally. Taking the constant out there would answer the - same shape `provided a > 0` and lose the negative half of the line by getting there first. +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 @@ -1316,9 +1322,12 @@ the rules for a root of an even power. |---|---|---| | `"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.2.0, 4.3.0 and 4.5.0: family 4 goes from 290 to 309 of 422, -with no row lost and two fewer timeouts. +`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 diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index f018ef04f..d1a5049e0 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -14654,6 +14654,93 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( return answer.InnerSimplified; } + /// + /// 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 @@ -14669,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; @@ -14690,14 +14780,19 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( 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 *= factor; + rest = rest == Number.Integer.One ? factor : rest * factor; restFactors++; single = factor; if (degree is not null) @@ -14717,6 +14812,18 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( } 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 @@ -14727,10 +14834,11 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( // 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. Taking the constant out here would answer the same - // shape `provided a > 0` and lose the negative half of the line, by getting - // there first. 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)`, + // 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; @@ -14739,14 +14847,26 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( if (evenness.IsEven && !evenness.IsZero) return node; var split = degree is not null && degree is not Number.Integer && degree.Evaled is not Number.Integer; - // Taking the constant out of a base of several factors rewrites the question - // without bringing its answer closer: `(c f g)^b` becomes `c^b (f g)^b` and the - // whole search runs again on a product that is no easier. It pays where the rest - // is a *single* factor, because `(c f)^b` then becomes a power of `f` that can - // meet another power of `f` beside it -- which is what answers - // `sqrt(b sec x)/sec(x)^(7/2)`. Over a product it only costs a descent: - // `(cos(x)^11 sin(x)^13)^(-1/4)` went from 10 to 47 seconds. - if (!split && (constant == Number.Integer.One || restFactors > 1)) + // 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 }) { @@ -14759,9 +14879,17 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( 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 a6521ecf4..001b6f5db 100644 --- a/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs @@ -163,15 +163,21 @@ public void APowerOfAMonomialIsDistributed() /// 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). Over a product of several factors it would - /// only rewrite the question into one no easier, at the price of a whole descent. - /// Rubi's 4.1.0, 4.2.0, 4.3.0 and 4.5.0. + /// 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);