[jmh530]

Partially addresses #11144. The broader discussion originated in #10318.

Add type-preserving float and double overloads for std.mathspecial.erf and erfc, using the corresponding Cephes kernels. Float calls use the native single-precision approximation. The existing real implementation is unchanged.

Implementation

Port the coefficients, approximation intervals, evaluation order, and underflow cutoffs from Cephes ndtr.c and ndtrf.c, with the double exponential-of-square calculation from expx2.c. Coefficients are reversed for Phobos' polynomial evaluator; implicit leading denominator coefficients are made explicit. Replace Cephes error reporting with Phobos return-value behavior and preserve @safe pure nothrow @nogc.

Source archives: https://netlib.org/cephes/cprob.tgz and https://netlib.org/cephes/single.tgz.

Compatibility and numerical limitations

Float/double arguments now return their own type, changing auto inference and available precision/range. Cast arguments to real to retain extended-precision behavior. Explicit integral overloads preserve real results for integral arguments and integer enums. Floating-point enums require an explicit cast to select an overload. A release note documents these changes.

The port deliberately retains the original Cephes limitations. In particular, erfc(10.0f) and erfc(27.0) return zero although their correctly rounded target-type results are nonzero subnormals. Improving tail handling is deferred to a follow-up.

The single-precision kernel also trades accuracy for speed compared with computing in real and narrowing. On the sampled positive normal-result grid with Windows LDC, maximum relative errors were approximately 1.23e-7 for float erf, 3.80e-6 for float erfc, 2.61e-16 for double erf, and 9.15e-16 for double erfc. These are measured sample maxima, not exhaustive bounds; zero results caused by the inherited cutoff are excluded from those relative-error figures.

Performance

Local microbenchmarks compare the new kernels against the original real implementation narrowed to the same target type. The table shows the range of median paired reductions in execution time for both functions, excluding early-underflow tail cases.

Compiler/platform	Float: inputs 0.5, 2, 5	Double: inputs 0.5, 2, 5, 10
DMD 2.113, Windows x64	69–80%	45–68%
LDC 1.43, Windows x64	62–72%	13–69%
LDC 1.36, Linux x64 (WSL)	88–97%	88–93%

The additional sampled tail inputs (float 9.5, 10, 11; double 27) were faster too, but those gains partly reflect returning zero early rather than computing subnormal erfc results. No sampled case was slower than the baseline.

Each case used seven rounds of 300,000 indirect calls, discarded the first round, rotated old/new function order, and consumed results. Flags were DMD -O -inline and LDC -O2, without fast-math. Host load varied; these results are not portable speed guarantees.

Validation

Focused optimized tests passed on Windows DMD x64 and x86, Windows LDC for the affected modules, and Linux LDC. Coverage includes return types, integral and enum calls, signed zero, infinities, NaNs, reference values, float approximation boundaries, and inherited underflow behavior.

The broad numerical harness samples 6,456 float/double outputs per target and separately checks that real results remain identical to the original implementation. Independent high-precision checks characterize numerical error and the inherited cutoff limitations.

The full repository suite and binary128 runtime have not been tested locally. The previously observed Windows LDC assertion in the gamma module also occurs on the unchanged baseline; the affected-module tests pass.

This PR was written entirely by an LLM/AI.

AI-Assisted-By: Codex (GPT-6 Astra Medium)

[ibuclaw]

Do we really need this?

And if so, can't this be done using constraints (optionally shorthand syntax)

real erf(I)(I x) if (isIntegral!I) => erf(cast(real) x);

[jmh530]

I checked this with Codex (GPT-6 Astra Medium) on DMD 2.113 and LDC 1.43. The compatibility overloads are needed to keep existing calls such as erf(1) compiling. Unfortunately, the constrained template still leaves that call ambiguous alongside the non-template floating-point overloads:

import std.traits : isIntegral;

real f(real x) { return x; }
double f(double x) { return x; }
float f(float x) { return x; }
real f(I)(I x) if (isIntegral!I) { return f(cast(real) x); }

void main()
{
    auto x = f(1); // Error: matches multiple overloads after implicit conversions
}

Explicitly calling f!int(1) works, but would require callers to change their code. Broadening the constraint to include bool and character types does not resolve the ambiguity either; those types are excluded by isIntegral but were accepted by the previous real-only API.

I would therefore retain the non-template compatibility overloads generated by static foreach. The existing tests cover the built-in types and their integer enums.

[ibuclaw]

Accepting bool and character types in a math function is a bit weird though

[jmh530]

I agree. Those were mainly kept to avoid breaking changes. Below compiles currently (or changing the bools to char). After making this change, it would be easier to deprecate that behavior in a future PR if that's the approach phobos wants to take.

import std.mathspecial;
import std.stdio: writeln;

void main() {
    bool x0 = 0;
    auto y0 = erf(x0);
    auto z0 = erfc(x0);
    writeln(y0); //prints 0
    writeln(z0); //prints 1
    bool x1 = 1;
    auto y1 = erf(x1);
    auto z1 = erfc(x1);
    writeln(y1); //prints ~0.842701
    writeln(z1); //prints ~0.157299
}

[jmh530]

Hmm, I guess if you had an int overload (or whatever), then it should also take a bool input.