When floating point contractions (like fma/fmadd instructions) are enabled, for certain parameter values of gsl_cdf_fdist_Q the function returns NaN instead of the actual outcome of something like 0.05166etc

Take this simple example:

```
#include <gsl/gsl_cdf.h>
#include <stdio.h>

int main() {
    printf("%.15g\n", gsl_cdf_fdist_Q(3.786820954867802, 1, 100000));
    return 0;
}
```

- When compiling normally on a 64-bit intel (i.e., `-O2`), all goes well and a float of ~0.05 get returned.
- When adding `-march=native` (or any modern CPU architecture that has fused multiply-add fp contractions; `-ffp-contract=fast` is default), NaN gets returned. As far as I can tell, this is because the call to beta_cont_frac does not converge within the maximum of 512 iteration.
- Adding `-ffp-contract=off` (i.e., flags `-O2 -march=native -ffp-contract=off`) returns the previous non-NaN output again.


I realize this might be inherent to this FP contractions, and I have just fixed the issue for myself by building with `-ffp-contract=off`. But given that it's not just a minor FP imprecision but the difference between an actual outcome and NaN, I thought I would report.