Re: Proposal: Trigonometric functions in degrees

Dean Rasheed <dean(dot)a(dot)rasheed(at)gmail(dot)com> writes:
> If I understand correctly there were 2 separate issues at play here:

> 1). On some platforms the compiler evaluates expressions like
> sin(constant) and comes up with a slightly different result than a
> runtime evaluation of the expression. The compiler-evaluated result
> is presumably a 64-bit IEEE float, but at runtime it may be using
> extended precision for intermediate results.

If I've not lost track of the bidding, both of the cases where we saw
that involved gcc on a platform where it's not the native (vendor
supplied) compiler. So my guess is that gcc was using a non-native
libm to do the pre-evaluation of sin(). It does not seem likely that
the gcc boys would intentionally do pre-evaluation differently from
run-time evaluation, but they could get burnt by what would arguably
be a build error in that particular copy of gcc. Cross-compilation
would be another way to hit the same type of hazard.

> That may well have been the sole contributing factor to the fact that
> sind(30) wasn't exactly 0.5.

Yeah. I am not sure whether the RADIANS_PER_DEGREE change fixed anything
or not, though I am not tempted to undo it; it may save us on some other
platform not currently represented in the buildfarm.

> 2). The compiler also sometimes rearranges expressions in ways that
> don't give the same result as evaluating in the order suggested by the
> parentheses. I think this actually explains the failure to get exactly
> 1 for tand(45). For x=45, this was being computed as
> cosd_0_to_60(90 - x) / cosd_0_to_60(x)
> so my guess is that it was inlining cosd_0_to_60(90 - x) and
> rearranging it to produce something different from cosd_0_to_60(x) for
> x=45.

Oh, interesting point. The inlining would have produced a subexpression
like

cos((90 - x) * RADIANS_PER_DEGREE)

For x=45, the result of 90-x would have been exact, so it's not obvious
where any change in results would have crept in --- but if the compiler
then tried to simplify to

cos((90 * RADIANS_PER_DEGREE) - (x * RADIANS_PER_DEGREE))

that could definitely change the roundoff behavior. OTOH, it's not
very clear why gcc would have done that; it saves no operations.
It'd be interesting to look at the produced assembly code on narwhal.

> I wonder if the same could have been achieved by disabling
> optimisation and inlining in those low-level functions, and also
> wrapping sin(x * RADIANS_PER_DEGREE) in a similar non-inlinable,
> non-optimised function to force it to be executed at runtime when
> passed a constant.

I considered that; in particular -ffloat-store would have helped with the
wider-intermediate-results problem, and indeed we might still be forced
into using that. I would rather not fall back to adding more
compiler-specific flags though. If we go that way we'll likely need a
custom fix for every new compiler we try to use.

Meanwhile, just when you thought it was safe to go back in the water,
cockatiel is still failing. It has the cos(60) != 0.5 problem, which
IIRC was exhibited by no other critter. Looking at the code,

cosd_0_to_60(double x)
{
return 1.0 - ((1.0 - cos(x * RADIANS_PER_DEGREE)) / one_minus_cos_60) / 2.0;
}

what seems likely is that the "1 - cos()" subtraction is being done in
a wider-than-double float register, which i686 does have, and producing
a different result than what was stored in one_minus_cos_60.

Perhaps we can fix this by rewriting as

float8 numerator = 1.0 - cos(x * RADIANS_PER_DEGREE);
return 1.0 - (numerator / one_minus_cos_60) / 2.0;

cockatiel's compiler does recognize -fexcess-precision=standard, and
my understanding of that is that the result put into "numerator" will
be rounded to double width, so that it should then match
"one_minus_cos_60".

Another idea would be to change the cache variable to just "cos_60" and
write "(1.0 - cos_60)" in the denominator --- but then we'd just be hoping
that the compiler does both subtractions the same way, which doesn't seem
necessarily guaranteed. Worse, I believe the 8087 has an FCOS instruction
which might deliver a wider-than-double result, so that maybe the problem
is not so much with the subtraction as with when rounding of the cos()
result happens. The code I show above seems more likely to match the
way one_minus_cos_60 is computed.

I'll go try it, though I guess we won't see results till tomorrow.

regards, tom lane