mesa
mesa copied to clipboard
MESAHub
Need a bit-for-bit safe error function in math
erf (e.g. https://gcc.gnu.org/onlinedocs/gfortran/ERF.html) does not produce bit-for-bit identical results across compilers, so is not currently appropriate for general usage in MESA in the occasional circumstance where one might want to use it as part of a fitting formula.
See discussion in #399, where an error function approximation is used in the fitting formula for the diffusion coefficients. That approximation is fine for the specific diffusion coefficient application there, but is not accurate to floating point precision, so is not appropriate as a general solution to use in math.
Could we just remake https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libgfortran/intrinsics/erfc_scaled.c;h=60982fb161c1ac81300d34ad960ce16cdf5997cd;hb=HEAD ? Doesn't look to hard to make and should be bit for bit once we use crlibm's exp.
Scrap that i see we would need gcc's one as well https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/erfq.c;h=82a65ab1e1e4267774248760f277ed6748c49e4c;hb=HEAD as that is hairy
If you have not already, you might consider a series expansion such as the one here, which we used in this paper.
Cool! Any chance you know how many terms are needed for accuracy to roughly floating point precision?
Hrm, after playing around with it a bit, it seems like none of them converge at a practical rate to reach that precision over the entire domain. Short of re-implementing the code Rob linked, using a degree-28 Gauss-Legendre quadrature to calculate erf as an integral seems to give accurate enough results up to large enough arguments that erf~1 to floating point precision.
I stumbled past a StackOverflow answer that refers to formula 7.1.26 in the Handbook of Mathematical Functions, which looks simple and is precise to roughly single precision. This just looks like a slightly higher-order version of what's in #399 and gets about two more s.f. over that implementation.
I still think that almost all user routines will be fine with the (non-standard) compiler implementations of erf, so this is a minor issue until we really need erf in MESA's core.
