Re: Proposal: Trigonometric functions in degrees

On Sat, Jan 23, 2016 at 12:08:40PM -0500, Tom Lane wrote:
> I wrote:
> > So the early returns from the buildfarm aren't very good:
> > * tern/sungazer isn't getting exactly 0.5 from sind(30).
>
> > The tern/sungazer result implies that this:
>
> > return (sin(x * (M_PI / 180.0)) / sin(30.0 * (M_PI / 180.0))) / 2.0;
>
> > is not producing exactly 0.5, which means that the two sin() calls aren't
> > producing identical results, which I suspect is best explained by the
> > theory that the compiler is rearranging 30.0 * (M_PI / 180.0) into
> > (30.0 * M_PI) / 180.0, and getting a slightly different number that way.
>
> > I think we could fix that by replacing (M_PI / 180.0) by a hard-wired
> > constant (computed to say 20 digits or so).
>
> So I pushed that, and tern/sungazer are still failing. Noah, could you
> trace through that and see exactly where it's going off the rails?

The second sin() is a constant, so gcc computes it immediately but sends the
first sin() to libm. The libm sin() is slightly more accurate. In %a
notation, AIX libm computes sin(30.0 * RADIANS_PER_DEGREE) as 0x1p-1 while gcc
computes it as 0x1.fffffffffffffp-2, a difference of one ULP. (Both "30.0 *
RADIANS_PER_DEGREE" and "30.0 * (M_PI / 180.0)" match the runtime computation
of 0x1.0c152382d7365p-1.)

To reliably produce exact answers, this code must delay all trigonometric
calculations to runtime. On sungazer, the float8 test happens to pass if I
rebuild float.c with -fno-builtin-sin; that leaves calls like acos(0.5) and
cos(60.0 * RADIANS_PER_DEGREE) unprotected, though.