[bug #64613] fp-contract=fast stops convergence in beta_inc_AXPY/beta_cont_frac

[Yannick Jadoul]

URL:
  <https://savannah.gnu.org/bugs/?64613>

                 Summary: fp-contract=fast stops convergence in
beta_inc_AXPY/beta_cont_frac
                   Group: GNU Scientific Library
               Submitter: 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
         Discussion Lock: Any
                 Release: 2.7.1


    _______________________________________________________

Follow-up Comments:


-------------------------------------------------------
Date: Wed 30 Aug 2023 04:58:52 PM UTC By: Yannick Jadoul <yannickjadoul>
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.







    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?64613>

_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/