The special functions compile, in double precision

[Rafael-SOWNet]

erf, erfc, erfi, Ei, li, Si, Ci, Shi and Chi had no compiled form. So Compile and Compile<TIn, TOut> threw UncompilableNodeException on them, and Solve's numerical search, which compiles, declined any equation with one in it (#1605). The interpreter's kernels are arbitrary-precision and take 0.7 to 5.5 ms a call at 20 digits. Newton's method makes about 3,000 calls an equation, so they cannot back a compiled form: #1607, suggested on #1605.

SpecialFunctions, in AngouriMath.Numerics beside the interval evaluation (moved there in the second commit, as suggested on #1607), gives each a double-precision routine, which both compilers call:

• Error functions: from the Faddeeva function, by Weideman's rational approximation with 40 terms (SIAM J. Numer. Anal. 31, 1994), and from erf's Taylor series near the origin.
• Ei: its series where the terms cancel little (|z| - Re z <= 3, or |z| <= 2), E1's continued fraction elsewhere out to |z| = 40, and its asymptotic series beyond.
• li, Si, Ci, Shi, Chi: from Ei and from their own series, by the relations the kernels state.

Each keeps the kernel's branch cuts and the kernel's values on the axes. An imaginary part of either sign of zero is read as the axis, as the kernels read it. Under Compile<TIn, TOut> the definitions are for Complex and double, and an integer or a float is widened to double first, as for the trigonometric functions. Over doubles, a value off the real line is NaN, as Math.Sqrt of a negative number is.

Solve now answers erf(x) = x/2 + 1/7 with its three real roots, each checked against mpmath, and more off the real line. An equation that isolates a special function, such as erf(x) = 1/2, is still left unsolved: the solver finds no inverse to write, so it never reaches the search. binomial(x, 3/2) = x/2 replaces erf as the test's example of what the search still declines.

Measured

• Accuracy. Compared with the kernels at 30 digits, through the compiled delegates, on about ten thousand points: a grid of [-14, 14]^2, both axes out to 40 and in to 1e-4, either side of the negative real axis at 1e-9 and 1e-3, and circles out to radius 300.
  • The integrals agree to within 6e-14 of the value.
  • The error functions agree to within 3e-13, which is the rounding of z^2 in their e^(-z^2) at large |z|.
  • Where the value overflows a double, the large part is infinite and the part the axis makes exact keeps its value, as the interpreter's does.
• Time: 0.5 to 2.4 µs a call over [-10, 10]^2, against the kernels' 0.7 to 5.5 ms.
• SpecialFunctionCompilationTest: 57 points through both compilers against the interpreter to 1e-12; real values and NaN over doubles; and five infinite values, whose exact part a complex halving had made NaN.
• NewtonFallbackTest: erf(x) = x/2 + 1/7 and Si(x) = x/2 + 1/7 are searched, each with its three real roots from mpmath, and every root it returns is checked. SolveNt still throws for what doesn't compile.
• BREAKING-CHANGES.md: the three entries that added the functions gain a row each, since 2.5.0 had none of them. The floors, the rounding and the extremes compile, and Solve's numerical search declines what it cannot represent #1605's entry names binomial in place of erf.
• The unit suite on 7290956 (net10.0): 14051 passed, 0 failed. The library builds for every target framework.
• Performance gate on 7290956: passed, with allocation within 0.6% of the baseline on all 19 gated benchmarks. FastExpression's interpreter gains nine cases.

Part of #1607. x!, gamma and phi under Compile<TIn, TOut> are the other half, next.

🤖 Generated with Claude Code