Terrible accuracy of DawsonF

Before DawsonF[] became built-in in Mathematica, I used the following method for (small to moderately-sized) real arguments:

dawson = With[{eps = $MachineEpsilon, e2 = $MachineEpsilon^2},
              Compile[{{z, _Real}},
                      Module[{a, b, c, d, f, h, w},
                             a = 2. z^2;
                             f = c = b = a + 1.;
                             a = w = -2. a; d = 0.;
                             If[c == 0., c = e2];
                             While[b += 2.;
                                   d = b + a d; If[d == 0., d = e2]; d = 1/d;
                                   c = b + a/c; If[c == 0., c = e2];
                                   a += w; f *= (h = c d);
                                   Abs[h - 1] > eps];
                             z/f], RuntimeAttributes -> {Listable}]];

This is based on using the Lentz-Thompson-Barnett algorithm to evaluate this CF representation for Dawson's integral. It is not the most efficient method, since power series will outperform the CF near the origin, and asymptotic series will do best for really large arguments. It is quite compact and respectably performant in its intended argument range, however.

Here is a plot of it compared against the built-in:

Plot[dawson[x] - DawsonF[x], {x, -20, 20}, Frame -> True, PlotRange -> All]

I submitted a bug report and this is the relevant line from what they replied:

Starting with Version 11.3 underflow in no longer trapped in machine arithmetic and Mathematica does not switch automatically to arbitrary precision. This provides a more efficient way to handle numerical calculations and brings Mathematica much more in line with the IEEE 754 standard for how floating point numbers are to be handled ( https://en.wikipedia.org/wiki/IEEE_754 )

This explains the origin of the bug. Apparently the bad machine precision algorithm was there all along, but in previous versions it fell back to arbitrary precision and thus went unnoticed (though it probably impacted performance).

Hope it gets fixed soon.

See @J.M. answer for a an algorithm that works in MachinePrecision.

Based on @EddieXiao's answer to Numerical underflow for a scaled error function, you could define your own DawsonF with:

dawsonF[x_] := -(I/2) E^-x^2 Sqrt[π]+I HermiteH[-1,I x]

For example:

Chop @ dawsonF[30.] //InputForm

General::munfl: Exp[-900.] is too small to represent as a normalized machine number; precision may be lost.

0.01667594140105917

vs.

DawsonF[30`19]

0.01667594140105918

This is faster than the built-in machine precision code for DawsonF, e.g.:

dawsonF[N @ Range[-10, 10]]; //AbsoluteTiming
DawsonF[N @ Range[-10, 10]]; //AbsoluteTiming

{0.005917, Null}

{0.008466, Null}

On the other hand, @JM's compiled code is about 4 times faster:

dawson[N @ Range[30, 40]] //AbsoluteTiming
Chop @ dawsonF[N @ Range[30, 40]] //AbsoluteTiming

{0.000078, {0.0166759, 0.0161374, 0.0156326, 0.0151585, 0.0147123, 0.0142916, 0.0138943, 0.0135185, 0.0131625, 0.0128247, 0.0125039}}

{0.000323, {0.0166759, 0.0161374, 0.0156326, 0.0151585, 0.0147123, 0.0142916, 0.0138943, 0.0135185, 0.0131625, 0.0128247, 0.0125039}}

Of course, @JM's compiled code is expecting real numbers, and so complex input doesn't use the compiled code:

dawson[3. + 2. I] //AbsoluteTiming
dawsonF[3. + 2. I] //AbsoluteTiming

CompiledFunction::cfsa: Argument 3. +2. I at position 1 should be a machine-size real number.

{0.001332, 0.110514 - 0.0771238 I}

{0.001298, 0.110514 - 0.0771238 I}

The non-compiled code is about the same speed as dawsonF. Finally, as @JM mentions, his compiled code is not meant to work for all possible arguments, e.g.:

dawson[3 + 5 I]
dawsonF[3. + 5. I]
N[DawsonF[3 + 5 I], 20]

CompiledFunction::cfsa: Argument 3+5 I at position 1 should be a machine-size real number.

14.7334 + 6.88742 I

-7.78086*10^6 + 1.21475*10^6 I

-7.7808580812920342136*10^6 + 1.2147471245770455984*10^6 I