diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 6aa4bdb85..9b9be3dcc 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1679,6 +1679,21 @@ where it was `0.1 * x ^ 2 / 2 + C` | `limit(sin(c*x)/x, x, 0)` | `NaN` | `c` | | `limit((1 - cos(c*x))/x^2, x, 0)` | `NaN` | `c ^ 2 / 2` | +### A half-odd power of a logarithm beside a power of `x` is integrated + +`x sqrt(ln(x))` was left unintegrated. `x^p F^n`, where `F` is `ln x` or `A + B ln(c x^r)` and +`n` is half an odd number, is now integrated onto the Gaussian's moments under `t = sqrt(F)`. +A polynomial beside it is summed a power at a time, and the logarithm of a power of a linear is +the same question under `u = d + e x` +([#1501](https://github.com/asc-community/AngouriMath/issues/1501)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x*sqrt(ln(x))".Integrate("x")` | `integral(x * sqrt(ln(x)), x)` | the antiderivative, with `erfi` | +| `"sqrt(ln(a*x^n))".Integrate("x")` | `integral(sqrt(ln(a * x ^ n)), x)` | the antiderivative, with `erfi` | +| `"x^2/ln(x)^(3/2)".Integrate("x")` | `integral(x ^ 2 / ln(x) ^ (3/2), x)` | the antiderivative, with `erfi` | +| `"(f + g*x)*(a + b*ln(c*(d + k*x)^n))^(3/2)".Integrate("x")` | `integral((f + g * x) * (a + b * ln(c * (d + k * x) ^ n)) ^ (3/2), x)` | the antiderivative | + ### `binomial(n, k)` is a function **Addition, not silent.** The binomial coefficient is a node, `Entity.Binomialf`, spelled diff --git a/Sources/.editorconfig b/Sources/.editorconfig index a8e90a293..3614f5364 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -636,3 +636,6 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Calculus/DowncastingOffCalculusTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + +[Tests/UnitTests/Calculus/HalfOddPowerOfALogarithmTest.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index b19c03a23..a3c203e66 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -8536,6 +8536,161 @@ private static Entity SquarefreePartAsWritten(Entity denominator, Entity.Variabl return sum.InnerSimplified; } + /// + /// A power of x times a half-odd power of a logarithm, x^p F^n with F = + /// ln x or A + B ln(c x^r) and 2n odd, onto the Gaussian's moments. + /// With t = sqrt(F) and F' = s/x, dx = 2 x t dt/s, and + /// x^(p + 1) e^(-(p + 1) F/s) is a constant, so the integral is that constant times + /// 2/s times the integral of t^(2n + 1) e^((p + 1) t^2/s): an even moment of + /// the Gaussian, which the table answers. At p = -1 it is + /// F^(n + 1)/(s (n + 1)). Rubi's 3.1.2, (d x)^m (a + b ln(c x^n))^p for a + /// half-odd p. + /// https://github.com/asc-community/AngouriMath/issues/1501 + /// + /// + /// The constant is written x^(p + 1) e^(-(p + 1) F/s), as Rubi writes it, rather + /// than with ln(c x^r) split into ln c + r ln x, which holds only where + /// c is positive. The logarithm to the right of a whole power of it is + /// 's. + /// + internal static Entity? SolveAPowerTimesAHalfOddPowerOfTheLogarithm(Entity expr, Entity.Variable x) + { + Entity? p = null; + Number.Rational? n = null; + Entity constant = Number.Integer.One; + Entity? logarithm = null; + Dictionary? polynomial = null; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (factor == x && p is null) + p = underneath ? Number.Integer.MinusOne : Number.Integer.One; + else if (factor is Powf(var @base, var exponent) && @base == x && !exponent.ContainsNode(x) && p is null) + p = underneath ? -exponent : exponent; + else if (factor is Powf(Mulf(var scale, var scaled), var scaledExponent) && scaled == x && !scale.ContainsNode(x) && !scaledExponent.ContainsNode(x) && p is null) + { + p = underneath ? -scaledExponent : scaledExponent; + var scalePower = MathS.Pow(scale, scaledExponent); + constant = underneath ? constant / scalePower : constant * scalePower; + } + else if (n is null && factor is Powf(var f, Number.Rational half) && half.ERational.Denominator.Equals(EInteger.FromInt32(2))) + { + if (f == MathS.Ln(x) || IsAffineInALogarithmOfAPowerOfX(f, x)) + (logarithm, n) = (f, underneath ? Number.Rational.Create(-half.ERational) : half); + // The logarithm of a power of a linear, Rubi's 3.3: the same question under + // u = d + e x. Declined, it was integrated by parts against the half-odd + // power's own antiderivative, which took (f + g x) (a + b ln(c (d + e x)^n))^(3/2) + // 37 s to give up. + else if (LinearUnderALogarithm(f, x) is { } linear) + return ByTheLinearUnderTheLogarithm(expr, x, linear); + else + return null; + } + // A polynomial beside it, which is how a linear outside the logarithm arrives once + // the linear inside it has been substituted: summed a power at a time. + else if (polynomial is null && !underneath && TreeAnalyzer.TryGetPolynomial(factor, x, out var monomials) + && monomials.Count > 1 && monomials.Keys.All(power => power.Sign >= 0 && power.CompareTo(EInteger.FromInt32(12)) <= 0)) + polynomial = monomials; + else + return null; + } + if (logarithm is null || n is null) + return null; + p = (p ?? Number.Integer.Zero).InnerSimplified; + if (p.Evaled is Number.Complex and not Number.Real) + return null; + Entity slopeOfF = Number.Integer.One; + if (logarithm != MathS.Ln(x)) + { + slopeOfF = Functions.PartialFractions.Bare((logarithm.Differentiate(x) * x).Simplify()); + if (slopeOfF.ContainsNode(x) || slopeOfF.Evaled is Number.Complex { IsZero: true }) + return null; + } + if (polynomial is null) + return HalfOddLogarithmicTerm(p, logarithm, n, slopeOfF, x) is { } single ? (constant * single).InnerSimplified : null; + Entity sum = Number.Integer.Zero; + foreach (var monomial in polynomial) + { + if (HalfOddLogarithmicTerm((p + Number.Integer.Create(monomial.Key)).InnerSimplified, logarithm, n, slopeOfF, x) is not { } term) + return null; + sum += monomial.Value * term; + } + return (constant * sum).InnerSimplified; + } + + /// + /// The linear d + e x, with d not zero, whose powers and constant multiples + /// make up the argument of the one logarithm in , or + /// . + /// + private static Entity? LinearUnderALogarithm(Entity expr, Entity.Variable x) + { + Entity? argument = null; + foreach (var node in expr.Nodes) + if (node is Logf(var @base, var antilogarithm) && @base == MathS.e && antilogarithm.ContainsNode(x)) + { + if (argument is not null && argument != antilogarithm) + return null; + argument = antilogarithm; + } + if (argument is null) + return null; + foreach (var node in argument.Nodes) + if (node is Sumf or Minusf && TreeAnalyzer.TryGetPolyLinear(node, x, out var slope, out var intercept) + && !TreeAnalyzer.IsZero(slope) && !TreeAnalyzer.IsZero(intercept)) + { + var placeholder = Variable.CreateUnique(argument, "u_lin"); + var inThePlaceholder = argument.Replace(child => child == node ? placeholder : child); + return !inThePlaceholder.ContainsNode(x) && IsAProductOfPowersOfTheVariable(inThePlaceholder, placeholder) ? node : null; + } + return null; + } + + /// + /// under u = d + e x, with the + /// d + e x under its logarithm, answered as a power or a polynomial of u + /// beside a half-odd power of A + B ln(c u^r) and written back in x. + /// + private static Entity? ByTheLinearUnderTheLogarithm(Entity expr, Entity.Variable x, Entity linear) + { + if (!TreeAnalyzer.TryGetPolyLinear(linear, x, out var slope, out var intercept) || slope is null || intercept is null) + return null; + var u = Variable.CreateUnique(expr, "u"); + var inU = (expr.Replace(node => node == linear ? u : node).Substitute(x, (u - intercept) / slope) / slope).InnerSimplified; + if (inU.ContainsNode(x) || LinearUnderALogarithm(inU, u) is not null) + return null; + return SolveAPowerTimesAHalfOddPowerOfTheLogarithm(inU, u)?.Substitute(u, linear).InnerSimplified; + } + + /// + /// The integral of x^q F^n for a half-odd n and F' = s/x: see + /// . + /// + private static Entity? HalfOddLogarithmicTerm(Entity q, Entity logarithm, Number.Rational n, Entity slopeOfF, Entity.Variable x) + { + var qPlusOne = (q + 1).InnerSimplified; + if (qPlusOne.Evaled is Number.Complex { IsZero: true }) + return MathS.Pow(logarithm, n + 1) / ((n + 1) * slopeOfF); + // t^(2n + 1) e^((q + 1) t^2/s) in a variable of its own, for the table's moments. + var t = Variable.CreateUnique(logarithm, "t"); + var power = (2 * n.ERational + ERational.One).ToEInteger(); + if (!power.CanFitInInt32() || power.Abs().CompareTo(EInteger.FromInt32(32)) > 0) + return null; + var gaussian = MathS.Pow(MathS.e, (qPlusOne / slopeOfF).InnerSimplified * MathS.Sqr(t)); + var inT = power.IsZero ? gaussian : MathS.Pow(t, Number.Integer.Create(power)) * gaussian; + if (IntegralPatterns.TryStandardIntegrals(inT, t) is not { } answerInT) + return null; + // x^(q + 1) e^(-(q + 1) ln(x)) is 1 wherever ln(x) is real, so for ln(x) itself it is left out. + Entity theConstant = logarithm == MathS.Ln(x) + ? Number.Integer.One + : MathS.Pow(x, qPlusOne) * MathS.Pow(MathS.e, -qPlusOne * logarithm / slopeOfF); + return 2 / slopeOfF * theConstant * answerInT.Substitute(t, MathS.Sqrt(logarithm)); + } + /// /// Whether is A + B ln(c x^r) with A, B, /// c and r free of and the logarithm present: one diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index a3343ebcc..715a162fd 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -569,6 +569,8 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // A power of x times a power of its logarithm, by the closed reduction: exact, // and by parts n times where the exponent is a symbol was not taken. if ((answer = IndefiniteIntegralSolver.SolveAPowerTimesAPowerOfTheLogarithm(expr, x)) is { }) return answer; + // And a half-odd power of it, onto the Gaussian's moments by t = sqrt(F). + if ((answer = IndefiniteIntegralSolver.SolveAPowerTimesAHalfOddPowerOfTheLogarithm(expr, x)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveByFlatteningAPowerOfAnExponential(expr, x, integrateByParts)) is { }) return answer; // An exponential of a multiple of a logarithm is a power of the argument, which is // how every inverse hyperbolic function under an exponential arrives. diff --git a/Sources/Tests/UnitTests/Calculus/HalfOddPowerOfALogarithmTest.cs b/Sources/Tests/UnitTests/Calculus/HalfOddPowerOfALogarithmTest.cs new file mode 100644 index 000000000..a8e8ed9a4 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/HalfOddPowerOfALogarithmTest.cs @@ -0,0 +1,86 @@ +// +// 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 System.Linq; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A power of x times a half-odd power of a logarithm, x^p (A + B ln(c x^r))^n, onto the + /// Gaussian's moments by t = sqrt(F): Rubi's 3.1.2, (d x)^m (a + b ln(c x^n))^p with + /// p half an odd number, every one of which was left unevaluated. + /// https://github.com/asc-community/AngouriMath/issues/1501 + /// + [Trait("Area", "Calculus")] + public sealed class HalfOddPowerOfALogarithmTest + { + /// Past 1, where every logarithm here is real and positive for the pins below. + private static readonly double[] Points = { 1.3, 1.9, 2.6, 3.4 }; + private static readonly (string Name, string Value)[] Pins = + { ("a", "3/2"), ("b", "2/3"), ("c", "5/4"), ("n", "3"), ("m", "1/2"), ("d", "1/2"), ("k", "2"), ("f", "3/4"), ("g", "2"), ("h", "1/3"), ("p", "2"), ("q", "3/2") }; + + private static void DifferentiatesBack(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Entity Pinned(Entity e) => Pins.Aggregate(e, (current, pin) => current.Substitute(pin.Name, pin.Value.ToEntity())); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in Points) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + 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-9, $"d/dx of {integral} is {got} at x = {at}, where {integrand} is {want}"); + } + } + + [Theory] + [InlineData("sqrt(ln(x))")] + [InlineData("x*sqrt(ln(x))")] + [InlineData("x^2/sqrt(ln(x))")] + [InlineData("sqrt(ln(x))/x^2")] + [InlineData("x^3*ln(x)^(3/2)")] + [InlineData("1/(x^3*ln(x)^(3/2))")] + [InlineData("x/ln(x)^(5/2)")] + public void OfTheLogarithmItself(string integrand) => DifferentiatesBack(integrand); + + [Theory] + [InlineData("sqrt(ln(x^2))")] + [InlineData("sqrt(ln(a*x^n))")] + [InlineData("x^3*sqrt(ln(a*x^n))")] + [InlineData("x^m/sqrt(ln(a*x^n))")] + [InlineData("(a + b*ln(c*x^n))^(1/2)")] + [InlineData("x^3*ln(a*x^n)^(3/2)")] + [InlineData("x^2/ln(a*x^n)^(3/2)")] + [InlineData("x^m/ln(a*x^n)^(5/2)")] + public void OfALogarithmOfAPower(string integrand) => DifferentiatesBack(integrand); + + /// + /// The logarithm of a power of a linear, with a polynomial beside it, under + /// u = d + e x: Rubi's 3.3. Declined, these went by parts against the half-odd + /// power's own antiderivative and took up to 55 s to give up. + /// + [Theory] + [InlineData("sqrt(ln(2 + 3*x))")] + [InlineData("(f + g*x)*(a + b*ln(c*(d + k*x)^n))^(3/2)")] + [InlineData("(f + g*x)^2*sqrt(a + b*ln(c*(d + k*x)^n))")] + [InlineData("(f + g*x)^3/sqrt(a + b*ln(c*(d + k*x)^n))")] + [InlineData("(f + g*x)/(a + b*ln(c*(d + k*x)^n))^(5/2)")] + [InlineData("(g + h*x)^2/sqrt(a + b*ln(c*(d*(k + f*x)^p)^q))")] + public void OfALogarithmOfAPowerOfALinear(string integrand) => DifferentiatesBack(integrand); + + /// Beside 1/x it is elementary, F^(n + 1)/(s (n + 1)). + [Fact] + public void BesideTheReciprocalItIsAPower() + => Assert.Equal("2 * sqrt(ln(x)) + C".ToEntity(), "1/(x*sqrt(ln(x)))".ToEntity().Integrate("x")); + } +}