Dear GSL developers, I'm working on a project where we use the gsl_sf_log_erfc function to calculate the logarithm of the complementary error function to high precision. After a while, we noticed that for large numbers, gsl_sf_log_erfc incorrectly returns -NAN. (The complementary error function of a large number is an extremely small positive number, so its logarithm should be rounded to negative infinity.) The bug can be easily demonstrated by printing the result of gsl_sf_log_erfc(HUGE_VAL). Interestingly, log(erfc(HUGE_VAL)) correctly returns -HUGE_VAL, but the loss of precision in other cases is unacceptable for us.
A similar problem occurs with gsl_sf_log_erfc(-HUGE_VAL), which returns -NAN instead of 0.693 (the natural log of 2). Once again, log(erfc(-HUGE_VAL)) gives the correct answer. Both of these problems with gsl_sf_log_erfc occur with any input of magnitude approximately 1e62 or greater. Is this something that can be fixed in GSL? It'd be nice to not have to detect and correct for erroneous results. -Alex
