PATCH: [BZ #14803] Different ULPs depending on size of long int in GCC
Joseph S. Myers
joseph@codesourcery.com
Mon Nov 5 13:30:00 GMT 2012
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
More information about the Libc-alpha
mailing list