A power of the variable times a sine or cosine of a logarithm, and a power of a monomial

BREAKING-CHANGES.md

@@ -1256,6 +1256,31 @@ here that writes one does: left inside, a rule below differentiated it, and
 | `"sqrt(a^2+2*a*b*x+b^2*x^2)*sqrt(c+pe*x+d*x^2)".Integrate("x")` | left unevaluated | the antiderivative |
 | `"(a^2+2*a*b*x+b^2*x^2)^(5/2)".Integrate("x")` | the antiderivative, through the substitution search | `sgn(a + b x) (a + b x)^6 |b|^5/(6 b)` up to the form |
 
+### A power of the variable times a sine or cosine of a logarithm, and a power of a monomial
+
+`x^2 sin(a + b ln(c x^n))` and `(c x^n)^b` were both left as written. Two rules, each exact:
+
+**The trigonometric of a logarithm** is a closed form rather than a search. Two rounds of parts
+close on the integrand -- the sine's remainder is the cosine's integral and the cosine's is the
+sine's -- so the pair is *solved*: `int x^m sin(L) dx` is
+`x^(m + 1)((m + 1) sin(L) - B cos(L))/((m + 1)^2 + B^2)` wherever `L' = B/x` for a constant `B`,
+which `a + b ln(c x^n)` is with `B = b n`; the cosine's is the same with
+`(m + 1) cos(L) + B sin(L)`. The hyperbolic twin needs no rule -- the library writes `sinh` as
+exponentials, which fold against the logarithm -- and the sine and cosine are nodes that fold
+against nothing.
+
+**A power of a monomial** is distributed: `(c x^n)^b` is `c^b x^(n b)`, for a positive `c`, and
+only where `n` is not whole. With a whole `n` the integrand is real at a negative `x` too, and
+there `(c x^n)^b` is a power of `|x|`, not of `x`: `(2 u^3)^(3/2)` is `2^(3/2) sgn(u) u^(9/2)`,
+which the rules for a root of an even power already answer with its sign
+([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"x^2*sin(a+b*ln(c*x^n))".Integrate("x")` | left unevaluated | `x^3(3 sin(L) - b n cos(L))/(9 + b^2 n^2)` |
+| `"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` |
+
 ### `binomial(n, k)` is a function
 
 **Addition, not silent.** The binomial coefficient is a node, `Entity.Binomialf`, spelled

Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs

@@ -14583,6 +14583,154 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum(
             return back.Nodes.Any(node => node == MathS.NaN) ? null : back;
         }
 
+        /// <summary>
+        /// A power of the variable times a sine or cosine of a logarithm, in closed form:
+        /// <c>int x^m sin(L) dx</c> is <c>x^(m + 1)((m + 1) sin(L) - B cos(L))/((m + 1)^2 + B^2)</c>
+        /// whenever <c>L' = B/x</c> for a constant <c>B</c>, which <c>a + b ln(c x^n)</c> is with
+        /// <c>B = b n</c>; the cosine's is the same with <c>(m + 1) cos(L) + B sin(L)</c>.
+        /// </summary>
+        /// <remarks>
+        /// Two rounds of parts close on the integrand -- the sine's remainder is the cosine's
+        /// integral and the cosine's is the sine's -- and the pair is solved rather than
+        /// iterated. Differentiating the answer is the whole proof: the cross terms cancel and
+        /// what is left is <c>x^m sin(L)((m + 1)^2 + B^2)</c>. The hyperbolic twin needs no rule,
+        /// the library writing <c>sinh</c> as exponentials that fold against the logarithm;
+        /// the sine and cosine are nodes and fold against nothing. Rubi's 4.7.5.
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </remarks>
+        internal static Entity? SolveAPowerTimesATrigonometricOfALogarithm(Entity expr, Entity.Variable x, bool integrateByParts)
+        {
+            Entity coefficient = Number.Integer.One;
+            Entity? degree = null;
+            Entity? argument = null;
+            var isSine = false;
+            foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr))
+            {
+                if (!factor.ContainsNode(x))
+                {
+                    coefficient = underneath ? coefficient / factor : coefficient * factor;
+                    continue;
+                }
+                switch (factor)
+                {
+                    case Sinf(var inner) when argument is null && !underneath:
+                        (argument, isSine) = (inner, true);
+                        continue;
+                    case Cosf(var inner) when argument is null && !underneath:
+                        (argument, isSine) = (inner, false);
+                        continue;
+                    case Variable bare when bare == x && degree is null:
+                        degree = underneath ? Number.Integer.Create(-1) : Number.Integer.One;
+                        continue;
+                    case Powf(Variable bare, var power) when bare == x && degree is null && !power.ContainsNode(x):
+                        degree = underneath ? (-power).InnerSimplified : power;
+                        continue;
+                    default:
+                        return null;
+                }
+            }
+            if (argument is null || !argument.ContainsNode(x))
+                return null;
+            // L' x is the constant B, which is what makes the pair close.
+            var b = (argument.Differentiate(x) * x).InnerSimplified;
+            b = Functions.PartialFractions.Bare(b.Simplify());
+            if (b.ContainsNode(x) || b.Evaled is Number.Complex { IsZero: true })
+                return null;
+            var m = degree ?? Number.Integer.Zero;
+            var mPlusOne = (m + Number.Integer.One).InnerSimplified;
+            // `(m + 1)^2 + B^2`, which the answer divides by. Symbolically as well as
+            // numerically: for `B = n sqrt(-9/n^2)` and `m = 2` it is `9 + (n sqrt(-9/n^2))^2`,
+            // which `InnerSimplified` leaves standing and `Simplify` folds to zero -- an
+            // imaginary `B` of the right size makes the pair of parts circular rather than
+            // closing, and the answer divided by nothing at all.
+            var scale = (MathS.Sqr(mPlusOne) + MathS.Sqr(b)).InnerSimplified;
+            if (scale.Evaled is Number.Complex { IsZero: true }
+                || scale.Vars.Any() && Functions.PartialFractions.Bare(scale.Simplify()).Evaled is Number.Complex { IsZero: true })
+                return null;
+            var sine = MathS.Sin(argument);
+            var cosine = MathS.Cos(argument);
+            var bracket = isSine ? mPlusOne * sine - b * cosine : mPlusOne * cosine + b * sine;
+            var answer = coefficient * MathS.Pow(x, mPlusOne) * bracket / scale;
+            return answer.InnerSimplified;
+        }
+
+        /// <summary>
+        /// A power of a monomial with the power distributed: <c>(c x^n)^b</c> is
+        /// <c>c^b x^(n b)</c>, which the power rule answers where the monomial as written is read
+        /// by nothing -- <c>x^2 sin(a + b ln(c x^n))</c> folds to a power of <c>c x^n</c> and
+        /// stops there.
+        /// </summary>
+        /// <remarks>
+        /// Exact on the domain the integrand has. With a symbolic <c>n</c>, <c>x^n</c> is real
+        /// only for a positive <c>x</c> -- at a negative one it is a complex number or nothing --
+        /// so the integrand is defined there and nowhere else, and on <c>x &gt; 0</c> the
+        /// identity holds for a positive <c>c</c>. That condition travels with the answer as
+        /// <c>provided c &gt; 0</c>, unless <c>c</c> is a positive number already; the
+        /// <c>(e x)^(n - 1)</c> of Rubi's 6.5.2 carries the same one. A whole exponent is
+        /// <see cref="SolveByDistributingWholePowersOfProducts"/>'s, and needs no condition.
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </remarks>
+        internal static Entity? SolveByDistributingAPowerOfAMonomial(Entity expr, Entity.Variable x, bool integrateByParts)
+        {
+            Entity assumed = Entity.Boolean.True;
+            var changed = false;
+            var written = expr.Replace(node =>
+            {
+                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.
+                Entity constant = Number.Integer.One;
+                Entity? degree = null;
+                foreach (var factor in Mulf.LinearChildren(@base))
+                {
+                    if (!factor.ContainsNode(x))
+                    {
+                        constant *= factor;
+                        continue;
+                    }
+                    if (degree is not null)
+                        return node;
+                    switch (factor)
+                    {
+                        case Variable bare when bare == x:
+                            degree = Number.Integer.One;
+                            break;
+                        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)
+                    return node;
+                if (constant != Number.Integer.One && constant.Evaled is not Number.Real { IsPositive: true })
+                {
+                    if (constant.Vars.Any() is false)
+                        return node;
+                    var positive = new Greaterf(constant, Number.Integer.Zero);
+                    assumed = assumed == Entity.Boolean.True ? positive : assumed & positive;
+                }
+                changed = true;
+                var distributed = MathS.Pow(x, (degree * exponent).InnerSimplified);
+                return constant == Number.Integer.One ? distributed : MathS.Pow(constant, exponent) * distributed;
+            });
+            if (!changed || written == expr)
+                return null;
+            if (Integration.ComputeAsTheSameQuestion(written.InnerSimplified, x, integrateByParts) is not { } answer)
+                return null;
+            return assumed == Entity.Boolean.True ? answer : answer.Provided(assumed);
+        }
+
         /// <summary>
         /// A pair of hyperbolic powers two apart whose coefficients kill the reduction's
         /// residual: <c>cosh(y)^p - (p - 1)/p cosh(y)^(p - 2)</c> is

Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs

@@ -511,6 +511,12 @@ private static Entity Normalized(Entity expr, Entity.Variable x) =>
             // A whole power of a product of a constant and the variable, as the product of
             // the powers, which is how the inverse hyperbolic secant and cosecant arrive.
             if ((answer = IndefiniteIntegralSolver.SolveByDistributingWholePowersOfProducts(expr, x, integrateByParts)) is { }) return answer;
+            // And a fractional or symbolic power of a monomial: `(c x^n)^b` is `c^b x^(n b)`
+            // for a positive `c`, on the `x > 0` where a symbolic `n` leaves the integrand real.
+            if ((answer = IndefiniteIntegralSolver.SolveByDistributingAPowerOfAMonomial(expr, x, integrateByParts)) is { }) return answer;
+            // A power of the variable times a sine or cosine of a logarithm, where two rounds of
+            // parts close on the integrand: `int x^m sin(a + b ln(c x^n))` in closed form.
+            if ((answer = IndefiniteIntegralSolver.SolveAPowerTimesATrigonometricOfALogarithm(expr, x, integrateByParts)) is { }) return answer;
             // A logarithm of a quotient that cancels with the functions in it as
             // indeterminates, which is how an inverse hyperbolic function of a hyperbolic
             // one arrives: atanh(tanh(u)) is 1/2 ln(e^(2u)) once its quotient is cancelled.

Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs

@@ -90,5 +90,70 @@ private static void DifferentiatesBack(string integrand)
         [InlineData("e^(1 + ln(x + 2)/2)")]
         [InlineData("tanh(ln(x))")]
         public void AHyperbolicFunctionOfALogarithm(string integrand) => DifferentiatesBack(integrand);
+
+        /// <summary>
+        /// A power of the variable times a sine or cosine of a logarithm, in closed form: two
+        /// rounds of parts close on the integrand, since the sine's remainder is the cosine's
+        /// integral and the cosine's is the sine's, and the pair is solved rather than iterated.
+        /// <c>int x^m sin(L)</c> is <c>x^(m+1)((m+1) sin(L) - B cos(L))/((m+1)^2 + B^2)</c>
+        /// wherever <c>L' = B/x</c>. The hyperbolic twin needs no rule, `sinh` being written as
+        /// exponentials that fold against the logarithm; the sine and cosine are nodes.
+        /// Rubi's 4.7.5.
+        /// </summary>
+        [Theory]
+        [InlineData("x^2*sin(a + b*ln(c*x^n))", "a=0.4,b=1.3,c=1.7,n=2.1")]
+        [InlineData("cos(a + b*ln(c*x^n))", "a=0.4,b=1.3,c=1.7,n=2.1")]
+        [InlineData("x^m*cos(a + b*ln(c*x^n))", "a=0.4,b=1.3,c=1.7,n=2.1,m=1.4")]
+        [InlineData("sin(2 + 3*ln(x))", "")]
+        [InlineData("x^2*sin(a + b*ln(x))/x", "a=0.4,b=1.3")]
+        public void APowerTimesATrigonometricOfALogarithm(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}");
+            }
+        }
+
+        /// <summary>
+        /// A fractional or symbolic power of a monomial, distributed: <c>(c x^n)^b</c> is
+        /// <c>c^b x^(n b)</c> on the <c>x &gt; 0</c> where a symbolic <c>n</c> leaves the
+        /// integrand real, for a positive <c>c</c> -- which the answer says.
+        /// </summary>
+        [Fact]
+        public void APowerOfAMonomialIsDistributed()
+        {
+            var integral = "(c*x^n)^b".ToEntity().Integrate("x").Substitute("C", 0);
+            Assert.DoesNotContain("integral(", integral.Stringize());
+            Assert.Contains(integral.Nodes, node => node is Entity.Providedf(_, var predicate) && predicate == "c > 0".ToEntity());
+            var pinned = integral.Substitute("c", 1.7).Substitute("n", 2.1).Substitute("b", 0.6);
+            var derivative = pinned.Differentiate("x");
+            var original = "(c*x^n)^b".ToEntity().Substitute("c", 1.7).Substitute("n", 2.1).Substitute("b", 0.6);
+            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.True(Math.Abs((double)(got - want).RealPart) < 1e-9, $"at x = {at}: {got} for {want}");
+            }
+        }
     }
 }