[Bug math/27121] precision error on M_2_SQRTPI

[adhemerval.zanella at linaro dot org]

https://sourceware.org/bugzilla/show_bug.cgi?id=27121

Adhemerval Zanella <adhemerval.zanella at linaro dot org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |adhemerval.zanella at linaro dot o
                   |                            |rg
         Resolution|---                         |INVALID
             Status|UNCONFIRMED                 |RESOLVED

--- Comment #1 from Adhemerval Zanella <adhemerval.zanella at linaro dot org> ---
The M_2_SQRTPI is correctly rounded for Binary64/double precision, to get the
correctly rounded 2/sqrt(pi) for float128 you need to use M_2_SQRTPIf128 (it is
a gnu extension so you will need to build with _GNU_SOURCE) added on glibc
2.26.

The following example shows that glibc exported value has the correct precision
for each exported type:

---

$ cat m_2_sqrtpi.c 
#include <mpfr.h>
#include <stdio.h>
#include <math.h>
#include <quadmath.h>

int main (int argc, char *argv[])
{
  const int prec = 4096;
  const mpfr_rnd_t rnd = MPFR_RNDN;

  mpfr_t pi;
  mpfr_init2 (pi, prec);
  mpfr_const_pi (pi, rnd);

  mpfr_t sqrtpi;
  mpfr_init2 (sqrtpi, prec);
  mpfr_sqrt (sqrtpi, pi, rnd);

  mpfr_t two;
  mpfr_init2 (two, prec);
  mpfr_set_d (two, 2.0, rnd);

  mpfr_t ret;
  mpfr_init2 (ret, prec);
  mpfr_div (ret, two, sqrtpi, rnd);

  mpfr_printf ("M_2_SQRTPI (mpfr)      = %.128Ra\n", ret);
  double ret_d = mpfr_get_d (ret, rnd);

  printf ("M_2_SQRTPI (mpfr)      = %a\n", ret_d);
  printf ("M_2_SQRTPI (glibc)     = %a\n", M_2_SQRTPI);

  {
    __float128 ret_f128 = mpfr_get_float128 (ret, rnd);
    char buf[128];
    quadmath_snprintf (buf, sizeof buf, "%Qa", ret_f128);
    printf ("M_2_SQRTPIf128 (mpfr)  = %s\n", buf);
  }

  {
    __float128 ret_f128 = M_2_SQRTPIf128;
    char buf[128];
    quadmath_snprintf (buf, sizeof buf, "%Qa", ret_f128);
    printf ("M_2_SQRTPIf128 (glibc) = %s\n", buf);
  }

  return 0;
}
$ gcc -Wall m_2_sqrtpi.c -D_GNU_SOURCE -DMPFR_WANT_FLOAT128=1 -lquadmath -lmpfr
-o m_2_sqrtpi && ./m_2_sqrtpi 
M_2_SQRTPI (mpfr)      =
0x1.20dd750429b6d11ae3a914fed7fd8688281341d7587cea2e7342b06199cc416180eb39f0b24e1e2281806c12d98f35d77a3e9ddc91c394f0e9eedf0efffd84a3p+0
M_2_SQRTPI (mpfr)      = 0x1.20dd750429b6dp+0
M_2_SQRTPI (glibc)     = 0x1.20dd750429b6dp+0
M_2_SQRTPIf128 (mpfr)  = 0x1.20dd750429b6d11ae3a914fed7fep+0
M_2_SQRTPIf128 (glibc) = 0x1.20dd750429b6d11ae3a914fed7fep+0

-- 
You are receiving this mail because:
You are on the CC list for the bug.