diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 050721360..a9b52ffd3 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -816,6 +816,26 @@ answered too: an odd root is the real root | `"1/sqrt((a + b*x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 100 characters | | `"1/sqrt(1 - coth(x)^2)".ToEntity().Integrate("x")` | `integral(...)` | 301 characters | +### A half-odd power of a negative multiple keeps its sign through the long division + +**Silent: a wrong answer, and answers where there were none.** `sqrt(-x) (x + 1)/x` was integrated as +`2i sqrt(x)^3/3 + 2i sqrt(x)` and simplified to `i (sqrt(x) + x^(-1/2))`, each the negative of the right value for +every `x < 0`, where the integrand is real. The long division read `(c x)^r` as `c^r x^r` for any constant `c` and a +power `r` that is not whole, which holds for `c > 0` and not for `c < 0`: there `c x > 0` and `c^r` and `x^r` are both +imaginary. Such a power now comes apart over a positive number only. `(b x)^(5/2) (x^2 - 1)/x` and Rubi's 4.3.0 +`(b sec(x))^m tan(x)^3` and `cot(x)^3 (b csc(x))^m` at a half-odd `m` were declined, and are answered with `b` under +the power, right for either sign of `b` +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sqrt(-x)*(x+1)/x".ToEntity().Integrate("x")` | `2i * (sqrt(x)) ^ 3 / 3 + 2i * sqrt(x) + C`, wrong for x < 0 | `(sqrt(-x) * x / (3/2) + sqrt(-x) / (1/2) provided not x = 0) + C` | +| `"(-3*x)^(5/2)*(x+1)/x".ToEntity().Integrate("x")` | `((-3) ^ (5/2) * x ^ (7/2) / (7/2) + (-3) ^ (5/2) * x ^ (5/2) / (5/2) provided not x = 0) + C`, wrong for x < 0 | `(((-3) * x) ^ (5/2) * x / (7/2) + ((-3) * x) ^ (5/2) / (5/2) provided not x = 0) + C` | +| `"sqrt(-x)*(x+1)/x".ToEntity().Simplify()` | `i * (sqrt(x) + x ^ (-1/2))`, wrong for x < 0 | `sqrt(-x) * (x + 1) / x` | +| `"(b*x)^(5/2)*(x^2-1)/x".ToEntity().Integrate("x")` | `integral(...)` | 86 characters | +| `"(b*sec(x))^(5/2)*tan(x)^3".ToEntity().Integrate("x")` | `integral(...)` | 96 characters | +| `"cot(x)^3*(b*csc(x))^(3/2)".ToEntity().Integrate("x")` | `integral(...)` | 99 characters | + ### A power of x past the twelfth beside a block is split off **Answers where there were none.** `(c + d x^3 + g x^6 + f x^9)/(x^14 (a + b x^3))` was declined, and so were diff --git a/Sources/AngouriMath/Functions/Continuous/Limits/Solvers/Limit.Solvers.cs b/Sources/AngouriMath/Functions/Continuous/Limits/Solvers/Limit.Solvers.cs index 94912b687..06a5ab1db 100644 --- a/Sources/AngouriMath/Functions/Continuous/Limits/Solvers/Limit.Solvers.cs +++ b/Sources/AngouriMath/Functions/Continuous/Limits/Solvers/Limit.Solvers.cs @@ -25,7 +25,7 @@ internal static class LimitSolvers return null; var monomials = AnalyticalSolving.PolynomialSolver.GatherMonomialInformation - (children, x); + (children, x, aPositiveVariable: true); if (monomials is null) return null; var filteredDictionary = new Dictionary(); foreach (var monomial in monomials) diff --git a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs index b0a51b193..48ff84912 100644 --- a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/PolynomialSolver.cs @@ -262,15 +262,19 @@ Entity GetMonomialByPower(EInteger power) => } /// Finds all terms of a polynomial + /// The terms of the sum. + /// The variable the polynomial is in. + /// Whether the caller reads the terms only where the + /// variable is positive: see . /// if polynomial is bad - internal static Dictionary? GatherMonomialInformation(IEnumerable terms, Variable subtree) + internal static Dictionary? GatherMonomialInformation(IEnumerable terms, Variable subtree, bool aPositiveVariable = false) where TPrimitive : TreeAnalyzer.IPrimitive, new() where T : notnull { var monomialsByPower = new Dictionary(); // here we fill the dictionary with information about monomials' coefficiants foreach (var child in terms) - if (ParseMonomial(subtree, child) is var (free, q)) + if (ParseMonomial(subtree, child, aPositiveVariable) is var (free, q)) monomialsByPower[q.Value] = monomialsByPower.TryGetValue(q.Value, out var power) ? power + free : free; else return null; @@ -294,6 +298,11 @@ internal static (Dictionary poly, Entity rem) GatherMonomialInformati } + /// Whether a factor is a number greater than zero, which a power of a product + /// that is not whole may take out of it. + private static bool IsAPositiveNumber(Entity factor) => + !factor.Vars.Any() && factor.Evaled is Real { IsPositive: true }; + /// Whether a monomial's power, integer, rational or decimal, is one. private static bool IsOne(T power) => power switch { @@ -303,7 +312,25 @@ internal static (Dictionary poly, Entity rem) GatherMonomialInformati _ => false, }; - internal static (Entity FreeMono, TPrimitive Power)? ParseMonomial(Variable aVar, Entity expr) + /// + /// A term read as a coefficient free of times a power of it, or + /// where it is not one. + /// + /// + /// A power of a product that is not whole is the product of the powers only where the + /// factor taken out is positive: (b x)^(5/2) is b^(5/2) x^(5/2) for + /// b > 0 or x > 0, and its negative at b < 0, x < 0, + /// where b x is positive and each of the two roots is imaginary. Read so, the long + /// division made b^(5/2) x^(5/2) + b^(5/2) x^(3/2) of (b x)^(5/2) (x + 1)/x, and + /// the integrator answered Rubi's (b sec(x))^(5/2) tan(x)^3 wrongly wherever + /// b sec(x) > 0 with b < 0; a negative number likewise, and + /// sqrt(-x) (x + 1)/x was simplified to i (sqrt(x) + x^(-1/2)), its negative + /// for every x < 0. A positive number comes out, and any factor where + /// says the variable is positive, as the limits at + /// +oo read it. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static (Entity FreeMono, TPrimitive Power)? ParseMonomial(Variable aVar, Entity expr, bool aPositiveVariable = false) where TPrimitive : TreeAnalyzer.IPrimitive, new() { var power = new TPrimitive(); @@ -332,10 +359,12 @@ internal static (Entity FreeMono, TPrimitive Power)? ParseMonomial(aVar, @base) is var (tmpFree, q)) + else if (ParseMonomial(aVar, @base, aPositiveVariable) is var (tmpFree, q)) { if (value is not Integer && !IsOne(q.Value)) return null; + if (value is not Integer && !aPositiveVariable && !IsAPositiveNumber(tmpFree)) + return null; freeMono *= MathS.Pow(tmpFree, value); power.AddMp(q.Value, value); } diff --git a/Sources/Tests/UnitTests/Calculus/AHalfOddPowerOfANegativeMultipleIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AHalfOddPowerOfANegativeMultipleIntegralTest.cs new file mode 100644 index 000000000..69481100d --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AHalfOddPowerOfANegativeMultipleIntegralTest.cs @@ -0,0 +1,117 @@ +// +// 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 +{ + /// + /// A power of a product that is not whole is not the product of the powers: (b u)^(5/2) + /// is -b^(5/2) u^(5/2) at b < 0, u < 0. The long division read it so, + /// and the integrator answered (b u)^(5/2) (u^2 - 1)/u, and through it Rubi's 4.3.0 + /// (b sec(x))^m tan(x)^3 at a half-odd m, wrongly wherever b sec(x) > 0 with + /// b < 0. A negative number did the same: sqrt(-x) (x + 1)/x was integrated, and + /// simplified, as i (sqrt(x) + x^(-1/2)), the negative of it for every x < 0. + /// Checked at both signs of b and of the variable. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + [Trait("Area", "Calculus")] + public sealed class AHalfOddPowerOfANegativeMultipleIntegralTest + { + private static readonly double[] AlgebraicPoints = { -2.1, -1.3, -0.6, 0.5, 1.2, 2.3 }; + // Both signs of the secant and of the cosecant. + private static readonly double[] TrigonometricPoints = { 0.4, 0.9, 2.0, 2.5, 3.0, 3.5, 4.0, 4.4 }; + + [Theory] + [InlineData("(b*x)^(5/2)*(x^2 - 1)/x", -0.4)] + [InlineData("(b*x)^(5/2)*(x^2 - 1)/x", 0.4)] + [InlineData("(b*x)^(5/2)*(x + 1)/x", -0.4)] + [InlineData("(b*x)^(5/2)*(x + 1)/x", 0.4)] + [InlineData("(b*x)^(3/2)*(x^2 + 1)/x", -0.4)] + [InlineData("(b*x)^(5/2)*x", -0.4)] + [InlineData("sqrt(b*x)*(1 + x)", -0.4)] + [InlineData("sqrt(-x)*(x + 1)/x", 0.0)] + [InlineData("(-x)^(3/2)*(x^2 - 1)/x", 0.0)] + [InlineData("(-3*x)^(5/2)*(x + 1)/x", 0.0)] + public void APowerOfAMultipleOfTheVariable(string integrand, double b) + => DifferentiatesBack(integrand, b, AlgebraicPoints); + + [Theory] + [InlineData("(b*sec(x))^(5/2)*tan(x)^3", -0.4)] + [InlineData("(b*sec(x))^(5/2)*tan(x)^3", 0.4)] + [InlineData("(b*sec(x))^(3/2)*tan(x)^3", -0.4)] + [InlineData("(b*sec(x))^(1/2)*tan(x)^5", -0.4)] + [InlineData("(b*sec(x))^(-1/2)*tan(x)^3", -0.4)] + [InlineData("cot(x)^3*(b*csc(x))^(3/2)", -0.4)] + [InlineData("cot(x)^3*(b*csc(x))^(3/2)", 0.4)] + [InlineData("cot(x)^5*(b*csc(x))^(1/2)", -0.4)] + public void APowerOfAMultipleOfTheSecantOrTheCosecant(string integrand, double b) + => DifferentiatesBack(integrand, b, TrigonometricPoints); + + [Theory] + [InlineData("sqrt(-x)*(x + 1)/x")] + [InlineData("(b*x)^(5/2)*(x + 1)/x")] + [InlineData("(b*x)^(3/2)*(x^2 + 1)/x")] + public void TheSimplificationKeepsTheValue(string written) + { + var expression = written.ToEntity(); + var simplified = expression.Simplify(); + foreach (var at in AlgebraicPoints) + { + Entity Pinned(Entity e) => e.Substitute("b", -0.4).Substitute("x", at); + var want = Pinned(expression).EvalNumerical(); + var got = Pinned(simplified).EvalNumerical(); + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"{written} simplified to {simplified} is {got} at x = {at}, b = -0.4, where it is {want}"); + } + } + + [Fact] + public void TheLongDivisionKeepsTheMultipleUnderThePower() + { + var u = MathS.Var("u"); + var dividend = "(b*u)^(5/2)*(u + 1)".ToEntity(); + var divisor = (Entity)u; + if (Functions.TreeAnalyzer.PolynomialLongDivision(dividend, divisor, genericCase: true, inTermsOf: u) is not var (quotient, remainder)) + return; + foreach (var (b, at) in new[] { (-0.4, -1.3), (-0.4, 0.7), (0.4, -1.3), (0.4, 0.7) }) + { + Entity Pinned(Entity e) => e.Substitute("b", b).Substitute(u, at); + var got = Pinned(quotient + remainder).EvalNumerical(); + var want = Pinned(dividend / divisor).EvalNumerical(); + Assert.True(((got - want).Abs() as Entity.Number.Real)!.EDecimal.ToDouble() < 1e-9 * Math.Max(1, (want.Abs() as Entity.Number.Real)!.EDecimal.ToDouble()), + $"the quotient and remainder are {got} at b = {b}, u = {at}, where the division is {want}"); + } + } + + private static void DifferentiatesBack(string integrand, double b, double[] points) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Entity Pinned(Entity e) => e.Substitute("b", b); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in points) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of {text}, the antiderivative of {integrand}, is {got} at x = {at}, b = {b}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +}