From a918d8f5674552bf9a0909abe1d6244532567c61 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Sat, 26 Sep 2026 23:12:03 +0000 Subject: [PATCH 1/9] The integrator's sampled-point checks evaluate in intervals, which read no setting HoldsAtSampledPoints and IsZeroAtPinnedSymbols compared hundred-digit decimals from EvalNumerical under DowncastingEnabled off. They now evaluate in intervals: doubles with outward rounding first, then decimal intervals of 40 and 80 digits where cancellation leaves the doubles too wide to tell. The accuracy is an argument and no global setting is read (#1019, item 10). The evaluator is internal, in AngouriMath.Numerics. The one scope left is around the simplification of the pinned answer in DerivativeHoldsAtSampledPoints, which folds number nodes: downcast, a pinned 1.37 is 137/100 and the answer to u^5/((u^2 - 1)^5 (2 a u + b (u^2 + 1))) took 28 s to simplify, against 3 ms in decimals. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- Sources/.editorconfig | 6 + .../Algebra/Polynomials/PartialFractions.cs | 76 +-- .../Numerics/IntervalEvaluation.cs | 523 ++++++++++++++++++ .../AngouriMath/Numerics/PreciseEvaluation.cs | 482 ++++++++++++++++ .../UnitTests/Core/IntervalEvaluationTest.cs | 131 +++++ 5 files changed, 1185 insertions(+), 33 deletions(-) create mode 100644 Sources/AngouriMath/Numerics/IntervalEvaluation.cs create mode 100644 Sources/AngouriMath/Numerics/PreciseEvaluation.cs create mode 100644 Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 98e221048..92187b349 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -606,3 +606,9 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Calculus/InverseHyperbolicSubstitutionTest.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/Numerics/*.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/IntervalEvaluationTest.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/Functions/Algebra/Polynomials/PartialFractions.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs index f05c9514f..84a4229b0 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs @@ -5,6 +5,7 @@ // Website: https://am.angouri.org. // +using AngouriMath.Numerics; using PeterO.Numbers; using System.Diagnostics.CodeAnalysis; using static AngouriMath.Entity; @@ -730,13 +731,12 @@ private static Entity WithThePrimitiveDenominator(Entity quotient) } /// - /// Whether , free of every variable but symbols, is zero at - /// two sets of values for them, off the integers and distinct; a value that is not - /// a number at either counts as not zero. + /// Whether , free of every variable but symbols, is zero -- within + /// 1e-30 -- at two sets of values for them, off the integers and distinct; a value + /// that is not a number at either, or that the evaluation cannot place, counts as not zero. /// private static bool IsZeroAtPinnedSymbols(Entity expr) { - using var _ = MathS.Settings.DowncastingEnabled.Set(false); var symbols = expr.Vars.ToList(); foreach (var seed in new[] { 0, 1 }) { @@ -748,7 +748,7 @@ private static bool IsZeroAtPinnedSymbols(Entity expr) pinned = pinned.Substitute(symbol, value); i++; } - if (pinned.Evaled is not Complex value_ || value_.IsNaN || value_.Abs().EDecimal.CompareTo(EDecimal.FromString("1e-30")) > 0) + if (Agree(pinned, Integer.Zero, 1e-30) is not true) return false; } return true; @@ -1864,6 +1864,10 @@ internal static Entity Bare(Entity e) /// internal static bool DerivativeHoldsAtSampledPoints(Entity antiderivative, Entity integrand, Variable x) { + // Off the downcasting for the simplification of the pinned answer, which folds its + // numbers: downcast, a pinned 1.37 is 137/100 and every constant of the answer is + // folded exactly -- 28 s for the answer to u^5/((u^2 - 1)^5 (2 a u + b (u^2 + 1))), + // where in decimals it is 3 ms. The comparison itself reads no setting. using var _ = MathS.Settings.DowncastingEnabled.Set(false); var (pinnedAntiderivative, pinnedIntegrand) = Pinned(antiderivative, integrand, x); return HoldsAtSampledPoints(Bare(pinnedAntiderivative).Differentiate(x), pinnedIntegrand, x); @@ -1892,46 +1896,52 @@ private static (Entity, Entity) Pinned(Entity left, Entity right, Variable x) /// /// Whether and agree, numerically, at a /// few points in with every other symbol pinned to a fixed value. - /// A point where either is undefined, or cannot be evaluated -- a derivative left - /// unevaluated, of a sign whose argument is not shown real -- is skipped, and at least - /// two must compare. The points may be the caller's, for an identity that holds on a - /// real domain only -- sqrt(1 - L^2) = sech(artanh(L)) inside (-1, 1) and - /// off by a sign outside it -- where the default set would compare where it does not - /// hold. + /// A point where either has no value -- a pole, a cut, a condition that fails -- or where + /// the evaluation cannot tell is skipped, and at least two must compare. The points may + /// be the caller's, for an identity that holds on a real domain only -- + /// sqrt(1 - L^2) = sech(artanh(L)) inside (-1, 1) and off by a sign outside + /// it -- where the default set would compare where it does not hold. /// + /// + /// In intervals, which read no setting (https://github.com/asc-community/AngouriMath/issues/1019, + /// item 10): doubles first, then forty and eighty digits where cancellation leaves the + /// doubles too wide to tell. + /// internal static bool HoldsAtSampledPoints(Entity left, Entity right, Variable x, string[]? points = null) { - // In decimals: a pinned value is a small rational, which the downcasting keeps - // exact through every operation, and a symbolic answer with powers in it - // evaluated so was a minute of gcd on integers of thousands of digits -- Rubi's - // `(A + B ln(e ((a + b x)/(c + d x))^n))/(a + b x)^3` never returned from its - // check. Off the downcasting, a decimal stays a decimal of a hundred digits. - using var _ = MathS.Settings.DowncastingEnabled.Set(false); (left, right) = Pinned(left, right, x); var compared = 0; foreach (var at in points ?? new[] { "0.29", "1.43", "3.17", "-0.61" }) { var point = Real.Create(EDecimal.FromString(at)); - Number.Complex l; - Number.Complex r; - try + switch (Agree(left.Substitute(x, point), right.Substitute(x, point), 1e-9)) { - l = left.Substitute(x, point).EvalNumerical(); - r = right.Substitute(x, point).EvalNumerical(); - } - catch (Core.Exceptions.CannotEvalException) - { - continue; + case false: + return false; + case true: + compared++; + break; } - if (l.IsNaN || r.IsNaN) - continue; - var difference = (l - r).Abs().EDecimal; - var scale = EDecimal.Max(EDecimal.One, l.Abs().EDecimal); - if (difference.CompareTo(scale.Multiply(EDecimal.FromString("1e-9"))) > 0) - return false; - compared++; } return compared >= 2; } + + /// + /// Whether the two numbers agree to of their size, or + /// of one where both are smaller: in double intervals, and where those cannot tell, in + /// decimal intervals of forty and then eighty digits; null where either has no value or + /// none can tell. + /// + private static bool? Agree(Entity left, Entity right, double relativeTolerance) + { + if (IntervalEvaluation.Of(left) is { } leftInterval && IntervalEvaluation.Of(right) is { } rightInterval + && IntervalEvaluation.Agree(leftInterval, rightInterval, relativeTolerance) is { } inDoubles) + return inDoubles; + foreach (var digits in new[] { 40, 80 }) + if (PreciseEvaluation.Of(left, digits) is { } leftPrecise && PreciseEvaluation.Of(right, digits) is { } rightPrecise + && PreciseEvaluation.Agree(leftPrecise, rightPrecise, relativeTolerance) is { } inDecimals) + return inDecimals; + return null; + } } } diff --git a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs new file mode 100644 index 000000000..0c57405d6 --- /dev/null +++ b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs @@ -0,0 +1,523 @@ +// +// 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; + +namespace AngouriMath.Numerics +{ + /// + /// A closed interval of doubles, [Low, High], that holds a real value: every operation + /// rounds its bounds outward, so the exact value of what was evaluated is inside. + /// + internal readonly struct Interval + { + internal readonly double Low; + internal readonly double High; + + internal Interval(double low, double high) + { + Low = low; + High = high; + } + + /// The one value, exactly. + internal static Interval Exactly(double value) => new(value, value); + + /// A value known to within a unit in the last place either way. + internal static Interval Around(double value) => new(Down(value), Up(value)); + + internal bool IsFinite => IsFiniteDouble(Low) && IsFiniteDouble(High); + internal bool ContainsZero => Low <= 0 && High >= 0; + internal double Magnitude => Math.Max(Math.Abs(Low), Math.Abs(High)); + /// The distance of the interval from zero: zero where it contains zero. + internal double Mignitude => ContainsZero ? 0 : Math.Min(Math.Abs(Low), Math.Abs(High)); + internal double Middle => Low / 2 + High / 2; + + /// Exactly zero, which arithmetic keeps exact: the imaginary part of a real value stays zero. + internal bool IsZero => Low == 0 && High == 0; + + public static Interval operator +(Interval a, Interval b) + => a.IsZero ? b : b.IsZero ? a : new(Down(a.Low + b.Low), Up(a.High + b.High)); + public static Interval operator -(Interval a, Interval b) + => b.IsZero ? a : a.IsZero ? -b : new(Down(a.Low - b.High), Up(a.High - b.Low)); + public static Interval operator -(Interval a) => new(-a.High, -a.Low); + + public static Interval operator *(Interval a, Interval b) + { + if (a.IsZero || b.IsZero) + return Exactly(0); + var p1 = a.Low * b.Low; + var p2 = a.Low * b.High; + var p3 = a.High * b.Low; + var p4 = a.High * b.High; + return new(Down(Math.Min(Math.Min(p1, p2), Math.Min(p3, p4))), Up(Math.Max(Math.Max(p1, p2), Math.Max(p3, p4)))); + } + + /// The quotient, or an interval that is not finite where the divisor may be zero. + public static Interval operator /(Interval a, Interval b) + { + if (a.IsZero && !b.ContainsZero) + return Exactly(0); + if (b.ContainsZero) + return new(double.NegativeInfinity, double.PositiveInfinity); + var q1 = a.Low / b.Low; + var q2 = a.Low / b.High; + var q3 = a.High / b.Low; + var q4 = a.High / b.High; + return new(Down(Math.Min(Math.Min(q1, q2), Math.Min(q3, q4))), Up(Math.Max(Math.Max(q1, q2), Math.Max(q3, q4)))); + } + + internal Interval Square() + { + if (ContainsZero) + return new(0, Up(Math.Max(Low * Low, High * High))); + var (small, large) = Low > 0 ? (Low, High) : (-High, -Low); + return new(Down(small * small), Up(large * large)); + } + + internal Interval Exp() => new(Wider(Math.Exp(Low), -1), Wider(Math.Exp(High), 1)); + + /// The natural logarithm of a positive interval; not finite otherwise. + internal Interval Log() + => Low > 0 ? new(Wider(Math.Log(Low), -1), Wider(Math.Log(High), 1)) : new(double.NaN, double.NaN); + + internal Interval Sqrt() + => Low >= 0 ? new(Math.Max(0, Wider(Math.Sqrt(Low), -1)), Wider(Math.Sqrt(High), 1)) : new(double.NaN, double.NaN); + + // The sine and cosine have a derivative no larger than one, so over an interval of + // half-width r they are within r of their value at the middle. + internal Interval Sin() => AroundTheMiddle(Math.Sin, 1); + internal Interval Cos() => AroundTheMiddle(Math.Cos, 1); + + internal Interval Atan() => new(Wider(Math.Atan(Low), -1), Wider(Math.Atan(High), 1)); + + internal Interval Asin() + => Low >= -1 && High <= 1 ? new(Wider(Math.Asin(Low), -1), Wider(Math.Asin(High), 1)) : new(double.NaN, double.NaN); + + internal Interval Acos() + => Low >= -1 && High <= 1 ? new(Wider(Math.Acos(High), -1), Wider(Math.Acos(Low), 1)) : new(double.NaN, double.NaN); + + internal Interval Abs() + => ContainsZero ? new(0, Magnitude) : Low > 0 ? this : -this; + + private Interval AroundTheMiddle(Func function, double slope) + { + var middle = Middle; + var value = function(middle); + var reach = Up(Math.Max(High - middle, middle - Low) * slope); + return new(Math.Max(-1, Wider(value - reach, -1)), Math.Min(1, Wider(value + reach, 1))); + } + + /// + /// A library function's result moved four units in the last place outward, and a little + /// more near zero: the platform's functions are within one of the true value, and the + /// slack keeps the claim true where they are not quite. + /// + private static double Wider(double value, int direction) + { + for (var i = 0; i < 4; i++) + value = direction > 0 ? Up(value) : Down(value); + return value + direction * 1e-300; + } + + internal static bool IsFiniteDouble(double value) => !double.IsNaN(value) && !double.IsInfinity(value); + + internal static double Up(double value) + { + if (double.IsNaN(value) || double.IsPositiveInfinity(value)) + return value; + if (value == 0) + return double.Epsilon; + var bits = BitConverter.DoubleToInt64Bits(value); + return BitConverter.Int64BitsToDouble(value > 0 ? bits + 1 : bits - 1); + } + + internal static double Down(double value) => -Up(-value); + } + + /// A rectangle of two s holding a complex value. + internal readonly struct ComplexInterval + { + internal readonly Interval Re; + internal readonly Interval Im; + + internal ComplexInterval(Interval re, Interval im) + { + Re = re; + Im = im; + } + + internal static ComplexInterval Real(Interval re) => new(re, Interval.Exactly(0)); + + internal bool IsFinite => Re.IsFinite && Im.IsFinite; + internal bool IsReal => Im.Low == 0 && Im.High == 0; + internal double Magnitude => Math.Sqrt(Re.Magnitude * Re.Magnitude + Im.Magnitude * Im.Magnitude); + + public static ComplexInterval operator +(ComplexInterval a, ComplexInterval b) => new(a.Re + b.Re, a.Im + b.Im); + public static ComplexInterval operator -(ComplexInterval a, ComplexInterval b) => new(a.Re - b.Re, a.Im - b.Im); + public static ComplexInterval operator -(ComplexInterval a) => new(-a.Re, -a.Im); + + public static ComplexInterval operator *(ComplexInterval a, ComplexInterval b) + => a.IsReal && b.IsReal ? Real(a.Re * b.Re) : new(a.Re * b.Re - a.Im * b.Im, a.Re * b.Im + a.Im * b.Re); + + public static ComplexInterval operator /(ComplexInterval a, ComplexInterval b) + { + if (b.IsReal) + return new(a.Re / b.Re, a.IsReal ? Interval.Exactly(0) : a.Im / b.Re); + var denominator = b.Re.Square() + b.Im.Square(); + return new((a.Re * b.Re + a.Im * b.Im) / denominator, (a.Im * b.Re - a.Re * b.Im) / denominator); + } + + internal ComplexInterval Exp() + { + var modulus = Re.Exp(); + return IsReal ? Real(modulus) : new(modulus * Im.Cos(), modulus * Im.Sin()); + } + + /// + /// The principal logarithm, where the rectangle neither holds zero nor straddles the + /// negative real axis, the cut; not finite otherwise. Off the cut the argument is + /// continuous on the rectangle and takes its extremes at the corners. + /// + internal ComplexInterval Log() + { + if (IsReal && Re.Low > 0) + return Real(Re.Log()); + // Exactly on the negative real axis the principal value is ln|x| + i pi, the cut's + // upper side, as the decimal evaluation has it: sqrt(-3) is i sqrt(3). Only a + // rectangle that straddles the cut is undecided. + if (IsReal && Re.High < 0) + return new((-Re).Log(), Interval.Around(Math.PI)); + if (Re.High <= 0 && Im.ContainsZero || Re.ContainsZero && Im.ContainsZero) + return new(new(double.NaN, double.NaN), new(double.NaN, double.NaN)); + var modulusSquared = Re.Square() + Im.Square(); + var logModulus = modulusSquared.Log(); + var a1 = Math.Atan2(Im.Low, Re.Low); + var a2 = Math.Atan2(Im.Low, Re.High); + var a3 = Math.Atan2(Im.High, Re.Low); + var a4 = Math.Atan2(Im.High, Re.High); + var argument = new Interval( + Interval.Down(Interval.Down(Math.Min(Math.Min(a1, a2), Math.Min(a3, a4)))), + Interval.Up(Interval.Up(Math.Max(Math.Max(a1, a2), Math.Max(a3, a4))))); + return new(new(logModulus.Low / 2, logModulus.High / 2), argument); + } + + /// A whole power by repeated squaring, a negative one as the reciprocal. + internal ComplexInterval Pow(long exponent) + { + if (exponent < 0) + return Real(Interval.Exactly(1)) / Pow(-exponent); + var result = Real(Interval.Exactly(1)); + var power = this; + while (exponent > 0) + { + if ((exponent & 1) == 1) + result *= power; + exponent >>= 1; + if (exponent > 0) + power *= power; + } + return result; + } + + /// The principal power, exp(w log z); not finite where the logarithm is not. + internal ComplexInterval Pow(ComplexInterval exponent) + { + // A real square root of a real that is not negative, without the logarithm's detour. + if (IsReal && Re.Low >= 0 && exponent.IsReal && exponent.Re.Low == 0.5 && exponent.Re.High == 0.5) + return Real(Re.Sqrt()); + if (IsReal && Re.Low > 0 && exponent.IsReal) + return Real((exponent.Re * Re.Log()).Exp()); + return (exponent * Log()).Exp(); + } + + internal ComplexInterval Sin() + { + if (IsReal) + return Real(Re.Sin()); + // sin(a + b i) = sin a cosh b + i cos a sinh b + var (cosh, sinh) = CoshSinh(Im); + return new(Re.Sin() * cosh, Re.Cos() * sinh); + } + + internal ComplexInterval Cos() + { + if (IsReal) + return Real(Re.Cos()); + // cos(a + b i) = cos a cosh b - i sin a sinh b + var (cosh, sinh) = CoshSinh(Im); + return new(Re.Cos() * cosh, -(Re.Sin() * sinh)); + } + + private static (Interval Cosh, Interval Sinh) CoshSinh(Interval b) + { + var up = b.Exp(); + var down = (-b).Exp(); + var half = Interval.Exactly(0.5); + return ((up + down) * half, (up - down) * half); + } + } + + /// + /// A number an expression stands for, as a that holds it: in + /// double arithmetic rounded outward, with no settings read -- neither the precision of the + /// decimal evaluation nor the downcasting -- so that a question needing a dozen digits is + /// answered in a dozen and gets the same answer however the caller's settings are set. The + /// pilot of #1019's item 10, evaluation to a requested accuracy, for the checks inside the + /// library; stays as it is beside it. + /// + /// + /// An interval is an answer with its error attached, so a comparison has three outcomes, not + /// two: the values agree to the accuracy asked, they differ by more, or the intervals are too + /// wide to tell -- where cancellation has eaten the digits, or a node is not one this reads, + /// and the caller asks the decimal evaluation instead. + /// https://github.com/asc-community/AngouriMath/issues/1019 + /// + internal static class IntervalEvaluation + { + /// + /// The interval evaluates to, or null where it holds a variable, a + /// node this does not read, or a value that is not finite -- a pole, a logarithm of zero, + /// a point on a branch cut. + /// + internal static ComplexInterval? Of(Entity expr) + { + var value = Evaluate(expr); + return value is { IsFinite: true } finite ? finite : null; + } + + /// + /// 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 + /// it, false where the whole of it is beyond, null where the intervals cannot tell. + /// + internal static bool? Agree(ComplexInterval left, ComplexInterval right, double relativeTolerance) + { + var difference = left - right; + var scale = Math.Max(1.0, Math.Max(left.Magnitude, right.Magnitude)); + var allowed = relativeTolerance * scale; + var largest = Math.Sqrt(difference.Re.Magnitude * difference.Re.Magnitude + difference.Im.Magnitude * difference.Im.Magnitude); + if (largest <= allowed) + return true; + var smallest = Math.Sqrt(difference.Re.Mignitude * difference.Re.Mignitude + difference.Im.Mignitude * difference.Im.Mignitude); + if (smallest > allowed) + return false; + return null; + } + + [AngouriMath.Core.ConstantField] private static readonly ComplexInterval Undefined = new(new(double.NaN, double.NaN), new(double.NaN, double.NaN)); + + private static ComplexInterval Evaluate(Entity expr) + { + switch (expr) + { + case Number.Complex number: + return FromNumber(number); + case Variable variable when variable == MathS.pi: + return ComplexInterval.Real(Interval.Around(Math.PI)); + case Variable variable when variable == MathS.e: + return ComplexInterval.Real(Interval.Around(Math.E)); + case Sumf(var left, var right): + return Evaluate(left) + Evaluate(right); + case Minusf(var left, var right): + return Evaluate(left) - Evaluate(right); + case Mulf(var left, var right): + return Evaluate(left) * Evaluate(right); + case Divf(var left, var right): + return Evaluate(left) / Evaluate(right); + case Powf(var @base, Number.Integer power) when power.EInteger.CanFitInInt64(): + return Evaluate(@base).Pow(power.EInteger.ToInt64Unchecked()); + case Powf(var @base, var exponent): + if (@base == MathS.e) + return Evaluate(exponent).Exp(); + return Evaluate(@base).Pow(Evaluate(exponent)); + case Sinf(var argument): + return Evaluate(argument).Sin(); + case Cosf(var argument): + return Evaluate(argument).Cos(); + case Tanf(var argument): + { + var inner = Evaluate(argument); + return inner.Sin() / inner.Cos(); + } + case Cotanf(var argument): + { + var inner = Evaluate(argument); + return inner.Cos() / inner.Sin(); + } + case Secantf(var argument): + return ComplexInterval.Real(Interval.Exactly(1)) / Evaluate(argument).Cos(); + case Cosecantf(var argument): + return ComplexInterval.Real(Interval.Exactly(1)) / Evaluate(argument).Sin(); + case Logf(var @base, var antilogarithm): + return Evaluate(antilogarithm).Log() / Evaluate(@base).Log(); + case Arctanf(var argument): + return RealOnly(Evaluate(argument), static x => x.Atan()); + case Arcsinf(var argument): + return RealOnly(Evaluate(argument), static x => x.Asin()); + case Arccosf(var argument): + return RealOnly(Evaluate(argument), static x => x.Acos()); + case Absf(var argument): + { + var inner = Evaluate(argument); + return inner.IsReal ? ComplexInterval.Real(inner.Re.Abs()) : ComplexInterval.Real((inner.Re.Square() + inner.Im.Square()).Sqrt()); + } + case Signumf(var argument): + { + var inner = Evaluate(argument); + if (!inner.IsReal || inner.Re.ContainsZero || !inner.Re.IsFinite) + return Undefined; + return ComplexInterval.Real(Interval.Exactly(inner.Re.Low > 0 ? 1 : -1)); + } + case Providedf(var inner, var predicate): + // The value where the condition holds, none where it fails, as the decimal + // evaluation has it; a condition the intervals cannot settle leaves the value + // unsettled too. + return Decide(predicate) switch + { + true => Evaluate(inner), + false => Undefined, + null => Undefined, + }; + case Piecewise piecewise: + // The first case whose condition holds; one that cannot be settled before it + // leaves the whole unsettled. + foreach (var @case in piecewise.Cases) + switch (Decide(@case.Predicate)) + { + case true: + return Evaluate(@case.Expression); + case null: + return Undefined; + } + return Undefined; + default: + return Undefined; + } + } + + /// + /// A condition at the point the expression was evaluated at: true, false, or null where + /// the intervals cannot tell -- an equality between two values their intervals do not + /// separate, a comparison of values that are not real, a node this does not read. + /// + private static bool? Decide(Entity predicate) + { + switch (predicate) + { + case Entity.Boolean value: + return value == Entity.Boolean.True; + case Notf(var argument): + return !Decide(argument); + case Andf(var left, var right): + { + var l = Decide(left); + if (l == false) return false; + var r = Decide(right); + if (r == false) return false; + return l == true && r == true ? true : null; + } + case Orf(var left, var right): + { + var l = Decide(left); + if (l == true) return true; + var r = Decide(right); + if (r == true) return true; + return l == false && r == false ? false : null; + } + case Xorf(var left, var right): + return Decide(left) is { } first && Decide(right) is { } second ? first != second : null; + case Impliesf(var assumption, var conclusion): + { + var a = Decide(assumption); + if (a == false) return true; + var c = Decide(conclusion); + if (c == true) return true; + return a == true && c == false ? false : null; + } + case Equalsf(var left, var right): + { + var difference = Evaluate(left) - Evaluate(right); + if (!difference.IsFinite) + return null; + if (difference.Re.Mignitude > 0 || difference.Im.Mignitude > 0) + return false; + return difference.Re.IsZero && difference.Im.IsZero ? true : null; + } + case Greaterf(var left, var right): + return Compare(left, right) is { } greater ? greater > 0 : null; + case GreaterOrEqualf(var left, var right): + return Compare(left, right, orEqual: true) is { } atLeast ? atLeast >= 0 : null; + case Lessf(var left, var right): + return Compare(right, left) is { } less ? less > 0 : null; + case LessOrEqualf(var left, var right): + return Compare(right, left, orEqual: true) is { } atMost ? atMost >= 0 : null; + default: + return null; + } + } + + /// + /// The sign of left - right, both real: 1 or -1 where the intervals are apart, 0 + /// where both are the same exact value (asked only for a comparison that allows it), and + /// null otherwise. + /// + private static int? Compare(Entity left, Entity right, bool orEqual = false) + { + var l = Evaluate(left); + var r = Evaluate(right); + if (!l.IsFinite || !r.IsFinite || !l.IsReal || !r.IsReal) + return null; + if (l.Re.Low > r.Re.High) + return 1; + if (l.Re.High < r.Re.Low) + return -1; + if (orEqual && l.Re.Low == l.Re.High && r.Re.Low == r.Re.High && l.Re.Low == r.Re.Low) + return 0; + return null; + } + + /// A function read on real arguments only, in this pilot. + private static ComplexInterval RealOnly(ComplexInterval argument, Func function) + => argument.IsReal && argument.Re.IsFinite ? ComplexInterval.Real(function(argument.Re)) : Undefined; + + private static ComplexInterval FromNumber(Number.Complex number) + { + if (number is Number.Rational rational) + return ComplexInterval.Real(FromRational(rational.ERational)); + var re = FromDecimal(number.RealPart.EDecimal); + var im = number.ImaginaryPart.EDecimal.IsZero ? Interval.Exactly(0) : FromDecimal(number.ImaginaryPart.EDecimal); + return new(re, im); + } + + /// + /// A rational as an interval: the double it rounds to, exactly where that is the rational + /// itself, and a unit in the last place either way otherwise -- from the exact value, not + /// from the hundred-digit decimal that stands for a third. + /// + private static Interval FromRational(ERational value) + { + var converted = value.ToDouble(); + if (!Interval.IsFiniteDouble(converted)) + return new(double.NaN, double.NaN); + return ERational.FromDouble(converted).CompareTo(value) == 0 ? Interval.Exactly(converted) : Interval.Around(converted); + } + + /// The decimal as an interval, the same way. + private static Interval FromDecimal(EDecimal value) + { + if (!value.IsFinite) + return new(double.NaN, double.NaN); + var converted = value.ToDouble(); + if (!Interval.IsFiniteDouble(converted)) + return new(double.NaN, double.NaN); + return EDecimal.FromDouble(converted).CompareTo(value) == 0 ? Interval.Exactly(converted) : Interval.Around(converted); + } + } +} diff --git a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs new file mode 100644 index 000000000..31ad49f4c --- /dev/null +++ b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs @@ -0,0 +1,482 @@ +// +// 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; + +namespace AngouriMath.Numerics +{ + /// + /// An interval of decimals at a working precision the evaluation carries: every operation + /// rounds its lower bound down and its upper bound up in that many digits, so the exact value + /// is inside. 's second tier, for where doubles cannot tell. + /// + internal readonly struct PreciseInterval + { + internal readonly EDecimal Low; + internal readonly EDecimal High; + + internal PreciseInterval(EDecimal low, EDecimal high) + { + Low = low; + High = high; + } + + internal static PreciseInterval Exactly(EDecimal value) => new(value, value); + + internal bool IsFinite => Low.IsFinite && High.IsFinite; + internal bool IsZero => Low.IsZero && High.IsZero; + internal bool ContainsZero => Low.Sign <= 0 && High.Sign >= 0; + internal EDecimal Magnitude => EDecimal.Max(Low.Abs(), High.Abs()); + internal EDecimal Mignitude => ContainsZero ? EDecimal.Zero : EDecimal.Min(Low.Abs(), High.Abs()); + + internal PreciseInterval Negate() => new(High.Negate(), Low.Negate()); + } + + /// A rectangle of two s. + internal readonly struct PreciseComplexInterval + { + internal readonly PreciseInterval Re; + internal readonly PreciseInterval Im; + + internal PreciseComplexInterval(PreciseInterval re, PreciseInterval im) + { + Re = re; + Im = im; + } + + internal bool IsFinite => Re.IsFinite && Im.IsFinite; + internal bool IsReal => Im.IsZero; + } + + /// + /// at a requested number of significant digits, for the + /// points where double intervals are too wide to tell -- the digits lost to cancellation -- + /// or land on a cut by a rounding. The precision is an argument, never a setting. + /// + internal sealed class PreciseEvaluation + { + private readonly int digits; + private readonly EContext down; + private readonly EContext up; + private readonly EContext near; + private readonly EDecimal unit; + private readonly PreciseComplexInterval undefined = + new(new(EDecimal.NaN, EDecimal.NaN), new(EDecimal.NaN, EDecimal.NaN)); + + private PreciseEvaluation(int digits) + { + this.digits = digits; + down = EContext.ForPrecisionAndRounding(digits, ERounding.Floor); + up = EContext.ForPrecisionAndRounding(digits, ERounding.Ceiling); + near = EContext.ForPrecisionAndRounding(digits + 10, ERounding.HalfEven); + unit = EDecimal.FromString("1E-" + (digits - 2)); + } + + /// The interval evaluates to in digits, or null. + internal static PreciseComplexInterval? Of(Entity expr, int digits) + { + var value = new PreciseEvaluation(digits).Evaluate(expr); + return value.IsFinite ? value : null; + } + + /// The three-way agreement of , in decimals. + internal static bool? Agree(PreciseComplexInterval left, PreciseComplexInterval right, double relativeTolerance) + { + var context = EContext.ForPrecision(30); + var reDifference = new PreciseInterval(left.Re.Low.Subtract(right.Re.High, context), left.Re.High.Subtract(right.Re.Low, context)); + var imDifference = new PreciseInterval(left.Im.Low.Subtract(right.Im.High, context), left.Im.High.Subtract(right.Im.Low, context)); + var scale = EDecimal.Max(EDecimal.One, EDecimal.Max( + EDecimal.Max(left.Re.Magnitude, left.Im.Magnitude), EDecimal.Max(right.Re.Magnitude, right.Im.Magnitude))); + var allowed = scale.Multiply(EDecimal.FromDouble(relativeTolerance), context); + var largest = EDecimal.Max(reDifference.Magnitude, imDifference.Magnitude); + if (largest.CompareTo(allowed) <= 0) + return true; + var smallest = EDecimal.Max(reDifference.Mignitude, imDifference.Mignitude); + if (smallest.CompareTo(allowed) > 0) + return false; + return null; + } + + // ---- real intervals ------------------------------------------------------------- + + private PreciseInterval Add(PreciseInterval a, PreciseInterval b) + => a.IsZero ? b : b.IsZero ? a : new(a.Low.Add(b.Low, down), a.High.Add(b.High, up)); + + private PreciseInterval Subtract(PreciseInterval a, PreciseInterval b) => Add(a, b.Negate()); + + private PreciseInterval Multiply(PreciseInterval a, PreciseInterval b) + { + if (a.IsZero || b.IsZero) + return PreciseInterval.Exactly(EDecimal.Zero); + var lows = new[] { a.Low.Multiply(b.Low, down), a.Low.Multiply(b.High, down), a.High.Multiply(b.Low, down), a.High.Multiply(b.High, down) }; + var highs = new[] { a.Low.Multiply(b.Low, up), a.Low.Multiply(b.High, up), a.High.Multiply(b.Low, up), a.High.Multiply(b.High, up) }; + return new(Min(lows), Max(highs)); + } + + private PreciseInterval Divide(PreciseInterval a, PreciseInterval b) + { + if (b.ContainsZero) + return new(EDecimal.NaN, EDecimal.NaN); + if (a.IsZero) + return PreciseInterval.Exactly(EDecimal.Zero); + var lows = new[] { a.Low.Divide(b.Low, down), a.Low.Divide(b.High, down), a.High.Divide(b.Low, down), a.High.Divide(b.High, down) }; + var highs = new[] { a.Low.Divide(b.Low, up), a.Low.Divide(b.High, up), a.High.Divide(b.Low, up), a.High.Divide(b.High, up) }; + return new(Min(lows), Max(highs)); + } + + private PreciseInterval Square(PreciseInterval a) + { + if (a.ContainsZero) + return new(EDecimal.Zero, EDecimal.Max(a.Low.Multiply(a.Low, up), a.High.Multiply(a.High, up))); + var small = a.Mignitude; + var large = a.Magnitude; + return new(small.Multiply(small, down), large.Multiply(large, up)); + } + + /// + /// A function's value over the interval, from its value at the middle, computed with ten + /// digits to spare, and a bound on its slope there: the value moves no further than the + /// slope times the half-width, and the evaluation's own error is a unit or two in the last + /// place of the working precision. + /// + private PreciseInterval AroundTheMiddle(PreciseInterval a, Func function, EDecimal slope) + { + var middle = a.Low.Add(a.High, near).Divide(EDecimal.FromInt32(2), near); + var halfWidth = a.High.Subtract(middle, up).Abs(); + var value = function(middle); + if (!value.IsFinite) + return new(EDecimal.NaN, EDecimal.NaN); + var reach = slope.Multiply(halfWidth, up).Add(value.Abs().Multiply(unit, up), up) + .Add(EDecimal.FromString("1E-" + (digits + 30)), up); + return new(value.Subtract(reach, down), value.Add(reach, up)); + } + + private PreciseInterval Exp(PreciseInterval a) + => AroundTheMiddle(a, x => x.Exp(near), a.High.Exp(up).Abs()); + + private PreciseInterval Log(PreciseInterval a) + => a.Low.Sign > 0 ? AroundTheMiddle(a, x => x.Log(near), EDecimal.One.Divide(a.Low, up)) : new(EDecimal.NaN, EDecimal.NaN); + + private PreciseInterval Sqrt(PreciseInterval a) + => a.Low.Sign > 0 ? AroundTheMiddle(a, x => x.Sqrt(near), EDecimal.One.Divide(a.Low.Sqrt(down).Multiply(EDecimal.FromInt32(2), down), up)) + : a.IsZero ? a : new(EDecimal.NaN, EDecimal.NaN); + + private PreciseInterval Sin(PreciseInterval a) => AroundTheMiddle(a, x => x.Sin(near), EDecimal.One); + private PreciseInterval Cos(PreciseInterval a) => AroundTheMiddle(a, x => x.Cos(near), EDecimal.One); + private PreciseInterval Atan(PreciseInterval a) => AroundTheMiddle(a, x => x.Arctan(near), EDecimal.One); + + private PreciseInterval Asin(PreciseInterval a) + => InsideTheUnitInterval(a) is { } slope ? AroundTheMiddle(a, x => x.Arcsin(near), slope) : new(EDecimal.NaN, EDecimal.NaN); + + private PreciseInterval Acos(PreciseInterval a) + => InsideTheUnitInterval(a) is { } slope ? AroundTheMiddle(a, x => x.Acos(near), slope) : new(EDecimal.NaN, EDecimal.NaN); + + /// The slope bound of the arcsine and arccosine, 1/sqrt(1 - x^2), strictly inside (-1, 1). + private EDecimal? InsideTheUnitInterval(PreciseInterval a) + { + var largest = a.Magnitude; + if (largest.CompareTo(EDecimal.One) >= 0) + return null; + return EDecimal.One.Divide(EDecimal.One.Subtract(largest.Multiply(largest, up), down).Sqrt(down), up); + } + + private static EDecimal Min(EDecimal[] values) + { + var result = values[0]; + foreach (var value in values) + if (value.IsNaN() || value.CompareTo(result) < 0) result = value; + return result; + } + + private static EDecimal Max(EDecimal[] values) + { + var result = values[0]; + foreach (var value in values) + if (value.IsNaN() || value.CompareTo(result) > 0) result = value; + return result; + } + + // ---- complex rectangles --------------------------------------------------------- + + private PreciseComplexInterval Real(PreciseInterval re) => new(re, PreciseInterval.Exactly(EDecimal.Zero)); + + private PreciseComplexInterval Add(PreciseComplexInterval a, PreciseComplexInterval b) => new(Add(a.Re, b.Re), Add(a.Im, b.Im)); + private PreciseComplexInterval Subtract(PreciseComplexInterval a, PreciseComplexInterval b) => new(Subtract(a.Re, b.Re), Subtract(a.Im, b.Im)); + + private PreciseComplexInterval Multiply(PreciseComplexInterval a, PreciseComplexInterval b) + => a.IsReal && b.IsReal ? Real(Multiply(a.Re, b.Re)) + : new(Subtract(Multiply(a.Re, b.Re), Multiply(a.Im, b.Im)), Add(Multiply(a.Re, b.Im), Multiply(a.Im, b.Re))); + + private PreciseComplexInterval Divide(PreciseComplexInterval a, PreciseComplexInterval b) + { + if (b.IsReal) + return new(Divide(a.Re, b.Re), a.IsReal ? PreciseInterval.Exactly(EDecimal.Zero) : Divide(a.Im, b.Re)); + var denominator = Add(Square(b.Re), Square(b.Im)); + return new(Divide(Add(Multiply(a.Re, b.Re), Multiply(a.Im, b.Im)), denominator), + Divide(Subtract(Multiply(a.Im, b.Re), Multiply(a.Re, b.Im)), denominator)); + } + + private PreciseComplexInterval Exp(PreciseComplexInterval a) + { + var modulus = Exp(a.Re); + return a.IsReal ? Real(modulus) : new(Multiply(modulus, Cos(a.Im)), Multiply(modulus, Sin(a.Im))); + } + + private PreciseComplexInterval Log(PreciseComplexInterval a) + { + if (a.IsReal && a.Re.Low.Sign > 0) + return Real(Log(a.Re)); + if (a.IsReal && a.Re.High.Sign < 0) + return new(Log(a.Re.Negate()), Pi()); + if (a.Re.High.Sign <= 0 && a.Im.ContainsZero || a.Re.ContainsZero && a.Im.ContainsZero) + return undefined; + var logModulus = Log(Add(Square(a.Re), Square(a.Im))); + var half = EDecimal.FromString("0.5"); + var re = new PreciseInterval(logModulus.Low.Multiply(half, down), logModulus.High.Multiply(half, up)); + // The argument over a rectangle off the cut takes its extremes at the corners. + var corners = new[] + { + a.Im.Low.Arctan2(a.Re.Low, near), a.Im.Low.Arctan2(a.Re.High, near), + a.Im.High.Arctan2(a.Re.Low, near), a.Im.High.Arctan2(a.Re.High, near), + }; + var low = Min(corners); + var high = Max(corners); + var slack = EDecimal.Max(low.Abs(), high.Abs()).Multiply(unit, up).Add(EDecimal.FromString("1E-" + (digits + 30)), up); + return new(re, new(low.Subtract(slack, down), high.Add(slack, up))); + } + + private PreciseComplexInterval Pow(PreciseComplexInterval a, long exponent) + { + if (exponent < 0) + return Divide(Real(PreciseInterval.Exactly(EDecimal.One)), Pow(a, -exponent)); + var result = Real(PreciseInterval.Exactly(EDecimal.One)); + var power = a; + while (exponent > 0) + { + if ((exponent & 1) == 1) + result = Multiply(result, power); + exponent >>= 1; + if (exponent > 0) + power = Multiply(power, power); + } + return result; + } + + private PreciseComplexInterval Pow(PreciseComplexInterval a, PreciseComplexInterval exponent) + { + if (a.IsReal && a.Re.Low.Sign >= 0 && exponent.IsReal && exponent.Re.Low.CompareTo(EDecimal.FromString("0.5")) == 0 + && exponent.Re.High.CompareTo(EDecimal.FromString("0.5")) == 0) + return Real(Sqrt(a.Re)); + if (a.IsReal && a.Re.Low.Sign > 0 && exponent.IsReal) + return Real(Exp(Multiply(exponent.Re, Log(a.Re)))); + return Exp(Multiply(exponent, Log(a))); + } + + private PreciseComplexInterval Sin(PreciseComplexInterval a) + { + if (a.IsReal) + return Real(Sin(a.Re)); + var (cosh, sinh) = CoshSinh(a.Im); + return new(Multiply(Sin(a.Re), cosh), Multiply(Cos(a.Re), sinh)); + } + + private PreciseComplexInterval Cos(PreciseComplexInterval a) + { + if (a.IsReal) + return Real(Cos(a.Re)); + var (cosh, sinh) = CoshSinh(a.Im); + return new(Multiply(Cos(a.Re), cosh), Multiply(Sin(a.Re), sinh).Negate()); + } + + private (PreciseInterval Cosh, PreciseInterval Sinh) CoshSinh(PreciseInterval b) + { + var expUp = Exp(b); + var expDown = Exp(b.Negate()); + var half = PreciseInterval.Exactly(EDecimal.FromString("0.5")); + return (Multiply(Add(expUp, expDown), half), Multiply(Subtract(expUp, expDown), half)); + } + + private PreciseInterval Pi() + { + var pi = EDecimal.FromString( + "3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196"); + return new(pi.RoundToPrecision(down), pi.RoundToPrecision(up)); + } + + private PreciseInterval E() + { + var e = EDecimal.FromString( + "2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901"); + return new(e.RoundToPrecision(down), e.RoundToPrecision(up)); + } + + private PreciseComplexInterval Evaluate(Entity expr) + { + switch (expr) + { + case Number.Rational rational: + { + var value = rational.ERational; + return Real(new(value.ToEDecimal(down), value.ToEDecimal(up))); + } + case Number.Complex number: + { + var re = number.RealPart.EDecimal; + var im = number.ImaginaryPart.EDecimal; + return new(new(re.RoundToPrecision(down), re.RoundToPrecision(up)), new(im.RoundToPrecision(down), im.RoundToPrecision(up))); + } + case Variable variable when variable == MathS.pi: + return Real(Pi()); + case Variable variable when variable == MathS.e: + return Real(E()); + case Sumf(var left, var right): + return Add(Evaluate(left), Evaluate(right)); + case Minusf(var left, var right): + return Subtract(Evaluate(left), Evaluate(right)); + case Mulf(var left, var right): + return Multiply(Evaluate(left), Evaluate(right)); + case Divf(var left, var right): + return Divide(Evaluate(left), Evaluate(right)); + case Powf(var @base, Number.Integer power) when power.EInteger.CanFitInInt64(): + return Pow(Evaluate(@base), power.EInteger.ToInt64Unchecked()); + case Powf(var @base, var exponent): + return @base == MathS.e ? Exp(Evaluate(exponent)) : Pow(Evaluate(@base), Evaluate(exponent)); + case Sinf(var argument): + return Sin(Evaluate(argument)); + case Cosf(var argument): + return Cos(Evaluate(argument)); + case Tanf(var argument): + { + var inner = Evaluate(argument); + return Divide(Sin(inner), Cos(inner)); + } + case Cotanf(var argument): + { + var inner = Evaluate(argument); + return Divide(Cos(inner), Sin(inner)); + } + case Secantf(var argument): + return Divide(Real(PreciseInterval.Exactly(EDecimal.One)), Cos(Evaluate(argument))); + case Cosecantf(var argument): + return Divide(Real(PreciseInterval.Exactly(EDecimal.One)), Sin(Evaluate(argument))); + case Logf(var @base, var antilogarithm): + return Divide(Log(Evaluate(antilogarithm)), Log(Evaluate(@base))); + case Arctanf(var argument): + return RealOnly(Evaluate(argument), Atan); + case Arcsinf(var argument): + return RealOnly(Evaluate(argument), Asin); + case Arccosf(var argument): + return RealOnly(Evaluate(argument), Acos); + case Absf(var argument): + { + var inner = Evaluate(argument); + if (!inner.IsReal) + return Real(Sqrt(Add(Square(inner.Re), Square(inner.Im)))); + var re = inner.Re; + return Real(re.ContainsZero ? new(EDecimal.Zero, re.Magnitude) : re.Low.Sign > 0 ? re : re.Negate()); + } + case Signumf(var argument): + { + var inner = Evaluate(argument); + if (!inner.IsFinite || !inner.IsReal || inner.Re.ContainsZero) + return undefined; + return Real(PreciseInterval.Exactly(inner.Re.Low.Sign > 0 ? EDecimal.One : EDecimal.One.Negate())); + } + case Providedf(var inner, var predicate): + return Decide(predicate) == true ? Evaluate(inner) : undefined; + case Piecewise piecewise: + foreach (var @case in piecewise.Cases) + switch (Decide(@case.Predicate)) + { + case true: + return Evaluate(@case.Expression); + case null: + return undefined; + } + return undefined; + default: + return undefined; + } + } + + private PreciseComplexInterval RealOnly(PreciseComplexInterval argument, Func function) + => argument.IsReal && argument.Re.IsFinite ? Real(function(argument.Re)) : undefined; + + private bool? Decide(Entity predicate) + { + switch (predicate) + { + case Entity.Boolean value: + return value == Entity.Boolean.True; + case Notf(var argument): + return !Decide(argument); + case Andf(var left, var right): + { + var l = Decide(left); + if (l == false) return false; + var r = Decide(right); + if (r == false) return false; + return l == true && r == true ? true : null; + } + case Orf(var left, var right): + { + var l = Decide(left); + if (l == true) return true; + var r = Decide(right); + if (r == true) return true; + return l == false && r == false ? false : null; + } + case Xorf(var left, var right): + return Decide(left) is { } first && Decide(right) is { } second ? first != second : null; + case Impliesf(var assumption, var conclusion): + { + var a = Decide(assumption); + if (a == false) return true; + var c = Decide(conclusion); + if (c == true) return true; + return a == true && c == false ? false : null; + } + case Equalsf(var left, var right): + { + var difference = Subtract(Evaluate(left), Evaluate(right)); + if (!difference.IsFinite) + return null; + if (difference.Re.Mignitude.Sign > 0 || difference.Im.Mignitude.Sign > 0) + return false; + return difference.Re.IsZero && difference.Im.IsZero ? true : null; + } + case Greaterf(var left, var right): + return Compare(left, right) is { } greater ? greater > 0 : null; + case GreaterOrEqualf(var left, var right): + return Compare(left, right, orEqual: true) is { } atLeast ? atLeast >= 0 : null; + case Lessf(var left, var right): + return Compare(right, left) is { } less ? less > 0 : null; + case LessOrEqualf(var left, var right): + return Compare(right, left, orEqual: true) is { } atMost ? atMost >= 0 : null; + default: + return null; + } + } + + private int? Compare(Entity left, Entity right, bool orEqual = false) + { + var l = Evaluate(left); + var r = Evaluate(right); + if (!l.IsFinite || !r.IsFinite || !l.IsReal || !r.IsReal) + return null; + if (l.Re.Low.CompareTo(r.Re.High) > 0) + return 1; + if (l.Re.High.CompareTo(r.Re.Low) < 0) + return -1; + if (orEqual && l.Re.Low.CompareTo(l.Re.High) == 0 && r.Re.Low.CompareTo(r.Re.High) == 0 && l.Re.Low.CompareTo(r.Re.Low) == 0) + return 0; + return null; + } + } +} diff --git a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs new file mode 100644 index 000000000..24cf4bfa4 --- /dev/null +++ b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs @@ -0,0 +1,131 @@ +// +// 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.Extensions; +using AngouriMath.Numerics; +using PeterO.Numbers; +using Xunit; + +namespace AngouriMath.Tests.Core +{ + /// + /// The interval evaluation holds the value it stands for -- the hundred-digit decimal + /// evaluation's, inside its bounds -- and is narrow where nothing cancels, declines what it + /// does not read, and compares in three outcomes. + /// #1019 + /// + public sealed class IntervalEvaluationTest + { + [Theory] + [InlineData("sqrt(2)")] + [InlineData("pi")] + [InlineData("e^2")] + [InlineData("sin(1)")] + [InlineData("cos(1/3)")] + [InlineData("tan(3/10)")] + [InlineData("sec(1) + csc(2) + cot(3)")] + [InlineData("ln(7)")] + [InlineData("log(2, 10)")] + [InlineData("arctan(2)")] + [InlineData("arcsin(1/3)")] + [InlineData("arccos(-2/3)")] + [InlineData("(2/3)^(1/3)")] + [InlineData("(1 + 2i)^3")] + [InlineData("1/(3 + 4i)")] + [InlineData("e^(i * pi/4)")] + [InlineData("sqrt(-3)")] + [InlineData("sin(1 + i)")] + [InlineData("cos(2 - 3i)")] + [InlineData("ln(-2 + i)")] + [InlineData("ln(-2)")] + [InlineData("abs(-3 + 4i)")] + [InlineData("(sin(1)^2 + cos(1)^2) * 7/3")] + [InlineData("sqrt(2) provided 1 < 2")] + [InlineData("piecewise(1 provided 2 < 1, sin(1) provided 3 > 2)")] + public void TheIntervalHoldsTheValue(string expression) + { + var expr = expression.ToEntity(); + var interval = IntervalEvaluation.Of(expr); + Assert.True(interval.HasValue, $"{expression} was not evaluated"); + var value = expr.EvalNumerical(); + Holds(interval!.Value.Re, value.RealPart.EDecimal, expression, "real"); + Holds(interval.Value.Im, value.ImaginaryPart.EDecimal, expression, "imaginary"); + } + + private static void Holds(Interval interval, EDecimal value, string expression, string part) + { + Assert.True(EDecimal.FromDouble(interval.Low).CompareTo(value) <= 0 && EDecimal.FromDouble(interval.High).CompareTo(value) >= 0, + $"the {part} part of {expression}, {value}, is outside [{interval.Low:R}, {interval.High:R}]"); + var scale = System.Math.Max(1.0, interval.Magnitude); + Assert.True(interval.High - interval.Low <= 1e-12 * scale, + $"the {part} part of {expression} is [{interval.Low:R}, {interval.High:R}], wider than it has to be"); + } + + /// + /// Declined: a variable, a pole, the logarithm of zero, a node the pilot does not read, + /// and a value whose condition fails. + /// + [Theory] + [InlineData("x + 1")] + [InlineData("1/0")] + [InlineData("ln(0)")] + [InlineData("tan(pi/2)")] + [InlineData("3!")] + [InlineData("1 provided 3 < 2")] + public void UndecidedWhereItShouldBe(string expression) + => Assert.Null(IntervalEvaluation.Of(expression.ToEntity())); + + /// + /// Where cancellation leaves double intervals too wide to tell, forty digits tell: the + /// hundred million cancels and leaves the sine, which the doubles hold only to within a + /// few hundred-millionths. + /// + [Fact] + public void MoreDigitsSettleWhatDoublesCannot() + { + var cancelled = "(10^8 + sin(1)) - 10^8".ToEntity(); + var sine = "sin(1)".ToEntity(); + Assert.Null(IntervalEvaluation.Agree(IntervalEvaluation.Of(cancelled)!.Value, IntervalEvaluation.Of(sine)!.Value, 1e-9)); + Assert.True(PreciseEvaluation.Agree(PreciseEvaluation.Of(cancelled, 40)!.Value, PreciseEvaluation.Of(sine, 40)!.Value, 1e-9)); + var nearby = "sin(1) + 1/1000000".ToEntity(); + Assert.False(PreciseEvaluation.Agree(PreciseEvaluation.Of(nearby, 40)!.Value, PreciseEvaluation.Of(sine, 40)!.Value, 1e-9)); + } + + [Theory] + [InlineData("sqrt(2)")] + [InlineData("e^2 - pi")] + [InlineData("sin(1 + i) * cos(2 - 3i)")] + [InlineData("ln(-2 + i) + arctan(2) - arcsin(1/3)")] + [InlineData("(2/3)^(1/3) + sec(1) + csc(2)")] + [InlineData("sqrt(2) provided 1 < 2")] + public void ThePreciseIntervalHoldsTheValue(string expression) + { + var expr = expression.ToEntity(); + var interval = PreciseEvaluation.Of(expr, 40); + Assert.True(interval.HasValue, $"{expression} was not evaluated"); + var value = expr.EvalNumerical(); + Assert.True(interval!.Value.Re.Low.CompareTo(value.RealPart.EDecimal) <= 0 && interval.Value.Re.High.CompareTo(value.RealPart.EDecimal) >= 0, + $"the real part of {expression}, {value.RealPart}, is outside [{interval.Value.Re.Low}, {interval.Value.Re.High}]"); + Assert.True(interval.Value.Im.Low.CompareTo(value.ImaginaryPart.EDecimal) <= 0 && interval.Value.Im.High.CompareTo(value.ImaginaryPart.EDecimal) >= 0, + $"the imaginary part of {expression}, {value.ImaginaryPart}, is outside [{interval.Value.Im.Low}, {interval.Value.Im.High}]"); + Assert.True(interval.Value.Re.High.Subtract(interval.Value.Re.Low).Abs().CompareTo(EDecimal.FromString("1E-30")) <= 0, + $"the real part of {expression} is wider than forty digits should leave it"); + } + + [Fact] + public void AgreementHasThreeOutcomes() + { + var one = IntervalEvaluation.Of("1".ToEntity())!.Value; + var identity = IntervalEvaluation.Of("sin(1)^2 + cos(1)^2".ToEntity())!.Value; + var nearby = IntervalEvaluation.Of("1 + 1/1000000".ToEntity())!.Value; + Assert.True(IntervalEvaluation.Agree(identity, one, 1e-9)); + Assert.False(IntervalEvaluation.Agree(nearby, one, 1e-9)); + var wide = ComplexInterval.Real(new Interval(0.9, 1.1)); + Assert.Null(IntervalEvaluation.Agree(wide, one, 1e-9)); + } + } +} From 062cb492ba666a9e18418cb225c039ddc50dd20f Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Sat, 26 Sep 2026 23:22:24 +0000 Subject: [PATCH 2/9] The evaluation reads arccot, arcsec and arccsc, and its bounds cover each function's own error The three are read on a real argument as the decimal evaluation defines them: arctan(1/x) with pi/2 at zero, arccos(1/x) and arcsin(1/x). Without them, a check over an integrand holding arccot could compare no point, and would turn every answer away. In doubles, the sine and cosine widen the function's value before adding the reach, and the complex logarithm's argument takes the same slack as every other function. In decimals, the sine, cosine and inverse functions carry the absolute error of the fixed point they are computed in, the sine's growing with the argument. A whole power is read directly only in the range of an int: 2^(-2^63) negated its exponent into itself. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Numerics/IntervalEvaluation.cs | 28 ++++++--- .../AngouriMath/Numerics/PreciseEvaluation.cs | 61 +++++++++++++------ .../UnitTests/Core/IntervalEvaluationTest.cs | 20 ++++++ 3 files changed, 84 insertions(+), 25 deletions(-) diff --git a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs index 0c57405d6..6354f3668 100644 --- a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs +++ b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs @@ -110,8 +110,11 @@ private Interval AroundTheMiddle(Func function, double slope) { var middle = Middle; var value = function(middle); - var reach = Up(Math.Max(High - middle, middle - Low) * slope); - return new(Math.Max(-1, Wider(value - reach, -1)), Math.Min(1, Wider(value + reach, 1))); + var reach = Up(Math.Max(Up(High - middle), Up(middle - Low)) * slope); + // The function's error is relative to its value, so the value is widened before the + // reach is added: widened after, the slack is a few units of the sum, which is far + // smaller than a unit of the value where the reach nearly cancels it. + return new(Math.Max(-1, Down(Wider(value, -1) - reach)), Math.Min(1, Up(Wider(value, 1) + reach))); } /// @@ -119,7 +122,7 @@ private Interval AroundTheMiddle(Func function, double slope) /// more near zero: the platform's functions are within one of the true value, and the /// slack keeps the claim true where they are not quite. /// - private static double Wider(double value, int direction) + internal static double Wider(double value, int direction) { for (var i = 0; i < 4; i++) value = direction > 0 ? Up(value) : Down(value); @@ -203,8 +206,8 @@ internal ComplexInterval Log() var a3 = Math.Atan2(Im.High, Re.Low); var a4 = Math.Atan2(Im.High, Re.High); var argument = new Interval( - Interval.Down(Interval.Down(Math.Min(Math.Min(a1, a2), Math.Min(a3, a4)))), - Interval.Up(Interval.Up(Math.Max(Math.Max(a1, a2), Math.Max(a3, a4))))); + Interval.Wider(Math.Min(Math.Min(a1, a2), Math.Min(a3, a4)), -1), + Interval.Wider(Math.Max(Math.Max(a1, a2), Math.Max(a3, a4)), 1)); return new(new(logModulus.Low / 2, logModulus.High / 2), argument); } @@ -331,8 +334,9 @@ private static ComplexInterval Evaluate(Entity expr) return Evaluate(left) * Evaluate(right); case Divf(var left, var right): return Evaluate(left) / Evaluate(right); - case Powf(var @base, Number.Integer power) when power.EInteger.CanFitInInt64(): - return Evaluate(@base).Pow(power.EInteger.ToInt64Unchecked()); + // In the range of an int, so that the reciprocal's negation cannot overflow. + case Powf(var @base, Number.Integer power) when power.EInteger.CanFitInInt32(): + return Evaluate(@base).Pow(power.EInteger.ToInt32Checked()); case Powf(var @base, var exponent): if (@base == MathS.e) return Evaluate(exponent).Exp(); @@ -356,13 +360,21 @@ private static ComplexInterval Evaluate(Entity expr) case Cosecantf(var argument): return ComplexInterval.Real(Interval.Exactly(1)) / Evaluate(argument).Sin(); case Logf(var @base, var antilogarithm): - return Evaluate(antilogarithm).Log() / Evaluate(@base).Log(); + return @base == MathS.e ? Evaluate(antilogarithm).Log() : Evaluate(antilogarithm).Log() / Evaluate(@base).Log(); case Arctanf(var argument): return RealOnly(Evaluate(argument), static x => x.Atan()); case Arcsinf(var argument): return RealOnly(Evaluate(argument), static x => x.Asin()); case Arccosf(var argument): return RealOnly(Evaluate(argument), static x => x.Acos()); + // As the decimal evaluation defines them: arctan(1/x), with pi/2 at zero, where + // the function jumps; arccos(1/x); arcsin(1/x). + case Arccotanf(var argument): + return RealOnly(Evaluate(argument), static x => x.IsZero ? Interval.Around(Math.PI / 2) : (Interval.Exactly(1) / x).Atan()); + case Arcsecantf(var argument): + return RealOnly(Evaluate(argument), static x => (Interval.Exactly(1) / x).Acos()); + case Arccosecantf(var argument): + return RealOnly(Evaluate(argument), static x => (Interval.Exactly(1) / x).Asin()); case Absf(var argument): { var inner = Evaluate(argument); diff --git a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs index 31ad49f4c..b27a1e1ca 100644 --- a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs +++ b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs @@ -66,6 +66,8 @@ internal sealed class PreciseEvaluation private readonly EContext up; private readonly EContext near; private readonly EDecimal unit; + private readonly EDecimal fixedPointError; + private readonly EDecimal least; private readonly PreciseComplexInterval undefined = new(new(EDecimal.NaN, EDecimal.NaN), new(EDecimal.NaN, EDecimal.NaN)); @@ -76,6 +78,8 @@ private PreciseEvaluation(int digits) up = EContext.ForPrecisionAndRounding(digits, ERounding.Ceiling); near = EContext.ForPrecisionAndRounding(digits + 10, ERounding.HalfEven); unit = EDecimal.FromString("1E-" + (digits - 2)); + fixedPointError = EDecimal.FromString("1E-" + (digits + 5)); + least = EDecimal.FromString("1E-" + (digits + 30)); } /// The interval evaluates to in digits, or null. @@ -142,40 +146,50 @@ private PreciseInterval Square(PreciseInterval a) /// /// A function's value over the interval, from its value at the middle, computed with ten /// digits to spare, and a bound on its slope there: the value moves no further than the - /// slope times the half-width, and the evaluation's own error is a unit or two in the last - /// place of the working precision. + /// slope times the half-width. The function's own error is a unit or two in the last place + /// of the working precision, relative to the value -- and, for the functions this library + /// computes in fixed point, besides, which is what is + /// left of it where the value is small. /// - private PreciseInterval AroundTheMiddle(PreciseInterval a, Func function, EDecimal slope) + private PreciseInterval AroundTheMiddle(PreciseInterval a, Func function, EDecimal slope, EDecimal absoluteError) { var middle = a.Low.Add(a.High, near).Divide(EDecimal.FromInt32(2), near); - var halfWidth = a.High.Subtract(middle, up).Abs(); + var halfWidth = EDecimal.Max(a.High.Subtract(middle, up), middle.Subtract(a.Low, up)); var value = function(middle); if (!value.IsFinite) return new(EDecimal.NaN, EDecimal.NaN); var reach = slope.Multiply(halfWidth, up).Add(value.Abs().Multiply(unit, up), up) - .Add(EDecimal.FromString("1E-" + (digits + 30)), up); + .Add(EDecimal.Max(absoluteError, least), up); return new(value.Subtract(reach, down), value.Add(reach, up)); } + // PeterO's exponential, logarithm and square root are within a unit of the value; the + // sine, the cosine and the inverse functions are this library's, in fixed point, and + // their error is absolute -- the sine's grows with the argument, whose reduction by + // 2 pi uses pi to the same fixed point. + private PreciseInterval Exp(PreciseInterval a) - => AroundTheMiddle(a, x => x.Exp(near), a.High.Exp(up).Abs()); + => AroundTheMiddle(a, x => x.Exp(near), a.High.Exp(up).Abs(), EDecimal.Zero); private PreciseInterval Log(PreciseInterval a) - => a.Low.Sign > 0 ? AroundTheMiddle(a, x => x.Log(near), EDecimal.One.Divide(a.Low, up)) : new(EDecimal.NaN, EDecimal.NaN); + => a.Low.Sign > 0 ? AroundTheMiddle(a, x => x.Log(near), EDecimal.One.Divide(a.Low, up), EDecimal.Zero) : new(EDecimal.NaN, EDecimal.NaN); private PreciseInterval Sqrt(PreciseInterval a) - => a.Low.Sign > 0 ? AroundTheMiddle(a, x => x.Sqrt(near), EDecimal.One.Divide(a.Low.Sqrt(down).Multiply(EDecimal.FromInt32(2), down), up)) + => a.Low.Sign > 0 ? AroundTheMiddle(a, x => x.Sqrt(near), EDecimal.One.Divide(a.Low.Sqrt(down).Multiply(EDecimal.FromInt32(2), down), up), EDecimal.Zero) : a.IsZero ? a : new(EDecimal.NaN, EDecimal.NaN); - private PreciseInterval Sin(PreciseInterval a) => AroundTheMiddle(a, x => x.Sin(near), EDecimal.One); - private PreciseInterval Cos(PreciseInterval a) => AroundTheMiddle(a, x => x.Cos(near), EDecimal.One); - private PreciseInterval Atan(PreciseInterval a) => AroundTheMiddle(a, x => x.Arctan(near), EDecimal.One); + private PreciseInterval Sin(PreciseInterval a) => AroundTheMiddle(a, x => x.Sin(near), EDecimal.One, ReducedArgumentError(a)); + private PreciseInterval Cos(PreciseInterval a) => AroundTheMiddle(a, x => x.Cos(near), EDecimal.One, ReducedArgumentError(a)); + private PreciseInterval Atan(PreciseInterval a) => AroundTheMiddle(a, x => x.Arctan(near), EDecimal.One, fixedPointError); + + private EDecimal ReducedArgumentError(PreciseInterval a) + => fixedPointError.Multiply(EDecimal.Max(EDecimal.One, a.Magnitude), up); private PreciseInterval Asin(PreciseInterval a) - => InsideTheUnitInterval(a) is { } slope ? AroundTheMiddle(a, x => x.Arcsin(near), slope) : new(EDecimal.NaN, EDecimal.NaN); + => InsideTheUnitInterval(a) is { } slope ? AroundTheMiddle(a, x => x.Arcsin(near), slope, fixedPointError) : new(EDecimal.NaN, EDecimal.NaN); private PreciseInterval Acos(PreciseInterval a) - => InsideTheUnitInterval(a) is { } slope ? AroundTheMiddle(a, x => x.Acos(near), slope) : new(EDecimal.NaN, EDecimal.NaN); + => InsideTheUnitInterval(a) is { } slope ? AroundTheMiddle(a, x => x.Acos(near), slope, fixedPointError) : new(EDecimal.NaN, EDecimal.NaN); /// The slope bound of the arcsine and arccosine, 1/sqrt(1 - x^2), strictly inside (-1, 1). private EDecimal? InsideTheUnitInterval(PreciseInterval a) @@ -247,7 +261,7 @@ private PreciseComplexInterval Log(PreciseComplexInterval a) }; var low = Min(corners); var high = Max(corners); - var slack = EDecimal.Max(low.Abs(), high.Abs()).Multiply(unit, up).Add(EDecimal.FromString("1E-" + (digits + 30)), up); + var slack = EDecimal.Max(low.Abs(), high.Abs()).Multiply(unit, up).Add(fixedPointError, up); return new(re, new(low.Subtract(slack, down), high.Add(slack, up))); } @@ -309,6 +323,13 @@ private PreciseInterval Pi() return new(pi.RoundToPrecision(down), pi.RoundToPrecision(up)); } + private PreciseInterval HalfPi() + { + var pi = Pi(); + var half = EDecimal.FromString("0.5"); + return new(pi.Low.Multiply(half, down), pi.High.Multiply(half, up)); + } + private PreciseInterval E() { var e = EDecimal.FromString( @@ -343,8 +364,8 @@ private PreciseComplexInterval Evaluate(Entity expr) return Multiply(Evaluate(left), Evaluate(right)); case Divf(var left, var right): return Divide(Evaluate(left), Evaluate(right)); - case Powf(var @base, Number.Integer power) when power.EInteger.CanFitInInt64(): - return Pow(Evaluate(@base), power.EInteger.ToInt64Unchecked()); + case Powf(var @base, Number.Integer power) when power.EInteger.CanFitInInt32(): + return Pow(Evaluate(@base), power.EInteger.ToInt32Checked()); case Powf(var @base, var exponent): return @base == MathS.e ? Exp(Evaluate(exponent)) : Pow(Evaluate(@base), Evaluate(exponent)); case Sinf(var argument): @@ -366,13 +387,19 @@ private PreciseComplexInterval Evaluate(Entity expr) case Cosecantf(var argument): return Divide(Real(PreciseInterval.Exactly(EDecimal.One)), Sin(Evaluate(argument))); case Logf(var @base, var antilogarithm): - return Divide(Log(Evaluate(antilogarithm)), Log(Evaluate(@base))); + return @base == MathS.e ? Log(Evaluate(antilogarithm)) : Divide(Log(Evaluate(antilogarithm)), Log(Evaluate(@base))); case Arctanf(var argument): return RealOnly(Evaluate(argument), Atan); case Arcsinf(var argument): return RealOnly(Evaluate(argument), Asin); case Arccosf(var argument): return RealOnly(Evaluate(argument), Acos); + case Arccotanf(var argument): + return RealOnly(Evaluate(argument), x => x.IsZero ? HalfPi() : Atan(Divide(PreciseInterval.Exactly(EDecimal.One), x))); + case Arcsecantf(var argument): + return RealOnly(Evaluate(argument), x => Acos(Divide(PreciseInterval.Exactly(EDecimal.One), x))); + case Arccosecantf(var argument): + return RealOnly(Evaluate(argument), x => Asin(Divide(PreciseInterval.Exactly(EDecimal.One), x))); case Absf(var argument): { var inner = Evaluate(argument); diff --git a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs index 24cf4bfa4..cae0fa263 100644 --- a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs +++ b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs @@ -33,6 +33,11 @@ public sealed class IntervalEvaluationTest [InlineData("arctan(2)")] [InlineData("arcsin(1/3)")] [InlineData("arccos(-2/3)")] + [InlineData("arccot(2)")] + [InlineData("arccot(-3/2)")] + [InlineData("arccot(0)")] + [InlineData("arcsec(3)")] + [InlineData("arccsc(-2)")] [InlineData("(2/3)^(1/3)")] [InlineData("(1 + 2i)^3")] [InlineData("1/(3 + 4i)")] @@ -74,6 +79,7 @@ private static void Holds(Interval interval, EDecimal value, string expression, [InlineData("1/0")] [InlineData("ln(0)")] [InlineData("tan(pi/2)")] + [InlineData("arcsec(1/2)")] [InlineData("3!")] [InlineData("1 provided 3 < 2")] public void UndecidedWhereItShouldBe(string expression) @@ -101,6 +107,8 @@ public void MoreDigitsSettleWhatDoublesCannot() [InlineData("sin(1 + i) * cos(2 - 3i)")] [InlineData("ln(-2 + i) + arctan(2) - arcsin(1/3)")] [InlineData("(2/3)^(1/3) + sec(1) + csc(2)")] + [InlineData("arccot(-2) + arcsec(3) + arccsc(-2) + arccot(0)")] + [InlineData("sin(10^6) + cos(10^6)")] [InlineData("sqrt(2) provided 1 < 2")] public void ThePreciseIntervalHoldsTheValue(string expression) { @@ -116,6 +124,18 @@ public void ThePreciseIntervalHoldsTheValue(string expression) $"the real part of {expression} is wider than forty digits should leave it"); } + /// + /// A whole power at the end of the range of a long is not negated into itself: the + /// reciprocal of 2^(2^63) was 2^(-2^63) again, and the recursion never ended. + /// + [Fact] + public void TheSmallestLongPowerFinishes() + { + var power = MathS.Pow(2, Entity.Number.Integer.Create(long.MinValue)); + Assert.Equal(0, IntervalEvaluation.Of(power)!.Value.Re.Mignitude); + Assert.Equal(0, PreciseEvaluation.Of(power, 40)!.Value.Re.Mignitude.Sign); + } + [Fact] public void AgreementHasThreeOutcomes() { From cb26ea5a9ba908da34d6616c93f2d454371b6cbd Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Sun, 27 Sep 2026 00:23:03 +0000 Subject: [PATCH 3/9] The evaluation decides membership in the reals and the complex numbers A value is real where its imaginary part is exactly zero, which real arithmetic keeps it, and not real where that part is certainly not zero. The pinned answer's derivative carries conditions of this kind on its constants, `-4 - 4 a^2 in RR`. Left undecided, they left the whole value undecided at every point, and the check rejected a correct answer: e^arsinh(a + b x)/x^3 came back unsolved, and four more of Rubi 7.1.4 and 7.1.5 timed out. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- Sources/AngouriMath/Numerics/IntervalEvaluation.cs | 12 ++++++++++++ Sources/AngouriMath/Numerics/PreciseEvaluation.cs | 9 +++++++++ .../Tests/UnitTests/Core/IntervalEvaluationTest.cs | 4 ++++ 3 files changed, 25 insertions(+) diff --git a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs index 6354f3668..b017b725a 100644 --- a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs +++ b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs @@ -462,6 +462,18 @@ private static ComplexInterval Evaluate(Entity expr) return false; return difference.Re.IsZero && difference.Im.IsZero ? true : null; } + // Real where the imaginary part is exactly zero, which real arithmetic keeps it; + // not real where it is certainly not zero. An answer's conditions ask this of its + // constants -- `-4 - 4 a^2 in RR` -- and legacy answers it. + case Set.Inf(var element, Set.SpecialSet.Reals): + { + var value = Evaluate(element); + if (!value.IsFinite) + return null; + return value.IsReal ? true : value.Im.Mignitude > 0 ? false : null; + } + case Set.Inf(var element, Set.SpecialSet.Complexes): + return Evaluate(element).IsFinite ? true : null; case Greaterf(var left, var right): return Compare(left, right) is { } greater ? greater > 0 : null; case GreaterOrEqualf(var left, var right): diff --git a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs index b27a1e1ca..a38ac1328 100644 --- a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs +++ b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs @@ -478,6 +478,15 @@ private PreciseComplexInterval RealOnly(PreciseComplexInterval argument, Func 0 ? false : null; + } + case Set.Inf(var element, Set.SpecialSet.Complexes): + return Evaluate(element).IsFinite ? true : null; case Greaterf(var left, var right): return Compare(left, right) is { } greater ? greater > 0 : null; case GreaterOrEqualf(var left, var right): diff --git a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs index cae0fa263..81c11c3f1 100644 --- a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs +++ b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs @@ -51,6 +51,8 @@ public sealed class IntervalEvaluationTest [InlineData("(sin(1)^2 + cos(1)^2) * 7/3")] [InlineData("sqrt(2) provided 1 < 2")] [InlineData("piecewise(1 provided 2 < 1, sin(1) provided 3 > 2)")] + [InlineData("sqrt(2) provided -4 - 4 * (137/100)^2 in RR")] + [InlineData("sqrt(2) provided sqrt(-2) in CC")] public void TheIntervalHoldsTheValue(string expression) { var expr = expression.ToEntity(); @@ -82,6 +84,7 @@ private static void Holds(Interval interval, EDecimal value, string expression, [InlineData("arcsec(1/2)")] [InlineData("3!")] [InlineData("1 provided 3 < 2")] + [InlineData("1 provided sqrt(-2) in RR")] public void UndecidedWhereItShouldBe(string expression) => Assert.Null(IntervalEvaluation.Of(expression.ToEntity())); @@ -110,6 +113,7 @@ public void MoreDigitsSettleWhatDoublesCannot() [InlineData("arccot(-2) + arcsec(3) + arccsc(-2) + arccot(0)")] [InlineData("sin(10^6) + cos(10^6)")] [InlineData("sqrt(2) provided 1 < 2")] + [InlineData("sqrt(2) provided -4 - 4 * (137/100)^2 in RR and sqrt(-2) in CC")] public void ThePreciseIntervalHoldsTheValue(string expression) { var expr = expression.ToEntity(); From c68e3ed131873fc96910b3c38c34c81b5f846f79 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Sun, 27 Sep 2026 02:22:12 +0000 Subject: [PATCH 4/9] The decimal evaluation reuses what it has worked out, and costs what the old one did On Rubi 6.6.7 the decimal tier made the check about ten times slower than the hundred-digit evaluation it replaced: 1/(a + b csch(c + d x)^2)^3 spent 3.9 s in five checks where the old evaluation spent 0.33 s. Four things, measured together to 4.9 s against 4.5 s for the whole integral: - The contexts are made once per precision. The library's constant cache and guard-digit contexts are keyed by the context instance, so pi and the logarithm tables were recomputed for every evaluation. - One evaluation serves every point of a check, remembering each node's value, so the subtrees free of the point -- which the substitutions share -- are worked out once. - The exponential and the logarithm are the library's fixed-point series, and the monotone functions are evaluated at the two ends rather than at the middle with a slope, which cost the exponential a second evaluation. - A product or quotient takes two directed operations where neither operand straddles zero, rather than eight. And an undefined operand stops the evaluation, and makes a product undefined even beside an exact zero: 0 * (1/0) was exactly 0. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Algebra/Polynomials/PartialFractions.cs | 26 ++-- .../Numerics/IntervalEvaluation.cs | 6 + .../AngouriMath/Numerics/PreciseEvaluation.cs | 131 +++++++++++++----- 3 files changed, 120 insertions(+), 43 deletions(-) diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs index 84a4229b0..53969e374 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs @@ -748,7 +748,7 @@ private static bool IsZeroAtPinnedSymbols(Entity expr) pinned = pinned.Substitute(symbol, value); i++; } - if (Agree(pinned, Integer.Zero, 1e-30) is not true) + if (Agree(pinned, Integer.Zero, 1e-30, new System.Lazy[] { new(() => PreciseEvaluation.In(40)), new(() => PreciseEvaluation.In(80)) }) is not true) return false; } return true; @@ -1896,8 +1896,9 @@ private static (Entity, Entity) Pinned(Entity left, Entity right, Variable x) /// /// Whether and agree, numerically, at a /// few points in with every other symbol pinned to a fixed value. - /// A point where either has no value -- a pole, a cut, a condition that fails -- or where - /// the evaluation cannot tell is skipped, and at least two must compare. The points may + /// A point where either has no value -- a pole, a cut, a condition that fails, a derivative + /// left unevaluated -- or where the evaluation cannot tell is skipped, and at least two must + /// compare. The points may /// be the caller's, for an identity that holds on a real domain only -- /// sqrt(1 - L^2) = sech(artanh(L)) inside (-1, 1) and off by a sign outside /// it -- where the default set would compare where it does not hold. @@ -1911,10 +1912,13 @@ internal static bool HoldsAtSampledPoints(Entity left, Entity right, Variable x, { (left, right) = Pinned(left, right, x); var compared = 0; + // The decimal evaluations are made once for all the points, so that the subtrees free + // of the point, which the substitutions share, are worked out once. + var inDecimals = new System.Lazy[] { new(() => PreciseEvaluation.In(40)), new(() => PreciseEvaluation.In(80)) }; foreach (var at in points ?? new[] { "0.29", "1.43", "3.17", "-0.61" }) { var point = Real.Create(EDecimal.FromString(at)); - switch (Agree(left.Substitute(x, point), right.Substitute(x, point), 1e-9)) + switch (Agree(left.Substitute(x, point), right.Substitute(x, point), 1e-9, inDecimals)) { case false: return false; @@ -1929,18 +1933,18 @@ internal static bool HoldsAtSampledPoints(Entity left, Entity right, Variable x, /// /// Whether the two numbers agree to of their size, or /// of one where both are smaller: in double intervals, and where those cannot tell, in - /// decimal intervals of forty and then eighty digits; null where either has no value or - /// none can tell. + /// the decimal evaluations of in turn -- forty and then + /// eighty digits; null where either has no value or none can tell. /// - private static bool? Agree(Entity left, Entity right, double relativeTolerance) + private static bool? Agree(Entity left, Entity right, double relativeTolerance, System.Lazy[] inDecimals) { if (IntervalEvaluation.Of(left) is { } leftInterval && IntervalEvaluation.Of(right) is { } rightInterval && IntervalEvaluation.Agree(leftInterval, rightInterval, relativeTolerance) is { } inDoubles) return inDoubles; - foreach (var digits in new[] { 40, 80 }) - if (PreciseEvaluation.Of(left, digits) is { } leftPrecise && PreciseEvaluation.Of(right, digits) is { } rightPrecise - && PreciseEvaluation.Agree(leftPrecise, rightPrecise, relativeTolerance) is { } inDecimals) - return inDecimals; + foreach (var precise in inDecimals) + if (precise.Value.ValueOf(left) is { } leftPrecise && precise.Value.ValueOf(right) is { } rightPrecise + && PreciseEvaluation.Agree(leftPrecise, rightPrecise, relativeTolerance) is { } decided) + return decided; return null; } } diff --git a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs index b017b725a..2e712d4d3 100644 --- a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs +++ b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs @@ -50,6 +50,10 @@ internal Interval(double low, double high) public static Interval operator *(Interval a, Interval b) { + // Unbounded before zero: nothing times no value has no value, and an interval that + // is not finite came from a divisor that may be zero. + if (!a.IsFinite || !b.IsFinite) + return new(double.NaN, double.NaN); if (a.IsZero || b.IsZero) return Exactly(0); var p1 = a.Low * b.Low; @@ -62,6 +66,8 @@ internal Interval(double low, double high) /// The quotient, or an interval that is not finite where the divisor may be zero. public static Interval operator /(Interval a, Interval b) { + if (double.IsNaN(a.Low) || double.IsNaN(a.High) || double.IsNaN(b.Low) || double.IsNaN(b.High)) + return new(double.NaN, double.NaN); if (a.IsZero && !b.ContainsZero) return Exactly(0); if (b.ContainsZero) diff --git a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs index a38ac1328..7ad8300b3 100644 --- a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs +++ b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs @@ -6,6 +6,7 @@ // using System; +using System.Collections.Generic; using PeterO.Numbers; using static AngouriMath.Entity; @@ -71,24 +72,51 @@ internal sealed class PreciseEvaluation private readonly PreciseComplexInterval undefined = new(new(EDecimal.NaN, EDecimal.NaN), new(EDecimal.NaN, EDecimal.NaN)); + // Values already worked out, by node: an expression substituted at several points shares + // every subtree free of the point, and each of those is worked out once. + private readonly Dictionary known = new(ByReference.Instance); + + // One set of contexts per precision, for the life of the process: the library's constant + // cache and its guard-digit contexts are keyed by the context instance, so contexts + // built afresh for every evaluation had pi and the logarithm tables recomputed each time. + [AngouriMath.Core.ConstantField] + private static readonly System.Collections.Concurrent.ConcurrentDictionary contexts = new(); + private PreciseEvaluation(int digits) { this.digits = digits; - down = EContext.ForPrecisionAndRounding(digits, ERounding.Floor); - up = EContext.ForPrecisionAndRounding(digits, ERounding.Ceiling); - near = EContext.ForPrecisionAndRounding(digits + 10, ERounding.HalfEven); + (down, up, near) = contexts.GetOrAdd(digits, static d => ( + EContext.ForPrecisionAndRounding(d, ERounding.Floor), + EContext.ForPrecisionAndRounding(d, ERounding.Ceiling), + EContext.ForPrecisionAndRounding(d + 10, ERounding.HalfEven))); unit = EDecimal.FromString("1E-" + (digits - 2)); fixedPointError = EDecimal.FromString("1E-" + (digits + 5)); least = EDecimal.FromString("1E-" + (digits + 30)); } + /// + /// An evaluation in digits that remembers what it has worked + /// out, for evaluating several expressions that share subtrees. + /// + internal static PreciseEvaluation In(int digits) => new(digits); + /// The interval evaluates to in digits, or null. - internal static PreciseComplexInterval? Of(Entity expr, int digits) + internal static PreciseComplexInterval? Of(Entity expr, int digits) => In(digits).ValueOf(expr); + + /// The interval evaluates to, or null. + internal PreciseComplexInterval? ValueOf(Entity expr) { - var value = new PreciseEvaluation(digits).Evaluate(expr); + var value = Evaluate(expr); return value.IsFinite ? value : null; } + private sealed class ByReference : IEqualityComparer + { + [AngouriMath.Core.ConstantField] internal static readonly ByReference Instance = new(); + public bool Equals(Entity? x, Entity? y) => ReferenceEquals(x, y); + public int GetHashCode(Entity obj) => System.Runtime.CompilerServices.RuntimeHelpers.GetHashCode(obj); + } + /// The three-way agreement of , in decimals. internal static bool? Agree(PreciseComplexInterval left, PreciseComplexInterval right, double relativeTolerance) { @@ -116,8 +144,20 @@ private PreciseInterval Add(PreciseInterval a, PreciseInterval b) private PreciseInterval Multiply(PreciseInterval a, PreciseInterval b) { + // Undefined before zero: nothing times no value has no value. + if (!a.IsFinite || !b.IsFinite) + return new(EDecimal.NaN, EDecimal.NaN); if (a.IsZero || b.IsZero) return PreciseInterval.Exactly(EDecimal.Zero); + // Where neither straddles zero, the ends are products of ends: two multiplications. + if (a.Low.Sign >= 0 && b.Low.Sign >= 0) + return new(a.Low.Multiply(b.Low, down), a.High.Multiply(b.High, up)); + if (a.High.Sign <= 0 && b.High.Sign <= 0) + return new(a.High.Multiply(b.High, down), a.Low.Multiply(b.Low, up)); + if (a.Low.Sign >= 0 && b.High.Sign <= 0) + return new(a.High.Multiply(b.Low, down), a.Low.Multiply(b.High, up)); + if (a.High.Sign <= 0 && b.Low.Sign >= 0) + return new(a.Low.Multiply(b.High, down), a.High.Multiply(b.Low, up)); var lows = new[] { a.Low.Multiply(b.Low, down), a.Low.Multiply(b.High, down), a.High.Multiply(b.Low, down), a.High.Multiply(b.High, down) }; var highs = new[] { a.Low.Multiply(b.Low, up), a.Low.Multiply(b.High, up), a.High.Multiply(b.Low, up), a.High.Multiply(b.High, up) }; return new(Min(lows), Max(highs)); @@ -125,13 +165,18 @@ private PreciseInterval Multiply(PreciseInterval a, PreciseInterval b) private PreciseInterval Divide(PreciseInterval a, PreciseInterval b) { - if (b.ContainsZero) + if (!a.IsFinite || !b.IsFinite || b.ContainsZero) return new(EDecimal.NaN, EDecimal.NaN); if (a.IsZero) return PreciseInterval.Exactly(EDecimal.Zero); - var lows = new[] { a.Low.Divide(b.Low, down), a.Low.Divide(b.High, down), a.High.Divide(b.Low, down), a.High.Divide(b.High, down) }; - var highs = new[] { a.Low.Divide(b.Low, up), a.Low.Divide(b.High, up), a.High.Divide(b.Low, up), a.High.Divide(b.High, up) }; - return new(Min(lows), Max(highs)); + // The divisor is of one sign, so the ends are quotients of ends: two divisions. + if (b.Low.Sign > 0) + return a.Low.Sign >= 0 ? new(a.Low.Divide(b.High, down), a.High.Divide(b.Low, up)) + : a.High.Sign <= 0 ? new(a.Low.Divide(b.Low, down), a.High.Divide(b.High, up)) + : new(a.Low.Divide(b.Low, down), a.High.Divide(b.Low, up)); + return a.Low.Sign >= 0 ? new(a.High.Divide(b.High, down), a.Low.Divide(b.Low, up)) + : a.High.Sign <= 0 ? new(a.High.Divide(b.Low, down), a.Low.Divide(b.High, up)) + : new(a.High.Divide(b.High, down), a.Low.Divide(b.High, up)); } private PreciseInterval Square(PreciseInterval a) @@ -163,42 +208,55 @@ private PreciseInterval AroundTheMiddle(PreciseInterval a, Func + /// A monotone function's value over the interval: its values at the two ends, computed + /// with ten digits to spare, each widened by the function's own error -- a unit or two in + /// the last place of the working precision, relative to the value, and, for a function + /// computed in fixed point, besides. + /// + private PreciseInterval AtTheEnds(PreciseInterval a, Func function, bool increasing, EDecimal absoluteError) + { + var atLow = function(a.Low); + var atHigh = a.Low.CompareTo(a.High) == 0 ? atLow : function(a.High); + if (!atLow.IsFinite || !atHigh.IsFinite) + return new(EDecimal.NaN, EDecimal.NaN); + var (low, high) = increasing ? (atLow, atHigh) : (atHigh, atLow); + var floor = EDecimal.Max(absoluteError, least); + return new(low.Subtract(low.Abs().Multiply(unit, up).Add(floor, up), down), + high.Add(high.Abs().Multiply(unit, up).Add(floor, up), up)); + } + + // The exponential and the logarithm are this library's fixed-point series -- a few + // microseconds, where PeterO's are hundreds -- within a unit of the value, the + // logarithm's error absolute near one. The square root is PeterO's. The sine, the cosine + // and the inverse functions are the library's too, in fixed point, and their error is + // absolute: the sine's grows with the argument, whose reduction by 2 pi uses pi to the + // same fixed point. private PreciseInterval Exp(PreciseInterval a) - => AroundTheMiddle(a, x => x.Exp(near), a.High.Exp(up).Abs(), EDecimal.Zero); + => AtTheEnds(a, x => x.Exponential(near), increasing: true, EDecimal.Zero); private PreciseInterval Log(PreciseInterval a) - => a.Low.Sign > 0 ? AroundTheMiddle(a, x => x.Log(near), EDecimal.One.Divide(a.Low, up), EDecimal.Zero) : new(EDecimal.NaN, EDecimal.NaN); + => a.Low.Sign > 0 ? AtTheEnds(a, x => x.NaturalLogarithm(near), increasing: true, fixedPointError) : new(EDecimal.NaN, EDecimal.NaN); private PreciseInterval Sqrt(PreciseInterval a) - => a.Low.Sign > 0 ? AroundTheMiddle(a, x => x.Sqrt(near), EDecimal.One.Divide(a.Low.Sqrt(down).Multiply(EDecimal.FromInt32(2), down), up), EDecimal.Zero) - : a.IsZero ? a : new(EDecimal.NaN, EDecimal.NaN); + => a.Low.Sign >= 0 ? AtTheEnds(a, x => x.Sqrt(near), increasing: true, EDecimal.Zero) : new(EDecimal.NaN, EDecimal.NaN); private PreciseInterval Sin(PreciseInterval a) => AroundTheMiddle(a, x => x.Sin(near), EDecimal.One, ReducedArgumentError(a)); private PreciseInterval Cos(PreciseInterval a) => AroundTheMiddle(a, x => x.Cos(near), EDecimal.One, ReducedArgumentError(a)); - private PreciseInterval Atan(PreciseInterval a) => AroundTheMiddle(a, x => x.Arctan(near), EDecimal.One, fixedPointError); + private PreciseInterval Atan(PreciseInterval a) => AtTheEnds(a, x => x.Arctan(near), increasing: true, fixedPointError); private EDecimal ReducedArgumentError(PreciseInterval a) => fixedPointError.Multiply(EDecimal.Max(EDecimal.One, a.Magnitude), up); private PreciseInterval Asin(PreciseInterval a) - => InsideTheUnitInterval(a) is { } slope ? AroundTheMiddle(a, x => x.Arcsin(near), slope, fixedPointError) : new(EDecimal.NaN, EDecimal.NaN); + => WithinTheUnitInterval(a) ? AtTheEnds(a, x => x.Arcsin(near), increasing: true, fixedPointError) : new(EDecimal.NaN, EDecimal.NaN); private PreciseInterval Acos(PreciseInterval a) - => InsideTheUnitInterval(a) is { } slope ? AroundTheMiddle(a, x => x.Acos(near), slope, fixedPointError) : new(EDecimal.NaN, EDecimal.NaN); + => WithinTheUnitInterval(a) ? AtTheEnds(a, x => x.Acos(near), increasing: false, fixedPointError) : new(EDecimal.NaN, EDecimal.NaN); - /// The slope bound of the arcsine and arccosine, 1/sqrt(1 - x^2), strictly inside (-1, 1). - private EDecimal? InsideTheUnitInterval(PreciseInterval a) - { - var largest = a.Magnitude; - if (largest.CompareTo(EDecimal.One) >= 0) - return null; - return EDecimal.One.Divide(EDecimal.One.Subtract(largest.Multiply(largest, up), down).Sqrt(down), up); - } + private static bool WithinTheUnitInterval(PreciseInterval a) + => a.Low.CompareTo(EDecimal.FromInt32(-1)) >= 0 && a.High.CompareTo(EDecimal.One) <= 0; private static EDecimal Min(EDecimal[] values) { @@ -338,6 +396,15 @@ private PreciseInterval E() } private PreciseComplexInterval Evaluate(Entity expr) + { + if (known.TryGetValue(expr, out var value)) + return value; + value = EvaluateNode(expr); + known[expr] = value; + return value; + } + + private PreciseComplexInterval EvaluateNode(Entity expr) { switch (expr) { @@ -357,13 +424,13 @@ private PreciseComplexInterval Evaluate(Entity expr) case Variable variable when variable == MathS.e: return Real(E()); case Sumf(var left, var right): - return Add(Evaluate(left), Evaluate(right)); + return Evaluate(left) is { IsFinite: true } lSumf ? Add(lSumf, Evaluate(right)) : undefined; case Minusf(var left, var right): - return Subtract(Evaluate(left), Evaluate(right)); + return Evaluate(left) is { IsFinite: true } lMinusf ? Subtract(lMinusf, Evaluate(right)) : undefined; case Mulf(var left, var right): - return Multiply(Evaluate(left), Evaluate(right)); + return Evaluate(left) is { IsFinite: true } lMulf ? Multiply(lMulf, Evaluate(right)) : undefined; case Divf(var left, var right): - return Divide(Evaluate(left), Evaluate(right)); + return Evaluate(left) is { IsFinite: true } lDivf ? Divide(lDivf, Evaluate(right)) : undefined; case Powf(var @base, Number.Integer power) when power.EInteger.CanFitInInt32(): return Pow(Evaluate(@base), power.EInteger.ToInt32Checked()); case Powf(var @base, var exponent): From 86b78cbe73ae74e08736b06bd7bbec8078194657 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Sun, 27 Sep 2026 02:23:39 +0000 Subject: [PATCH 5/9] The declined cases are declined in decimals too, and a pole beside a zero is one of them Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Tests/UnitTests/Core/IntervalEvaluationTest.cs | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs index 81c11c3f1..831b9fd17 100644 --- a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs +++ b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs @@ -73,8 +73,9 @@ private static void Holds(Interval interval, EDecimal value, string expression, } /// - /// Declined: a variable, a pole, the logarithm of zero, a node the pilot does not read, - /// and a value whose condition fails. + /// Declined, in doubles and in decimals: a variable, a pole, the logarithm of zero, a node + /// the evaluation does not read, a value whose condition fails -- and a product with a + /// pole in it, beside an exact zero or not. /// [Theory] [InlineData("x + 1")] @@ -85,8 +86,13 @@ private static void Holds(Interval interval, EDecimal value, string expression, [InlineData("3!")] [InlineData("1 provided 3 < 2")] [InlineData("1 provided sqrt(-2) in RR")] + [InlineData("0 * (1/0)")] + [InlineData("(1/0) * 0")] public void UndecidedWhereItShouldBe(string expression) - => Assert.Null(IntervalEvaluation.Of(expression.ToEntity())); + { + Assert.Null(IntervalEvaluation.Of(expression.ToEntity())); + Assert.Null(PreciseEvaluation.Of(expression.ToEntity(), 40)); + } /// /// Where cancellation leaves double intervals too wide to tell, forty digits tell: the From 8a3780200f5a673ccaea47f89d10ed35b9611059 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Sun, 27 Sep 2026 04:08:39 +0000 Subject: [PATCH 6/9] An exact zero factor makes a product zero again, whatever the other factor The decimal evaluation simplifies 0 x to 0 before it evaluates x, and an answer's derivative can carry a term whose coefficient is exactly zero at the pinned values beside a factor with no value there. Reading such a product as undefined, as the previous commit did, skipped points the old check compared: Rubi 6.6.7's 1/(a + b csch(c + d x)^2) came back unsolved, 6.4.7's (a + b coth(x)^2)^(3/2) tanh(x)^2 timed out, and 6.3.7's tanh(c + d x)^2/(a + b tanh(c + d x)^2) got an answer forty times as long. With the rule reverted all three are answered as on master. A product's operands are both evaluated again; a sum's, difference's and quotient's still stop at an undefined first operand. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Numerics/IntervalEvaluation.cs | 13 ++++++------ .../AngouriMath/Numerics/PreciseEvaluation.cs | 6 ++---- .../UnitTests/Core/IntervalEvaluationTest.cs | 21 +++++++++++++++---- 3 files changed, 26 insertions(+), 14 deletions(-) diff --git a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs index 2e712d4d3..97515e523 100644 --- a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs +++ b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs @@ -48,12 +48,15 @@ internal Interval(double low, double high) => b.IsZero ? a : a.IsZero ? -b : new(Down(a.Low - b.High), Up(a.High - b.Low)); public static Interval operator -(Interval a) => new(-a.High, -a.Low); + /// + /// The product, exactly zero where either factor is exactly zero, whatever the other -- + /// as the decimal evaluation has it, whose simplification takes 0 x to 0 + /// first. An answer's derivative can carry a term whose coefficient is exactly zero at + /// the pinned values beside a factor with no value there, and reading the term as having + /// none turned correct answers away: Rubi 6.6.7's 1/(a + b csch(c + d x)^2). + /// public static Interval operator *(Interval a, Interval b) { - // Unbounded before zero: nothing times no value has no value, and an interval that - // is not finite came from a divisor that may be zero. - if (!a.IsFinite || !b.IsFinite) - return new(double.NaN, double.NaN); if (a.IsZero || b.IsZero) return Exactly(0); var p1 = a.Low * b.Low; @@ -66,8 +69,6 @@ internal Interval(double low, double high) /// The quotient, or an interval that is not finite where the divisor may be zero. public static Interval operator /(Interval a, Interval b) { - if (double.IsNaN(a.Low) || double.IsNaN(a.High) || double.IsNaN(b.Low) || double.IsNaN(b.High)) - return new(double.NaN, double.NaN); if (a.IsZero && !b.ContainsZero) return Exactly(0); if (b.ContainsZero) diff --git a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs index 7ad8300b3..2e3e8d891 100644 --- a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs +++ b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs @@ -142,11 +142,9 @@ private PreciseInterval Add(PreciseInterval a, PreciseInterval b) private PreciseInterval Subtract(PreciseInterval a, PreciseInterval b) => Add(a, b.Negate()); + /// The product, exactly zero where either factor is, as 's is. private PreciseInterval Multiply(PreciseInterval a, PreciseInterval b) { - // Undefined before zero: nothing times no value has no value. - if (!a.IsFinite || !b.IsFinite) - return new(EDecimal.NaN, EDecimal.NaN); if (a.IsZero || b.IsZero) return PreciseInterval.Exactly(EDecimal.Zero); // Where neither straddles zero, the ends are products of ends: two multiplications. @@ -428,7 +426,7 @@ private PreciseComplexInterval EvaluateNode(Entity expr) case Minusf(var left, var right): return Evaluate(left) is { IsFinite: true } lMinusf ? Subtract(lMinusf, Evaluate(right)) : undefined; case Mulf(var left, var right): - return Evaluate(left) is { IsFinite: true } lMulf ? Multiply(lMulf, Evaluate(right)) : undefined; + return Multiply(Evaluate(left), Evaluate(right)); case Divf(var left, var right): return Evaluate(left) is { IsFinite: true } lDivf ? Divide(lDivf, Evaluate(right)) : undefined; case Powf(var @base, Number.Integer power) when power.EInteger.CanFitInInt32(): diff --git a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs index 831b9fd17..156955ddd 100644 --- a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs +++ b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs @@ -74,8 +74,7 @@ private static void Holds(Interval interval, EDecimal value, string expression, /// /// Declined, in doubles and in decimals: a variable, a pole, the logarithm of zero, a node - /// the evaluation does not read, a value whose condition fails -- and a product with a - /// pole in it, beside an exact zero or not. + /// the evaluation does not read, and a value whose condition fails. /// [Theory] [InlineData("x + 1")] @@ -86,8 +85,6 @@ private static void Holds(Interval interval, EDecimal value, string expression, [InlineData("3!")] [InlineData("1 provided 3 < 2")] [InlineData("1 provided sqrt(-2) in RR")] - [InlineData("0 * (1/0)")] - [InlineData("(1/0) * 0")] public void UndecidedWhereItShouldBe(string expression) { Assert.Null(IntervalEvaluation.Of(expression.ToEntity())); @@ -146,6 +143,22 @@ public void TheSmallestLongPowerFinishes() Assert.Equal(0, PreciseEvaluation.Of(power, 40)!.Value.Re.Mignitude.Sign); } + /// + /// An exact zero factor makes the product zero, whatever the other factor -- as the + /// decimal evaluation's simplification of 0 x has it. A derivative can hold a term + /// whose coefficient is exactly zero at the pinned values beside a factor with no value + /// there, and reading the product as undefined lost Rubi's + /// 1/(a + b csch(c + d x)^2). + /// + [Theory] + [InlineData("0 * (1/0)")] + [InlineData("(1/0) * 0")] + public void AnExactZeroFactorIsZero(string expression) + { + Assert.True(IntervalEvaluation.Of(expression.ToEntity()) is { Re.IsZero: true, Im.IsZero: true }); + Assert.True(PreciseEvaluation.Of(expression.ToEntity(), 40) is { Re.IsZero: true, Im.IsZero: true }); + } + [Fact] public void AgreementHasThreeOutcomes() { From c0890035a3e9f1fa2fc9bf45ccac2a8e6a117aff Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Sun, 27 Sep 2026 07:29:08 +0000 Subject: [PATCH 7/9] The changelog records the one answer the interval check changes Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 1d4a3c545..2b23ce021 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1501,6 +1501,25 @@ that came out in the expanded square before now come out in `x^2 + a/b`. Rubi's 1.2.2.2, 1.2.2.4, 1.2.2.7 and 1.2.3.2, the 574 problems that count with a perfect square written with symbols: 388 to 409, no row lost. +### The integrator's check evaluates in intervals + +The integrator accepts a candidate answer when its derivative agrees with the integrand at a few +points, with the symbols pinned to numbers. It evaluated both sides in hundred-digit decimals with +the downcasting off. There, a decimal below `1e-16` divides as zero and one below about `1e-199` +is zero, so a correct answer whose derivative passed through such a value was turned away and a +later route answered instead. The check now evaluates in intervals, which have neither limit and +read no setting ([#1019](https://github.com/asc-community/AngouriMath/issues/1019), item 10). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"e^(c*(a+b*x))*(cosh(a*c+b*c*x)^2)^(3/2)".Integrate("x")` | left unevaluated | the antiderivative, in `tanh(a c + b c x)` -- where `master` between the releases gave it in exponentials, the route the old check allowed | + +Rubi's suites against `master`: the 1774-problem suite, 4.7.2, 4.3.10 and 4.4.10, 7.1.4 to 7.2.5, +and samples of families 1, 4, 5, 6 and 7. That is 4056 problems, with one answer changed (the row +above). One problem is lost to the 5 s budget: `(c + d x^2)^(5/2)/(x^2 (a + b x^2)^2)`, which +`master` solved in 5.3 s, is now answered correctly by a slower route that the old check turned +away. + ### `binomial(n, k)` is a function **Addition, not silent.** The binomial coefficient is a node, `Entity.Binomialf`, spelled From 88ac2ae78f2e39763604272bfa19e17c36515129 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Mon, 28 Sep 2026 12:09:07 +0000 Subject: [PATCH 8/9] The interval check reads the special functions The error functions, the exponential and logarithmic integrals and the sine and cosine integrals are nodes on master now, and the evaluator had no case for them. An answer whose derivative keeps one -- `x erf(x)`, by parts, differentiates to `erf(x)` and more -- had no value at any point and would have been turned away. Each is read in the decimal tier: its value at the middle of the argument's rectangle, which the library computes with ten digits to spare, widened by a bound on the slope times the rectangle's reach from the middle. The slope is the derivative's magnitude in intervals over the rectangle, and all nine derivatives are elementary. Off the real line, a rectangle over the cut of Ei, Ci or Chi (left of 0) or of li (left of 1) is undefined, since the function jumps there. Co-Authored-By: Claude Opus 5.5 --- .../AngouriMath/Numerics/PreciseEvaluation.cs | 93 +++++++++++++++++++ .../UnitTests/Core/IntervalEvaluationTest.cs | 40 ++++++++ 2 files changed, 133 insertions(+) diff --git a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs index 2e3e8d891..b7b660495 100644 --- a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs +++ b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs @@ -492,11 +492,104 @@ private PreciseComplexInterval EvaluateNode(Entity expr) return undefined; } return undefined; + // The special functions, each through its derivative, which is elementary; Ei, Ci + // and Chi jump across the real line left of 0, and li left of 1. + case Erff(var argument): + return BySlope(Evaluate(argument), Number.Erf, z => Gaussian(z, negative: true), cutEnd: null); + case Erfcf(var argument): + return BySlope(Evaluate(argument), Number.Erfc, z => Gaussian(z, negative: true), cutEnd: null); + case Erfif(var argument): + return BySlope(Evaluate(argument), Number.Erfi, z => Gaussian(z, negative: false), cutEnd: null); + case Eif(var argument): + return BySlope(Evaluate(argument), Number.Ei, z => Divide(Exp(z), z), cutEnd: EDecimal.Zero); + case Lif(var argument): + return BySlope(Evaluate(argument), Number.Li, z => Divide(Real(PreciseInterval.Exactly(EDecimal.One)), Log(z)), cutEnd: EDecimal.One); + case Sif(var argument): + return BySlope(Evaluate(argument), Number.Si, z => Divide(Sin(z), z), cutEnd: null); + case Cif(var argument): + return BySlope(Evaluate(argument), Number.Ci, z => Divide(Cos(z), z), cutEnd: EDecimal.Zero); + case Shif(var argument): + return BySlope(Evaluate(argument), Number.Shi, z => Divide(Half(Subtract(Exp(z), Exp(Negate(z)))), z), cutEnd: null); + case Chif(var argument): + return BySlope(Evaluate(argument), Number.Chi, z => Divide(Half(Add(Exp(z), Exp(Negate(z)))), z), cutEnd: EDecimal.Zero); default: return undefined; } } + /// + /// A special function over the rectangle: its value at the middle, which the library + /// works out with ten digits to spare, widened by the most the function can move within + /// the rectangle, and by a unit in the last place of the working precision for its own + /// error. The most it can move is a bound on the slope times the farthest the rectangle + /// reaches from the middle, since along the segment from the middle to any point the + /// function changes by the integral of its derivative there. The bound is the + /// 's magnitude in intervals over the whole rectangle, and + /// where that has no value -- a singularity inside -- neither has the function. Off the + /// real line, a rectangle over the cut along the real line up to + /// is undefined as well: the function jumps there, and the + /// segment crosses the jump. On the real line itself there is no jump to cross, and the + /// value is the one the function takes on the cut. + /// + private PreciseComplexInterval BySlope(PreciseComplexInterval a, Func function, + Func derivative, EDecimal? cutEnd) + { + if (!a.IsFinite) + return undefined; + if (cutEnd is { } end && !a.IsReal && a.Im.ContainsZero && a.Re.Low.CompareTo(end) <= 0) + return undefined; + var two = EDecimal.FromInt32(2); + var middleRe = a.Re.Low.Add(a.Re.High, near).Divide(two, near); + var middleIm = a.Im.Low.Add(a.Im.High, near).Divide(two, near); + var halfRe = EDecimal.Max(a.Re.High.Subtract(middleRe, up), middleRe.Subtract(a.Re.Low, up)); + var halfIm = EDecimal.Max(a.Im.High.Subtract(middleIm, up), middleIm.Subtract(a.Im.Low, up)); + var reach = halfRe.Multiply(halfRe, up).Add(halfIm.Multiply(halfIm, up), up).Sqrt(up); + var moved = EDecimal.Zero; + // A point needs no slope, and Si's derivative at 0 has none written. + if (!reach.IsZero) + { + var slope = derivative(a); + if (!slope.IsFinite) + return undefined; + var steepest = slope.Re.Magnitude.Multiply(slope.Re.Magnitude, up) + .Add(slope.Im.Magnitude.Multiply(slope.Im.Magnitude, up), up).Sqrt(up); + moved = steepest.Multiply(reach, up); + } + // The library's functions read the working precision from the settings; they are + // given this evaluation's, and no downcasting, which would round the value. + Number.Complex value; + using (MathS.Settings.DecimalPrecisionContext.Set(near)) + using (MathS.Settings.DowncastingEnabled.Set(false)) + value = function(Number.Complex.Create(middleRe, middleIm)); + var (re, im) = (value.RealPart.EDecimal, value.ImaginaryPart.EDecimal); + if (!re.IsFinite || !im.IsFinite) + return undefined; + var floor = moved.Add(EDecimal.Max(fixedPointError, least), up); + var reReach = floor.Add(re.Abs().Multiply(unit, up), up); + var imReach = floor.Add(im.Abs().Multiply(unit, up), up); + // On the real line an imaginary part that is exactly zero at the middle is zero over + // the interval, which has no singularity in it: the function is real along the line + // there, or it would have jumped. + var imaginary = a.IsReal && im.IsZero ? PreciseInterval.Exactly(EDecimal.Zero) : new(im.Subtract(imReach, down), im.Add(imReach, up)); + return new(new(re.Subtract(reReach, down), re.Add(reReach, up)), imaginary); + } + + /// The error functions' derivative, 2/sqrt(pi) e^(-z^2), or e^(z^2) for erfi. + private PreciseComplexInterval Gaussian(PreciseComplexInterval z, bool negative) + { + var square = Multiply(z, z); + var twoOverRootPi = Divide(PreciseInterval.Exactly(EDecimal.FromInt32(2)), Sqrt(Pi())); + return Multiply(Real(twoOverRootPi), Exp(negative ? Negate(square) : square)); + } + + private PreciseComplexInterval Negate(PreciseComplexInterval a) => new(a.Re.Negate(), a.Im.Negate()); + + private PreciseComplexInterval Half(PreciseComplexInterval a) + { + var half = PreciseInterval.Exactly(EDecimal.FromString("0.5")); + return new(Multiply(a.Re, half), Multiply(a.Im, half)); + } + private PreciseComplexInterval RealOnly(PreciseComplexInterval argument, Func function) => argument.IsReal && argument.Re.IsFinite ? Real(function(argument.Re)) : undefined; diff --git a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs index 156955ddd..07680c75d 100644 --- a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs +++ b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs @@ -85,6 +85,10 @@ private static void Holds(Interval interval, EDecimal value, string expression, [InlineData("3!")] [InlineData("1 provided 3 < 2")] [InlineData("1 provided sqrt(-2) in RR")] + [InlineData("Ci(0)")] + [InlineData("li(1)")] + [InlineData("Ci(-2 + sin(pi) * i)")] + [InlineData("Ei(sin(pi))")] public void UndecidedWhereItShouldBe(string expression) { Assert.Null(IntervalEvaluation.Of(expression.ToEntity())); @@ -117,6 +121,19 @@ public void MoreDigitsSettleWhatDoublesCannot() [InlineData("sin(10^6) + cos(10^6)")] [InlineData("sqrt(2) provided 1 < 2")] [InlineData("sqrt(2) provided -4 - 4 * (137/100)^2 in RR and sqrt(-2) in CC")] + [InlineData("erf(1/2) + erfc(3) + erfi(-2/3)")] + [InlineData("erf(1 + i)")] + [InlineData("Ei(3/2) + Ei(-2)")] + [InlineData("Ei(1 - 2i)")] + [InlineData("li(1/2) + li(3)")] + [InlineData("li(-2)")] + [InlineData("Si(2) + Si(-3 + i)")] + [InlineData("Ci(3/10)")] + [InlineData("Ci(-3)")] + [InlineData("Ci(2 + 5i)")] + [InlineData("Ci(2 + sin(pi) * i)")] + [InlineData("Shi(3/2) + Chi(-1) + Chi(1 + 2i)")] + [InlineData("e^(-2) * Ei(sqrt(3)) + Si(pi/4)^2")] public void ThePreciseIntervalHoldsTheValue(string expression) { var expr = expression.ToEntity(); @@ -159,6 +176,29 @@ public void AnExactZeroFactorIsZero(string expression) Assert.True(PreciseEvaluation.Of(expression.ToEntity(), 40) is { Re.IsZero: true, Im.IsZero: true }); } + /// + /// The special functions are read, so an answer whose derivative keeps one -- by parts, + /// x erf(x) differentiates to erf(x) and more -- is checked, where reading + /// none of them left no point to compare at and the answer was turned away. On the real + /// line left of 0 the values are the ones on the cut: Chi(-0.61) is + /// Chi(0.61) + i pi, and li(-0.61) is complex. + /// + [Theory] + [InlineData("x * erf(x) + e^(-x^2)/sqrt(pi)", "erf(x)")] + [InlineData("x * erfi(x) - e^(x^2)/sqrt(pi)", "erfi(x)")] + [InlineData("x * Si(x) + cos(x)", "Si(x)")] + [InlineData("x * Ci(x) - sin(x)", "Ci(x)")] + [InlineData("x * Ei(x) - e^x", "Ei(x)")] + [InlineData("x * Chi(x) - sinh(x)", "Chi(x)")] + [InlineData("x * li(x) - Ei(2 * ln(x))", "li(x)")] + public void ADerivativeKeepingASpecialFunctionIsChecked(string antiderivative, string integrand) + { + var x = MathS.Var("x"); + var derivative = antiderivative.ToEntity().Differentiate(x); + Assert.True(Functions.PartialFractions.HoldsAtSampledPoints(derivative, integrand.ToEntity(), x)); + Assert.False(Functions.PartialFractions.HoldsAtSampledPoints(derivative, (integrand + " + 1/1000").ToEntity(), x)); + } + [Fact] public void AgreementHasThreeOutcomes() { From 539e91eea5a55e76a21f9ef759ae72b51e09f576 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 29 Sep 2026 02:53:41 +0000 Subject: [PATCH 9/9] The interval evaluation reads the Iverson bracket Master gained the Iverson bracket (#1567) after this branch last merged it, and neither tier had a case: an answer holding one had no value at any point. Both read it now, 1 where its statement is decided true over the interval and 0 where false, and leave it undecided otherwise, as a piecewise is. Co-Authored-By: Claude Opus 5.5 --- Sources/AngouriMath/Numerics/IntervalEvaluation.cs | 9 +++++++++ Sources/AngouriMath/Numerics/PreciseEvaluation.cs | 9 +++++++++ Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs | 3 +++ 3 files changed, 21 insertions(+) diff --git a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs index 97515e523..39448e488 100644 --- a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs +++ b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs @@ -416,6 +416,15 @@ private static ComplexInterval Evaluate(Entity expr) return Undefined; } return Undefined; + case Iversonf(var statement): + // 1 where the statement holds and 0 where it does not; one the intervals + // cannot settle leaves the value unsettled, as a piecewise does. + return Decide(statement) switch + { + true => ComplexInterval.Real(Interval.Exactly(1)), + false => ComplexInterval.Real(Interval.Exactly(0)), + null => Undefined, + }; default: return Undefined; } diff --git a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs index b7b660495..6261c61c7 100644 --- a/Sources/AngouriMath/Numerics/PreciseEvaluation.cs +++ b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs @@ -492,6 +492,15 @@ private PreciseComplexInterval EvaluateNode(Entity expr) return undefined; } return undefined; + // The Iverson bracket, 1 where its statement holds and 0 where it does not, and + // undefined where the statement is not decided over the rectangle, as a piecewise is. + case Iversonf(var statement): + return Decide(statement) switch + { + true => Real(PreciseInterval.Exactly(EDecimal.One)), + false => Real(PreciseInterval.Exactly(EDecimal.Zero)), + null => undefined, + }; // The special functions, each through its derivative, which is elementary; Ei, Ci // and Chi jump across the real line left of 0, and li left of 1. case Erff(var argument): diff --git a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs index 07680c75d..4e53787e4 100644 --- a/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs +++ b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs @@ -53,6 +53,7 @@ public sealed class IntervalEvaluationTest [InlineData("piecewise(1 provided 2 < 1, sin(1) provided 3 > 2)")] [InlineData("sqrt(2) provided -4 - 4 * (137/100)^2 in RR")] [InlineData("sqrt(2) provided sqrt(-2) in CC")] + [InlineData("iverson(2 > 1) * sqrt(2) + iverson(sin(1) < 0)")] public void TheIntervalHoldsTheValue(string expression) { var expr = expression.ToEntity(); @@ -89,6 +90,7 @@ private static void Holds(Interval interval, EDecimal value, string expression, [InlineData("li(1)")] [InlineData("Ci(-2 + sin(pi) * i)")] [InlineData("Ei(sin(pi))")] + [InlineData("iverson(sin(pi) > 0)")] public void UndecidedWhereItShouldBe(string expression) { Assert.Null(IntervalEvaluation.Of(expression.ToEntity())); @@ -134,6 +136,7 @@ public void MoreDigitsSettleWhatDoublesCannot() [InlineData("Ci(2 + sin(pi) * i)")] [InlineData("Shi(3/2) + Chi(-1) + Chi(1 + 2i)")] [InlineData("e^(-2) * Ei(sqrt(3)) + Si(pi/4)^2")] + [InlineData("iverson(pi > 3) * sqrt(2) + iverson(e < 2)")] public void ThePreciseIntervalHoldsTheValue(string expression) { var expr = expression.ToEntity();