diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 529a33eea..0f129d5d3 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -2121,6 +2121,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. + ### A quantifier over the whole numbers reads the bounds in its hypothesis, and residues written with `mod` `forall k in ZZ : (0 < k and k < 7) implies 7 divides binomial(7, k)` was left as written, because no diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 808e5fac1..2d7db79f0 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -708,5 +708,11 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Core/DowncastingOffComparisonTest.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 + [Tests/UnitTests/Calculus/HyperbolicIntegralIntegrationTest.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 cfa49a3f3..5eeb84532 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, new System.Lazy[] { new(() => PreciseEvaluation.In(40)), new(() => PreciseEvaluation.In(80)) }) 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); @@ -1905,48 +1909,58 @@ 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, 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. /// + /// + /// 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; + // 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 ?? DefaultSamplePoints) { 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, inDecimals)) { - 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 + /// 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, 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 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; + } + [ConstantField] private static readonly string[] DefaultSamplePoints = { "0.29", "1.43", "3.17", "-0.61" }; diff --git a/Sources/AngouriMath/Numerics/IntervalEvaluation.cs b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs new file mode 100644 index 000000000..39448e488 --- /dev/null +++ b/Sources/AngouriMath/Numerics/IntervalEvaluation.cs @@ -0,0 +1,563 @@ +// +// 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); + + /// + /// 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) + { + 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(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))); + } + + /// + /// 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. + /// + internal 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.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); + } + + /// 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); + // 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(); + 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 @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); + 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; + 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; + } + } + + /// + /// 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; + } + // 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): + 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..6261c61c7 --- /dev/null +++ b/Sources/AngouriMath/Numerics/PreciseEvaluation.cs @@ -0,0 +1,685 @@ +// +// 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 System.Collections.Generic; +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 EDecimal fixedPointError; + private readonly EDecimal least; + 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, 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) => In(digits).ValueOf(expr); + + /// The interval evaluates to, or null. + internal PreciseComplexInterval? ValueOf(Entity 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) + { + 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()); + + /// The product, exactly zero where either factor is, as 's is. + private PreciseInterval Multiply(PreciseInterval a, PreciseInterval b) + { + 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)); + } + + private PreciseInterval Divide(PreciseInterval a, PreciseInterval b) + { + if (!a.IsFinite || !b.IsFinite || b.ContainsZero) + return new(EDecimal.NaN, EDecimal.NaN); + if (a.IsZero) + return PreciseInterval.Exactly(EDecimal.Zero); + // 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) + { + 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. 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, EDecimal absoluteError) + { + var middle = a.Low.Add(a.High, near).Divide(EDecimal.FromInt32(2), near); + 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.Max(absoluteError, least), up); + return new(value.Subtract(reach, down), value.Add(reach, up)); + } + + /// + /// 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) + => AtTheEnds(a, x => x.Exponential(near), increasing: true, EDecimal.Zero); + + private PreciseInterval Log(PreciseInterval a) + => 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 ? 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) => 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) + => WithinTheUnitInterval(a) ? AtTheEnds(a, x => x.Arcsin(near), increasing: true, fixedPointError) : new(EDecimal.NaN, EDecimal.NaN); + + private PreciseInterval Acos(PreciseInterval a) + => WithinTheUnitInterval(a) ? AtTheEnds(a, x => x.Acos(near), increasing: false, fixedPointError) : new(EDecimal.NaN, EDecimal.NaN); + + 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) + { + 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(fixedPointError, 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 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( + "2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901"); + return new(e.RoundToPrecision(down), e.RoundToPrecision(up)); + } + + 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) + { + 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 Evaluate(left) is { IsFinite: true } lSumf ? Add(lSumf, Evaluate(right)) : undefined; + 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 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(): + 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): + 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 @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); + 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; + // 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): + 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; + + 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 Set.Inf(var element, Set.SpecialSet.Reals): + { + var value = Evaluate(element); + if (!value.IsFinite) + return null; + return value.IsReal ? true : value.Im.Mignitude.Sign > 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): + 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..4e53787e4 --- /dev/null +++ b/Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs @@ -0,0 +1,217 @@ +// +// 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("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)")] + [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)")] + [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(); + 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, in doubles and in decimals: a variable, a pole, the logarithm of zero, a node + /// the evaluation does not read, and a value whose condition fails. + /// + [Theory] + [InlineData("x + 1")] + [InlineData("1/0")] + [InlineData("ln(0)")] + [InlineData("tan(pi/2)")] + [InlineData("arcsec(1/2)")] + [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))")] + [InlineData("iverson(sin(pi) > 0)")] + public void UndecidedWhereItShouldBe(string expression) + { + 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 + /// 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("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")] + [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")] + [InlineData("iverson(pi > 3) * sqrt(2) + iverson(e < 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"); + } + + /// + /// 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); + } + + /// + /// 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 }); + } + + /// + /// 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() + { + 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)); + } + } +}