diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 66785bb64..4f6eb215f 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -308,6 +308,26 @@ partial fractions divide. It is taken out in front of the integral first | `"(b*d + 2*c*d*x)^9/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | past a minute | 76,249 characters, the length the content-free one already had | | `"(b*d + 2*c*d*x)^8/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | past a minute | 65,236 characters | +### A whole power is taken out of a power of a product, in front of a constant + +**Answers where there were none.** `(c sin(a + b x)^3)^(1/3)/x` was declined: nothing read a function under a root of +its own power. `(c R f^n)^r` is `K f^(n r) R^r` with `K = (c R f^n)^r/(f^(n r) R^r)`, constant wherever it is defined, +and only a factor of the integrand is rewritten. A constant stays under the root inside `K`, so a negative `c` is +answered too: an odd root is the real root +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(c*sin(a + b*x)^3)^(1/3)".ToEntity().Integrate("x")` | `integral(...)` | 35 characters | +| `"(c*sin(a + b*x)^3)^(2/3)".ToEntity().Integrate("x")` | `integral(...)` | 78 characters | +| `"(c*sin(a + b*x)^3)^(1/3)/x".ToEntity().Integrate("x")` | `integral(...)` | 105 characters | +| `"(c*sin(a + b*x)^3)^(2/3)/x^3".ToEntity().Integrate("x")` | `integral(...)` | 391 characters | +| `"(c*sin(a + b*x^n)^3)^(1/3)/x".ToEntity().Integrate("x")` | `integral(...)` | 129 characters | +| `"sqrt((x^2 - 1)^2/(x*(1 + x^2)))/(1 + x^2)".ToEntity().Integrate("x")` | `integral(...)` | 167 characters | +| `"1/sqrt((a^2 + x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 104 characters | +| `"1/sqrt((a + b*x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 100 characters | +| `"1/sqrt(1 - coth(x)^2)".ToEntity().Integrate("x")` | `integral(...)` | 301 characters | + ### A power of x past the twelfth beside a block is split off **Answers where there were none.** `(c + d x^3 + g x^6 + f x^9)/(x^14 (a + b x^3))` was declined, and so were diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 217d26362..c9d3c9d81 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -19031,6 +19031,95 @@ node is Powf(var @base, Number.Rational r) && byRadicand.TryGetValue(@base, out } } + /// + /// A power of a product with a factor under a whole power, that power taken out: + /// (c R f^n)^r is K f^(n r) R^r, where K = (c R f^n)^r/(f^(n r) R^r) and + /// c is free of x. The logarithmic derivatives of both sides are + /// r (R'/R + n f'/f), so K is + /// constant wherever it is defined, and it stands in front of the answer the way Rubi + /// writes (c sin(a + b x)^3)^(1/3) csc(a + b x) in front of its. + /// + /// + /// (c sin(a + b x)^3)^(1/3)/x, sqrt((x^2 - 1)^2/(x (1 + x^2)))/(1 + x^2) and + /// 1/sqrt((a^2 + x^2)^3) were declined: nothing reads a function under a root of its + /// own power. With the power taken out each is K times an integrand the chain + /// answers. After every other rule, and at the question asked only. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveByTakingTheWholePowersOutOfAPowerOfAProduct(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked || takingTheWholePowersOut) + return null; + var multiple = Variable.CreateUnique(expr, "k_out"); + Entity? constant = null; + var rewritten = AFactor(expr); + if (constant is null) + return null; + var simplified = rewritten.InnerSimplified; + Entity? answer; + takingTheWholePowersOut = true; + try + { + answer = Integration.ComputeAsTheSameQuestion(simplified, x, integrateByParts); + } + finally + { + takingTheWholePowersOut = false; + } + if (answer is null) + return null; + answer = answer.Substitute(multiple, constant); + return answer.ContainsNode(multiple) ? null : answer; + + // Only a factor of the integrand: inside a function, `ln(1/(c x) + sqrt(1/(c^2 x^2) - 1))` + // for `asech(c x)`, K is not a factor of anything, and the question it writes is no shorter. + Entity AFactor(Entity node) => node switch + { + Mulf(var left, var right) => AFactor(left) * AFactor(right), + Divf(var above, var below) => AFactor(above) / AFactor(below), + _ => TakeOut(node), + }; + + Entity TakeOut(Entity node) + { + if (constant is not null || node is not Powf(var @base, Number.Rational exponent) || exponent is Number.Integer + || !@base.ContainsNode(x)) + return node; + var (above, below) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(@base)); + Entity outside = Number.Integer.One; + Entity inside = Number.Integer.One; + var taken = false; + foreach (var (side, sign) in new[] { (above, 1), (below, -1) }) + foreach (var factor in Mulf.LinearChildren(side)) + { + var (factorBase, factorExponent) = factor switch + { + Powf(var b, Number.Integer n) => (b, n), + _ => (factor, Number.Integer.One), + }; + if (!factorBase.ContainsNode(x)) + continue; + if (factorExponent.EInteger.Abs().CompareTo(EInteger.FromInt32(2)) >= 0) + { + outside *= MathS.Pow(factorBase, (sign * factorExponent * exponent).InnerSimplified); + taken = true; + } + else + inside = sign > 0 ? inside * factor : inside / factor; + } + if (!taken) + return node; + // A factor free of x stays under the power, in K: taken out, `c^(1/3)` is only + // `c`'s cube root for `c > 0`. + var written = inside == Number.Integer.One ? outside : outside * MathS.Pow(inside, exponent); + constant = node / written; + return multiple * written; + } + } + + /// Set while asks its question. + [System.ThreadStatic] private static bool takingTheWholePowersOut; + /// /// A rational function of x and one square root of a quadratic in x, /// rationalised by an Euler substitution and handed to the rational integrator. @@ -28451,7 +28540,11 @@ Entity Rewritten(Entity side, Number.Integer sign) 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 }) + // An odd denominator is the real root, `(-8)^(1/3) = -2`, which splits across a + // product for any real constant: `(c sin(x)^3)^(1/3)` is `c^(1/3) sin(x)` for a + // negative `c` as well, and asking `c > 0` left it undefined there. + var aRealRoot = exponent is Number.Rational { ERational: var fraction } && !fraction.Denominator.IsEven; + if (constant != Number.Integer.One && constant.Evaled is not Number.Real { IsPositive: true } && !aRealRoot) { if (constant.Vars.Any() is false) return node; diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 6ad575041..f7b5b5a8b 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -1200,6 +1200,9 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // And a sum with the root below the bar, over its conjugate: after every rule above, // which answer the numeric spellings, `1/(a c + b c x^2 + d sqrt(a + b x^2))` among them. if ((answer = IndefiniteIntegralSolver.SolveByTheConjugateOfASumWithARootOfAQuadraticBelowTheBar(expr, x, integrateByParts)) is { }) return answer; + // And a power of a product with a factor under a whole power, that power taken out in + // front of a constant: `(c sin(a + b x)^3)^(1/3)` is `K sin(a + b x)`. + if ((answer = IndefiniteIntegralSolver.SolveByTakingTheWholePowersOutOfAPowerOfAProduct(expr, x, integrateByParts)) is { }) return answer; // Product-to-sum among the rewrites rather than before them, because a product of // trigonometric functions of *equal* arguments is a power and wants a different tool; // this only fires where the arguments differ, which is exactly what every substitution diff --git a/Sources/Tests/UnitTests/Calculus/AWholePowerTakenOutOfAPowerOfAProductIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AWholePowerTakenOutOfAPowerOfAProductIntegralTest.cs new file mode 100644 index 000000000..64f096a96 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AWholePowerTakenOutOfAPowerOfAProductIntegralTest.cs @@ -0,0 +1,80 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A power of a product with a factor under a whole power, that power taken out in front of + /// a constant: (c sin(a + b x)^3)^(1/3) is K sin(a + b x). + /// #718 + /// + /// + /// Rubi's 4.1.12 and 1.3.2 were declined. Each answer is differentiated back with the constant + /// under the root on both signs: an odd root is the real root, and a negative c used to + /// be answered provided c > 0. + /// + [Trait("Area", "Calculus")] + public sealed class AWholePowerTakenOutOfAPowerOfAProductIntegralTest + { + private static void DifferentiatesBack(string integrand, double[] points, params (string, double)[] pins) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Assert.DoesNotContain("NaN", 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)); + Assert.True(difference / scale < 1e-8, + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, " + + $"where the integrand is {want}"); + } + Assert.True(compared >= 3, + $"only {compared} of {points.Length} points were comparable for {integrand}"); + } + + [Theory] + [InlineData("(c*sin(a+b*x)^3)^(1/3)")] + [InlineData("(c*sin(a+b*x)^3)^(2/3)")] + [InlineData("(c*sin(a+b*x)^3)^(1/3)/x")] + [InlineData("(c*sin(a+b*x)^3)^(2/3)/x^3")] + [InlineData("(c*sin(a+b*x^2)^3)^(1/3)/x^3")] + public void AConstantOfEitherSignUnderAnOddRoot(string integrand) + { + var points = new[] { -1.4, -0.3, 0.8, 2.25 }; + DifferentiatesBack(integrand, points, ("a", 1.0 / 3), ("b", 2), ("c", 2)); + DifferentiatesBack(integrand, points, ("a", 1.0 / 3), ("b", 2), ("c", -3)); + } + + [Theory] + [InlineData("sqrt((-1+x^2)^2/(x*(1+x^2)))/(1+x^2)", new[] { 0.5, 0.8, 2.0, 3.0 })] + [InlineData("1/sqrt((a^2+x^2)^3)", new[] { -2.0, 0.3, 2.5 })] + [InlineData("1/sqrt((a+b*x^2)^3)", new[] { -1.3, 0.3, 1.1 })] + public void ARootOfAPowerOfAPolynomial(string integrand, double[] points) + => DifferentiatesBack(integrand, points, ("a", 2), ("b", -1)); + } +}