diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 7b1c16722..c8d388420 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1444,6 +1444,25 @@ multiple of it, `(d - c^2 d x^2)^(3/2)`, is written over it the same way, with Rubi's 5.1.4, 5.1.5, 5.2.4 and 5.2.5, all 504 problems that count: 405 to 463, no row lost, 77 timeouts to 19. +### A perfect square in a power of the variable is read as one + +`sqrt(a^2 + 2 a b x^2 + b^2 x^4) sqrt(c + e x + d x^2)` was left unevaluated. The radicand is +`(a + b x^2)^2`, so its root is the modulus `|a + b x^2|`, but the rule that reads a root of a +perfect square read only a quadratic in `x`. It now reads the same quadratic in a power of `x`, +`a x^(2k) + b x^k + c`, and puts `sgn(x^k + h)` in front of the integral, as it does for `k = 1`. +There is no sign where `k` is even and `h` is a positive number, since `x^k + h` is then +positive. The shift `h` is simplified where it holds a symbol: `2 a b/(2 b^2)` is `a/b`, for +`k = 1` as well. + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sqrt(c+pe*x+d*x^2)*sqrt(a^2+2*a*b*x^2+b^2*x^4)".Integrate("x")` | left unevaluated | the antiderivative, with `sgn(x^2 + a/b) sqrt(b^2)` in front | +| `"x/sqrt(a^2+2*a*b*x^3+b^2*x^6)".Integrate("x")` | left unevaluated | logarithms and an arctangent in `(a/b)^(1/3)`, with `sgn(x^3 + a/b)` in front | +| `"(1+2*x^2+x^4)^(3/2)".Integrate("x")` | left unevaluated | `x + 3 x^3/3 + 3 x^5/5 + x^7/7`, with no sign | + +Rubi's 1.2.2.2, 1.2.2.4, 1.2.2.7 and 1.2.3.2, the 574 problems that count with a perfect square +written with symbols: 182 to 388, no row lost, 19 timeouts to 9. + ### `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 9dd36ac1a..f2d62e9e1 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -10444,7 +10444,10 @@ private static bool IsAPerfectSquareDiscriminant(Entity a, Entity b, Entity c) /// A square root of a perfect square in is the modulus: /// sqrt(x^2) is |x|, which for a real x is sgn(x) x, and /// sqrt(a (x + h)^2) is sqrt(a) sgn(x + h) (x + h) for a positive number - /// a. Every rule that reads a root of a quadratic assumes it is not one: + /// a, and a square in a power of x the same way: + /// sqrt(a^2 + 2 a b x^2 + b^2 x^4) is sqrt(b^2) sgn(x^2 + a/b) (x^2 + a/b), + /// the square Rubi's 1.2.2 and 1.2.3 write. Every rule that reads a root of a quadratic + /// assumes it is not one: /// 1/sqrt(1 + csch(x)^2) under u = tanh(x) is sqrt(u^2)/(u^2 - 1), /// and the table rule for sqrt(a w^2 + b w + c) beside a linear, reached after the /// partial fractions and a reciprocal, answered it with a logarithm of zero. @@ -10467,7 +10470,8 @@ private static bool IsAPerfectSquareDiscriminant(Entity a, Entity b, Entity c) internal static Entity? SolveByTakingARootOfAPerfectSquare(Entity expr, Entity.Variable x, bool integrateByParts) { // Asked of every integrand at every depth, so the radicand is read as a polynomial - // only where it is small enough to be a written square: x^2, or a x^2 + b x + c. + // only where it is small enough to be a written square: x^2, a x^2 + b x + c, or + // the same quadratic in a power of x, a x^(2k) + b x^k + c. // At the top, any half-odd power of a perfect square is the power of the modulus; // below it the square root only, since a sign written for a substitution's variable // is a factor the rest of the search carries -- `(1 + 2u + u^2)^(5/2)` under @@ -10493,12 +10497,20 @@ bool MayBeARootOfASquare(Entity node, Entity.Variable x) var r = (Number.Rational)exponent; if (!TreeAnalyzer.TryGetPolynomial(radicand, x, out var monomials) || monomials.Count == 0) return node; + // A quadratic in w = x^k: a w^2 + b w + c, the square being of w + h. Beyond + // k = 1 only with the constant term written, since `a x^(2k)` alone is a power + // of a monomial, which another rule distributes. var degree = monomials.Keys.Max()!; - if (!degree.Equals(EInteger.FromInt32(2)) || monomials.Keys.Any(k => k.Sign < 0)) + if (degree.IsZero || !degree.IsEven || monomials.Keys.Any(k => k.Sign < 0)) return node; - var a = monomials.TryGetValue(EInteger.FromInt32(2), out var a2) ? a2 : Number.Integer.Zero; - var b = monomials.TryGetValue(EInteger.One, out var b1) ? b1 : Number.Integer.Zero; + var half = degree / EInteger.FromInt32(2); + if (monomials.Keys.Any(k => !k.IsZero && !k.Equals(half) && !k.Equals(degree)) + || !half.Equals(EInteger.One) && !monomials.ContainsKey(EInteger.Zero)) + return node; + var a = monomials[degree]; + var b = monomials.TryGetValue(half, out var b1) ? b1 : Number.Integer.Zero; var c = monomials.TryGetValue(EInteger.Zero, out var c0) ? c0 : Number.Integer.Zero; + Entity w = half.Equals(EInteger.One) ? x : MathS.Pow(x, Number.Integer.Create(half)); // A positive leading coefficient, since `sqrt(a (x + h)^2)` is `sqrt(a) |x + h|`: // a positive number, or one positive for a real parameter -- `b^2` is, and // `a^2 + 2 a b x + b^2 x^2` is the square Rubi writes -- whose condition travels @@ -10517,17 +10529,22 @@ bool MayBeARootOfASquare(Entity node, Entity.Variable x) if (b.Evaled is Number.Complex and not Number.Real || c.Evaled is Number.Complex and not Number.Real || !IsAPerfectSquareDiscriminant(a, b, c)) return node; - // a (x + h)^2 with h = b/(2a): (a (x + h)^2)^r is a^r |x + h|^(2r), and with 2r odd - // (a rational is held in lowest terms) that is a^r sgn(x + h) (x + h)^(2r). + // a (w + h)^2 with h = b/(2a): (a (w + h)^2)^r is a^r |w + h|^(2r), and with 2r odd + // (a rational is held in lowest terms) that is a^r sgn(w + h) (w + h)^(2r). var h = (b / (Number.Integer.Create(2) * a)).InnerSimplified; - var linear = h == Number.Integer.Zero ? x : (x + h).InnerSimplified; + // `2 a b/(2 b^2)` is `a/b`, which the normalisation does not cancel. + if (h.Vars.Any()) + h = Functions.PartialFractions.Bare(h.Simplify()); + var linear = h == Number.Integer.Zero ? w : (w + h).InnerSimplified; var twoR = Number.Integer.Create(r.ERational.Numerator); Entity power = twoR == Number.Integer.One ? linear : MathS.Pow(linear, twoR); // The sign in front of the integral rather than inside it: it is constant // between the linear factor's zeros, and a rule below that differentiates the // integrand cannot evaluate `derivative(sgn(...))` -- which is the exception // `sqrt(a^2 + 2abx + b^2x^2) sqrt(c + ex + dx^2)` threw with it left in place. - signs = signs * MathS.Signum(linear); + // None where it is plainly one: an even power of x plus a positive number. + if (!(half.IsEven && h.Evaled is Number.Real { IsPositive: true })) + signs = signs * MathS.Signum(linear); changed = true; return leading == Number.Integer.One ? power : MathS.Pow(leading, r) * power; }); diff --git a/Sources/Tests/UnitTests/Calculus/RootOfAPerfectSquareIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/RootOfAPerfectSquareIntegralTest.cs index b664da26a..834bd354e 100644 --- a/Sources/Tests/UnitTests/Calculus/RootOfAPerfectSquareIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/RootOfAPerfectSquareIntegralTest.cs @@ -92,6 +92,42 @@ Entity Pin(Entity e) => e.Substitute("a", 0.9).Substitute("b", 1.7).Substitute(" Assert.True(compared >= 5, $"only {compared} points were comparable for {integrand}"); } + /// + /// The same square in a power of the variable: a^2 + 2 a b x^2 + b^2 x^4 is + /// (a + b x^2)^2, and its root is sqrt(b^2) sgn(x^2 + a/b) (x^2 + a/b). Pinned + /// with a and b of opposite signs, so that the sign changes among the points. + /// Rubi's 1.2.2.7 and 1.2.3.2. + /// #718 + /// + [Theory] + [InlineData("sqrt(c + pe*x + d*x^2)*sqrt(a^2 + 2*a*b*x^2 + b^2*x^4)")] + [InlineData("x*sqrt(c + pe*x + d*x^2)*sqrt(a^2 + 2*a*b*x^2 + b^2*x^4)")] + [InlineData("x/sqrt(a^2 + 2*a*b*x^3 + b^2*x^6)")] + public void ASquareInAPowerOfTheVariable(string integrand) + { + var integral = integrand.ToEntity().Integrate("x").Substitute("C", 0); + Assert.DoesNotContain("integral(", integral.Stringize()); + Assert.DoesNotContain("NaN", integral.Stringize()); + Entity Pin(Entity e) => e.Substitute("a", -0.9).Substitute("b", 1.7).Substitute("c", 2.1) + .Substitute("pe", 0.3).Substitute("d", 1.3); + var derivative = Pin(integral).Differentiate("x"); + var original = Pin(integrand.ToEntity()); + 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) + 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}"); + } + Assert.True(compared >= 5, $"only {compared} points were comparable for {integrand}"); + } + [Fact] public void TheSignIsTheSignOfTheLinearFactor() {