GNU Scientific Library - Bugs: bug #64613, fp-contract=fast stops convergence...
You are not allowed to post comments on this tracker with your current authentication level.
bug #64613: fp-contract=fast stops convergence in beta_inc_AXPY/beta_cont_frac
Submitter: | Yannick Jadoul <yannickjadoul> | ||
Submitted: | Wed 30 Aug 2023 04:58:52 PM UTC | ||
Category: | Accuracy problem | Severity: | 3 - Normal |
Operating System: | Linux 5.15.0-82-generic #91-Ubuntu SMP Mon Aug 14 14:14:14 UTC 2023 x86_64 x86_64 x86_64 GNU/Linux | Status: | None |
Assigned to: | None | Open/Closed: | Open |
Release: | 2.7.1 |
No files currently attached
Depends on the following items: None found
Items that depend on this one: None found
Carbon-Copy List
No changes have been made to this item
Powered by Savane 3.13-02a9.
Corresponding source code
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.