The interval evaluation reads the special functions (#1621)

The interval evaluation, the fast tier the integrator's checks run on (#1497), read no special
function. A check with erf, Ei or Si in it was undecided there and went on to the decimal
evaluation, which calls the arbitrary-precision kernels at 40 and then 80 digits. It now takes the
double-precision routines in AngouriMath.Numerics, as suggested on #1607.

Each value is enclosed as its routine's value at the rectangle's middle, widened by the routine's
error and by how far the function can move across the rectangle. The routines agree with the
kernels to within 3e-13 of the value, measured and not proved, and are allowed 1e-10 of it, and
1e-10 more near a zero. The movement is the largest the derivative is on the rectangle times its
half diagonal; the derivatives are elementary, so the intervals bound them themselves. A rectangle
that straddles a cut of Ei, li, Ci or Chi, or where the value or the derivative overflows a double,
has no enclosure and stays undecided, as before. A real value stays real where the function is.

Every enclosure was checked against the kernels at 30 digits on #1613's grid, about 17,800 per
function, each point as an exact double and nudged by a relative 1e-20: no kernel value falls
outside its enclosure. IntervalEvaluationTest gains 17 special-function rows.

Over Rubi's problems in families 0-7 whose antiderivative uses one of the nine functions, 2,354
in all, master solves 1975 and this 1976, with 0 wrong and 0 errors on both. The one move,
6.1.3:32 (sinh(a + b x^2)^3/x^3, answered in Chi and Shi), is a timeout on master. Run on its own
it takes 27 to 29 s on both builds, so the change does not bring it inside the budget by itself.

Part of #1607.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Author: Rafael-SOWNet

Sources/AngouriMath/Numerics/IntervalEvaluation.cs

@@ -8,6 +8,7 @@
 using System;
 using PeterO.Numbers;
 using static AngouriMath.Entity;
+using NumericsComplex = System.Numerics.Complex;
 
 namespace AngouriMath.Numerics
 {
@@ -382,6 +383,27 @@ private static ComplexInterval Evaluate(Entity expr)
                     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());
+                // The special functions, from their double-precision routines.
+                // https://github.com/asc-community/AngouriMath/issues/1607
+                case Erff(var argument):
+                    return Special(Evaluate(argument), SpecialFunctions.Erf, static z => TwoOverSqrtPi * (-(z * z)).Exp(), static _ => false, static _ => true);
+                case Erfcf(var argument):
+                    return Special(Evaluate(argument), SpecialFunctions.Erfc, static z => TwoOverSqrtPi * (-(z * z)).Exp(), static _ => false, static _ => true);
+                case Erfif(var argument):
+                    return Special(Evaluate(argument), SpecialFunctions.Erfi, static z => TwoOverSqrtPi * (z * z).Exp(), static _ => false, static _ => true);
+                case Eif(var argument):
+                    return Special(Evaluate(argument), SpecialFunctions.Ei, static z => z.Exp() / z, static z => StraddlesTheAxisLeftOf(z, 0), static _ => true);
+                case Lif(var argument):
+                    // Its cut is where the logarithm's is, and where Ei's is under it: 0 < x < 1.
+                    return Special(Evaluate(argument), SpecialFunctions.Li, static z => One / z.Log(), static z => StraddlesTheAxisLeftOf(z, 1), static x => x.Low >= 0);
+                case Sif(var argument):
+                    return Special(Evaluate(argument), SpecialFunctions.Si, static z => z.Sin() / z, static _ => false, static _ => true);
+                case Cif(var argument):
+                    return Special(Evaluate(argument), SpecialFunctions.Ci, static z => z.Cos() / z, static z => StraddlesTheAxisLeftOf(z, 0), static x => x.Low > 0);
+                case Shif(var argument):
+                    return Special(Evaluate(argument), SpecialFunctions.Shi, static z => Sinh(z) / z, static _ => false, static _ => true);
+                case Chif(var argument):
+                    return Special(Evaluate(argument), SpecialFunctions.Chi, static z => Cosh(z) / z, static z => StraddlesTheAxisLeftOf(z, 0), static x => x.Low > 0);
                 case Absf(var argument):
                 {
                     var inner = Evaluate(argument);
@@ -523,6 +545,57 @@ private static ComplexInterval Evaluate(Entity expr)
             return null;
         }
 
+        [AngouriMath.Core.ConstantField] private static readonly ComplexInterval One = ComplexInterval.Real(Interval.Exactly(1));
+        [AngouriMath.Core.ConstantField] private static readonly ComplexInterval TwoOverSqrtPi = ComplexInterval.Real(Interval.Around(2 / Math.Sqrt(Math.PI)));
+
+        private static ComplexInterval Sinh(ComplexInterval z) => (z.Exp() - (-z).Exp()) * ComplexInterval.Real(Interval.Exactly(0.5));
+        private static ComplexInterval Cosh(ComplexInterval z) => (z.Exp() + (-z).Exp()) * ComplexInterval.Real(Interval.Exactly(0.5));
+
+        /// <summary>Whether the rectangle crosses the real axis left of <paramref name="end"/>, where a cut runs.</summary>
+        private static bool StraddlesTheAxisLeftOf(ComplexInterval z, double end)
+            => !z.IsReal && z.Im.ContainsZero && z.Re.Low < end;
+
+        /// <summary>
+        /// A special function over <paramref name="argument"/>. The value is the double-precision
+        /// routine's at the rectangle's middle. It is widened by the routine's error, and by the most
+        /// the function can move across the rectangle: the largest its derivative is there, which
+        /// these intervals bound, times the rectangle's half diagonal. The value is real where the
+        /// argument is real and <paramref name="realOn"/> says the function is real there.
+        /// </summary>
+        /// <remarks>
+        /// The routines agree with the interpreter's kernels to within <c>3e-13</c> of the value,
+        /// measured on some ten thousand points rather than proved. They are allowed <c>1e-10</c> of
+        /// the value, and <c>1e-10</c> more for a value near a zero, where a relative error means
+        /// nothing. A rectangle that straddles a cut, where the function jumps, has no enclosure
+        /// this way. It is left undecided, and the caller asks the decimal evaluation.
+        /// </remarks>
+        private static ComplexInterval Special(ComplexInterval argument, Func<NumericsComplex, NumericsComplex> function,
+            Func<ComplexInterval, ComplexInterval> derivative, Func<ComplexInterval, bool> straddlesACut, Func<Interval, bool> realOn)
+        {
+            if (!argument.IsFinite || straddlesACut(argument))
+                return Undefined;
+            var middle = new NumericsComplex(argument.Re.Middle, argument.Im.Middle);
+            var value = function(middle);
+            if (!Interval.IsFiniteDouble(value.Real) || !Interval.IsFiniteDouble(value.Imaginary))
+                return Undefined;
+            var reach = Math.Max(argument.Re.High - middle.Real, middle.Real - argument.Re.Low);
+            var imReach = Math.Max(argument.Im.High - middle.Imaginary, middle.Imaginary - argument.Im.Low);
+            var radius = Interval.Up(Math.Sqrt(reach * reach + imReach * imReach));
+            var spread = 0.0;
+            if (radius > 0)
+            {
+                var slope = derivative(argument);
+                if (!slope.IsFinite)
+                    return Undefined;
+                spread = Interval.Up(Interval.Up(slope.Magnitude) * radius);
+            }
+            var error = Interval.Up(1e-10 * (Math.Sqrt(value.Real * value.Real + value.Imaginary * value.Imaginary) + 1) + spread);
+            var re = new Interval(Interval.Down(value.Real - error), Interval.Up(value.Real + error));
+            if (argument.IsReal && realOn(argument.Re))
+                return ComplexInterval.Real(re);
+            return new(re, new Interval(Interval.Down(value.Imaginary - error), Interval.Up(value.Imaginary + error)));
+        }
+
         /// <summary>A function read on real arguments only, in this pilot.</summary>
         private static ComplexInterval RealOnly(ComplexInterval argument, Func<Interval, Interval> function)
             => argument.IsReal && argument.Re.IsFinite ? ComplexInterval.Real(function(argument.Re)) : Undefined;

Sources/Tests/UnitTests/Core/IntervalEvaluationTest.cs

@@ -73,6 +73,52 @@ private static void Holds(Interval interval, EDecimal value, string expression,
                 $"the {part} part of {expression} is [{interval.Low:R}, {interval.High:R}], wider than it has to be");
         }
 
+        /// <summary>
+        /// The special functions, from their double-precision routines. Each is allowed <c>1e-10</c>
+        /// of its value, and <c>1e-10</c> more near a zero, against a measured error of at most
+        /// <c>3e-13</c>. An argument that is itself an interval, <c>sqrt(2)</c> or <c>pi</c>, widens the
+        /// value by the most the function moves across it. A real value stays real where the
+        /// function is real, and <c>Ci(-2)</c> is <c>Ci(2) + i pi</c>.
+        /// <a href="https://github.com/asc-community/AngouriMath/issues/1607">#1607</a>
+        /// </summary>
+        [Theory]
+        [InlineData("erf(1/2)")]
+        [InlineData("erf(sqrt(2))")]
+        [InlineData("erfc(3)")]
+        [InlineData("erfi(1 + i)")]
+        [InlineData("erf(2 - 3i)")]
+        [InlineData("Ei(2)")]
+        [InlineData("Ei(-3/2)")]
+        [InlineData("Ei(1 + 2i)")]
+        [InlineData("li(5)")]
+        [InlineData("li(1/3)")]
+        [InlineData("Si(pi)")]
+        [InlineData("Si(7)")]
+        [InlineData("Ci(3)")]
+        [InlineData("Ci(-2)")]
+        [InlineData("Shi(-1)")]
+        [InlineData("Chi(4)")]
+        [InlineData("erf(1/2) + Si(2) * Ei(1)")]
+        public void TheIntervalHoldsASpecialFunctionsValue(string expression)
+        {
+            var expr = expression.ToEntity();
+            var interval = IntervalEvaluation.Of(expr);
+            Assert.True(interval.HasValue, $"{expression} was not evaluated");
+            var value = expr.EvalNumerical();
+            HoldsWithin(interval!.Value.Re, value.RealPart.EDecimal, 1e-9, expression, "real");
+            HoldsWithin(interval.Value.Im, value.ImaginaryPart.EDecimal, 1e-9, expression, "imaginary");
+            if (value.ImaginaryPart.EDecimal.IsZero)
+                Assert.True(interval.Value.IsReal, $"{expression} is real, and its interval is not");
+        }
+
+        private static void HoldsWithin(Interval interval, EDecimal value, double width, 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}]");
+            Assert.True(interval.High - interval.Low <= width * System.Math.Max(1.0, interval.Magnitude),
+                $"the {part} part of {expression} is [{interval.Low:R}, {interval.High:R}], wider than it has to be");
+        }
+
         /// <summary>
         /// 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.