diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 9ad713119..c75e5a311 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -747,6 +747,35 @@ they are written ([#718](https://github.com/asc-community/AngouriMath/issues/718 | `"1/((a + b*x^2)^(9/2)*(1 + x^2))".ToEntity().Integrate("x")` | `integral(...)` | the same, in 1 s | | `"x^2/((a + b*x^2)^(7/2)*(c + d*x^2))".ToEntity().Integrate("x")` | `integral(...)` | the same by the sign of `c (a d - b c)`, in 1 s | +### A binomial differential with symbols in it, or a power of `x` that is not whole, is integrated + +**Answers where there were none.** Chebyshev's theorem says when `x^m (a + b x^n)^(p/q)` has an +elementary antiderivative, and it is about the exponents only. The rule for it read whole `m` and +`n` and rational numbers for `a` and `b`, and declined the rest. It reads rational `m` and `n` and +anything free of `x` for `a` and `b` now, and a power of a multiple of `x`, `(c x)^(5/2)`, as +`x^(5/2)` times `(c x)^(5/2)/x^(5/2)`, which is constant on either side of zero. Where `m + 1 + +n (p/q + 1)` is zero the answer is the one product of powers it is. With a symbol in the +coefficients the rule is asked after the rules for a root of a quadratic and for a rational +function of `x^n` beside the root of its binomial, which answer what they share with it more +shortly. Rubi's 1.1.2.2 and 1.1.3.2, `(c x)^m (a + b x^n)^p` +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x^2/(a + b*x^4)^(3/4)".ToEntity().Integrate("x")` | `integral(...)` | logarithms and an arctangent of `(a + b x^4)^(1/4)/x` | +| `"x^6*(a + b*x^4)^(1/4)".ToEntity().Integrate("x")` | `integral(...)` | the same, beside a rational function of it | +| `"x^(7/3)*(a + b*x^2)^(1/3)".ToEntity().Integrate("x")` | `integral(...)` | a rational function, two logarithms and an arctangent of `(a + b x^2)^(1/3)/x^(2/3)`, written in `x^(1/3)` | +| `"(c*x)^(5/2)/(a - b*x^2)^(3/4)".ToEntity().Integrate("x")` | `integral(...)` | `(c x)^(5/2)/x^(5/2)` times the same in `(a - b x^2)^(1/4)/x^(1/2)` | +| `"(a - b*x^2)^(1/4)/(c*x)^(15/2)".ToEntity().Integrate("x")` | `integral(...)` | powers of `(a - b x^2)^(1/4)/(c x)^(1/2)` | + +Of Rubi's 3,163 problems in those two files with an answer in functions the library has, 96 more +are answered than on the unreleased master, and none fewer. Nine that the unreleased master answered +are written otherwise: a rule ahead of this one takes the integrand apart and asks a sub-integral +that this one answers now, so the rule ahead finishes where it used to decline. Four of the nine come +out shorter and five longer -- `1/((c x)^(5/3) (a + b x^2)^(2/3))` was +`-3/(2 a) x (c x)^(-5/3) (a + b x^2)^(1/3)` there, and is the same function written in `x^(1/3)`. +2.5.0 declined all nine. + ### A power of a multiple of a quadratic's derivative beside a power of the quadratic is a binomial **Improvement, not silent.** `(b d + 2 c d x)^m (a + b x + c x^2)^p`, with a power that is not whole, diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 14787aa3f..2ac4d83e4 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -4468,7 +4468,7 @@ Entity ConstantNegative() /// /// The third case comes down to the second. Where s + p/q is whole instead, /// x = 1/y turns x^m (a + b x^n)^(p/q) dx into - /// -y^m' (b + a y^n)^(p/q) dy with m' = -m - 2 - n p/q, a whole number, and + /// -y^m' (b + a y^n)^(p/q) dy with m' = -m - 2 - n p/q and /// (m' + 1)/n = -(s + p/q), whole — so it is the second case in y, with the /// roles of a and b exchanged. Its u, whose q-th power is /// b + a/x^n, is put back as (a + b x^n)^(1/q)/x^(n/q), which has that power @@ -4482,6 +4482,31 @@ Entity ConstantNegative() /// antiderivative at all, which is worth knowing before anyone goes looking. /// /// + /// The theorem is about the exponents, and m and n are read as rationals and + /// a and b as anything free of the variable, since neither substitution + /// reads the coefficients: x^(7/3) (a + b x^2)^(1/3) is the third case. A power of + /// a multiple of the variable, (c x)^(5/2), is read as x^(5/2) times + /// (c x)^(5/2)/x^(5/2), whose derivative is zero: it is c^(5/2) for a + /// positive c x and constant on either side of zero whatever the signs, so it + /// stands outside the integral as written. Rubi's 1.1.2.2 and 1.1.3.2, + /// (c x)^m (a + b x^n)^p, are written so. + /// + /// + /// Where m + 1 + n (p/q + 1) is zero, a case of the third, the answer is one product + /// of powers, x^(m+1) (a + b x^n)^(p/q+1)/(a (m + 1)), and is given as that. + /// finds it too, but is asked at the + /// top and one level below only, and below that the third case wrote it out through the + /// rational integrator, many terms where one will do. + /// + /// + /// Here with numbers only. With a symbol in a, b or the multiple it is + /// , asked after the rules written for + /// a root of a quadratic and for a rational function of x^n beside the root of its + /// binomial, which answer what they share with this more shortly: asked first, + /// 1/(a - b x^4)^(1/4) came out as two logarithms and two arctangents where + /// gives two arctangents. + /// + /// /// Closed in one step: the polynomial case is a sum of powers, and the rational case /// goes to the rational integrator directly, not back into the chain -- which is what /// lets it be volunteered at any depth rather than asked at the top only @@ -4492,29 +4517,46 @@ Entity ConstantNegative() /// https://github.com/asc-community/AngouriMath/issues/718 /// internal static Entity? SolveABinomialDifferential(Entity expr, Entity.Variable x) + => SolveABinomialDifferential(expr, x, withSymbols: false); + + /// + /// where a symbol stands + /// in a, b or the multiple, and nowhere else: what the first asking left + /// for the rules after it, and they declined. + /// + internal static Entity? SolveABinomialDifferentialWithSymbols(Entity expr, Entity.Variable x) + => SolveABinomialDifferential(expr, x, withSymbols: true); + + private static Entity? SolveABinomialDifferential(Entity expr, Entity.Variable x, bool withSymbols) { if (!TryReadABinomialDifferential(expr, x, out var power, out var exponent, out var inner, out var free, out var leading, out var factor)) return null; + if (withSymbols != (free.Vars.Any() || leading.Vars.Any() || factor.Vars.Any(symbol => symbol != x))) + return null; var u = Variable.CreateUnique(expr, "u_binom"); - var bracket = (Number.Rational.Create(free) - + Number.Rational.Create(leading) * MathS.Pow(x, Number.Integer.Create(inner))).InnerSimplified; + var bracket = (free + leading * MathS.Pow(x, Number.Rational.Create(inner))).InnerSimplified; var root = MathS.Pow(bracket, Number.Rational.Create(EInteger.One, exponent.Denominator)); + // m + 1 + n (p/q + 1) = 0, a case of the third with one term for its answer: + // (x^(m+1) (a + b x^n)^(p/q+1))' is (m + 1) a x^m (a + b x^n)^(p/q) exactly there. + var raisedPower = power.Add(ERational.One); + if (!raisedPower.IsZero && raisedPower.Add(inner.Multiply(exponent.Add(ERational.One))).IsZero) + return (factor * MathS.Pow(x, Number.Rational.Create(raisedPower)) + * MathS.Pow(bracket, Number.Rational.Create(exponent.Add(ERational.One))) + / (free * Number.Rational.Create(raisedPower))).InnerSimplified; + // The second case: s = (m + 1)/n whole. - if ((power + 1) % inner == 0) + var s = raisedPower.Divide(inner); + if (s.IsInteger()) return IntegrateABinomialDifferentialInTheSecondCase( power, inner, exponent, free, leading, u, root, factor); // The third: s + p/q whole, taken to the second by x = 1/y. - var sPlusP = ERational.Create(power + 1, inner).Add(exponent); - if (!sPlusP.IsInteger()) + if (!s.Add(exponent).IsInteger()) return null; - var nTimesP = exponent.Multiply(EInteger.FromInt32(inner)); - if (!nTimesP.IsInteger() || !nTimesP.Numerator.CanFitInInt32()) - return null; - var reflectedPower = -power - 2 - nTimesP.ToLowestTerms().Numerator.ToInt32Unchecked(); - var back = root / MathS.Pow(x, Number.Rational.Create(ERational.Create(EInteger.FromInt32(inner), exponent.Denominator))); + var reflectedPower = power.Negate().Subtract(ERational.FromInt32(2)).Subtract(inner.Multiply(exponent)); + var back = root / MathS.Pow(x, Number.Rational.Create(inner.Divide(ERational.FromEInteger(exponent.Denominator)))); if (IntegrateABinomialDifferentialInTheSecondCase( reflectedPower, inner, exponent, leading, free, u, back, Number.Integer.MinusOne) is not { } reflected) return null; @@ -4529,18 +4571,21 @@ Entity ConstantNegative() /// s <= 0. /// private static Entity? IntegrateABinomialDifferentialInTheSecondCase( - int power, int inner, ERational exponent, ERational free, ERational leading, + ERational power, ERational inner, ERational exponent, Entity free, Entity leading, Entity.Variable u, Entity back, Entity factor) { - if ((power + 1) % inner != 0) + var whole = power.Add(ERational.One).Divide(inner); + if (!whole.IsInteger()) + return null; + var wholeInteger = whole.ToEIntegerIfExact(); + if (!wholeInteger.CanFitInInt32()) return null; - var s = (power + 1) / inner; + var s = wholeInteger.ToInt32Checked(); var q = exponent.Denominator.ToInt32Checked(); var p = exponent.Numerator.ToInt32Checked(); - var a = Number.Rational.Create(free); - var b = Number.Rational.Create(leading); - var outside = Number.Integer.Create(q) - / (Number.Integer.Create(inner) * MathS.Pow(b, Number.Integer.Create(s))); + var a = free; + var b = leading; + var outside = Number.Integer.Create(q) / (Number.Rational.Create(inner) * ToThe(b, s)); if (s < 1) { @@ -4570,16 +4615,23 @@ Entity ConstantNegative() if (raised == 0) return null; // the power rule would divide by zero; not a polynomial after all total += Number.Integer.Create(binomial) - * MathS.Pow(-a, Number.Integer.Create(s - 1 - i)) + * ToThe(-a, s - 1 - i) * MathS.Pow(u, Number.Integer.Create(raised)) / Number.Integer.Create(raised); binomial = binomial * (s - 1 - i) / (i + 1); } return (factor * outside * total.Substitute(u, back)).InnerSimplified; + + // The zeroth power as the one it is: with a symbol in the base, `b^0` is simplified + // to 1 provided `b` is not zero, a condition the answer has no use for. + static Entity ToThe(Entity @base, int power) + => power == 0 ? Number.Integer.One : MathS.Pow(@base, Number.Integer.Create(power)); } /// - /// Reads as a rational multiple of x^m (a + b x^n)^(p/q), - /// with q above one and n at least two. + /// Reads as a multiple of x^m (a + b x^n)^(p/q), with + /// m and n rational, q above one, n neither zero nor one, and + /// the multiple, a and b free of the variable -- or the multiple + /// (c x^k)^r/x^(k r), whose derivative is zero, of a power of a monomial. /// /// /// A whole exponent on the bracket is a polynomial and wants expanding rather than this; @@ -4588,19 +4640,19 @@ Entity ConstantNegative() /// shortly. /// private static bool TryReadABinomialDifferential( - Entity expr, Entity.Variable x, out int power, out ERational exponent, - out int inner, out ERational free, out ERational leading, out Entity factor) + Entity expr, Entity.Variable x, out ERational power, out ERational exponent, + out ERational inner, out Entity free, out Entity leading, out Entity factor) { - power = 0; + power = ERational.Zero; exponent = ERational.Zero; - inner = 0; - free = ERational.Zero; - leading = ERational.Zero; + inner = ERational.Zero; + free = 0; + leading = 0; factor = 1; Entity? bracket = null; var exponentFound = ERational.Zero; - var powerFound = 0; + var powerFound = ERational.Zero; Entity constantFactor = 1; if (!Read(expr, 1) || bracket is null) return false; @@ -4622,18 +4674,18 @@ private static bool TryReadABinomialDifferential( return false; if (constantPart is null || monomial is null) return false; - if (!TryReadAMonomial(monomial, x, out var degree, out var coefficient) || degree < 2) + if (!TryReadAMonomial(monomial, x, out var degree, out var coefficient) + || degree.IsZero || degree.CompareTo(ERational.One) == 0) return false; - if (constantPart.Evaled is not Number.Rational constantValue - || coefficient.Evaled is not Number.Rational coefficientValue - || coefficientValue.ERational.IsZero) + if (VanishesIdentically(constantPart) || VanishesIdentically(coefficient)) return false; power = powerFound; exponent = lowest; inner = degree; - free = constantValue.ERational; - leading = coefficientValue.ERational; + // A number as the number it is, so that the answer is written as it was for numbers. + free = constantPart.Evaled is Number.Rational freeValue ? freeValue : constantPart.InnerSimplified; + leading = coefficient.Evaled is Number.Rational leadingValue ? leadingValue : coefficient.InnerSimplified; factor = constantFactor; return true; @@ -4641,27 +4693,45 @@ bool Read(Entity node, int multiplicity) { switch (node) { + case var constant when !constant.ContainsNode(x): + constantFactor = multiplicity > 0 + ? constantFactor * MathS.Pow(constant, multiplicity) + : constantFactor / MathS.Pow(constant, -multiplicity); + return true; case Variable v when v == x: - powerFound += multiplicity; + powerFound = powerFound.Add(ERational.FromInt32(multiplicity)); return true; case Mulf(var left, var right): return Read(left, multiplicity) && Read(right, multiplicity); case Divf(var above, var below): return Read(above, multiplicity) && Read(below, -multiplicity); - case Powf(var @base, Number.Integer whole) - when @base == x && whole.EInteger.CanFitInInt32(): - powerFound += multiplicity * whole.EInteger.ToInt32Checked(); + case Powf(var @base, Number.Rational raised) + when TryReadAMonomial(@base, x, out var ofTheBase, out _): + // x^r, or (c x^k)^r as x^(k r) times (c x^k)^r/x^(k r): that is + // c^r for a positive c x^k, and constant on either side of zero whatever + // the signs, so it stands outside as written. + var raisedHere = raised.ERational.Multiply(ERational.FromInt32(multiplicity)); + var ofTheVariable = ofTheBase.Multiply(raised.ERational); + powerFound = powerFound.Add(ofTheBase.Multiply(raisedHere)); + if (@base != x) + { + var standingOutside = node / MathS.Pow(x, Number.Rational.Create(ofTheVariable)); + constantFactor = multiplicity > 0 + ? constantFactor * MathS.Pow(standingOutside, multiplicity) + : constantFactor / MathS.Pow(standingOutside, -multiplicity); + } return true; - case Powf(var @base, Number.Rational raised) when @base.ContainsNode(x): + case Powf(var @base, Number.Rational raised): if (bracket is not null && bracket != @base) return false; bracket = @base; - exponentFound += ERational.FromInt32(multiplicity) * raised.ERational; + exponentFound = exponentFound.Add(ERational.FromInt32(multiplicity).Multiply(raised.ERational)); return true; - case Number.Rational rational when node is not Powf: - constantFactor = multiplicity > 0 - ? constantFactor * MathS.Pow(rational, multiplicity) - : constantFactor / MathS.Pow(rational, -multiplicity); + case Sumf or Minusf: + if (bracket is not null && bracket != node) + return false; + bracket = node; + exponentFound = exponentFound.Add(ERational.FromInt32(multiplicity)); return true; default: return false; @@ -4669,18 +4739,21 @@ bool Read(Entity node, int multiplicity) } } - /// Reads c * x^n, giving the degree and the coefficient. - private static bool TryReadAMonomial(Entity term, Entity.Variable x, out int degree, out Entity coefficient) + /// + /// Reads c * x^n, giving the degree, any rational, and the coefficient, anything + /// free of the variable. + /// + private static bool TryReadAMonomial(Entity term, Entity.Variable x, out ERational degree, out Entity coefficient) { - degree = 0; + degree = ERational.Zero; coefficient = 1; switch (term) { case Variable v when v == x: - degree = 1; + degree = ERational.One; return true; - case Powf(var @base, Number.Integer whole) when @base == x && whole.EInteger.CanFitInInt32(): - degree = whole.EInteger.ToInt32Checked(); + case Powf(var @base, Number.Rational raised) when @base == x: + degree = raised.ERational; return true; case Mulf(var left, var right) when !left.ContainsNode(x): if (!TryReadAMonomial(right, x, out degree, out var fromRight)) @@ -4692,6 +4765,17 @@ private static bool TryReadAMonomial(Entity term, Entity.Variable x, out int deg return false; coefficient = right * fromLeft; return true; + case Divf(var above, var below) when !below.ContainsNode(x): + if (!TryReadAMonomial(above, x, out degree, out var fromAbove)) + return false; + coefficient = fromAbove / below; + return true; + case Divf(var above, var below) when !above.ContainsNode(x): + if (!TryReadAMonomial(below, x, out var belowDegree, out var fromBelow)) + return false; + degree = belowDegree.Negate(); + coefficient = above / fromBelow; + return true; default: return false; } diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 8cae2ee5f..3086fd9c7 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -1028,6 +1028,10 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // same reciprocal: `1/(x^2 sqrt(Q))` is a polynomial over the root of the reversed // quadratic, which the rules for those answer. if ((answer = IndefiniteIntegralSolver.SolveByTheReciprocalBesideARootOfAQuadratic(expr, x, integrateByParts)) is { }) return answer; + // The binomial differential with a symbol in its coefficients, asked here rather than + // beside the one with numbers: the rules since then answer what they share with it + // more shortly, `1/(a - b x^4)^(1/4)` by two arctangents where this gives four terms. + if ((answer = IndefiniteIntegralSolver.SolveABinomialDifferentialWithSymbols(expr, x)) is { }) return answer; // A rational function of x and one cube root of a polynomial: no substitution // rationalises it, and the elementary ones are logarithms of `L - y` for linear L // whose cube agrees with the polynomial at the poles, found by an ansatz. diff --git a/Sources/Tests/UnitTests/Calculus/SymbolicBinomialDifferentialTest.cs b/Sources/Tests/UnitTests/Calculus/SymbolicBinomialDifferentialTest.cs new file mode 100644 index 000000000..58633629d --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/SymbolicBinomialDifferentialTest.cs @@ -0,0 +1,85 @@ +// +// 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 +{ + /// + /// The binomial differential x^m (a + b x^n)^(p/q) with symbols for a and + /// b, with a power of x that is not whole, or with a power of a multiple of + /// x in place of the power of x. Chebyshev's theorem is about the exponents, + /// and each of these was declined. + /// #718 + /// + /// + /// Checked by differentiating back with a = 2.3, b = 0.7 and the given + /// c, at points where the integrand is real -- on both sides of zero wherever it is + /// real on both, an odd root of a negative being real here. + /// + [Trait("Area", "Calculus")] + public sealed class SymbolicBinomialDifferentialTest + { + private static void DifferentiatesBack(string integrand, double c, double[] points) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Assert.DoesNotContain("NaN", integral.Stringize()); + Entity Pinned(Entity e) => e.Substitute("a", 2.3).Substitute("b", 0.7).Substitute("c", c); + 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(); + Assert.True(Math.Abs((double)want.ImaginaryPart) < 1e-12, $"the integrand at x = {at} is {want}, not real"); + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + } + + /// + /// A whole power of x beside a root of a binomial with symbols in it, in the third + /// case: x^2/(a + b x^4)^(3/4) is (m + 1)/n + p/q = 0. + /// + [Theory] + [InlineData("x^2/(a + b*x^4)^(3/4)")] + [InlineData("(a + b*x^4)^(1/4)/x^2")] + [InlineData("x^6*(a + b*x^4)^(1/4)")] + public void SymbolsInTheBinomial(string integrand) + => DifferentiatesBack(integrand, 1.3, new[] { -1.4, -0.6, 0.5, 1.2 }); + + /// + /// A power of x that is not whole: x^(7/3) (a + b x^2)^(1/3) is the third + /// case, (m + 1)/n + p/q = 2. + /// + [Theory] + [InlineData("x^(7/3)*(a + b*x^2)^(1/3)")] + [InlineData("x^(1/3)*(a + b*x^2)^(1/3)")] + public void APowerOfXThatIsNotWhole(string integrand) + => DifferentiatesBack(integrand, 1.3, new[] { -1.4, -0.6, 0.3, 0.7, 1.9 }); + + /// + /// A power of a multiple of x, read as the power of x times + /// (c x)^r/x^r, which stands outside: checked for a positive c on the + /// positive side, and for a negative one on the negative side, where c x is + /// positive again. + /// + [Theory] + [InlineData("(c*x)^(7/3)/(a + b*x^2)^(2/3)", 1.3, new[] { 0.3, 0.7, 1.2, 1.6 })] + [InlineData("(c*x)^(7/3)/(a + b*x^2)^(2/3)", -1.3, new[] { -1.6, -1.2, -0.7, -0.3 })] + [InlineData("(c*x)^(5/2)/(a - b*x^2)^(3/4)", 1.3, new[] { 0.3, 0.7, 1.2, 1.6 })] + [InlineData("(c*x)^(5/2)/(a - b*x^2)^(3/4)", -1.3, new[] { -1.6, -1.2, -0.7, -0.3 })] + [InlineData("(a - b*x^2)^(1/4)/(c*x)^(15/2)", 1.3, new[] { 0.3, 0.7, 1.2, 1.6 })] + [InlineData("(a - b*x^2)^(1/4)/(c*x)^(15/2)", -1.3, new[] { -1.6, -1.2, -0.7, -0.3 })] + public void APowerOfAMultipleOfX(string integrand, double c, double[] points) + => DifferentiatesBack(integrand, c, points); + } +}