From 949133dac7782eef9f003fd04b956888009fd2db Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 2 Oct 2026 09:22:14 +0000 Subject: [PATCH] A root of a square in a fractional power of x is the modulus as well A sum of powers of x that are whole multiples of one fractional k is read as a polynomial in w = x^k, which is exact for every x, and the root of a square in it is answered for a positive x, where w is real, as the modulus of a linear in w. Rubi's 1.2.3.2 (a^2 + b^2/x^(2/5) + 2 a b/x^(1/5))^(5/2) is a square in x^(-1/5), and was declined: 0.5 s. A root of an even order plus a positive number, sqrt(x) + 1, is not given a sign. Part of #718. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 15 ++++ .../Integration/IndefiniteIntegralSolver.cs | 75 ++++++++++++++++++- .../RootOfAPerfectSquareIntegralTest.cs | 31 ++++++++ 3 files changed, 117 insertions(+), 4 deletions(-) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 8564b93e5..16c6a39e8 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -2210,6 +2210,21 @@ coefficient of either sign ([#1670](https://github.com/asc-community/AngouriMath | `"1/sqrt(-a^2-2*a*b*x-b^2*x^2)".Integrate("x")` | a piecewise whose arcsine divides by `sqrt(0)`, `NaN` for every real `b` but zero | a logarithm over `sqrt(-b^2)`, `provided b^2 > 0` | | `"1/(x*sqrt(a^2+2*a*b*x+b^2*x^2))".Integrate("x")`, the same with `-a^2 - 2 a b x - b^2 x^2`, or with `d + e x` for `x` | left unevaluated | the antiderivative, `provided b^2 > 0` | +### A root of a square in a fractional power of `x` is the modulus of a linear in that power + +**Answers where there were none.** A sum of powers of `x` that are whole multiples of one +fractional `k`, such as `a^2 + 2 a b x^(1/3) + b^2 x^(2/3)`, is read as a polynomial in `w = x^k`, +and the root of a square in it as the modulus of a linear in `w`, as the root of a square +quadratic is. `w` is real for a positive `x`, and the answer says so, `provided x > 0`. Rubi's +1.2.3.2 ([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sqrt(a^2+2*a*b*x^(1/3)+b^2*x^(2/3))".Integrate("x")` | left unevaluated | `sgn(x^(1/3) + a/b) sqrt(b^2) (3/4 x^(4/3) + a/b x) provided x > 0 and b^2 > 0` | +| `"1/sqrt(a^2+b^2*x+2*a*b*sqrt(x))".Integrate("x")` | left unevaluated | `sgn(sqrt(x) + a/b) (2 sqrt(x) - 2 a/b ln(a/b + sqrt(x)))/sqrt(b^2)`, under the same condition | +| `"(a^2+b^2/x^(2/5)+2*a*b/x^(1/5))^(5/2)".Integrate("x")`, Rubi 1.2.3.2 #659 | left unevaluated | an antiderivative in `x^(-1/5)`, under the same condition | +| `"(1+2*x^(1/2)+x)^(3/2)".Integrate("x")` | left unevaluated | `2 (x/2 + x^(3/2) + 3 x^2/4 + x^(5/2)/5) provided x > 0` | + ### 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: diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index a801299a2..89f23ce05 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -13779,8 +13779,19 @@ bool MayBeARootOfASquare(Entity node, Entity.Variable x) return node; var (radicand, exponent) = (Powf)node; var r = (Number.Rational)exponent; + // A quadratic in a fractional power of x, `a^2 + 2 a b x^(-1/5) + b^2 x^(-2/5)`, is + // read in `w = x^k`, exact for every x since each power of x is a whole power of w; + // and real only for a positive x, which the answer then says. + Entity wBase = x; if (!TreeAnalyzer.TryGetPolynomial(radicand, x, out var monomials) || monomials.Count == 0) - return node; + { + if (!atTheTop || InAFractionalPowerOfX(radicand, x) is not var (inW, k)) + return node; + monomials = inW; + wBase = MathS.Pow(x, k); + var positive = new Greaterf(x, Number.Integer.Zero); + assumed = assumed == Entity.Boolean.True ? positive : assumed & positive; + } // 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. @@ -13799,7 +13810,7 @@ bool MayBeARootOfASquare(Entity node, Entity.Variable x) 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)); + Entity w = half.Equals(EInteger.One) ? wBase : MathS.Pow(wBase, Number.Integer.Create(half)); // A leading coefficient of known sign, since `sqrt(a (x + h)^2)` is // `sqrt(a) |x + h|` for either: a number, or one positive or negative for a real // parameter -- `b^2` is, and `a^2 + 2 a b x + b^2 x^2` is the square Rubi writes, @@ -13831,8 +13842,10 @@ bool MayBeARootOfASquare(Entity node, Entity.Variable x) // 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. - // None where it is plainly one: an even power of x plus a positive number. - if (!whole && !(half.IsEven && h.Evaled is Number.Real { IsPositive: true })) + // None where it is plainly one: an even power of x plus a positive number, or a + // root of x of an even order plus one, `sqrt(x) + 1`. + if (!whole && !((half.IsEven || wBase is Powf(_, Number.Rational { ERational.Denominator.IsEven: true })) + && 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; @@ -13844,6 +13857,60 @@ bool MayBeARootOfASquare(Entity node, Entity.Variable x) return assumed == Entity.Boolean.True ? answer : answer.Provided(assumed); } + /// + /// as a polynomial in w = x^k for a rational k that is + /// not whole, each of its terms a constant times a power of that is a + /// whole, non-negative multiple of k; otherwise. + /// + private static (Dictionary Monomials, Number.Rational K)? InAFractionalPowerOfX(Entity sum, Entity.Variable x) + { + var terms = new List<(ERational Power, Entity Coefficient)>(); + foreach (var term in Sumf.LinearChildren(sum)) + { + var power = ERational.Zero; + Entity coefficient = Number.Integer.One; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(term)) + { + if (!factor.ContainsNode(x)) + { + coefficient = underneath ? coefficient / factor : coefficient * factor; + continue; + } + var exponent = factor == x ? ERational.One + : factor is Powf(var @base, Number.Rational q) && @base == x ? q.ERational + : (ERational?)null; + if (exponent is null) + return null; + power = underneath ? power.Subtract(exponent) : power.Add(exponent); + } + terms.Add((power, coefficient)); + } + var powers = terms.Select(term => term.Power).Where(q => !q.IsZero).Distinct().ToList(); + if (powers.Count == 0 || powers.All(q => q.IsInteger()) || powers.Any(q => q.Sign > 0) && powers.Any(q => q.Sign < 0)) + return null; + // k the greatest common divisor of the powers, with their sign. + var numerators = EInteger.Zero; + var denominators = EInteger.One; + foreach (var q in powers) + { + numerators = numerators.Gcd(q.Numerator.Abs()); + denominators = denominators.Multiply(q.Denominator).Divide(denominators.Gcd(q.Denominator)); + } + var k = ERational.Create(numerators, denominators); + if (powers[0].Sign < 0) + k = k.Negate(); + var monomials = new Dictionary(); + foreach (var (power, coefficient) in terms) + { + var degree = power.Divide(k); + if (!degree.IsInteger() || degree.Sign < 0) + return null; + var key = degree.ToEInteger(); + monomials[key] = monomials.TryGetValue(key, out var so) ? so + coefficient : coefficient; + } + return (monomials, Number.Rational.Create(k)); + } + /// /// A whole power of a product in which the variable stands beside a constant, written /// as the product of the powers: (c x)^2 is c^2 x^2, exact for a whole diff --git a/Sources/Tests/UnitTests/Calculus/RootOfAPerfectSquareIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/RootOfAPerfectSquareIntegralTest.cs index 80eea18dd..6494308f7 100644 --- a/Sources/Tests/UnitTests/Calculus/RootOfAPerfectSquareIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/RootOfAPerfectSquareIntegralTest.cs @@ -228,6 +228,37 @@ Entity Pin(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute(" Assert.True(compared >= 6, $"only {compared} points were comparable for {integrand}"); } + /// + /// A square in a fractional power of x, Rubi's 1.2.3.2 (a^2 + b^2/x^(2/5) + 2 a b/x^(1/5))^(5/2): + /// read in w = x^k, every power of x a whole power of w, and real for a positive x, + /// which the answer says. Differentiated back where it says, with a b of either + /// sign, so that w + a/b changes sign among the points for the first two. + /// + [Theory] + [InlineData("(a^2 + b^2/x^(2/5) + 2*a*b/x^(1/5))^(5/2)")] + [InlineData("sqrt(a^2 + 2*a*b*x^(1/3) + b^2*x^(2/3))")] + [InlineData("(1 + 2*x^(1/2) + x)^(3/2)")] + [InlineData("sqrt(4 + 4*x^(1/4) + x^(1/2))/x")] + public void ASquareInAFractionalPowerOfTheVariable(string integrand) + { + var integral = integrand.ToEntity().Integrate("x").Substitute("C", 0); + Assert.DoesNotContain("integral(", integral.Stringize()); + foreach (var b in new[] { 0.7, -0.7 }) + { + Entity Pin(Entity e) => e.Substitute("a", 1.3).Substitute("b", b); + var derivative = Pin(integral).Differentiate("x"); + var original = Pin(integrand.ToEntity()); + foreach (var at in new[] { 0.01, 0.31, 0.83, 2.41, 9.0 }) + { + 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); + Assert.True(difference / Math.Max(1.0, Math.Abs((double)want.RealPart)) < 1e-9, + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, b = {b}, where the integrand is {want}: {integral}"); + } + } + } + [Fact] public void TheSignIsTheSignOfTheLinearFactor() {