DawsonF[30.]
returns 0. The correct value is 0.016676...
At least it prints a warning message,
General::munfl: Exp[-900.] is too small to represent as a normalized
machine number; precision may be lost.
I am on Mathematica 11.3.
The problem is that DawsonF[x]
is being computed as Exp(-x^2) * Erfi[x]
(times constant factors), which in this case is a product of a very small quantity times a very large one, resulting in under/overflow. This is a VERY bad algorithm. The point of having a DawsonF
in the first place is to bypass this multiplication and return the result without under/overflow (see the section on Numerical Recipes book, for example).
I know I can use N[DawsonF[30], 20]
to obtain an accurate result, but this can be slower, and there is no reason why DawsonF
could not work in machine precision.
I will submit a bug report to Wolfram, but I wanted to post this here to get some feedback before. If the community agrees please tag this as a bug.
Are there other examples like this in Mathematica of special functions that just don't work in machine precision?
Update: I submitted a bug report and they replied (see my answer). See J.M.'s answer for a workaround.
Answer
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.
Comments
Post a Comment