Re: PATCH: [BZ #14803] Different ULPs depending on size of long intin GCC

["Joseph S. Myers" <joseph at codesourcery dot com>]

On Mon, 5 Nov 2012, H.J. Lu wrote:

> I don't know how to compute  the result of rounding pi/2 - pio2_hi, to
> nearest in C. I simply replaced

Using MPFR,

#include <stdio.h>
#include <mpfr.h>

int
main (void)
{
  mpfr_t pi2_64, pi2_lo, pi2_lo_64, pi4_64;
  mpfr_init2 (pi2_64, 64);
  mpfr_init2 (pi2_lo, 300);
  mpfr_init2 (pi2_lo_64, 64);
  mpfr_init2 (pi4_64, 64);
  mpfr_const_pi (pi2_64, MPFR_RNDN);
  mpfr_div_ui (pi2_64, pi2_64, 2, MPFR_RNDN);
  mpfr_div_ui (pi4_64, pi2_64, 2, MPFR_RNDN);
  mpfr_const_pi (pi2_lo, MPFR_RNDN);
  mpfr_div_ui (pi2_lo, pi2_lo, 2, MPFR_RNDN);
  mpfr_sub (pi2_lo_64, pi2_lo, pi2_64, MPFR_RNDN);
  mpfr_printf ("%Ra %Ra %Ra\n", pi2_64, pi2_lo_64, pi4_64);
  return 0;
}

gives:

0x1.921fb54442d1846ap+0 -0x7.6733ae8fe47c65d8p-68 0xc.90fdaa22168c235p-4

> pio2_lo = 2.9127320560933561582586004641843300502121E-20L,

FWIW, it looks like this is the low part you get if you take the high part 
as being pi/2 rounded to 65 bits towards zero (not to nearest), then write 
pi/2 minus that high part to 40 decimal places after the '.' (as opposed 
to writing a 64-bit or 65-bit approximation of that difference to 40 
decimal places).

I see no apparent reason in this code to need pio2_lo to be positive, so 
think rounding to nearest (giving the values I give above) is appropriate 
here.

> +  pio2_lo = 0x8.98cc51701b839a2p-68,

In any case, missing 'L' suffix here; all three constants need such a 
suffix.

-- 
Joseph S. Myers
joseph@codesourcery.com