From 2f3b0d1c1ecb9d7f3757aaa9b77862d4f5508e3c Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 23 Oct 2026 22:52:32 +0000 Subject: [PATCH 1/2] Write an expression with floors and ceilings over a bounded range as the piecewise it is there MathS.ExperimentalFeatures.PiecewiseOver(expr, x, from, to) gives expr over [from, to] as a piecewise with each floor and ceiling of x replaced by its value on each piece: the pieces partition the range, their ends are exact, and an end belongs to the piece whose value it has, which the closedness of each piece's interval says. floor(x^2) over [0, 2] is 0 on [0, 1), 1 on [1, sqrt(2)), 2 on [sqrt(2), sqrt(3)), 3 on [sqrt(3), 2) and 4 at 2. The ends are the roots of g(x) = k for each whole k in a step argument g's range there, from the solver, each checked to be one exactly. That they are all the ends is proved: between two consecutive ends g reaches no whole number, by bisection, each part settled by intervals that hold every value there -- g inside one (k, k + 1), or g strictly monotone with no whole number between its values at the part's ends, or, beside an end where g = k, the first derivative that does not vanish there of one sign throughout, so that g - k has one sign beside the end by Taylor's theorem. The last is what the turning points of a sine need: floor(sin(pi x (4 - x)/2)) over [0, 4] reaches 1, -1 and 0 where its argument's derivative is 0. Where the proof does not go through -- a root the solver does not give, an argument with a pole in the range or an infinite slope at an end, a step inside a step, another jump in x -- the answer is null, and never a partition with a piece missing. Experimental, because the issue asks for the contract to be designed and none has been agreed. The interval evaluator takes a variable bound to a range, which every one of its operations already encloses. Part of https://github.com/asc-community/AngouriMath/issues/1807. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- Sources/.editorconfig | 3 + .../Experimental/MathS.Experimental.cs | 26 ++ .../Functions/Continuous/StepPieces.cs | 407 ++++++++++++++++++ .../Numerics/IntervalEvaluation.cs | 18 + .../UnitTests/Calculus/PiecewiseOverTest.cs | 91 ++++ Sources/Tests/UnitTests/Common/PublicApi.txt | 1 + 6 files changed, 546 insertions(+) create mode 100644 Sources/AngouriMath/Functions/Continuous/StepPieces.cs create mode 100644 Sources/Tests/UnitTests/Calculus/PiecewiseOverTest.cs diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 43723eecc..be901a7f0 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -39,6 +39,9 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [AngouriMath/Functions/TreeAnalyzer/PerfectPower.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/StepPieces.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/Common/LogarithmOfAPerfectPowerTest.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/Experimental/MathS.Experimental.cs b/Sources/AngouriMath/Convenience/Experimental/MathS.Experimental.cs index f7795a0de..c45667c64 100644 --- a/Sources/AngouriMath/Convenience/Experimental/MathS.Experimental.cs +++ b/Sources/AngouriMath/Convenience/Experimental/MathS.Experimental.cs @@ -187,6 +187,32 @@ public static Entity ExpandCosineOfSum(IReadOnlyList<(Entity SinX, Entity CosX)> /// public static Entity? SymbolicFormOfCosine(Entity angle) => TrigonometricAngleExpansion.SymbolicFormOfCosine(angle)?.InnerSimplified; + + /// + /// over [, ] as the piecewise it is + /// there, with each floor and ceiling of replaced by its value on each piece: the pieces + /// partition the range, their ends are exact, and an end belongs to the piece whose value it has, which the + /// closedness of each piece's interval says. floor(x^2) over [0, 2] is 0 for + /// x in [0; 1), 1 for x in [1; sqrt(2)), 2 for x in [sqrt(2); sqrt(3)), + /// 3 for x in [sqrt(3); 2), and 4 for x = 2. + /// + /// + /// The ends are the points where a floor's or ceiling's argument is a whole number, from the solver; that + /// there are no others is proved, by intervals that hold every value the argument takes between two ends. + /// Where that cannot be proved -- a root the solver does not give, an argument with a pole in the range, a + /// floor inside another, or any other jump in -- the answer is , + /// and never a partition with a piece missing. + /// https://github.com/asc-community/AngouriMath/issues/1807 + /// + /// The expression, with floors and ceilings of in it + /// The variable + /// The start of the range, a finite real number + /// The end of the range, a finite real number above + /// + /// The piecewise, or where its pieces cannot be proved to be all of them + /// + public static Entity.Piecewise? PiecewiseOver(Entity expr, Entity.Variable x, Entity from, Entity to) + => Functions.Algebra.StepPieces.Of(expr, x, from, to); } } } \ No newline at end of file diff --git a/Sources/AngouriMath/Functions/Continuous/StepPieces.cs b/Sources/AngouriMath/Functions/Continuous/StepPieces.cs new file mode 100644 index 000000000..c91c8e1b6 --- /dev/null +++ b/Sources/AngouriMath/Functions/Continuous/StepPieces.cs @@ -0,0 +1,407 @@ +// +// 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.Collections.Generic; +using System.Linq; +using AngouriMath.Numerics; +using PeterO.Numbers; +using static AngouriMath.Entity; +using static AngouriMath.Entity.Number; + +namespace AngouriMath.Functions.Algebra +{ + /// + /// An expression with floors and ceilings of x in it, written over a bounded range as the piecewise it is + /// there: each floor and ceiling replaced by its value on each piece, the pieces' ends exact, and each end owned by + /// the piece whose value it has. floor(x^2) over [0, 2] is 0 on [0, 1), 1 on + /// [1, sqrt(2)), 2 on [sqrt(2), sqrt(3)), 3 on [sqrt(3), 2) and 4 at 2. + /// https://github.com/asc-community/AngouriMath/issues/1807 + /// + /// + /// + /// The ends are where a step's argument g is a whole number: the roots of g(x) = k for each whole + /// k in g's range there, as the solver gives them, each checked to be one exactly. What makes the + /// answer a partition rather than a list of the roots the solver happened to find is the proof that between two + /// consecutive ends g reaches no whole number: by bisection, each part settled by intervals that hold every + /// value there (see ) -- g inside one (k, k + 1), or g + /// strictly monotone with no whole number between its values at the part's ends, or, next to an end where + /// g = k, the first derivative of g that does not vanish there of one sign throughout, so that by + /// Taylor's theorem g - k has one sign beside the end. A root the solver missed is a part none of these + /// settles, and the whole is declined: , not an answer with a piece missing. + /// + /// + /// The last criterion is for the ends where g only touches a whole number, which are most of them for a + /// sine: floor(sin(pi x (4 - x)/2)) over [0, 4] reaches 1, -1 and 0 at the turning points of its + /// argument, where g' = 0 and no monotone part has the end in it. + /// + /// + internal static class StepPieces + { + /// The most whole numbers one argument may pass through over the range. + private const int MostWholeNumbers = 256; + /// How far either way a solution's integer parameter is tried. + private const int ParameterWindow = 64; + /// The most intervals evaluated to settle one gap. + private const int MostEvaluationsPerGap = 4096; + /// The highest derivative asked at an end where the argument touches a whole number. + private const int HighestDerivative = 4; + + internal static Piecewise? Of(Entity expr, Variable x, Entity from, Entity to) + { + if (from.Evaled is not Real lower || !lower.IsFinite || to.Evaled is not Real upper || !upper.IsFinite || lower >= upper) + return null; + var steps = expr.Nodes.Where(node => node is Floorf or Ceilf && node.ContainsNode(x)).Distinct().ToList(); + // Any other jump in x is one this does not read, and a step inside a step's argument is a jump in it. + if (expr.Nodes.Any(node => node is Piecewise or Signumf or Iversonf or Providedf && node.ContainsNode(x))) + return null; + var arguments = new List(); + foreach (var step in steps) + { + var argument = ArgumentOf(step); + if (argument.Nodes.Any(node => node is Floorf or Ceilf) || argument.Vars.Any(variable => variable != x)) + return null; + if (!arguments.Contains(argument)) + arguments.Add(argument); + } + + // The ends: the range's own, and every point where an argument is a whole number. + var range = new Interval(AsADouble(lower.EDecimal).Low, AsADouble(upper.EDecimal).High); + var ends = new List<(Entity At, EDecimal Value)> { (from, lower.EDecimal), (to, upper.EDecimal) }; + var wholeAt = new List>(); + foreach (var argument in arguments) + { + if (WholeNumbersReached(argument, x, lower, upper, range) is not { } reached) + return null; + wholeAt.Add(reached.ToDictionary(root => root.At, root => root.Whole)); + ends.AddRange(reached.Select(root => (root.At, root.Value))); + } + if (Ordered(ends) is not { } ordered) + return null; + + // Each argument's value at each end and on each gap between two. + var points = new Value[arguments.Count, ordered.Count]; + var gaps = new EInteger[arguments.Count, ordered.Count - 1]; + for (var a = 0; a < arguments.Count; a++) + { + var argument = arguments[a]; + var derivatives = new Derivatives(argument, x); + for (var i = 0; i < ordered.Count; i++) + if (AtAnEnd(argument, x, ordered[i], wholeAt[a]) is { } value) + points[a, i] = value; + else + return null; + for (var i = 0; i + 1 < ordered.Count; i++) + if (OnAGap(argument, x, derivatives, (ordered[i], points[a, i]), (ordered[i + 1], points[a, i + 1])) is { } between) + gaps[a, i] = between; + else + return null; + } + + // The expression on each end and gap in turn, and the runs of one expression joined into one piece. + var regions = new List(); + for (var i = 0; i < ordered.Count; i++) + { + regions.Add(Resolved(expr, steps, arguments, step => points[arguments.IndexOf(ArgumentOf(step)), i].Of(step))); + if (i + 1 < ordered.Count) + regions.Add(Resolved(expr, steps, arguments, step => Value.Between(gaps[arguments.IndexOf(ArgumentOf(step)), i]).Of(step))); + } + var cases = new List(); + for (var start = 0; start < regions.Count;) + { + var end = start; + while (end + 1 < regions.Count && regions[end + 1] == regions[start]) + end++; + // Region 2i is the end i, region 2i + 1 the gap after it. + var (left, leftClosed) = (ordered[start / 2].At, start % 2 == 0); + var (right, rightClosed) = end % 2 == 0 ? (ordered[end / 2].At, true) : (ordered[end / 2 + 1].At, false); + Entity condition = start == end && start % 2 == 0 + ? MathS.Equality(x, left) + : MathS.Sets.ElementInSet(x, MathS.Interval(left, leftClosed, right, rightClosed)); + cases.Add(new Providedf(regions[start], condition)); + start = end + 1; + } + return new Piecewise(cases); + } + + private static Entity ArgumentOf(Entity step) => step switch + { + Floorf(var argument) => argument, + Ceilf(var argument) => argument, + _ => throw new Core.Exceptions.AngouriBugException($"{step} is not a step") + }; + + private static Entity Resolved(Entity expr, List steps, List arguments, System.Func valueOf) + { + var resolved = expr; + foreach (var step in steps) + resolved = resolved.Substitute(step, valueOf(step)); + return resolved.InnerSimplified; + } + + /// + /// An argument's value at a point: a whole number exactly, or strictly between Floor and Floor + 1. + /// + private readonly struct Value + { + internal readonly EInteger Floor; + internal readonly bool IsWhole; + private Value(EInteger floor, bool isWhole) => (Floor, IsWhole) = (floor, isWhole); + internal static Value Whole(EInteger value) => new(value, true); + internal static Value Between(EInteger floor) => new(floor, false); + internal Entity Of(Entity step) => Integer.Create(step is Ceilf && !IsWhole ? Floor.Add(EInteger.One) : Floor); + } + + /// + /// The points of the range where is a whole number, with that number; null where + /// the argument is not real and finite over the range, passes too many whole numbers, or the solver gives + /// something other than a finite set for one of them. + /// + private static List<(Entity At, EDecimal Value, EInteger Whole)>? WholeNumbersReached(Entity argument, Variable x, Real lower, Real upper, Interval range) + { + if (IntervalEvaluation.Over(argument, x, range) is not { IsReal: true } over) + return null; + var least = EInteger.FromInt64((long)System.Math.Floor(over.Re.Low)); + var most = EInteger.FromInt64((long)System.Math.Ceiling(over.Re.High)); + if (most.Subtract(least).CompareTo(EInteger.FromInt32(MostWholeNumbers)) > 0) + return null; + var found = new List<(Entity At, EDecimal Value, EInteger Whole)>(); + for (var k = least; k.CompareTo(most) <= 0; k = k.Add(EInteger.One)) + { + var whole = Integer.Create(k); + if ((argument - whole).SolveEquation(x) is not Set.FiniteSet solutions) + return null; + foreach (var solution in solutions) + { + var parameters = solution.Vars.Where(variable => variable != x && !argument.Vars.Contains(variable)).Distinct().ToList(); + if (parameters.Count > 1) + return null; + var candidates = parameters.Count == 0 + ? new[] { solution } + : Enumerable.Range(-ParameterWindow, 2 * ParameterWindow + 1).Select(n => solution.Substitute(parameters[0], n)); + foreach (var candidate in candidates) + { + // The intervals first, which turn most candidates away for the cost of a few doubles: one + // outside the range, or off the real line. An interval that cannot tell -- the root of a + // discriminant that is exactly zero, at a double root -- leaves it to the decimals. + if (IntervalEvaluation.Of(candidate) is { } enclosure + && (enclosure.IsReal ? enclosure.Re.High < range.Low || enclosure.Re.Low > range.High : !enclosure.Im.ContainsZero)) + continue; + if (candidate.Evaled is not Real { IsFinite: true } at || at < lower || at > upper) + continue; + var root = Bare(candidate.Simplify()); + // A root exactly, or none: its end is owned by the value there. + if (!TreeAnalyzer.IsZero(Bare((argument.Substitute(x, root) - whole).Simplify()))) + return null; + // Two solutions can be one root: sin(u) = 0 is u = 2 pi n and u = pi + 2 pi n. + if (!found.Any(known => known.At == root)) + found.Add((root, at.EDecimal, k)); + } + } + } + return found; + } + + /// + /// The ends in increasing order, each once: two that are one point by their values are one point exactly, or + /// the order is not known and this is null. + /// + private static List<(Entity At, EDecimal Value)>? Ordered(List<(Entity At, EDecimal Value)> ends) + { + var sorted = ends.OrderBy(end => end.Value).ToList(); + var ordered = new List<(Entity At, EDecimal Value)>(); + foreach (var end in sorted) + { + if (ordered.Count > 0 && ordered[^1] is var last && AlikeAt(last.Value, end.Value)) + { + if (last.At == end.At) + continue; + if (!TreeAnalyzer.IsZero(Bare((last.At - end.At).Simplify()))) + return null; + if (end.At.Complexity < last.At.Complexity) + ordered[^1] = end; + continue; + } + ordered.Add(end); + } + return ordered; + + static bool AlikeAt(EDecimal a, EDecimal b) + => a.Subtract(b).Abs().CompareTo(EDecimal.FromString("1E-40").Multiply(EDecimal.Max(EDecimal.One, a.Abs()))) <= 0; + } + + /// + /// 's value at an end: the whole number it was found to be there, or one its + /// interval at the point holds strictly between two, or a whole number the simplification shows it to be. + /// + private static Value? AtAnEnd(Entity argument, Variable x, (Entity At, EDecimal Value) end, Dictionary wholeAt) + { + if (wholeAt.TryGetValue(end.At, out var known)) + return Value.Whole(known); + var there = argument.Substitute(x, end.At); + if (IntervalEvaluation.Of(there) is not { IsReal: true } enclosure) + return null; + var floor = System.Math.Floor(enclosure.Re.Low); + if (enclosure.Re.Low > floor && enclosure.Re.High < floor + 1) + return Value.Between(EInteger.FromInt64((long)floor)); + // An end found for another argument, or the range's own, where this one may be whole too -- or so near a + // whole number that the intervals cannot tell, and neither can this. + if (there.Evaled is not Real { IsFinite: true } value) + return null; + var nearest = value.EDecimal.RoundToIntegerNoRoundedFlag(EContext.CliDecimal.WithRounding(ERounding.HalfEven)).ToEInteger(); + return TreeAnalyzer.IsZero(Bare((there - Integer.Create(nearest)).Simplify())) ? Value.Whole(nearest) : null; + } + + /// The derivatives of an argument, taken as they are asked for. + private sealed class Derivatives + { + private readonly List taken; + private readonly Variable x; + private readonly Dictionary<(int Order, Entity At), bool> vanishing = new(); + internal Derivatives(Entity argument, Variable x) => (taken, this.x) = (new List { argument }, x); + + internal Entity this[int order] + { + get + { + while (taken.Count <= order) + taken.Add(Bare(taken[^1].Differentiate(x)).InnerSimplified); + return taken[order]; + } + } + + /// Whether the derivative of this order is exactly zero at the point, by simplification. + internal bool VanishesAt(int order, Entity at) + { + if (!vanishing.TryGetValue((order, at), out var vanishes)) + vanishing[(order, at)] = vanishes = TreeAnalyzer.IsZero(Bare(this[order].Substitute(x, at).Simplify())); + return vanishes; + } + } + + /// + /// The whole number j with strictly between j and j + 1 on the + /// whole open gap between two consecutive ends, or null where the bisection does not settle every part of it. + /// + private static EInteger? OnAGap(Entity argument, Variable x, Derivatives derivatives, + ((Entity At, EDecimal Value) End, Value Value) left, ((Entity At, EDecimal Value) End, Value Value) right) + { + if (IntervalEvaluation.Of(left.End.At) is not { IsReal: true } leftAt || IntervalEvaluation.Of(right.End.At) is not { IsReal: true } rightAt) + return null; + var budget = MostEvaluationsPerGap; + EInteger? settled = null; + var pending = new Stack<(Side Left, Side Right)>(); + pending.Push((Side.AnEnd(leftAt.Re, left.Value, left.End.At), Side.AnEnd(rightAt.Re, right.Value, right.End.At))); + while (pending.Count > 0) + { + var (l, r) = pending.Pop(); + if ((budget -= 4) < 0 || !(l.Range.High < r.Range.Low)) + return null; + if (Settle(l, r) is { } j) + { + if (settled is { } before && !before.Equals(j)) + return null; + settled = j; + continue; + } + var middle = l.Range.High / 2 + r.Range.Low / 2; + if (!(l.Range.High < middle && middle < r.Range.Low) + || IntervalEvaluation.Over(argument, x, Interval.Exactly(middle)) is not { IsReal: true } atTheMiddle) + return null; + var floor = System.Math.Floor(atTheMiddle.Re.Low); + // A part that meets a whole number at a split point is not settled by splitting it there. + if (!(atTheMiddle.Re.Low > floor && atTheMiddle.Re.High < floor + 1)) + return null; + var inside = Side.Inside(middle, Value.Between(EInteger.FromInt64((long)floor))); + pending.Push((inside, r)); + pending.Push((l, inside)); + } + return settled; + + EInteger? Settle(Side l, Side r) + { + var over = new Interval(l.Range.Low, r.Range.High); + if (IntervalEvaluation.Over(argument, x, over) is not { IsReal: true } values) + return null; + // Inside one (j, j + 1) throughout. + var floor = System.Math.Floor(values.Re.Low); + if (values.Re.Low > floor && values.Re.High < floor + 1) + return EInteger.FromInt64((long)floor); + // Strictly monotone, with no whole number strictly between its values at the two sides. + if (IntervalEvaluation.Over(derivatives[1], x, over) is { IsReal: true } slope && !slope.Re.ContainsZero + && Between(l.Value, r.Value) is { } monotone) + return monotone; + // Beside an end where it is a whole number k, the first derivative that does not vanish there of one + // sign throughout: g - k then has one sign beside the end, by Taylor's theorem. + foreach (var (end, beside, onTheRight) in new[] { (l, r, true), (r, l, false) }) + { + if (end.At is not { } at || !end.Value.IsWhole || beside.At is not null) + continue; + for (var order = 1; order <= HighestDerivative; order++) + { + if (IntervalEvaluation.Over(derivatives[order], x, over) is { IsReal: true } derivative && !derivative.Re.ContainsZero) + { + var above = derivative.Re.Low > 0 == (onTheRight || order % 2 == 0); + var k = end.Value.Floor; + // And no other whole number reached on the part. + if (above ? values.Re.High < (double)k.ToInt64Checked() + 1 : values.Re.Low > (double)k.ToInt64Checked() - 1) + return above ? k : k.Subtract(EInteger.One); + break; + } + if (!derivatives.VanishesAt(order, at)) + break; + } + } + return null; + } + } + + /// + /// The whole number j with every value strictly between the two in (j, j + 1), for a function + /// strictly monotone between two points whose values are those; null where a whole number is strictly between. + /// + private static EInteger? Between(Value a, Value b) + { + if (a.IsWhole && b.IsWhole) + { + var difference = a.Floor.Subtract(b.Floor).Abs(); + return difference.Equals(EInteger.One) ? EInteger.Min(a.Floor, b.Floor) : null; + } + if (a.IsWhole || b.IsWhole) + { + var (whole, other) = a.IsWhole ? (a, b) : (b, a); + return other.Floor.Equals(whole.Floor) || other.Floor.Equals(whole.Floor.Subtract(EInteger.One)) ? other.Floor : null; + } + return a.Floor.Equals(b.Floor) ? a.Floor : null; + } + + /// + /// One side of a part being settled: an end of the gap, with its exact point and value, or a split point inside + /// it, a double, with the argument's value there. + /// + private readonly struct Side + { + internal readonly Interval Range; + internal readonly Value Value; + internal readonly Entity? At; + private Side(Interval range, Value value, Entity? at) => (Range, Value, At) = (range, value, at); + internal static Side AnEnd(Interval range, Value value, Entity at) => new(range, value, at); + internal static Side Inside(double point, Value value) => new(Interval.Exactly(point), value, null); + } + + private static Entity Bare(Entity expr) => expr is Providedf(var inner, _) ? inner : expr; + + /// + /// The double a bound is, exactly, where it is one -- 0, a half, a whole number -- and the interval of one ulp + /// either way otherwise. Widened, the range of sqrt(x) from 0 held a negative and had no value. + /// + private static Interval AsADouble(EDecimal value) + { + var nearest = value.ToDouble(); + return EDecimal.FromDouble(nearest).CompareTo(value) == 0 ? Interval.Exactly(nearest) : Interval.Around(nearest); + } + } +} diff --git a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs index aba50555d..abd6e8ac5 100644 --- a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs +++ b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs @@ -317,6 +317,22 @@ internal static class IntervalEvaluation return value is { IsFinite: true } finite ? finite : null; } + /// + /// An interval that holds every value takes for in + /// , or null as has it. Every operation is an enclosure over + /// an interval as well as at a point -- the sine and cosine by their slope about the middle -- so + /// the result is one, if a loose one: x x over [-1, 2] is [-2, 4]. + /// + internal static ComplexInterval? Over(Entity expr, Variable x, Interval range) + { + var (outerVariable, outerRange) = (boundVariable, boundRange); + (boundVariable, boundRange) = (x, range); + try { return Of(expr); } + finally { (boundVariable, boundRange) = (outerVariable, outerRange); } + } + [System.ThreadStatic] private static Variable? boundVariable; + [System.ThreadStatic] private static Interval boundRange; + /// /// Whether the two agree to of the larger's size (or /// of one, where both are smaller): true where the whole of their difference is within @@ -396,6 +412,8 @@ private static ComplexInterval EvaluateOtherThanArithmetic(Entity expr) return ComplexInterval.Real(Interval.Around(Math.PI)); case Variable variable when variable == MathS.e: return ComplexInterval.Real(Interval.Around(Math.E)); + case Variable variable when boundVariable is not null && variable == boundVariable: + return ComplexInterval.Real(boundRange); case Powf(var @base, var exponent): { if (@base == MathS.e) diff --git a/Sources/Tests/UnitTests/Calculus/PiecewiseOverTest.cs b/Sources/Tests/UnitTests/Calculus/PiecewiseOverTest.cs new file mode 100644 index 000000000..b37a6026b --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/PiecewiseOverTest.cs @@ -0,0 +1,91 @@ +// +// 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; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// An expression with floors and ceilings in it over a bounded range, as the piecewise it is there: the pieces + /// partition the range, their ends are exact, and an end belongs to the piece whose value it has -- or null, where + /// that cannot be proved. The expressions are the ones + /// #1807 names. + /// + [Trait("Area", "Calculus")] + public sealed class PiecewiseOverTest + { + private static Entity.Piecewise? Pieces(string expr, string from, string to) + => MathS.ExperimentalFeatures.PiecewiseOver(expr.ToEntity(), "x", from.ToEntity(), to.ToEntity()); + + /// Each piece as value on interval, with the interval as the library writes it. + private static string Written(Entity.Piecewise pieces) + => string.Join(", ", pieces.Cases.Select(@case => $"{@case.Expression.Stringize()} for {@case.Predicate.Stringize()}")); + + [Theory] + [InlineData("floor(x^2)", "0", "2", + "0 for x in [0; 1), 1 for x in [1; sqrt(2)), 2 for x in [sqrt(2); sqrt(3)), 3 for x in [sqrt(3); 2), 4 for x = 2")] + [InlineData("floor(x^2/2)", "0", "3", + "0 for x in [0; sqrt(2)), 1 for x in [sqrt(2); 2), 2 for x in [2; sqrt(6)), 3 for x in [sqrt(6); 2 * sqrt(2)), 4 for x in [2 * sqrt(2); 3]")] + // The argument touches 1, -1 and 0 at its turning points, where its derivative is 0: the value there is the + // one the point has, and the pieces either side of it are one value. + [InlineData("floor(sin(pi*x*(4 - x)/2))", "0", "4", + "0 for x in [0; 2 - sqrt(3)), 1 for x = 2 - sqrt(3), 0 for x in (2 - sqrt(3); 2 - sqrt(2)], -1 for x in (2 - sqrt(2); 2), " + + "0 for x = 2, -1 for x in (2; 2 + sqrt(2)), 0 for x in [2 + sqrt(2); 2 + sqrt(3)), 1 for x = 2 + sqrt(3), 0 for x in (2 + sqrt(3); 4]")] + [InlineData("ceil(sin(pi*x*(4 - x)/2))", "0", "4", + "0 for x = 0, 1 for x in (0; 2 - sqrt(2)), 0 for x in [2 - sqrt(2); 1), -1 for x = 1, 0 for x in (1; 3), -1 for x = 3, " + + "0 for x in (3; 2 + sqrt(2)], 1 for x in (2 + sqrt(2); 4), 0 for x = 4")] + [InlineData("floor(sin(pi*x/2))", "0", "4", + "0 for x in [0; 1), 1 for x = 1, 0 for x in (1; 2], -1 for x in (2; 4), 0 for x = 4")] + // Two steps, the pieces their common refinement, and the runs of one value joined. + [InlineData("floor(x^2) + ceil(x)", "-1", "2", + "0 for x in [-1; 0], 1 for x in (0; 1), 2 for x = 1, 3 for x in (1; sqrt(2)), 4 for x in [sqrt(2); sqrt(3)), 5 for x in [sqrt(3); 2), 6 for x = 2")] + // What is beside the step is carried into each piece. + [InlineData("x*floor(x^2)", "0", "2", + "0 for x in [0; 1), x for x in [1; sqrt(2)), x * 2 for x in [sqrt(2); sqrt(3)), x * 3 for x in [sqrt(3); 2), x * 4 for x = 2")] + [InlineData("floor(e^x)", "0", "1", "1 for x in [0; ln(2)), 2 for x in [ln(2); 1]")] + public void IsThePiecewiseItIsOverTheRange(string expr, string from, string to, string expected) + { + var pieces = Pieces(expr, from, to); + Assert.NotNull(pieces); + Assert.Equal(expected, Written(pieces!)); + } + + /// + /// The piecewise has the expression's value at every point of the range, a point inside a piece and each end + /// alike: exact substitution, which is what floor(sin(pi x/2)) at 4 got wrong in 2.5.0. + /// + [Theory] + [InlineData("floor(sin(pi*x*(4 - x)/2))", "0", "4", new[] { "0", "1/10", "2 - sqrt(3)", "1/2", "1", "2", "5/2", "3", "2 + sqrt(2)", "4" })] + [InlineData("floor(x^2) + ceil(x)", "-1", "2", new[] { "-1", "-1/2", "0", "1/2", "1", "sqrt(2)", "3/2", "sqrt(3)", "2" })] + public void HasTheExpressionsValueAtEachPoint(string expr, string from, string to, string[] points) + { + var pieces = Pieces(expr, from, to); + Assert.NotNull(pieces); + foreach (var point in points) + { + var want = expr.ToEntity().Substitute("x", point.ToEntity()).Simplify(); + var got = pieces!.Substitute("x", point.ToEntity()).Simplify(); + Assert.Equal(want, got); + } + } + + /// + /// Declined, not answered with a piece missing: a pole of the argument in the range, a floor inside a floor, a + /// range that is not one, and an argument with another symbol in it. + /// + [Theory] + [InlineData("floor(1/x)", "-1", "1")] + [InlineData("floor(floor(x)/2 + x)", "0", "3")] + [InlineData("floor(x^2)", "2", "0")] + [InlineData("floor(x^2)", "0", "a")] + [InlineData("floor(a*x)", "0", "1")] + public void IsNullWhereThePiecesCannotBeProvedComplete(string expr, string from, string to) + => Assert.Null(Pieces(expr, from, to)); + } +} diff --git a/Sources/Tests/UnitTests/Common/PublicApi.txt b/Sources/Tests/UnitTests/Common/PublicApi.txt index f42b990c5..eefad7424 100644 --- a/Sources/Tests/UnitTests/Common/PublicApi.txt +++ b/Sources/Tests/UnitTests/Common/PublicApi.txt @@ -3092,6 +3092,7 @@ AngouriMath.MathS+ExperimentalFeatures.ExpandSineArgumentMultiplied(AngouriMath. AngouriMath.MathS+ExperimentalFeatures.ExpandSineOfSum(System.Collections.Generic.IReadOnlyList>) : AngouriMath.Entity AngouriMath.MathS+ExperimentalFeatures.GetCosineOfHalvedAngle(AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity AngouriMath.MathS+ExperimentalFeatures.GetSineOfHalvedAngle(AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity +AngouriMath.MathS+ExperimentalFeatures.PiecewiseOver(AngouriMath.Entity, AngouriMath.Entity+Variable, AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity+Piecewise AngouriMath.MathS+ExperimentalFeatures.SolveDiophantineEquation(AngouriMath.Entity+Number+Integer, AngouriMath.Entity+Number+Integer, AngouriMath.Entity+Number+Integer) : System.Nullable> AngouriMath.MathS+ExperimentalFeatures.SymbolicFormOfCosine(AngouriMath.Entity) : AngouriMath.Entity AngouriMath.MathS+ExperimentalFeatures.SymbolicFormOfSine(AngouriMath.Entity) : AngouriMath.Entity From 783521c5b52717659bfdec67c5551d400d65878a Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 23 Oct 2026 22:57:47 +0000 Subject: [PATCH 2/2] Integrate a floor or a ceiling of something in x over its pieces integral(floor(x^2), x, 0, 2) stayed as written: a floor or a ceiling was a break only where its argument was x itself. Between two numeric bounds it is now the sum of the integrals of the pieces StepPieces proves, each with the step a constant, which is 5 - sqrt(2) - sqrt(3); a piece that is a point adds nothing. Where the pieces are not proved complete, or a bound is symbolic, the integral is left as written as before. WhatCannotBeSplitIsLeftAsWritten pinned floor(x^2) from 0 to 2 as one that could not be split, which is what changed. Part of https://github.com/asc-community/AngouriMath/issues/1807. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 17 +++++++ .../Integration/BreakpointIntegration.cs | 46 +++++++++++++++++-- .../Calculus/BreakpointIntegrationTest.cs | 35 ++++++++++++-- 3 files changed, 92 insertions(+), 6 deletions(-) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index cc161e3d1..2be02df9b 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -110,6 +110,23 @@ as `a - b = 0`, as a statement's own solve reads it: `"x + 1 = 3".Solve("x")` wa | `MathS.Equations("x + y = 3", "x - y = 1").Solve("x", "y")` | `null` | `[[2, 1]]` | | `MathS.Equations("x + y = 3", "x - y = 1").ToString()` | `x + y = 3 = 0` and `x - y = 1 = 0` | `x + y - 3 = 0` and `x - y - 1 = 0` | +### A definite integral of a floor or a ceiling of something in x is the sum over its pieces + +**Silent: values where the integral was left as written.** A floor or a ceiling was split only where its argument was +`x` itself, and `integral(floor(x^2), x, 0, 2)` stayed as written. Between two bounds that are numbers it is now the +sum of the integrals of the pieces `MathS.ExperimentalFeatures.PiecewiseOver` writes the integrand in, each with the +step a constant there: the points where the argument is a whole number, from the solver, with a proof that there are +no others. Where that proof does not go through, or a bound is symbolic, the integral is left as written, as before +([#1807](https://github.com/asc-community/AngouriMath/issues/1807)). + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `integral(floor(x^2), x, 0, 2)` | `integral(floor(x ^ 2), x, 0, 2)` | `5 - sqrt(2) - sqrt(3)` | +| `integral(floor(2*x), x, 0, 3)` | `integral(floor(2 * x), x, 0, 3)` | `15/2` | +| `integral(floor(sin(pi*x*(4-x)/2)), x, 0, 4)` | as written | `(-2) * sqrt(2)` | +| `integral(floor(1/x), x, 1/4, 1)` | as written | `13/12` | +| `integral(floor(1/x), x, -1, 1)` | as written | as written: a pole inside the range | + ### The polygamma functions, and the derivative of the factorial **Silent: a new reading, and an answer where there was none.** diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/BreakpointIntegration.cs b/Sources/AngouriMath/Functions/Continuous/Integration/BreakpointIntegration.cs index 144a349f7..720214237 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/BreakpointIntegration.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/BreakpointIntegration.cs @@ -54,6 +54,12 @@ namespace AngouriMath.Functions.Algebra /// (#1174). /// /// + /// Where something in the variable is a whole number -- a floor or a ceiling of + /// x^2, of sin(pi x/2) -- the breaks are those points, which + /// finds and proves to be all of them over a range with numeric + /// bounds; each piece is then integrated with its step a constant. + /// + /// /// Offered only where every piece then resolves; an integral rewritten as an unevaluated sum /// of unevaluated integrals is not an answer. Question I.2 of /// #1212, and the @@ -65,7 +71,7 @@ internal static class BreakpointIntegration { /// /// Whether the integrand can jump in : a piecewise whose conditions - /// mention it, or a floor or ceiling of it. + /// mention it, or a floor or ceiling of something in it. /// internal static bool HasABreak(Entity expr, Variable x) => expr.Nodes.Any(node => IsABreak(node, x)); @@ -73,8 +79,8 @@ internal static bool HasABreak(Entity expr, Variable x) private static bool IsABreak(Entity node, Variable x) => node switch { - Floorf(var argument) => argument == x, - Ceilf(var argument) => argument == x, + Floorf(var argument) => argument.ContainsNode(x), + Ceilf(var argument) => argument.ContainsNode(x), Piecewise piecewise => piecewise.Cases.Any(c => c.Predicate.ContainsNode(x)), // A condition on the variable is a one-case piecewise whose other case is undefined. Providedf provided => provided.Predicate.ContainsNode(x), @@ -85,11 +91,45 @@ private static bool IsABreak(Entity node, Variable x) { if (expr.Nodes.Any(node => node is Piecewise or Providedf && IsABreak(node, x))) return OverCaseBoundaries(expr, x, from, to); + if (expr.Nodes.Any(node => node is Floorf(var f) && f != x && f.ContainsNode(x) || node is Ceilf(var c) && c != x && c.ContainsNode(x))) + return OverStepPieces(expr, x, from, to); if (HasAStep(expr, x)) return OverUnitIntervals(expr, x, from, to); return null; } + // ---- a floor or a ceiling of something in the variable: the pieces StepPieces proves ------- + + /// + /// The integral as the sum of the integrals of the pieces writes the integrand in over + /// the range, each with its steps replaced by their values there, so with no jump in it: floor(x^2) from 0 + /// to 2 is (sqrt(2) - 1) + 2 (sqrt(3) - sqrt(2)) + 3 (2 - sqrt(3)), which is 5 - sqrt(2) - sqrt(3). + /// A piece that is a point adds nothing. Where the pieces are not proved complete, or a piece does not integrate, + /// the integral is left as written. https://github.com/asc-community/AngouriMath/issues/1807 + /// + private static Entity? OverStepPieces(Entity expr, Variable x, Entity from, Entity to) + { + if (from.Evaled is not Real lower || !lower.IsFinite || to.Evaled is not Real upper || !upper.IsFinite) + return null; + if (lower == upper) + return Integer.Zero; + if (upper < lower) + return OverStepPieces(expr, x, to, from) is { } reversed ? (-reversed).InnerSimplified : null; + if (StepPieces.Of(expr, x, from, to) is not { } pieces) + return null; + Entity total = Integer.Zero; + foreach (var @case in pieces.Cases) + { + if (@case.Predicate is not Set.Inf(_, Set.Interval interval)) + continue; + var piece = @case.Expression.Integrate(x, interval.Left, interval.Right); + if (piece.Nodes.Any(node => node is Integralf)) + return null; + total += piece; + } + return total.InnerSimplified; + } + // ---- finitely many breakpoints: a piecewise -------------------------------------------- private static Entity? OverCaseBoundaries(Entity expr, Variable x, Entity from, Entity to) diff --git a/Sources/Tests/UnitTests/Calculus/BreakpointIntegrationTest.cs b/Sources/Tests/UnitTests/Calculus/BreakpointIntegrationTest.cs index 85e8febd5..85e76b572 100644 --- a/Sources/Tests/UnitTests/Calculus/BreakpointIntegrationTest.cs +++ b/Sources/Tests/UnitTests/Calculus/BreakpointIntegrationTest.cs @@ -93,6 +93,35 @@ public void AnIntegralOverWholeBoundsIsSummedOverUnitIntervals(string integral, public void ABoundThatIsNotWholeContributesItsOwnPiece(string integral, string expected) => Assert.Equal(expected.ToEntity(), integral.ToEntity().Simplify()); + /// + /// A floor or a ceiling of something in x, integrated over the pieces + /// proves, each with the step a constant: + /// floor(x^2) from 0 to 2 is (sqrt(2) - 1) + 2 (sqrt(3) - sqrt(2)) + 3 (2 - sqrt(3)). + /// #1807 + /// + [Theory] + [InlineData("integral(floor(x^2), x, 0, 2)", "5 - sqrt(2) - sqrt(3)")] + [InlineData("integral(floor(2x), x, 0, 3)", "15/2")] + [InlineData("integral(x * floor(x^2), x, 0, 2)", "3")] + [InlineData("integral(floor(sin(pi * x * (4 - x)/2)), x, 0, 4)", "-2 sqrt(2)")] + [InlineData("integral(ceil(3x - 1), x, 0, 1)", "1")] + [InlineData("integral(floor(1/x), x, 1/4, 1)", "13/12")] + [InlineData("integral(floor(x^2), x, 2, 0)", "sqrt(2) + sqrt(3) - 5")] + public void AStepOfSomethingInXIsIntegratedOverItsPieces(string integral, string expected) + { + var result = integral.ToEntity().Simplify(); + Assert.DoesNotContain(result.Nodes, node => node is Entity.Integralf); + var difference = (result - expected.ToEntity()).EvalNumerical(); + Assert.True(difference.Abs() < 1e-30, $"{integral} is {result}, not {expected}"); + } + + /// A symbolic bound, and a pole of the step's argument inside the range: left as written. + [Theory] + [InlineData("integral(floor(x^2), x, 0, a)")] + [InlineData("integral(floor(1/x), x, -1, 1)")] + public void AStepWhosePiecesAreNotProvedIsLeftAsWritten(string integral) + => Assert.IsType(integral.ToEntity().Simplify()); + // A symbolic bound is declined rather than integrated as if the step were constant, // which used to answer (n - floor(n))^3 / 3 for this one; with the bound a number the // same integral splits. @@ -116,10 +145,10 @@ public void TheSumTheSplitLeavesIsWhatAnswersTheIntegral() "integral(floor(x)^floor(x), x, 1, +oo)".ToEntity().Simplify()); } - // Not split and not integrated either: a floor of something other than the variable - // is not a break this reads, and an integrand the unit interval cannot integrate stays. + // Split and not integrated: an integrand the unit interval cannot integrate stays. A floor of + // something other than the variable was here too, floor(x^2) from 0 to 2, and is now split + // at its pieces: see AStepOfSomethingInXIsIntegratedOverItsPieces. [Theory] - [InlineData("integral(floor(x^2), x, 0, 2)")] [InlineData("integral(floor(x) * e^(x^3 + x) , x, 0, 2)")] public void WhatCannotBeSplitIsLeftAsWritten(string integral) => Assert.IsType(integral.ToEntity().Simplify());