From 7b7f5e1108117d79b8580290151b0680f8a08ca8 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 30 Sep 2026 09:37:15 +0000 Subject: [PATCH 1/2] A sum over a set: sum(f, x in S) sum took a range and nothing else, so a sum over a set had no spelling: sum(f, x in S) was a parse error asking for the range form's four arguments. SumOverSetf reads the second argument as `variable in set` and binds the variable over the term and the set, as max(f, t in S) does; it prints \sum_{x \in S} f, and MathS.Sum(f, x, S) builds it. The range form is untouched. A sum over the roots of a polynomial is a sum over its set of roots, sum(f(w), w in { w : p(w) = 0 }), with no index and no order. The terms are added up where the set is finite and its members are known: a listed set of numbers, distinct by value, and the roots of a polynomial where the solver writes all of them -- counted against the degree of the square-free part, so a root the solver missed leaves the sum as written rather than a term short. The solver is asked only where every irreducible factor has degree at most four or is a binomial. Evaluated to a number, the roots it cannot write are found to the working precision by Durand-Kerner, one square-free factor at a time, all of them or none; the polynomial is read off the tree, so that holds with downcasting off. A set whose members may coincide, and one that is not finite, are left as written. Termwise differentiation where the set does not depend on the variable, zero in the bound name; SymPy export as a sum of substitutions or a RootSum of the square-free part. Part of #1285. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 28 +++ Sources/.editorconfig | 6 + Sources/AngouriMath/Convenience/MathS.cs | 32 +++ Sources/AngouriMath/Core/Antlr/AngouriMath.g | 2 +- .../Core/Antlr/AngouriMathParser.cs | 2 +- Sources/AngouriMath/Core/Domains.Classes.cs | 8 + .../Entity.Continuous.Calculus.Classes.cs | 23 ++ .../Core/Entity/Entity.Definition.cs | 3 + .../Serialization/Entity.Serialization.cs | 1 + Sources/AngouriMath/Docs/Usage/Syntax.md | 15 +- .../Functions/Algebra/SumOverSet.cs | 220 ++++++++++++++++++ .../Functions/Continuous/Differentiation.cs | 16 ++ .../EquationSolver/InvertNode.Classes.cs | 7 + .../Solvers/NumericalSolving/DurandKerner.cs | 136 +++++++++++ .../Evaluation.Continuous.Calculus.Classes.cs | 19 ++ .../Output/Latex/Latex.Calculus.Classes.cs | 8 + .../ToString/ToString.Calculus.Classes.cs | 9 + .../ToSympy/ToSympy.Calculus.Classes.cs | 14 ++ Sources/AngouriMath/Functions/Substitute.cs | 14 ++ .../Functions/TreeAnalyzer/Sort.Classes.cs | 5 + .../UnitTests/Calculus/SumOverSetTest.cs | 134 +++++++++++ Sources/Tests/UnitTests/Common/PublicApi.txt | 21 ++ 22 files changed, 719 insertions(+), 4 deletions(-) create mode 100644 Sources/AngouriMath/Functions/Algebra/SumOverSet.cs create mode 100644 Sources/AngouriMath/Functions/Continuous/Solvers/NumericalSolving/DurandKerner.cs create mode 100644 Sources/Tests/UnitTests/Calculus/SumOverSetTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 03e1be289..c3d2e3b62 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1608,6 +1608,34 @@ the ends are about, `n >= 0`. With two names it is left as it was, `n >= m - 1` |---|---|---| | `"sum(2 k + 1, k, 0, n - 1)".Simplify()` | `piecewise((2 * n + (n - 1) ^ 2 - 1) provided (n - 1 >= -1), 0 provided True)` | `piecewise((n ^ 2) provided (n >= 0), 0 provided True)` | +### A sum over a set: `sum(f, x in S)` + +`sum` took a range and nothing else, so a sum over a set had no spelling: `sum(f, x in S)` was a +parse error asking for the four arguments of the range form. `SumOverSetf` reads the second +argument as `variable in set` and binds the variable over the term and the set, as `max(f, t in S)` +does: spelt `sum(f, x in S)`, printed `\sum_{x \in S} f` in LaTeX, and built by +`MathS.Sum(f, x, S)`. The range form `sum(f, k, from, to)` is untouched. A sum over the roots of a +polynomial is a sum over its set of roots, `sum(f(w), w in { w : p(w) = 0 })`, with no index and no +order. + +The terms are added up where the set is finite and its members are known: a listed set of numbers, +distinct by value, and the roots of a polynomial where the solver writes all of them. That is +checked by counting them against the degree of the square-free part, so a root the solver missed +leaves the sum as written rather than one term short. Evaluated to a number, the roots it cannot +write are found to the working precision by the Durand–Kerner iteration, all of them or none. A set +with a symbol among its members, whose members may coincide, and a set that is not finite are left +as written. [#1285](https://github.com/asc-community/AngouriMath/issues/1285). Both columns +measured on a build, `v2.5.0` against this change. + +| | Was (2.5.0) | Is | +|---|---|---| +| `"sum(x^2, x in {1, 2, 3})".ToEntity()` | `FunctionArgumentCountException`: *sum should have exactly 4 arguments but 2 arguments are provided* | `14` once simplified | +| `"sum(w^2, w in { w : w^3 - 2w + 1 = 0 })".ToEntity().Simplify()` | the same exception | `4` | +| `"sum(w^5, w in { w : w^5 + w + 3 = 0 })".ToEntity().Evaled` | the same exception | `-15`, from the five roots found numerically | +| `"sum(k, k in ZZ intersect [1; 10])".ToEntity().Simplify()` | the same exception | `55`, as `sum(k, k, 1, 10)` is | +| `"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 + 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/.editorconfig b/Sources/.editorconfig index 9a89de0e1..d964644f0 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -279,6 +279,12 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [AngouriMath/Functions/Continuous/ExtremumOverSet.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n +[AngouriMath/Functions/Algebra/SumOverSet.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + +[AngouriMath/Functions/Continuous/Solvers/NumericalSolving/DurandKerner.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + [Tests/UnitTests/Core/Transformations/LogarithmSpellingsTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Convenience/MathS.cs b/Sources/AngouriMath/Convenience/MathS.cs index 4ba860f7d..f96334374 100644 --- a/Sources/AngouriMath/Convenience/MathS.cs +++ b/Sources/AngouriMath/Convenience/MathS.cs @@ -7476,6 +7476,38 @@ public static Entity Limit(Entity expr, Entity var, Entity dest, ApproachFrom ap public static Entity Sum(Entity expr, Entity var, Entity from, Entity to) => new Summationf(expr, var, from, to); + /// + /// The sum of over the members of as + /// ranges over it, each member counted once: written + /// sum(expr, var in over). A sum over the roots of a polynomial is a sum over its + /// set of roots, sum(f(w), w in { w : p(w) = 0 }), with no index and no order. + /// + /// The term, as a function of . + /// The name that ranges over the set, bound by the sum. + /// The set. + /// + /// The terms added up where the set is finite and its members are known: a listed set of + /// numbers, or the roots of a polynomial that the solver writes in full. Evaluated to a + /// number, the roots it cannot write are found numerically, all of them or none. Left as + /// written otherwise: a set with a symbol among its members, which may coincide, or a set + /// that is not finite. + /// + /// #1285 + /// + /// + /// using System; + /// using static AngouriMath.MathS; + /// + /// Console.WriteLine(Sum("w^2", "w", "{ w : w^3 - 2w + 1 = 0 }").Simplify()); + /// + /// Prints + /// + /// 4 + /// + /// + public static Entity Sum(Entity expr, Entity var, Entity over) + => new SumOverSetf(expr, var, over); + /// /// The largest value takes as ranges over /// , written max(expr, var in over). Answered over a finite diff --git a/Sources/AngouriMath/Core/Antlr/AngouriMath.g b/Sources/AngouriMath/Core/Antlr/AngouriMath.g index 26183810f..e4576d85c 100644 --- a/Sources/AngouriMath/Core/Antlr/AngouriMath.g +++ b/Sources/AngouriMath/Core/Antlr/AngouriMath.g @@ -571,7 +571,7 @@ atom returns[Entity value] | 'limit(' args = function_arguments ')' { Assert("limit", 3, $args.list.Count); $value = MathS.Limit($args.list[0], $args.list[1], $args.list[2]); } | 'limitleft(' args = function_arguments ')' { Assert("limitleft", 3, $args.list.Count); $value = MathS.Limit($args.list[0], $args.list[1], $args.list[2], AngouriMath.Core.ApproachFrom.Left); } | 'limitright(' args = function_arguments ')' { Assert("limitright", 3, $args.list.Count); $value = MathS.Limit($args.list[0], $args.list[1], $args.list[2], AngouriMath.Core.ApproachFrom.Right); } - | 'sum(' args = function_arguments ')' { Assert("sum", 4, $args.list.Count); $value = MathS.Sum($args.list[0], $args.list[1], $args.list[2], $args.list[3]); } + | 'sum(' args = function_arguments ')' { if ($args.list.Count == 2 && AngouriMath.Functions.ExtremumOverSet.AsRange($args.list[1]) is var (sumVar, sumOver)) $value = MathS.Sum($args.list[0], sumVar, sumOver); else { Assert("sum", 4, $args.list.Count); $value = MathS.Sum($args.list[0], $args.list[1], $args.list[2], $args.list[3]); } } | 'product(' args = function_arguments ')' { Assert("product", 4, $args.list.Count); $value = MathS.Product($args.list[0], $args.list[1], $args.list[2], $args.list[3]); } | 'signum(' args = function_arguments ')' { Assert("signum", 1, $args.list.Count); $value = MathS.Signum($args.list[0]); } | 'sgn(' args = function_arguments ')' { Assert("sgn", 1, $args.list.Count); $value = MathS.Signum($args.list[0]); } diff --git a/Sources/AngouriMath/Core/Antlr/AngouriMathParser.cs b/Sources/AngouriMath/Core/Antlr/AngouriMathParser.cs index 7fd0ac2a4..1e11d2dc1 100644 --- a/Sources/AngouriMath/Core/Antlr/AngouriMathParser.cs +++ b/Sources/AngouriMath/Core/Antlr/AngouriMathParser.cs @@ -4449,7 +4449,7 @@ public AtomContext atom() { _localctx.args = function_arguments(); State = 916; Match(T__30); - Assert("sum", 4, _localctx.args.list.Count); _localctx.value = MathS.Sum(_localctx.args.list[0], _localctx.args.list[1], _localctx.args.list[2], _localctx.args.list[3]); + if (_localctx.args.list.Count == 2 && AngouriMath.Functions.ExtremumOverSet.AsRange(_localctx.args.list[1]) is var (sumVar, sumOver)) _localctx.value = MathS.Sum(_localctx.args.list[0], sumVar, sumOver); else { Assert("sum", 4, _localctx.args.list.Count); _localctx.value = MathS.Sum(_localctx.args.list[0], _localctx.args.list[1], _localctx.args.list[2], _localctx.args.list[3]); } } break; case 100: diff --git a/Sources/AngouriMath/Core/Domains.Classes.cs b/Sources/AngouriMath/Core/Domains.Classes.cs index adc24a380..1da51eaa3 100644 --- a/Sources/AngouriMath/Core/Domains.Classes.cs +++ b/Sources/AngouriMath/Core/Domains.Classes.cs @@ -280,6 +280,14 @@ partial record Productf internal override Domain DefaultCodomain => Domain.Complex; } + partial record SumOverSetf + { + /// + public override Domain Codomain { get; protected init; } = Domain.Complex; + /// + internal override Domain DefaultCodomain => Domain.Complex; + } + partial record Maximumf { /// diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Calculus.Classes.cs b/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Calculus.Classes.cs index 0bd9c7fbc..bf85c70b2 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Calculus.Classes.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Calculus.Classes.cs @@ -119,6 +119,29 @@ public override Entity Replace(Func func) => protected override Entity[] InitDirectChildren() => new[] { Expression, Var, From, To }; } + /// + /// A sum over a set: sum(f(x), x in S) is the sum of f(x) over the members + /// x of S, each counted once. A binder, like , whose + /// range is a set rather than two bounds, and so has no index and no order: the sum over + /// the roots of a polynomial is sum(f(w), w in { w : p(w) = 0 }), which needs no + /// first root. + /// https://github.com/asc-community/AngouriMath/issues/1285 + /// + public sealed partial record SumOverSetf(Entity Expression, Entity Var, Entity Over) : CalculusOperator(Expression, Var) + { + /// The set the variable ranges over. + public Entity Over { get; init; } = Binding.Of(Var).In(Over); + + private SumOverSetf New(Entity expression, Entity var, Entity over) => + ReferenceEquals(Expression, expression) && ReferenceEquals(Var, var) && ReferenceEquals(Over, over) + ? this : new(expression, var, over) { Codomain = Codomain }; + /// + public override Entity Replace(Func func) => + func(New(Expression.Replace(func), Var, Over.Replace(func))); + /// + protected override Entity[] InitDirectChildren() => new[] { Expression, Var, Over }; + } + /// /// A node of limit /// diff --git a/Sources/AngouriMath/Core/Entity/Entity.Definition.cs b/Sources/AngouriMath/Core/Entity/Entity.Definition.cs index 21ad2c884..ff904fe04 100644 --- a/Sources/AngouriMath/Core/Entity/Entity.Definition.cs +++ b/Sources/AngouriMath/Core/Entity/Entity.Definition.cs @@ -694,6 +694,9 @@ public IReadOnlyList Vars => BoundBy(index, body, from, to), Productf(var body, var index, var from, var to) => BoundBy(index, body, from, to), + // A sum over a set binds the name that ranges over the set: + // sum(w * a, w in { w : w^2 = 2 }) is a function of a alone. + SumOverSetf(var body, var bound, var over) => BoundBy(bound, body, over), // An extremum over a set binds the variable that ranges over the set: // max(t^2 + a, t in [0; 1]) is a function of a alone. Maximumf(var body, var bound, var over) => BoundBy(bound, body, over), diff --git a/Sources/AngouriMath/Core/Serialization/Entity.Serialization.cs b/Sources/AngouriMath/Core/Serialization/Entity.Serialization.cs index 6cbef6aa1..c0f7a6f05 100644 --- a/Sources/AngouriMath/Core/Serialization/Entity.Serialization.cs +++ b/Sources/AngouriMath/Core/Serialization/Entity.Serialization.cs @@ -135,6 +135,7 @@ [EntityJsonConverter] partial record SpecialSet [EntityJsonConverter] partial record Sumf; [EntityJsonConverter] partial record Summationf; [EntityJsonConverter] partial record Maximumf; + [EntityJsonConverter] partial record SumOverSetf; [EntityJsonConverter] partial record Minimumf; [EntityJsonConverter] partial record Argmaxf; [EntityJsonConverter] partial record Argminf; diff --git a/Sources/AngouriMath/Docs/Usage/Syntax.md b/Sources/AngouriMath/Docs/Usage/Syntax.md index 0674a0f60..e81844e15 100644 --- a/Sources/AngouriMath/Docs/Usage/Syntax.md +++ b/Sources/AngouriMath/Docs/Usage/Syntax.md @@ -328,8 +328,9 @@ undefined at `0`; `Si(+oo)` is `pi / 2` and `Ci(+oo)` is `0`. **Calculus** — `derivative(expr, var, order)`, `integral(expr, var)`, `integral(expr, var, from, to)`, `limit(expr, var, dest)`, `limitleft(...)`, `limitright(...)`; -`max(expr, var in set)` and `min` for the extremum of an expression over a set, `argmax` and -`argmin` for the set of points where it is taken — `max(a, b)` of two values is still the larger. +`sum(expr, var in set)` for a sum over a set; `max(expr, var in set)` and `min` for the extremum of +an expression over a set, `argmax` and `argmin` for the set of points where it is taken — +`max(a, b)` of two values is still the larger. `min(S)` with one argument that is a set is its least member, `min(x, x in S)`, and the range may carry conditions, `min(x, x in PP and x > 14)` being the least of `{ x in PP : x > 14 }`; decided for `PP`, `ZZ+`, `ZZ*` and their cuts by an interval or bounds (`min(ZZ+)` is `1`, @@ -487,6 +488,16 @@ at `c = 0`, where the empty product is `1`. Where the index appears in the body must be a concrete integer of at least one, `factorial` having no value at the negative integers; `product(k, k, 0, n)` therefore stays as written. +`sum(body, name in set)` is a sum **over a set**: the body at each member of the set, each member +counted once, added up. `sum(x^2, x in {1, 2, 3})` is `14`. A sum over the roots of a polynomial is a +sum over its set of roots, `sum(w^2, w in { w : w^3 - 2w + 1 = 0 })`, which is `4`, and it needs no +index and no order, a set having neither. The terms are added up where the set is finite and its +members are known: a listed set of numbers, or roots the solver writes in full, which is checked by +counting them against the polynomial's number of distinct roots. A set whose members might coincide, +such as `{ a, b }`, and a set that is not finite are left as written. Evaluated to a number, roots +the solver cannot write are found numerically, all of them to the working precision or none, so +`sum(w^5, w in { w : w^5 + w + 3 = 0 })` evaluates to `-15`. + **The declared name means the name, and only inside the operator that declares it.** Two names show this, because both mean something else in this language: diff --git a/Sources/AngouriMath/Functions/Algebra/SumOverSet.cs b/Sources/AngouriMath/Functions/Algebra/SumOverSet.cs new file mode 100644 index 000000000..8171af6b0 --- /dev/null +++ b/Sources/AngouriMath/Functions/Algebra/SumOverSet.cs @@ -0,0 +1,220 @@ +// +// 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 PeterO.Numbers; +using static AngouriMath.Entity; +using static AngouriMath.Entity.Number; +using static AngouriMath.Entity.Set; + +namespace AngouriMath.Functions +{ + /// + /// The value of a sum over a set, sum(f(x), x in S): the terms written out and added + /// where the members of S are known and there are finitely many of them, and the sum + /// left as written otherwise. + /// + /// + /// + /// A sum over a set counts each member once, which is what makes it well defined without an + /// order: a set has no first member, and addition does not care. So the members have to be + /// known to be distinct before they are added, and a listed set is read only where its + /// elements are numbers. sum(x, x in { a, b }) is a + b when a and + /// b differ and a when they are equal, and nothing here can tell which, so it is + /// left as written. + /// + /// + /// The sum over the roots of a polynomial is the case the node exists for: + /// sum(f(w), w in { w : p(w) = 0 }). Where the solver writes every root in closed form + /// the roots are added; every is checked, by counting the distinct roots it gave + /// against the degree of the square-free part of p, which is how many distinct roots + /// p has. A root set the solver answers only in part is left as written, since a sum + /// that misses a term is a wrong answer. Evaluated to a number, the roots the solver cannot + /// write are found to the working precision, all of them or none, by + /// . And a sum over a set that is not finite -- an + /// interval, the integers -- is a series, which is not this node's to answer. + /// #1285 + /// + /// + internal static class SumOverSet + { + /// + /// The sum, or where the members of the set are not known or not + /// finitely many. + /// + internal static Entity? Value(Entity expression, Entity var, Entity over, bool isExact) + { + if (var is not Variable x || over.InnerSimplified is not Set set) + return null; + if ((Members(set) ?? (isExact ? null : NumericMembers(set))) is not { } members + || members.Count > Summationf.MaxExpandedTerms) + return null; + Entity? sum = null; + foreach (var member in members) + { + var term = expression.Substitute(x, member); + sum = sum is null ? term : sum + term; + } + sum ??= Integer.Zero; + return isExact ? sum.InnerSimplified : sum.Evaled; + } + + /// + /// The members of , each once, or where they + /// are not known to be distinct or there are not finitely many. + /// + private static IReadOnlyList? Members(Set set) + => set switch + { + FiniteSet finite => Distinct(finite.Elements), + // { w in S : p(w) = 0 }, which the set builder holds as w in S and p(w) = 0: the + // roots that are members of S. + ConditionalSet { Var: Variable w } builder when builder.DeclaredMembership is (Set universe, Equalsf(var left, var right)) + => Roots(left - right, w) is { } roots ? Within(roots, universe) : null, + ConditionalSet { Var: Variable w, Predicate: Equalsf(var left, var right) } => Roots(left - right, w), + _ => null, + }; + + /// + /// The members of that are in , or + /// where that is not decided for one of them. + /// + private static IReadOnlyList? Within(IReadOnlyList members, Set universe) + { + var inside = new List(); + foreach (var member in members) + { + if (!universe.TryContains(member, out var contains)) + return null; + if (contains) + inside.Add(member); + } + return inside; + } + + /// + /// The elements, where each is a number and no two are equal; + /// otherwise. Two numbers written differently can be the same number, so they are + /// compared by value, and a pair that agrees to the working precision is taken for one + /// that may be equal and the set is not read. + /// + private static IReadOnlyList? Distinct(IEnumerable elements) + { + var members = new List(); + var values = new List(); + foreach (var element in elements) + { + if (element.Evaled is not Complex value || !value.IsFinite) + return null; + foreach (var seen in values) + if (seen == value) + return null; + members.Add(element); + values.Add(value); + } + return members; + } + + /// + /// The roots of the polynomial in , each + /// once, where the solver writes all of them; otherwise. + /// + /// + /// The solver is asked only where it can write every root -- each irreducible factor of + /// degree at most four, or a binomial -- since past that there is nothing it could write, + /// and asking would cost a search every time a simplification meets the sum. + /// + private static IReadOnlyList? Roots(Entity polynomial, Variable w) + { + if (SquareFreeParts(polynomial, w) is not { } parts) + return null; + var expected = 0; + foreach (var part in parts) + expected += part.Factor.Degree; + if (expected > 4 && !(PolynomialFactorization.FactorComplete(polynomial, w) is { } factored + && factored.Parts.All(static part => part.Factor.Degree <= 4 || IsBinomial(part.Factor)))) + return null; + if (polynomial.SolveEquation(w) is not FiniteSet solutions + || Distinct(solutions.Elements.Select(static root => root.InnerSimplified)) is not { } roots) + return null; + return roots.Count == expected ? roots : null; + } + + /// + /// The members of a set of roots to the working precision, where the solver cannot write + /// them: each square-free part of the polynomial by , + /// all its roots or none. The parts share no root, so the roots are distinct. + /// + private static IReadOnlyList? NumericMembers(Set set) + => set switch + { + ConditionalSet { Var: Variable w } builder when builder.DeclaredMembership is (Set universe, Equalsf(var left, var right)) + => NumericRoots(left - right, w) is { } roots ? Within(roots, universe) : null, + ConditionalSet { Var: Variable w, Predicate: Equalsf(var left, var right) } => NumericRoots(left - right, w), + _ => null, + }; + + private static IReadOnlyList? NumericRoots(Entity polynomial, Variable w) + { + if (SquareFreeParts(polynomial, w) is not { } parts) + return null; + var context = MathS.Settings.DecimalPrecisionContext.Value; + var roots = new List(); + foreach (var part in parts) + { + if (part.Factor.Degree < 1) + continue; + if (Algebra.NumericalSolving.DurandKerner.Roots(part.Factor, context) is not { } found) + return null; + roots.AddRange(found); + } + return roots; + } + + /// + /// The square-free decomposition of in , + /// where it is a polynomial with rational coefficients of degree at least one; its parts' + /// degrees add up to the number of distinct roots. otherwise. + /// + /// + /// Read as reads a polynomial, off the tree, so the + /// coefficients are the rationals written rather than values evaluated: with downcasting + /// off an evaluated 3 is a decimal and not the rational 3, and a reading + /// through evaluation would refuse the polynomial -- in the setting numerical checks + /// run in. + /// + private static IReadOnlyList? SquareFreeParts(Entity polynomial, Variable w) + { + if (MultivariatePolynomial.TryParse(polynomial, new Dictionary { [w] = 0 }) is not { } parsed) + return null; + var degree = parsed.DegreeIn(0); + if (degree < 1 || degree > IntegerPolynomial.MaxDegree) + return null; + var rational = new ERational[degree + 1]; + for (var i = 0; i < rational.Length; i++) + rational[i] = ERational.Zero; + foreach (var term in parsed.Terms) + rational[MultivariatePolynomial.PowerOf(term.Key, 0)] = term.Value; + var denominator = EInteger.One; + foreach (var coefficient in rational) + denominator = denominator.Divide(denominator.Gcd(coefficient.Denominator)).Multiply(coefficient.Denominator); + var whole = new EInteger[rational.Length]; + for (var i = 0; i < whole.Length; i++) + whole[i] = rational[i].Numerator.Multiply(denominator.Divide(rational[i].Denominator)); + var poly = IntegerPolynomial.Create(whole); + return poly.Degree < 1 ? null : SquareFreeDecomposition.Decompose(poly.PrimitivePart()); + } + + private static bool IsBinomial(IntegerPolynomial poly) + { + var terms = 0; + for (var power = 0; power <= poly.Degree; power++) + if (!poly[power].IsZero) + terms++; + return terms == 2; + } + } +} diff --git a/Sources/AngouriMath/Functions/Continuous/Differentiation.cs b/Sources/AngouriMath/Functions/Continuous/Differentiation.cs index 5a6a72fd7..9d16d3b6d 100644 --- a/Sources/AngouriMath/Functions/Continuous/Differentiation.cs +++ b/Sources/AngouriMath/Functions/Continuous/Differentiation.cs @@ -476,6 +476,22 @@ protected override Entity InnerDifferentiate(Variable variable) => } #pragma warning restore IDE0054 // Use compound assignment + partial record SumOverSetf + { + /// + /// + /// Term by term, where the set does not depend on the variable: the derivative of a + /// sum of finitely many terms is the sum of their derivatives, and a sum over a set + /// that is not finite is not evaluated in the first place. The name the sum binds is + /// not a free variable of it, so nothing varies with it. + /// https://github.com/asc-community/AngouriMath/issues/1285 + /// + protected override Entity InnerDifferentiate(Variable variable) => + Var == variable ? Integer.Zero + : Over.ContainsNode(variable) ? base.InnerDifferentiate(variable) + : New(Expression.InnerDifferentiate(variable), Var, Over); + } + partial record Limitf { /// diff --git a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/InvertNode.Classes.cs b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/InvertNode.Classes.cs index 3b1a9a60d..9f84d96cf 100644 --- a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/InvertNode.Classes.cs +++ b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/InvertNode.Classes.cs @@ -329,6 +329,13 @@ partial record Integralf : Enumerable.Empty(); } + partial record SumOverSetf + { + // The unknown sits under a binder; see Summationf below. + private protected override IEnumerable? InvertNode(Entity value, Entity x) => + null; + } + partial record Maximumf { // The unknown sits under a binder; see Summationf below. diff --git a/Sources/AngouriMath/Functions/Continuous/Solvers/NumericalSolving/DurandKerner.cs b/Sources/AngouriMath/Functions/Continuous/Solvers/NumericalSolving/DurandKerner.cs new file mode 100644 index 000000000..e11a1f7be --- /dev/null +++ b/Sources/AngouriMath/Functions/Continuous/Solvers/NumericalSolving/DurandKerner.cs @@ -0,0 +1,136 @@ +// +// 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 PeterO.Numbers; +using static AngouriMath.Entity.Number; + +namespace AngouriMath.Functions.Algebra.NumericalSolving +{ + /// + /// Every root of a square-free polynomial with whole coefficients, to the working precision, + /// by the Durand–Kerner iteration: each approximation moves by p(z_i) / prod_{j != i} + /// (z_i - z_j), which is Newton's step for the root it is nearest once the others are + /// close to theirs, so the iteration converges to all the roots at once and quadratically at + /// the end. + /// + /// + /// + /// All the roots or none: a sum over the roots of a polynomial that is missing one is a wrong + /// answer, so where the iteration does not settle, or settles on two approximations too close + /// to tell apart, there is no answer. A square-free polynomial has no repeated root, so two + /// approximations that close mean the iteration has not separated them. + /// + /// + /// The coefficients are real, so the non-real roots come in conjugate pairs, and a root whose + /// imaginary part is smaller than half the separation that was checked is real: were it not, + /// it and its conjugate would be closer together than that. It is written as a real number + /// then, rather than as one with an imaginary part of 1e-110 that no root has. + /// + /// #1285 + /// + internal static class DurandKerner + { + /// Digits carried beyond the working precision, and dropped at the end. + private const int GuardDigits = 10; + + /// How many sweeps over all the roots before the iteration is taken not to settle. + private const int MaxSweeps = 1000; + + /// + /// The roots of , which has degree at least one and no repeated + /// root, at the precision of ; where they + /// are not settled. + /// + internal static IReadOnlyList? Roots(IntegerPolynomial squareFree, EContext context) + { + var degree = squareFree.Degree; + if (degree < 1) + return null; + var digits = context.Precision.ToInt32Checked(); + var work = context.WithPrecision(digits + GuardDigits); + var lead = EDecimal.FromEInteger(squareFree[degree]); + var monic = new EDecimal[degree]; + for (var i = 0; i < degree; i++) + monic[i] = EDecimal.FromEInteger(squareFree[i]).Divide(lead, work); + + // Cauchy's bound: every root is within 1 + max |a_i / a_n| of zero. The approximations + // start spread round that circle, turned off the axes so that no two start as a + // conjugate pair or on a line of symmetry the iteration cannot leave. + var radius = EDecimal.One; + foreach (var coefficient in monic) + radius = EDecimal.Max(radius, EDecimal.One.Add(coefficient.Abs(), work)); + var r = radius.ToDouble(); + var re = new EDecimal[degree]; + var im = new EDecimal[degree]; + for (var i = 0; i < degree; i++) + { + var angle = 2 * Math.PI * i / degree + 0.4; + re[i] = EDecimal.FromDouble(r * Math.Cos(angle)); + im[i] = EDecimal.FromDouble(r * Math.Sin(angle)); + } + + var tolerance = EDecimal.Create(EInteger.One, EInteger.FromInt32(-(digits + 2))).Multiply(radius, work); + var settled = false; + for (var sweep = 0; sweep < MaxSweeps && !settled; sweep++) + { + var largest = EDecimal.Zero; + for (var i = 0; i < degree; i++) + { + // p(z_i) by Horner's scheme, on the monic polynomial. + var (pRe, pIm) = (EDecimal.One, EDecimal.Zero); + for (var k = degree - 1; k >= 0; k--) + { + (pRe, pIm) = Times(pRe, pIm, re[i], im[i], work); + pRe = pRe.Add(monic[k], work); + } + var (qRe, qIm) = (EDecimal.One, EDecimal.Zero); + for (var j = 0; j < degree; j++) + if (j != i) + (qRe, qIm) = Times(qRe, qIm, re[i].Subtract(re[j], work), im[i].Subtract(im[j], work), work); + if (qRe.IsZero && qIm.IsZero) + return null; + var (stepRe, stepIm) = Over(pRe, pIm, qRe, qIm, work); + re[i] = re[i].Subtract(stepRe, work); + im[i] = im[i].Subtract(stepIm, work); + largest = EDecimal.Max(largest, stepRe.Abs().Add(stepIm.Abs(), work)); + } + settled = largest.CompareTo(tolerance) <= 0; + } + if (!settled) + return null; + + // Distinct to half the working precision, which a square-free polynomial's roots are + // unless the iteration has not separated two of them. + var separation = EDecimal.Create(EInteger.One, EInteger.FromInt32(-(digits / 2))); + for (var i = 0; i < degree; i++) + for (var j = i + 1; j < degree; j++) + if (re[i].Subtract(re[j], work).Abs().Add(im[i].Subtract(im[j], work).Abs(), work).CompareTo(separation) <= 0) + return null; + + var half = separation.Divide(EDecimal.FromInt32(2), work); + var roots = new List(degree); + for (var i = 0; i < degree; i++) + { + var imaginary = im[i].Abs().CompareTo(half) < 0 ? EDecimal.Zero : im[i].RoundToPrecision(context); + roots.Add(Complex.Create(re[i].RoundToPrecision(context), imaginary)); + } + return roots; + } + + private static (EDecimal Re, EDecimal Im) Times(EDecimal aRe, EDecimal aIm, EDecimal bRe, EDecimal bIm, EContext work) + => (aRe.Multiply(bRe, work).Subtract(aIm.Multiply(bIm, work), work), + aRe.Multiply(bIm, work).Add(aIm.Multiply(bRe, work), work)); + + private static (EDecimal Re, EDecimal Im) Over(EDecimal aRe, EDecimal aIm, EDecimal bRe, EDecimal bIm, EContext work) + { + var norm = bRe.Multiply(bRe, work).Add(bIm.Multiply(bIm, work), work); + return (aRe.Multiply(bRe, work).Add(aIm.Multiply(bIm, work), work).Divide(norm, work), + aIm.Multiply(bRe, work).Subtract(aRe.Multiply(bIm, work), work).Divide(norm, work)); + } + } +} diff --git a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Calculus.Classes.cs b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Calculus.Classes.cs index dba589426..877c145a0 100644 --- a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Calculus.Classes.cs +++ b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Calculus.Classes.cs @@ -208,6 +208,25 @@ protected override Entity InnerSimplify(bool isExact) => Functions.ExtremumOverSet.Points(Expression, Var, Over, largest: false) ?? this; } + partial record SumOverSetf + { + /// + /// A sum over a set is a well-formed expression whatever the set, as + /// is whatever its bounds; adding its terms up is an + /// evaluation, not a claim about where it is defined. + /// + private protected override Entity IntrinsicCondition => Boolean.True; + /// + /// + /// Where the terms are not added up, the sum is still written in its simplest parts, so + /// that a derivative taken term by term reads as one: 2 y w and not + /// 0 * y ^ 2 + 2 * y ^ 1 * 1 * w. + /// + protected override Entity InnerSimplify(bool isExact) + => Functions.SumOverSet.Value(Expression, Var, Over, isExact) + ?? New(Expression.InnerSimplified, Var, Over.InnerSimplified); + } + // TODO: rewrite this part too public partial record Summationf { diff --git a/Sources/AngouriMath/Functions/Output/Latex/Latex.Calculus.Classes.cs b/Sources/AngouriMath/Functions/Output/Latex/Latex.Calculus.Classes.cs index 195cde5c7..2bd05f539 100644 --- a/Sources/AngouriMath/Functions/Output/Latex/Latex.Calculus.Classes.cs +++ b/Sources/AngouriMath/Functions/Output/Latex/Latex.Calculus.Classes.cs @@ -66,6 +66,14 @@ private protected override string LatexizeNode() } } + public partial record SumOverSetf + { + /// + private protected override string LatexizeNode() => + @"\sum_{" + Var.Latexize() + @" \in " + Over.Latexize() + "} " + + Expression.Latexize(Expression.Priority < Priority.Sum); + } + public partial record Maximumf { /// diff --git a/Sources/AngouriMath/Functions/Output/ToString/ToString.Calculus.Classes.cs b/Sources/AngouriMath/Functions/Output/ToString/ToString.Calculus.Classes.cs index 6f8ab28a0..5962c76d2 100644 --- a/Sources/AngouriMath/Functions/Output/ToString/ToString.Calculus.Classes.cs +++ b/Sources/AngouriMath/Functions/Output/ToString/ToString.Calculus.Classes.cs @@ -36,6 +36,15 @@ Range is var (from, to) public override string ToString() => Stringize(); } + public partial record SumOverSetf + { + /// + private protected override string StringizeNode() => + $"sum({Expression.Stringize()}, {Var.Stringize()} in {Over.Stringize()})"; + /// + public override string ToString() => Stringize(); + } + public partial record Maximumf { /// diff --git a/Sources/AngouriMath/Functions/Output/ToSympy/ToSympy.Calculus.Classes.cs b/Sources/AngouriMath/Functions/Output/ToSympy/ToSympy.Calculus.Classes.cs index f6025d275..c98eb5f8e 100644 --- a/Sources/AngouriMath/Functions/Output/ToSympy/ToSympy.Calculus.Classes.cs +++ b/Sources/AngouriMath/Functions/Output/ToSympy/ToSympy.Calculus.Classes.cs @@ -21,6 +21,20 @@ public partial record Integralf internal override string ToSymPy() => $"sympy.integrate({Expression.ToSymPy()}, {(Range is var (from, to) ? $"({Var.ToSymPy()}, {from.ToSymPy()}, {to.ToSymPy()})" : Var.ToSymPy())})"; } + public partial record SumOverSetf + { + // SymPy's RootSum counts a repeated root as often as it is repeated, and a set counts + // it once, so the polynomial goes over as its square-free part. + internal override string ToSymPy() => Over switch + { + Set.FiniteSet finite => + $"sympy.Add(*[({Expression.ToSymPy()}).subs({Var.ToSymPy()}, member) for member in [{string.Join(", ", finite.Elements.Select(static e => e.ToSymPy()))}]])", + Set.ConditionalSet { Var: Variable w, Predicate: Equalsf(var left, var right) } when w == Var => + $"sympy.RootSum(sympy.sqf_part(sympy.Poly({(left - right).ToSymPy()}, {Var.ToSymPy()})), sympy.Lambda({Var.ToSymPy()}, {Expression.ToSymPy()}))", + _ => throw new NotSufficientlySupportedException("SymPy has no sum over a set; a listed set or the roots of a polynomial can be exported"), + }; + } + public partial record Maximumf { internal override string ToSymPy() => diff --git a/Sources/AngouriMath/Functions/Substitute.cs b/Sources/AngouriMath/Functions/Substitute.cs index a77cf72d5..af6b6af44 100644 --- a/Sources/AngouriMath/Functions/Substitute.cs +++ b/Sources/AngouriMath/Functions/Substitute.cs @@ -637,6 +637,20 @@ public override Entity Substitute(Entity x, Entity value) } } + partial record SumOverSetf + { + /// + // The variable is bound; renamed first, as a summation's index is. + public override Entity Substitute(Entity x, Entity value) + { + if (this == x) + return value; + var replacement = Variable.CreateTemp((x + value + Expression + Var).Vars); + var renamed = Expression.Substitute(Var, replacement).Substitute(x, value).Substitute(replacement, Var); + return New(renamed, Var, Over.Substitute(x, value)); + } + } + partial record Maximumf { /// diff --git a/Sources/AngouriMath/Functions/TreeAnalyzer/Sort.Classes.cs b/Sources/AngouriMath/Functions/TreeAnalyzer/Sort.Classes.cs index a8c5f4c5d..bd4346add 100644 --- a/Sources/AngouriMath/Functions/TreeAnalyzer/Sort.Classes.cs +++ b/Sources/AngouriMath/Functions/TreeAnalyzer/Sort.Classes.cs @@ -201,6 +201,11 @@ public partial record Productf private protected override string SortHashName(SortLevel level) => "productf_"; } + public partial record SumOverSetf + { + private protected override string SortHashName(SortLevel level) => "sumoversetf_"; + } + public partial record Maximumf { private protected override string SortHashName(SortLevel level) => "maximumf_"; diff --git a/Sources/Tests/UnitTests/Calculus/SumOverSetTest.cs b/Sources/Tests/UnitTests/Calculus/SumOverSetTest.cs new file mode 100644 index 000000000..c3f5d7118 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/SumOverSetTest.cs @@ -0,0 +1,134 @@ +// +// 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 AngouriMath; +using AngouriMath.Extensions; +using PeterO.Numbers; +using Xunit; +using static AngouriMath.Entity; +using static AngouriMath.Entity.Number; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// sum(f(x), x in S): the sum of an expression over the members of a set, each counted + /// once, as a binder. https://github.com/asc-community/AngouriMath/issues/1285 + /// + [Trait("Area", "Calculus")] + public sealed class SumOverSetTest + { + [Fact] + public void ItIsWrittenAsItIsRead() + { + var sum = "sum(x^2, x in {1, 2, 3})".ToEntity(); + Assert.IsType(sum); + Assert.Equal("sum(x ^ 2, x in { 1, 2, 3 })", sum.Stringize()); + Assert.Equal(sum, sum.Stringize().ToEntity()); + } + + /// The range form keeps its meaning, and two arguments that are not a range are the error they were. + [Fact] + public void TheRangeFormIsUnchanged() + { + Assert.IsType("sum(k, k, 1, 10)".ToEntity()); + Assert.Equal(Integer.Create(55), "sum(k, k, 1, 10)".ToEntity().Simplify()); + Assert.ThrowsAny(() => "sum(x, y)".ToEntity()); + } + + [Theory] + [InlineData("sum(x^2, x in {1, 2, 3})", "14")] + [InlineData("sum(x, x in {})", "0")] + [InlineData("sum(7, x in {1, 2, 3})", "21")] + // An integer range is listed, so the set form agrees with the range form. + [InlineData("sum(k, k in ZZ intersect [1; 10])", "55")] + [InlineData("sum(w^2, w in { w : w^3 - 2w + 1 = 0 })", "4")] + [InlineData("sum(w, w in { w : w^2 + 1 = 0 })", "0")] + [InlineData("sum(1/w, w in { w : w^2 - 3w + 2 = 0 })", "3/2")] + // Over the reals only the real roots count: w^4 - 1 has 1 and -1 there, and i and -i besides. + [InlineData("sum(w^2, w in { w in RR : w^4 - 1 = 0 })", "2")] + [InlineData("sum(w^2, w in { w : w^4 - 1 = 0 })", "0")] + public void TheTermsAreAddedUp(string sum, string value) + => Assert.Equal(value.ToEntity(), sum.ToEntity().Simplify()); + + /// A set counts a member once, so a double root is one term. + [Fact] + public void ARepeatedRootIsOneMember() + => Assert.Equal(Integer.Create(3), "sum(w, w in { w : (w - 1)^2 (w - 2) = 0 })".ToEntity().InnerSimplified); + + /// + /// Left as written: members that may coincide, a set that is not finite, and roots the + /// solver cannot write, which are not added up symbolically. + /// + [Theory] + [InlineData("sum(x, x in {a, b})")] + [InlineData("sum(x, x in [0; 1])")] + [InlineData("sum(k, k in ZZ)")] + [InlineData("sum(w, w in { w : w^5 + w + 3 = 0 })")] + public void WhatItCannotListIsLeftAsWritten(string sum) + => Assert.IsType(sum.ToEntity().InnerSimplified); + + /// + /// Evaluated to a number, the roots the solver cannot write are found numerically, all of + /// them, and the sum agrees with the exact value the coefficients give: for + /// w^5 + w + 3 the roots add to 0, their reciprocals to -1/3, and + /// their fifth powers to -15, since each root has w^5 = -w - 3. + /// + [Theory] + [InlineData("sum(w, w in { w : w^5 + w + 3 = 0 })", "0")] + [InlineData("sum(1/w, w in { w : w^5 + w + 3 = 0 })", "-1/3")] + [InlineData("sum(w^5, w in { w : w^5 + w + 3 = 0 })", "-15")] + public void NumericallyOverRootsTheSolverCannotWrite(string sum, string exact) + { + AssertClose(exact.ToEntity().Evaled, sum.ToEntity().Evaled); + // And evaluated with downcasting off, as a numerical check evaluates an answer that was + // built as usual: parsed afresh, since a cached entity would hand back the value it was + // evaluated to with downcasting on, and parsed outside the scope, since there a literal + // 3 is a decimal and nothing is a polynomial with rational coefficients. + var (fresh, reference) = (MathS.FromString(sum, useCache: false), MathS.FromString(exact, useCache: false)); + using var _ = MathS.Settings.DowncastingEnabled.Set(false); + AssertClose(reference.Evaled, fresh.Evaled); + } + + private static void AssertClose(Entity expected, Entity actual) + { + var e = Assert.IsAssignableFrom(expected); + var a = Assert.IsAssignableFrom(actual); + var tolerance = EDecimal.FromString("1E-50"); + Assert.True(e.RealPart.EDecimal.Subtract(a.RealPart.EDecimal).Abs().CompareTo(tolerance) < 0, $"{a} is not {e}"); + Assert.True(e.ImaginaryPart.EDecimal.Subtract(a.ImaginaryPart.EDecimal).Abs().CompareTo(tolerance) < 0, $"{a} is not {e}"); + } + + /// The name ranges over the set, so the sum is not a function of it. + [Fact] + public void TheNameIsBound() + { + var sum = "sum(w * a, w in { w : w^5 + w + 3 = 0 })".ToEntity(); + Assert.Equal(new[] { MathS.Var("a") }, sum.FreeVariables); + Assert.Equal(sum, sum.Substitute("w", 5)); + Assert.Equal("sum(w * 2, w in { w : w ^ 5 + w + 3 = 0 })".ToEntity(), sum.Substitute("a", 2)); + } + + /// + /// Term by term where the set does not depend on the variable, and zero in the bound name: + /// the derivative of y^2 times the fifth powers of the roots, which add to + /// -15, is a sum over the same set with the value -30 y. + /// + [Fact] + public void ItIsDifferentiatedTermByTerm() + { + var sum = "sum(w^5 * y^2, w in { w : w^5 + w + 3 = 0 })".ToEntity(); + var derivative = sum.Differentiate("y"); + Assert.IsType(derivative); + AssertClose(Integer.Create(-30), derivative.Substitute("y", 1).Evaled); + Assert.Equal(Integer.Zero, sum.Differentiate("w")); + } + + [Fact] + public void ItPrintsAsASumUnderItsRange() + => Assert.Equal(@"\sum_{x \in \left\{ 1, 2, 3 \right\}} {x}^{2}", "sum(x^2, x in {1, 2, 3})".ToEntity().Latexize()); + } +} diff --git a/Sources/Tests/UnitTests/Common/PublicApi.txt b/Sources/Tests/UnitTests/Common/PublicApi.txt index d7302cae4..a7e36a4ca 100644 --- a/Sources/Tests/UnitTests/Common/PublicApi.txt +++ b/Sources/Tests/UnitTests/Common/PublicApi.txt @@ -2618,6 +2618,25 @@ AngouriMath.Entity+Statement.ctor() AngouriMath.Entity+Statement.ctor(AngouriMath.Entity+Statement) AngouriMath.Entity+Statement.op_Equality(AngouriMath.Entity+Statement, AngouriMath.Entity+Statement) : System.Boolean AngouriMath.Entity+Statement.op_Inequality(AngouriMath.Entity+Statement, AngouriMath.Entity+Statement) : System.Boolean +AngouriMath.Entity+SumOverSetf.$() : AngouriMath.Entity+SumOverSetf +AngouriMath.Entity+SumOverSetf.Codomain { } : AngouriMath.Core.Domain +AngouriMath.Entity+SumOverSetf.Deconstruct(AngouriMath.Entity&, AngouriMath.Entity&, AngouriMath.Entity&) : System.Void +AngouriMath.Entity+SumOverSetf.EqualityContract { } : System.Type +AngouriMath.Entity+SumOverSetf.Equals(AngouriMath.Entity+CalculusOperator) : System.Boolean +AngouriMath.Entity+SumOverSetf.Equals(AngouriMath.Entity+SumOverSetf) : System.Boolean +AngouriMath.Entity+SumOverSetf.Equals(System.Object) : System.Boolean +AngouriMath.Entity+SumOverSetf.GetHashCode() : System.Int32 +AngouriMath.Entity+SumOverSetf.InitDirectChildren() : AngouriMath.Entity[] +AngouriMath.Entity+SumOverSetf.InnerDifferentiate(AngouriMath.Entity+Variable) : AngouriMath.Entity +AngouriMath.Entity+SumOverSetf.InnerSimplify(System.Boolean) : AngouriMath.Entity +AngouriMath.Entity+SumOverSetf.Over { } : AngouriMath.Entity +AngouriMath.Entity+SumOverSetf.PrintMembers(System.Text.StringBuilder) : System.Boolean +AngouriMath.Entity+SumOverSetf.Replace(System.Func) : AngouriMath.Entity +AngouriMath.Entity+SumOverSetf.Substitute(AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity +AngouriMath.Entity+SumOverSetf.ToString() : System.String +AngouriMath.Entity+SumOverSetf.ctor(AngouriMath.Entity, AngouriMath.Entity, AngouriMath.Entity) +AngouriMath.Entity+SumOverSetf.op_Equality(AngouriMath.Entity+SumOverSetf, AngouriMath.Entity+SumOverSetf) : System.Boolean +AngouriMath.Entity+SumOverSetf.op_Inequality(AngouriMath.Entity+SumOverSetf, AngouriMath.Entity+SumOverSetf) : System.Boolean AngouriMath.Entity+Sumf.$() : AngouriMath.Entity+Sumf AngouriMath.Entity+Sumf.Addend { } : AngouriMath.Entity AngouriMath.Entity+Sumf.Augend { } : AngouriMath.Entity @@ -3235,6 +3254,7 @@ AngouriMath.MathS.SolveEquation(AngouriMath.Entity, AngouriMath.Entity+Variable) AngouriMath.MathS.SolveOde(AngouriMath.Entity, AngouriMath.Entity+Variable, AngouriMath.Entity+Variable) : AngouriMath.Entity AngouriMath.MathS.Sqr(AngouriMath.Entity) : AngouriMath.Entity AngouriMath.MathS.Sqrt(AngouriMath.Entity) : AngouriMath.Entity +AngouriMath.MathS.Sum(AngouriMath.Entity, AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity AngouriMath.MathS.Sum(AngouriMath.Entity, AngouriMath.Entity, AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity AngouriMath.MathS.Tan(AngouriMath.Entity) : AngouriMath.Entity AngouriMath.MathS.ToBaseN(AngouriMath.Entity+Number+Real, System.Int32) : System.String @@ -3425,6 +3445,7 @@ class AngouriMath.Entity+Sif class AngouriMath.Entity+Signumf class AngouriMath.Entity+Sinf class AngouriMath.Entity+Statement +class AngouriMath.Entity+SumOverSetf class AngouriMath.Entity+Sumf class AngouriMath.Entity+Summationf class AngouriMath.Entity+Tanf From 6a55dbdeaf89cc0754f5323d1411f74ff7390d8d Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 30 Sep 2026 11:32:11 +0000 Subject: [PATCH 2/2] Roots are written out only where they read as simply as the sum Rational roots, the roots of a quadratic and the roots of a binomial are written and added; the roots of an irreducible cubic or quartic are not. In radicals they read no more simply than the sum over them, and for a real cubic with three real roots Cardano's formula writes them with complex radicals, so sum(w^2, w in { w : w^3 + w + 1 = 0 }) is kept as it is and evaluates to -2. A numeric value goes to the numeric roots at once rather than asking the solver first, since a check evaluates at many points and the roots do not depend on them. SquareFreeParts and NumericRoots are internal, for the interval evaluation that reads them next. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 8 ++-- Sources/AngouriMath/Convenience/MathS.cs | 4 +- Sources/AngouriMath/Docs/Usage/Syntax.md | 6 ++- .../Functions/Algebra/SumOverSet.cs | 38 ++++++++++++------- .../UnitTests/Calculus/SumOverSetTest.cs | 5 +++ 5 files changed, 42 insertions(+), 19 deletions(-) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index c3d2e3b62..9e8a0d10b 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1619,9 +1619,11 @@ polynomial is a sum over its set of roots, `sum(f(w), w in { w : p(w) = 0 })`, w order. The terms are added up where the set is finite and its members are known: a listed set of numbers, -distinct by value, and the roots of a polynomial where the solver writes all of them. That is -checked by counting them against the degree of the square-free part, so a root the solver missed -leaves the sum as written rather than one term short. Evaluated to a number, the roots it cannot +distinct by value, and the roots of a polynomial where every one is rational, a root of a quadratic +or a root of a binomial. That is checked by counting them against the degree of the square-free +part, so a root the solver missed leaves the sum as written rather than one term short. The roots of +an irreducible cubic or quartic are not written out in radicals, which read no more simply than the +sum. Evaluated to a number, the roots it cannot write are found to the working precision by the Durand–Kerner iteration, all of them or none. A set with a symbol among its members, whose members may coincide, and a set that is not finite are left as written. [#1285](https://github.com/asc-community/AngouriMath/issues/1285). Both columns diff --git a/Sources/AngouriMath/Convenience/MathS.cs b/Sources/AngouriMath/Convenience/MathS.cs index f96334374..6252e36a3 100644 --- a/Sources/AngouriMath/Convenience/MathS.cs +++ b/Sources/AngouriMath/Convenience/MathS.cs @@ -7487,7 +7487,9 @@ public static Entity Sum(Entity expr, Entity var, Entity from, Entity to) /// The set. /// /// The terms added up where the set is finite and its members are known: a listed set of - /// numbers, or the roots of a polynomial that the solver writes in full. Evaluated to a + /// numbers, or roots of a polynomial that are rational, roots of a quadratic or roots of a + /// binomial -- those of an irreducible cubic or quartic are kept as the sum, radicals + /// reading no more simply than it. Evaluated to a /// number, the roots it cannot write are found numerically, all of them or none. Left as /// written otherwise: a set with a symbol among its members, which may coincide, or a set /// that is not finite. diff --git a/Sources/AngouriMath/Docs/Usage/Syntax.md b/Sources/AngouriMath/Docs/Usage/Syntax.md index e81844e15..ab54245aa 100644 --- a/Sources/AngouriMath/Docs/Usage/Syntax.md +++ b/Sources/AngouriMath/Docs/Usage/Syntax.md @@ -492,8 +492,10 @@ must be a concrete integer of at least one, `factorial` having no value at the n counted once, added up. `sum(x^2, x in {1, 2, 3})` is `14`. A sum over the roots of a polynomial is a sum over its set of roots, `sum(w^2, w in { w : w^3 - 2w + 1 = 0 })`, which is `4`, and it needs no index and no order, a set having neither. The terms are added up where the set is finite and its -members are known: a listed set of numbers, or roots the solver writes in full, which is checked by -counting them against the polynomial's number of distinct roots. A set whose members might coincide, +members are known: a listed set of numbers, or roots that are rational, roots of a quadratic or roots +of a binomial, which is checked by counting them against the polynomial's number of distinct roots. +The roots of an irreducible cubic or quartic are not written out in radicals, which read no more +simply than the sum. A set whose members might coincide, such as `{ a, b }`, and a set that is not finite are left as written. Evaluated to a number, roots the solver cannot write are found numerically, all of them to the working precision or none, so `sum(w^5, w in { w : w^5 + w + 3 = 0 })` evaluates to `-15`. diff --git a/Sources/AngouriMath/Functions/Algebra/SumOverSet.cs b/Sources/AngouriMath/Functions/Algebra/SumOverSet.cs index 8171af6b0..ccb8764b1 100644 --- a/Sources/AngouriMath/Functions/Algebra/SumOverSet.cs +++ b/Sources/AngouriMath/Functions/Algebra/SumOverSet.cs @@ -28,10 +28,12 @@ namespace AngouriMath.Functions /// /// /// The sum over the roots of a polynomial is the case the node exists for: - /// sum(f(w), w in { w : p(w) = 0 }). Where the solver writes every root in closed form - /// the roots are added; every is checked, by counting the distinct roots it gave - /// against the degree of the square-free part of p, which is how many distinct roots - /// p has. A root set the solver answers only in part is left as written, since a sum + /// sum(f(w), w in { w : p(w) = 0 }). Where every root is rational, a root of a + /// quadratic or a root of a binomial, the solver writes them and they are added; every + /// is checked, by counting the distinct roots it gave against the degree of the square-free + /// part of p, which is how many distinct roots p has. The roots of an + /// irreducible cubic or quartic are not written out in radicals, which read no more simply + /// than the sum. A root set the solver answers only in part is left as written, since a sum /// that misses a term is a wrong answer. Evaluated to a number, the roots the solver cannot /// write are found to the working precision, all of them or none, by /// . And a sum over a set that is not finite -- an @@ -49,8 +51,10 @@ internal static class SumOverSet { if (var is not Variable x || over.InnerSimplified is not Set set) return null; - if ((Members(set) ?? (isExact ? null : NumericMembers(set))) is not { } members - || members.Count > Summationf.MaxExpandedTerms) + // A numeric value needs the roots as numbers, and those do not need the solver: asking + // it first would be a solve per evaluation, and a check evaluates at many points. + var members = isExact ? Members(set) : NumericMembers(set) ?? Members(set); + if (members is null || members.Count > Summationf.MaxExpandedTerms) return null; Entity? sum = null; foreach (var member in members) @@ -123,9 +127,13 @@ internal static class SumOverSet /// once, where the solver writes all of them; otherwise. /// /// - /// The solver is asked only where it can write every root -- each irreducible factor of - /// degree at most four, or a binomial -- since past that there is nothing it could write, - /// and asking would cost a search every time a simplification meets the sum. + /// Only roots that read as simply as the sum does are written out: rational ones, the + /// roots of a quadratic, and the roots of a binomial. The roots of an irreducible cubic or + /// quartic in radicals are no simpler than the sum over them, and for a real polynomial + /// with three real roots Cardano's formula writes them with complex radicals, so that sum + /// is kept as it is; its value is still found numerically. The solver is not asked for + /// anything past that, which also spares a search every time a simplification meets the + /// sum. /// private static IReadOnlyList? Roots(Entity polynomial, Variable w) { @@ -134,8 +142,8 @@ internal static class SumOverSet var expected = 0; foreach (var part in parts) expected += part.Factor.Degree; - if (expected > 4 && !(PolynomialFactorization.FactorComplete(polynomial, w) is { } factored - && factored.Parts.All(static part => part.Factor.Degree <= 4 || IsBinomial(part.Factor)))) + if (expected > 2 && !(PolynomialFactorization.FactorComplete(polynomial, w) is { } factored + && factored.Parts.All(static part => part.Factor.Degree <= 2 || IsBinomial(part.Factor)))) return null; if (polynomial.SolveEquation(w) is not FiniteSet solutions || Distinct(solutions.Elements.Select(static root => root.InnerSimplified)) is not { } roots) @@ -157,7 +165,11 @@ internal static class SumOverSet _ => null, }; - private static IReadOnlyList? NumericRoots(Entity polynomial, Variable w) + /// + /// The roots of in to the working + /// precision, each once, all of them or . + /// + internal static IReadOnlyList? NumericRoots(Entity polynomial, Variable w) { if (SquareFreeParts(polynomial, w) is not { } parts) return null; @@ -186,7 +198,7 @@ internal static class SumOverSet /// through evaluation would refuse the polynomial -- in the setting numerical checks /// run in. /// - private static IReadOnlyList? SquareFreeParts(Entity polynomial, Variable w) + internal static IReadOnlyList? SquareFreeParts(Entity polynomial, Variable w) { if (MultivariatePolynomial.TryParse(polynomial, new Dictionary { [w] = 0 }) is not { } parsed) return null; diff --git a/Sources/Tests/UnitTests/Calculus/SumOverSetTest.cs b/Sources/Tests/UnitTests/Calculus/SumOverSetTest.cs index c3f5d7118..c845783d1 100644 --- a/Sources/Tests/UnitTests/Calculus/SumOverSetTest.cs +++ b/Sources/Tests/UnitTests/Calculus/SumOverSetTest.cs @@ -68,6 +68,9 @@ public void ARepeatedRootIsOneMember() [InlineData("sum(x, x in [0; 1])")] [InlineData("sum(k, k in ZZ)")] [InlineData("sum(w, w in { w : w^5 + w + 3 = 0 })")] + // An irreducible cubic: Cardano would write its roots in complex radicals, which read no + // more simply than the sum. + [InlineData("sum(w^2, w in { w : w^3 + w + 1 = 0 })")] public void WhatItCannotListIsLeftAsWritten(string sum) => Assert.IsType(sum.ToEntity().InnerSimplified); @@ -81,6 +84,8 @@ public void WhatItCannotListIsLeftAsWritten(string sum) [InlineData("sum(w, w in { w : w^5 + w + 3 = 0 })", "0")] [InlineData("sum(1/w, w in { w : w^5 + w + 3 = 0 })", "-1/3")] [InlineData("sum(w^5, w in { w : w^5 + w + 3 = 0 })", "-15")] + // The cubic kept as a sum still has its value: the squares of its roots add to -2. + [InlineData("sum(w^2, w in { w : w^3 + w + 1 = 0 })", "-2")] public void NumericallyOverRootsTheSolverCannotWrite(string sum, string exact) { AssertClose(exact.ToEntity().Evaled, sum.ToEntity().Evaled);