Re: [COMMITTERS] pgsql: Add trigonometric functions that work in degrees. - Mailing list pgsql-hackers

On 19 April 2016 at 05:16, Noah Misch <noah@leadboat.com> wrote:
> On Mon, Apr 18, 2016 at 11:56:55PM -0400, Tom Lane wrote:
>> Noah Misch <noah@leadboat.com> writes:
>> > On Mon, Apr 18, 2016 at 10:22:59PM -0400, Tom Lane wrote:
>> >> We could alternatively set extra_float_digits to its max value and hope
>> >> that off-by-one-in-the-last-place values would get printed as something
>> >> visibly different from the exact result.  I'm not sure I want to trust
>> >> that that works reliably; but maybe it would be worth printing the
>> >> result both ways, just to provide additional info when there's a failure.
>>
>> > We'd have an independent problem if extra_float_digits=3 prints the same
>> > digits for distinguishable float values, so I wouldn't mind relying on it not
>> > to do that.  But can we expect the extra_float_digits=3 representation of
>> > those particular values to be the same for every implementation?
>>
>> Hm?  The expected answer is exact (30, 45, or whatever) in each case.
>> If we get some residual low-order digits then it's a failure, so we don't
>> need to worry about whether it's the same failure everywhere.
>
> Does something forbid snprintf implementations from printing '45'::float8 as
> 45.0000000000000001 under extra_float_digits=3?

I'm not sure it's really worth having the test output something like
45.0000000000000001 since that extra detail doesn't really seem
particularly useful beyond the fact that the result wasn't exactly 45.
Also you'd have to be careful how you modified the test, since it's
possible that 45.0000000000000001 might be printed as '45' even under
extra_float_digits=3 and so there'd be a risk of the test passing when
it ought to fail, unless you also printed something else out to
indicate exactness.

Thinking about the actual failure in this case (asind(0.5) not
returning exactly 30) a couple of theories spring to mind. One is that
the compiler is being more aggressive over function inlining, so
init_degree_constants() is being inlined, and then asin_0_5 is being
evaluated at compile time rather than runtime, giving a slightly
different result, as happened in the original version of this patch.
Another theory is that the compiler is performing an unsafe ordering
rearrangement and transforming (asin(x) / asin_0_5) * 30.0 into
asin(x) * (30.0 / asin_0_5). This might happen for example under
something like -funsafe-math-optimizations.

I did a search for other people encountering similar problems and I
came across the following quite interesting example in the Gnu
Scientific Library's implementation of log1p() --
https://github.com/ampl/gsl/blob/master/sys/log1p.c. The reason the
code is now written in that way is in response to this bug:
https://lists.gnu.org/archive/html/bug-gsl/2007-07/msg00008.html with
an overly aggressive compiler.

So using that technique, we might try using a volatile local variable
to enforce the desired evaluation order, e.g.:
   volatile double tmp;
   tmp = asin(x) / asin_0_5;   return tmp * 30.0;

A similar trick could be used to protect against
init_degree_constants() being inlined, by writing it as
   volatile double one_half = 0.5;
   asin_0_5 = asin(one_half);

since then the compiler wouldn't be able to assume that one_half was
still equal to 0.5, and wouldn't be able to optimise away the runtime
evaluation of asin() in favour of a compile-time constant.

These are both somewhat unconventional uses of volatile, but I think
they stand a better chance of being more portable, and also more
future-proof against compilers that might in the future make more
advanced code inlining and rearrangement choices.

Of course this is all pure speculation at this stage, but it seems
like it's worth a try.

Regards,
Dean