diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 5589f7269..1a872dc34 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1659,6 +1659,27 @@ measured on a build, `v2.5.0` against this change. | `"sum(x, x in {a, b})".ToEntity().Simplify()` | the same exception | `sum(x, x in { a, b })`, since `a` and `b` may be equal | | `"sum(k, k, 1, 10)".ToEntity().Simplify()` | `55` | the same | +### A condition on a name a binder binds stays in the binder + +**Fix.** `DomainCondition` conjoined every child's condition, a binder's body included, so a +summation's, a product's or a definite integral's condition named its index as though it were free, +and a sum over a set's did the same: the product rule leaves `0 * sum(...)`, which keeps it, and the +derivative of `2 sum(ln(x - w), w in { w : w^3 + w + 1 = 0 })` was `… provided not x - w = 0`, which +`EvalNumerical` could not decide. The condition is required at every value the name ranges over +now: a conjunction over a few numbers, a resultant over the roots of a polynomial, and `forall` +over the range otherwise. It is not left out, which would give `0 * sum(1/(k - x), k, 1, n)` a value +at `x = 1` ([#1632](https://github.com/asc-community/AngouriMath/issues/1632)). A definite integral's +condition is required between its limits, in either order. That is sufficient for a value and not +necessary: `ln(t)` over `[0; 1]` converges, to -1, and its condition reads it as undefined, since the +integrand has none at 0. Nothing is given a value it does not have that way. The integral itself, +where it is worked out, is not affected. + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sum(1/(k - x), k, 1, 3)".ToEntity().DomainCondition` | `not k - x = 0` | `not 1 - x = 0 and not 2 - x = 0 and not 3 - x = 0` | +| `"sum(1/(k - x), k, 1, n)".ToEntity().DomainCondition` | `not k - x = 0` | `forall k in { k : k in ZZ and 1 <= k and k <= n } : not k - x = 0` | +| `"integral(1/(t - x), t, 0, 1)".ToEntity().DomainCondition` | `not t - x = 0` | `forall t in [0; 1] : not t - x = 0` | + ### `(a + b asech(c x))/(d + e x)^2` is integrated, and a root written apart with `|x|` no longer needs a parity Rubi's 7.5.1 with a symbolic linear below the bar ran for ten minutes without an answer. By diff --git a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Definition.cs b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Definition.cs index 5dbfde217..7d111ef5d 100644 --- a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Definition.cs +++ b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Definition.cs @@ -7,6 +7,7 @@ using AngouriMath.Core.Transformations; using HonkSharp.Laziness; +using PeterO.Numbers; namespace AngouriMath { @@ -41,14 +42,252 @@ partial record Entity /// - Singularities and poles (points where a function is undefined) /// - Piecewise continuity (tracking where discontinuities occur) /// - public Entity DomainCondition => domainCondition.GetValue(static @this => @this.DirectChildren.Aggregate(@this.IntrinsicCondition, (accum, curr) => + public Entity DomainCondition => domainCondition.GetValue(static @this => ScopedToItsBinder(@this, @this.DirectChildren.Aggregate(@this.IntrinsicCondition, (accum, curr) => (accum, curr.DomainCondition) switch { (Boolean(true), Boolean(true)) => Boolean.True, (var l, Boolean(true)) => l, (Boolean(true), var r) => r, (var l, var r) => l & r, - }), this).InnerSimplified; + })), this).InnerSimplified; private LazyPropertyA domainCondition; + + /// + /// with each conjunct that mentions a name + /// binds read the way the binder reads it: required at every + /// value the name ranges over. What a sum's body says about its index is a condition + /// inside the sum, and outside the sum that name is free, or another expression's. + /// + /// + /// + /// 0 * f keeps the condition f is defined under, so the product rule carried + /// not x - w = 0 out of sum(ln(x - w), w in { w : w^3 + w + 1 = 0 }) with + /// w free in it, and the derivative could not be evaluated. Leaving the conjunct out + /// would have made 0 * sum(1/(k - x), k, 1, n) a plain 0 at x = 1, where the + /// sum has no value. + /// + /// + /// Over the roots of a polynomial p with rational coefficients, "q(w) != 0 at + /// every root" is Res_w(p, q) != 0, which mentions the other names alone -- here + /// not x^3 + x + 1 = 0, where the logarithms are singular -- since p's leading + /// coefficient is a number that is not zero. A finite range of numbers is the conjunction + /// over it. Anything else is kept as forall over the range: the set, the whole + /// numbers from the lower bound to the upper, the interval between a definite integral's + /// limits. A set builder admits the members its predicate is defined at, so its own + /// condition on its name is part of whom it admits rather than of where the set is defined. + /// A limit does not need its body defined throughout anything, and is left as it was. + /// https://github.com/asc-community/AngouriMath/issues/1632 + /// + /// + /// For a definite integral this is sufficient for a value and not necessary: an integral + /// whose integrand is undefined only at a point it still converges past, ln(t) over + /// [0; 1], is read as undefined here. That is the direction in which nothing is given + /// a value it does not have, and a pointwise condition cannot tell a convergent improper + /// integral from a divergent one. The integral itself, where it is worked out, is not + /// affected: the derivative of 2 integral(x ln(t), t, 0, 1) is -2. + /// + /// + private static Entity ScopedToItsBinder(Entity binder, Entity condition) + { + if (condition is Boolean) + return condition; + (Entity Name, Entity? Range)? scope = binder switch + { + SumOverSetf(_, var name, var over) => (name, over), + Maximumf(_, var name, var over) => (name, over), + Minimumf(_, var name, var over) => (name, over), + Argmaxf(_, var name, var over) => (name, over), + Argminf(_, var name, var over) => (name, over), + Quantifier quantifier => (quantifier.Var, quantifier.Over), + Summationf(_, var index, var from, var to) => (index, WholeNumbersBetween(index, from, to)), + Productf(_, var index, var from, var to) => (index, WholeNumbersBetween(index, from, to)), + Integralf { Range: { } limits } integral => (integral.Var, Between(integral.Var, limits.from, limits.to)), + Set.ConditionalSet(var name, _) => (name, null), + _ => null + }; + if (scope is not var (bound, range)) + return condition; + var names = bound.VarsAndConsts; + Entity kept = Boolean.True; + foreach (var conjunct in Conjuncts(condition)) + { + var scoped = !conjunct.FreeVariables.Any(names.Contains) ? conjunct + : range is null ? null + : OverTheRange(conjunct, bound, range); + if (scoped is not null) + kept = kept is Boolean(true) ? scoped : kept & scoped; + } + return kept; + } + + /// + /// The segment between and in either order: + /// an integral from 1 to 0 ranges over [0; 1], where [1; 0] is empty and would make + /// any condition over it true. Two numbers are put in order; otherwise it is written as the + /// set of values between the two, since [a; b] \/ [b; a] simplifies to { b }. + /// + private static Entity Between(Entity name, Entity from, Entity to) + { + if (from is Number.Real low && to is Number.Real high) + return low <= high ? MathS.Interval(low, high) : MathS.Interval(high, low); + return new Set.ConditionalSet(name, (from <= name) & (name <= to) | (to <= name) & (name <= from)); + } + + /// The whole numbers from to , listed where there are a few of them. + private static Entity WholeNumbersBetween(Entity index, Entity from, Entity to) + { + if (from is Number.Integer low && to is Number.Integer high + && high.EInteger.Subtract(low.EInteger).CompareTo(16) < 0) + { + var members = new List(); + for (var k = low.EInteger; k.CompareTo(high.EInteger) <= 0; k = k.Add(1)) + members.Add(Number.Integer.Create(k)); + return new Set.FiniteSet(members); + } + return new Set.ConditionalSet(index, index.In(MathS.Sets.Z) & (from <= index) & (index <= to)); + } + + /// required at every value takes in . + private static Entity OverTheRange(Entity conjunct, Entity bound, Entity range) + { + if (bound is Variable name) + { + if (range is Set.ConditionalSet { Var: Variable root, Predicate: Equalsf(var left, var right) } + && AtEveryRoot(conjunct, name, (left - right).Substitute(root, name)) is { } exact) + return exact; + if (range is Set.FiniteSet finite) + { + Entity each = Boolean.True; + foreach (var member in finite.Elements) + { + var at = conjunct.Substitute(name, member); + each = each is Boolean(true) ? at : each & at; + } + return each; + } + } + return new Forallf(bound, range, conjunct); + } + + /// + /// not Res_w(p, q) = 0 for a conjunct not g = h whose q = g - h is a + /// polynomial in : true exactly where q is not zero at any + /// root of , whose coefficients must be rational. Null for any + /// other conjunct or polynomial. + /// + private static Entity? AtEveryRoot(Entity conjunct, Variable name, Entity polynomial) + { + if (conjunct is not Notf(Equalsf(var g, var h)) + || Functions.SumOverSet.SquareFreeParts(polynomial, name) is not { } parts + || !Functions.TreeAnalyzer.TryGetPolynomial(g - h, name, out var terms)) + return null; + // The roots once each: the product of the square-free parts. + var p = Functions.IntegerPolynomial.Create(new[] { EInteger.One }); + foreach (var part in parts) + if (p.Multiply(part.Factor) is { } product) + p = product; + else + return null; + var n = p.Degree; + var m = 0; + foreach (var power in terms.Keys) + { + if (power.Sign < 0 || !power.CanFitInInt32()) + return null; + m = System.Math.Max(m, power.ToInt32Unchecked()); + } + if (m == 0 || n < 1) + return null; + Entity Q(int power) => terms.TryGetValue(EInteger.FromInt32(power), out var c) ? c : Number.Integer.Zero; + Entity P(int power) => Number.Integer.Create(p[power]); + Entity resultant; + if (m == 1) + { + // a^n p(-b/a), written without dividing by a: sum of p_k (-b)^k a^(n - k). The + // zeroth and first powers are written as themselves, since a power 0 of a base + // that might be zero is 1 only where it is not, and that condition is not this one. + var (a, b) = (Q(1), Q(0)); + static Entity Power(Entity @base, int exponent) + => exponent == 0 ? Number.Integer.One : exponent == 1 ? @base : MathS.Pow(@base, exponent); + resultant = Number.Integer.Zero; + for (var k = 0; k <= n; k++) + resultant += P(k) * Power(-b, k) * Power(a, n - k); + } + else + { + // Past a linear, only with rational coefficients, where the resultant is a number + // and only whether it is zero matters: the Sylvester matrix in whole numbers, its + // determinant by Bareiss's elimination. A symbolic q of higher degree is left to + // forall. The determinant of a matrix of entities would reach GenericTensor, whose + // operations native AOT cannot compile, from a property every node has. + var denominators = EInteger.One; + var rational = new ERational[m + 1]; + for (var power = 0; power <= m; power++) + { + if (Q(power).Evaled is not Number.Rational coefficient) + return null; + rational[power] = coefficient.ERational; + denominators = denominators.Divide(denominators.Gcd(rational[power].Denominator)).Multiply(rational[power].Denominator); + } + var size = n + m; + var sylvester = new EInteger[size, size]; + for (var row = 0; row < size; row++) + for (var column = 0; column < size; column++) + { + sylvester[row, column] = EInteger.Zero; + if (row < m && column - row >= 0 && column - row <= n) + sylvester[row, column] = p[n - (column - row)]; + else if (row >= m && column - (row - m) >= 0 && column - (row - m) <= m) + { + var q = rational[m - (column - (row - m))]; + sylvester[row, column] = q.Numerator.Multiply(denominators.Divide(q.Denominator)); + } + } + return WholeNumberDeterminant(sylvester).IsZero ? Boolean.False : Boolean.True; + } + return !resultant.Equalizes(Number.Integer.Zero); + } + + /// The determinant of a square matrix of whole numbers, by Bareiss's fraction-free elimination. + private static EInteger WholeNumberDeterminant(EInteger[,] matrix) + { + var size = matrix.GetLength(0); + var a = (EInteger[,])matrix.Clone(); + var negated = false; + var previous = EInteger.One; + for (var k = 0; k < size - 1; k++) + { + if (a[k, k].IsZero) + { + var swap = -1; + for (var i = k + 1; i < size && swap < 0; i++) + if (!a[i, k].IsZero) + swap = i; + if (swap < 0) + return EInteger.Zero; + for (var j = 0; j < size; j++) + (a[k, j], a[swap, j]) = (a[swap, j], a[k, j]); + negated = !negated; + } + for (var i = k + 1; i < size; i++) + for (var j = k + 1; j < size; j++) + a[i, j] = a[i, j].Multiply(a[k, k]).Subtract(a[i, k].Multiply(a[k, j])).Divide(previous); + previous = a[k, k]; + } + return negated ? a[size - 1, size - 1].Negate() : a[size - 1, size - 1]; + } + + private static IEnumerable Conjuncts(Entity condition) + { + if (condition is Andf(var left, var right)) + { + foreach (var conjunct in Conjuncts(left)) + yield return conjunct; + foreach (var conjunct in Conjuncts(right)) + yield return conjunct; + } + else + yield return condition; + } /// /// Returns the intrinsic condition under which this specific operation is defined, diff --git a/Sources/Tests/UnitTests/Core/BinderConditionTest.cs b/Sources/Tests/UnitTests/Core/BinderConditionTest.cs new file mode 100644 index 000000000..90b42daa7 --- /dev/null +++ b/Sources/Tests/UnitTests/Core/BinderConditionTest.cs @@ -0,0 +1,133 @@ +// +// 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.Linq; +using AngouriMath.Extensions; +using Xunit; +using static AngouriMath.Entity; + +namespace AngouriMath.Tests.Core +{ + /// + /// What a binder's body says about the name it binds stays inside the binder: it is required + /// at every value the name ranges over, and never left free outside, nor left out. + /// #1632 + /// + [Trait("Area", "Core")] + public sealed class BinderConditionTest + { + private const string LogarithmsOverRoots = "sum(ln(x - w), w in { w : w^3 + w + 1 = 0 })"; + + /// + /// The product rule leaves 0 * sum(...), which keeps the sum's condition: that was + /// not x - w = 0 with w free, and the derivative could not be evaluated. + /// + [Fact] + public void ADerivativeKeepsNoBoundName() + { + var derivative = ("2 * " + LogarithmsOverRoots).ToEntity().Differentiate("x"); + Assert.DoesNotContain(derivative.FreeVariables, v => v.Name == "w"); + Assert.True(derivative.Substitute("x", 0.3).EvalNumerical() is Number.Complex); + } + + /// + /// Over the roots of p, "x - w != 0 at every root" is the resultant + /// p(x) != 0. With p = w^2 - 1 that is false at the roots 1 and -1 and true + /// elsewhere, and it mentions x alone. At 0 too: the zeroth power in the resultant + /// was written as a power, which is 1 only where its base is not zero. + /// + [Theory] + [InlineData("1", false)] + [InlineData("-1", false)] + [InlineData("2", true)] + [InlineData("1/2", true)] + [InlineData("0", true)] + public void OverRootsTheConditionIsTheResultant(string at, bool defined) + { + var condition = "sum(ln(x - w), w in { w : w^2 - 1 = 0 })".ToEntity().DomainCondition; + Assert.Equal(new[] { "x" }, condition.Vars.Select(v => v.Name).ToArray()); + Assert.Equal(defined ? MathS.Boolean.True : MathS.Boolean.False, condition.Substitute("x", at.ToEntity()).Simplify()); + } + + /// + /// The root polynomial itself as the condition: every term is undefined, so the sum is + /// defined nowhere, and the resultant of p with itself is zero. + /// + [Fact] + public void ASumOfTermsUndefinedAtEveryRootIsDefinedNowhere() + => Assert.Equal(MathS.Boolean.False, + "sum(1/(w^3 + w + 1), w in { w : w^3 + w + 1 = 0 })".ToEntity().DomainCondition); + + /// + /// Not left out: 0 * sum(1/(k - x), k, 1, n) is undefined wherever the sum is, and a + /// plain 0 would give it a value at x = 1. + /// + [Fact] + public void AZeroTimesASumKeepsTheSumsCondition() + { + var product = "0 * sum(1/(k - x), k, 1, n)".ToEntity().InnerSimplified; + Assert.IsType(product); + Assert.DoesNotContain(product.FreeVariables, v => v.Name == "k"); + Assert.True(product.Substitute("x", 1).Substitute("n", 3).Simplify().IsNaN); + Assert.Equal(0, product.Substitute("x", "1/2".ToEntity()).Substitute("n", 3).Simplify()); + } + + /// + /// A definite integral's condition holds between its limits in either order: from 1 to 0 + /// ranges over [0; 1] as well, where [1; 0] is empty and would make any condition + /// over it true. + /// + [Theory] + [InlineData("integral(1/(t - x), t, 0, 1)", "1/2", false)] + [InlineData("integral(1/(t - x), t, 0, 1)", "2", true)] + [InlineData("integral(1/(t - x), t, 1, 0)", "1/2", false)] + [InlineData("integral(1/(t - x), t, 1, 0)", "2", true)] + public void AnIntegralsConditionHoldsBetweenItsLimits(string integral, string at, bool defined) + { + var condition = integral.ToEntity().DomainCondition; + Assert.DoesNotContain(condition.FreeVariables, v => v.Name == "t"); + Assert.Equal(defined ? MathS.Boolean.True : MathS.Boolean.False, condition.Substitute("x", at.ToEntity()).Simplify()); + } + + /// Between symbolic limits, whichever is larger: no name of the integral's is left free. + [Fact] + public void BetweenSymbolicLimitsNoBoundNameIsLeft() + => Assert.DoesNotContain("integral(1/(t - x), t, a, b)".ToEntity().DomainCondition.FreeVariables, v => v.Name == "t"); + + /// + /// Sufficient for a value and not necessary: ln(t) over [0; 1] converges, to + /// -1, and is read as undefined, since the integrand has no value at 0; over [1; 2] + /// it is defined throughout. That is the direction in which nothing is given a value it does + /// not have. + /// + [Theory] + [InlineData("integral(ln(t), t, 0, 1)", false)] + [InlineData("integral(ln(t), t, 1, 2)", true)] + public void AnIntegralUndefinedAtAPointOfItsRangeIsReadAsUndefined(string integral, bool defined) + => Assert.Equal(defined ? MathS.Boolean.True : MathS.Boolean.False, integral.ToEntity().DomainCondition.Simplify()); + + /// A few whole numbers are the conjunction over them, and none is true. + [Theory] + [InlineData("sum(1/(k - x), k, 1, 3)", "2", false)] + [InlineData("sum(1/(k - x), k, 1, 3)", "5", true)] + [InlineData("sum(1/(k - x), k, 3, 1)", "2", true)] + public void AFiniteRangeIsTheConjunctionOverIt(string sum, string at, bool defined) + { + var condition = sum.ToEntity().DomainCondition; + Assert.DoesNotContain(condition.Vars, v => v.Name == "k"); + Assert.Equal(defined ? MathS.Boolean.True : MathS.Boolean.False, condition.Substitute("x", at.ToEntity()).Simplify()); + } + + /// + /// A set builder admits the members its predicate is defined at: its condition on its own + /// name is whom it admits, not where the set is defined. + /// + [Fact] + public void ASetBuildersOwnConditionIsMembership() + => Assert.DoesNotContain("{ w : 1/w > a }".ToEntity().DomainCondition.Vars, v => v.Name == "w"); + } +}