diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 5456a356d..06fff7bea 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1175,6 +1175,29 @@ power of the secant or the cosecant under a root goes on now too. The cosecant's | `"sqrt(1 + csc(x)^2)".ToEntity().Integrate("x")` | `integral(...)` | logarithms and an arctangent in `sqrt(2 tan(x)^2 + 1)`, times `sgn(tan(x))` | | `"1/sqrt(-1 + csc(x)^2)".ToEntity().Integrate("x")` | `integral(...)` | `sgn(tan(x)) ln(1 + tan(x)^2)/2` | +### The rational canonical form keeps the points where the expression has no value + +**Wrong answers, silent.** `CanonicalizeAsRationalFunction` and `Transformation.RationalCanonicalization` kept the condition a +cancelled factor owes, `x/x` being `1 provided not x = 0`, and dropped the one that dividing by a +quotient owes: `1/(1/x)` was `x`, which is `0` at zero, where `1/(1/x)` has no value. The +denominator a quotient turns over is now kept as a condition too, and so is the denominator under a +numerator that vanishes: `0/x` is `0 provided not x = 0`. The condition is written canonically, +because it is part of a form whose point is that equal trees mean equal functions: the square-free +part of what was excluded, without what the reduced denominator still excludes. So `x^2/x^2` and +`x/x` meet, and a cancelled factor that the denominator still has is not repeated. The form now +gathers through the first step of `AsSingleFraction`, so the two agree on what a single quotient is +and on where it has a value. [#1618](https://github.com/asc-community/AngouriMath/issues/1618). +Both columns measured on a build, `v2.5.0` against this change. + +| `….ToEntity().CanonicalizeAsRationalFunction()` of | Was (2.5.0) | Is | +|---|---|---| +| `"1/(1/x)"` | `x` | `x provided not x = 0` | +| `"1/(1 + 1/x)"` | `x / (x + 1)` | `x / (x + 1) provided not x = 0` | +| `"(a/b)/(c/d)"` | `a * d / (b * c)` | `a * d / (b * c) provided not d = 0` | +| `"0/x"` | `0` | `0 provided not x = 0` | +| `"x^2/x^2"` | `1 provided not x ^ 2 = 0` | `1 provided not x = 0` | +| `"(x^2 - 1)/(x^2 + 2x + 1) + 1/(x + 1)"` | `x / (x + 1) provided not x ^ 2 + 2 * x + 1 = 0` | `x / (x + 1)` | + ### A partial-fraction coefficient with symbols in it is in lowest terms, its rational content included **Improvement, not silent.** The decomposition over written factors with symbols among their diff --git a/Sources/AngouriMath/Core/Transformations/Transformation.Catalogue.cs b/Sources/AngouriMath/Core/Transformations/Transformation.Catalogue.cs index e3fbc7626..eac925e4c 100644 --- a/Sources/AngouriMath/Core/Transformations/Transformation.Catalogue.cs +++ b/Sources/AngouriMath/Core/Transformations/Transformation.Catalogue.cs @@ -352,16 +352,18 @@ private static class PolynomialFactorizationHolder /// expressions must not quietly hand back a normalisation that merely resembles one. /// /// - /// The expression is gathered into a single quotient — which is the part nothing else - /// in the library does, and without which 1/x + 1/y and (x + y)/(x*y) - /// could never meet — then reduced by the multivariate greatest common divisor and - /// scaled so the denominator is monic in the lexicographic monomial order. + /// The expression is written as a single fraction, as + /// writes it -- without which 1/x + 1/y and (x + y)/(x*y) could never + /// meet -- then reduced by the multivariate greatest common divisor and scaled so the + /// denominator is monic in the lexicographic monomial order. /// /// /// Cancelling carries its condition. x/x is not 1, so where a /// factor of positive degree comes out the answer says the factor is nonzero, as the - /// rest of the library already does. Gathering over a common denominator widens + /// rest of the library already does; and so does turning a quotient over, since + /// 1/(1/x) is not x either. Gathering over a common denominator widens /// nothing by itself: a sum is defined exactly where its terms are. + /// #1618 /// /// /// Nothing runs this by default. See diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/RationalFunction.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/RationalFunction.cs index fdc6cae63..1abc165db 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/RationalFunction.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/RationalFunction.cs @@ -30,8 +30,8 @@ namespace AngouriMath.Functions /// #934. /// /// - /// Four steps, and only the first is new. The expression is gathered into a single - /// quotient — nothing else in the library does that, and without it 1/x + 1/y and + /// Four steps. The expression is written as a single fraction, by the step + /// takes -- without it 1/x + 1/y and /// (x + y)/(x*y) could never meet. Then the numerator and denominator are divided /// by their multivariate greatest common divisor, which /// computes and verifies. Then both are scaled so that the @@ -40,12 +40,14 @@ namespace AngouriMath.Functions /// /// /// The domain is preserved rather than assumed away. Cancelling a common factor - /// widens the domain — x/x is not 1 — so where a factor of positive degree - /// is cancelled the answer carries the condition that it is nonzero, which is what the - /// library already does elsewhere and what keeps "equal trees means equal expressions" - /// true rather than nearly true. Gathering over a common denominator does not widen - /// anything: a sum is defined exactly where its terms are, and the product of the - /// denominators vanishes exactly where one of them does. + /// widens the domain — x/x is not 1 — and so does dividing by a quotient, + /// which moves its denominator into the numerator: 1/(1/x) is not x. So the + /// answer carries the condition that neither vanishes, which is what the library already + /// does elsewhere and what keeps "equal trees means equal expressions" true rather than + /// nearly true; the condition is written canonically too, square-free and without what + /// the denominator already excludes. A sum is defined exactly where its terms are, and the + /// common denominator vanishes exactly where one of theirs does, so gathering one widens + /// nothing. #1618 /// /// internal static class RationalFunction @@ -57,9 +59,10 @@ internal static class RationalFunction private const int MaxComplexity = 256; /// - /// Raising to a power is where a gathered quotient explodes, so the exponent is - /// bounded before is asked; that has its - /// own bound on the number of terms, which catches the rest. + /// Raising a quotient to a power is where a gathered quotient explodes, so a power of + /// anything but a polynomial is bounded before anything is multiplied out; a polynomial's + /// power has 's own bound on the number of + /// terms. /// private const int MaxExponent = 32; @@ -84,23 +87,24 @@ internal static bool TryCanonicalize(Entity expr, [NotNullWhen(true)] out Entity indices[variables[i]] = i; var variableCount = variables.Length; - if (!TryGather(expr, indices, variableCount, out var numerator, out var denominator)) + if (!TryGather(expr, indices, variableCount, out var numerator, out var denominator, out var excluded)) return false; - // A vanishing denominator is not a rational function, and a vanishing numerator - // is zero however it was written. + // A vanishing denominator is not a rational function. A vanishing numerator is zero + // wherever the expression has a value: everywhere its denominator does not vanish. if (denominator.IsZero) return false; + var order = new int[variableCount]; + for (var i = 0; i < order.Length; i++) + order[i] = i; if (numerator.IsZero) { - canonical = Integer.Create(0); + if (excluded.Multiply(denominator) is not { } wholeDenominator + || Conditioned(Integer.Create(0), wholeDenominator, MultivariatePolynomial.One(variableCount), order, variables) is not { } zero) + return false; + canonical = zero; return true; } - var order = new int[variableCount]; - for (var i = 0; i < order.Length; i++) - order[i] = i; - - var cancelled = MultivariatePolynomial.One(variableCount); if (variableCount > 0 && PolynomialGcd.Gcd(numerator, denominator, order, 0) is { } divisor && !divisor.IsConstant) @@ -115,9 +119,11 @@ internal static bool TryCanonicalize(Entity expr, [NotNullWhen(true)] out Entity || reducedBottom.Multiply(divisor) is not { } checkedBottom || !checkedBottom.SameAs(denominator)) return false; + if (excluded.Multiply(divisor) is not { } withTheCancelled) + return false; numerator = reducedTop; denominator = reducedBottom; - cancelled = divisor; + excluded = withTheCancelled; } // Scaled so the denominator leads with one, under the same lexicographic monomial @@ -138,125 +144,86 @@ internal static bool TryCanonicalize(Entity expr, [NotNullWhen(true)] out Entity var quotient = denominator.IsConstant ? numerator.ToEntity(variables) : numerator.ToEntity(variables) / denominator.ToEntity(variables); - - canonical = cancelled.IsConstant - ? quotient - : new Providedf(quotient, !cancelled.ToEntity(variables).EqualTo(0)); - return true; + canonical = Conditioned(quotient, excluded, denominator, order, variables); + return canonical is not null; } /// - /// as a single quotient of polynomials, gathering a sum of - /// quotients over a common denominator. The denominator is never zero and never - /// simplified away; it is 1 for a polynomial. + /// with the condition that it is not at a zero of + /// , where the expression it came from had no value, except + /// where vanishing already says so; + /// where a step of the arithmetic declined. /// /// - /// The common denominator is the product rather than the least common multiple. Both - /// are correct and the product is cheaper to build; what it costs is a larger - /// intermediate, which the greatest common divisor then removes — so the answer is the - /// same and only the work in between differs. + /// The condition is part of the form, so it is canonical too: two expressions with the same + /// value and the same points where they have none have to meet. So the polynomial written + /// is the square-free part of , which vanishes exactly where it + /// does -- x^2/x^2 and x/x are both 1 provided not x = 0 -- without the + /// factors it shares with the denominator, and with whole coprime coefficients and a + /// positive leading one. The square-free part is E / gcd(E, dE/dx_1, …, dE/dx_n): + /// over the rationals a factor to the power k divides each derivative to the power + /// k - 1, and one of them no further. /// - private static bool TryGather( - Entity expr, IReadOnlyDictionary indices, int variableCount, - [NotNullWhen(true)] out MultivariatePolynomial? numerator, - [NotNullWhen(true)] out MultivariatePolynomial? denominator) + private static Entity? Conditioned(Entity quotient, MultivariatePolynomial excluded, MultivariatePolynomial denominator, + IReadOnlyList order, IReadOnlyList variables) { - numerator = denominator = null; - - // A polynomial is its own numerator, and this is the common case, so it is tried - // before the expression is taken apart. - if (MultivariatePolynomial.TryParse(expr, indices) is { } whole) - { - numerator = whole; - denominator = MultivariatePolynomial.One(variableCount); - return true; - } - - switch (expr) - { - case Sumf(var left, var right): - return TryCombine(left, right, subtract: false, indices, variableCount, - out numerator, out denominator); - - case Minusf(var left, var right): - return TryCombine(left, right, subtract: true, indices, variableCount, - out numerator, out denominator); - - case Mulf(var left, var right): - { - if (!TryGather(left, indices, variableCount, out var leftTop, out var leftBottom) - || !TryGather(right, indices, variableCount, out var rightTop, out var rightBottom)) - return false; - if (leftTop.Multiply(rightTop) is not { } top - || leftBottom.Multiply(rightBottom) is not { } bottom) - return false; - numerator = top; - denominator = bottom; - return true; - } - - case Divf(var left, var right): + if (excluded.IsConstant) + return excluded.IsZero ? null : quotient; + var repeated = excluded; + for (var i = 0; i < excluded.VariableCount; i++) + if (excluded.DegreeIn(i) > 0) { - if (!TryGather(left, indices, variableCount, out var leftTop, out var leftBottom) - || !TryGather(right, indices, variableCount, out var rightTop, out var rightBottom)) - return false; - // Dividing by a quotient that is identically zero is not a rational - // function, and inverting it here would quietly produce one. - if (rightTop.IsZero) - return false; - if (leftTop.Multiply(rightBottom) is not { } top - || leftBottom.Multiply(rightTop) is not { } bottom) - return false; - numerator = top; - denominator = bottom; - return true; + if (PolynomialGcd.Gcd(repeated, excluded.DerivativeIn(i), order, 0) is not { } common) + return null; + repeated = common; } - - case Powf(var @base, Integer power): - { - var exponent = power.EInteger; - if (exponent.Abs().CompareTo(EInteger.FromInt32(MaxExponent)) > 0) - return false; - if (!TryGather(@base, indices, variableCount, out var baseTop, out var baseBottom)) - return false; - var magnitude = exponent.Abs().ToInt32Checked(); - var negative = exponent.Sign < 0; - if (negative && baseTop.IsZero) - return false; - var top = negative ? baseBottom : baseTop; - var bottom = negative ? baseTop : baseBottom; - if (top.Power(magnitude) is not { } raisedTop - || bottom.Power(magnitude) is not { } raisedBottom) - return false; - numerator = raisedTop; - denominator = raisedBottom; - return true; - } - - default: - return false; - } + if (excluded.DivideExact(repeated) is not { } squareFree + || PolynomialGcd.Gcd(squareFree, denominator, order, 0) is not { } shared + || squareFree.DivideExact(shared) is not { } left) + return null; + return left.IsConstant + ? quotient + : new Providedf(quotient, !left.Normalized().ToEntity(variables).EqualTo(0)); } /// - /// A sum or a difference of two quotients, over the product of their denominators. + /// as a single quotient of polynomials, and the polynomial whose + /// zeros are the other points where has no value: the product of + /// the denominators that dividing by a quotient turned over. The denominator is never + /// zero and never simplified away; it is 1 for a polynomial. /// - private static bool TryCombine( - Entity left, Entity right, bool subtract, - IReadOnlyDictionary indices, int variableCount, + /// + /// Gathered by , the step + /// takes, so that the two agree on what a single + /// quotient of an expression is and on where it has a value; this form is that quotient + /// reduced and normalised. + /// + private static bool TryGather( + Entity expr, IReadOnlyDictionary indices, int variableCount, [NotNullWhen(true)] out MultivariatePolynomial? numerator, - [NotNullWhen(true)] out MultivariatePolynomial? denominator) + [NotNullWhen(true)] out MultivariatePolynomial? denominator, + [NotNullWhen(true)] out MultivariatePolynomial? excluded) { - numerator = denominator = null; - if (!TryGather(left, indices, variableCount, out var leftTop, out var leftBottom) - || !TryGather(right, indices, variableCount, out var rightTop, out var rightBottom)) - return false; - if (leftTop.Multiply(rightBottom) is not { } crossLeft - || rightTop.Multiply(leftBottom) is not { } crossRight - || leftBottom.Multiply(rightBottom) is not { } bottom) + numerator = denominator = excluded = null; + foreach (var node in expr.Nodes) + if (node is Powf(var @base, Integer power) + && power.EInteger.Abs().CompareTo(EInteger.FromInt32(MaxExponent)) > 0 + && MultivariatePolynomial.TryParse(@base, indices) is null) + return false; + var carried = new List(); + var (top, bottom) = SingleQuotient.OverLeastCommonDenominator(expr, carried); + if (MultivariatePolynomial.TryParse(top, indices) is not { } parsedTop + || MultivariatePolynomial.TryParse(bottom, indices) is not { } parsedBottom) return false; - numerator = subtract ? crossLeft.Subtract(crossRight) : crossLeft.Add(crossRight); - denominator = bottom; + var product = MultivariatePolynomial.One(variableCount); + foreach (var turnedOver in carried) + { + if (MultivariatePolynomial.TryParse(turnedOver, indices) is not { } parsed + || product.Multiply(parsed) is not { } next) + return false; + product = next; + } + (numerator, denominator, excluded) = (parsedTop, parsedBottom, product); return true; } } diff --git a/Sources/Tests/UnitTests/Core/Transformations/RationalCanonicalFormTest.cs b/Sources/Tests/UnitTests/Core/Transformations/RationalCanonicalFormTest.cs index 79e40156f..bb6c38aac 100644 --- a/Sources/Tests/UnitTests/Core/Transformations/RationalCanonicalFormTest.cs +++ b/Sources/Tests/UnitTests/Core/Transformations/RationalCanonicalFormTest.cs @@ -50,7 +50,9 @@ private static Entity Required(string expression) [InlineData("2 * x / (4 * y)", "x / (2 * y)")] [InlineData("(x + 1) / (x + 2)", "(2 * x + 2) / (2 * x + 4)")] [InlineData("x / y + 1", "(x + y) / y")] - [InlineData("1 / (1 / x)", "x")] + [InlineData("1 / (1 / x)", "x ^ 2 / x")] + [InlineData("x ^ 2 / x ^ 2", "x / x")] + [InlineData("(a / b) / (c / d)", "(a * d ^ 2) / (b * c * d)")] [InlineData("x ^ (-2)", "1 / x ^ 2")] [InlineData("(a + b) / (a * b)", "1/b + 1/a")] public void OneFunctionHasOneForm(string left, string right) @@ -65,9 +67,52 @@ public void OneFunctionHasOneForm(string left, string right) [InlineData("x / y", "y / x")] [InlineData("(x + 1) / (x + 2)", "(x + 2) / (x + 1)")] [InlineData("x / (2 * y)", "x / (3 * y)")] + // Undefined at zero and zero there: not one function, however the value agrees elsewhere. + // https://github.com/asc-community/AngouriMath/issues/1618 + [InlineData("1 / (1 / x)", "x")] + [InlineData("1 / (1 + 1 / x)", "x / (x + 1)")] + [InlineData("0 / x", "0")] public void DifferentFunctionsDoNotCollide(string left, string right) => Assert.NotEqual(Required(left), Required(right)); + /// + /// Dividing by a quotient moves its denominator into the numerator, where it no longer + /// stops the form having a value; the expression has none there, so the form says the + /// denominator is nonzero. The same for a quotient raised to a negative power, and for a + /// numerator that vanishes, which is zero only where its denominator does not. + /// https://github.com/asc-community/AngouriMath/issues/1618 + /// + [Theory] + [InlineData("1 / (1 / x)", "x", "0")] + [InlineData("1 / (1 + 1 / x)", "x", "0")] + [InlineData("(a / b) / (c / d)", "d", "0")] + [InlineData("(1 / x) ^ (-1)", "x", "0")] + [InlineData("(x / y) ^ (-2)", "y", "0")] + [InlineData("x / (y / x)", "x", "0")] + [InlineData("0 / x", "x", "0")] + public void WhereTheExpressionHasNoValueNeitherHasTheForm(string expression, string variable, string at) + { + Entity original = expression; + var form = Required(expression); + foreach (var name in original.Vars) + { + Entity value = name.Name == variable ? at : "1"; + original = original.Substitute(name, value); + form = form.Substitute(name, value); + } + Assert.Equal(MathS.NaN, original.Evaled); + Assert.Equal(MathS.NaN, form.Evaled); + } + + /// + /// The condition is part of the form, so it is canonical too, and it says only what the + /// denominator does not: (x^2 - 1)/(x^2 + 2x + 1) + 1/(x + 1) cancels + /// (x + 1)^2, and x/(x + 1) is already undefined at -1. + /// + [Fact] + public void TheDenominatorsExclusionIsNotRepeated() + => Assert.Equal(Required("x / (x + 1)"), Required("(x ^ 2 - 1) / (x ^ 2 + 2 * x + 1) + 1 / (x + 1)")); + /// /// The case where "the same function" is a trap. `(x^2 - 1)/(x + 1)` is undefined at /// `x = -1` and `x - 1` is not, so they are *not* the same function and the form must