Thread: Re: [COMMITTERS] pgsql: Add trigonometric functions that work in degrees.

On 01/22/2016 03:46 PM, Tom Lane wrote:
> Add trigonometric functions that work in degrees.
>
> The implementations go to some lengths to deliver exact results for values
> where an exact result can be expected, such as sind(30) = 0.5 exactly.

I have a host here that is having regression test failures from this commit:

--- src/test/regress/expected/float8.out
+++ src/test/regress/results/float8.out
@@ -490,9 +490,9 @@   x   | asind | acosd | atand ------+-------+-------+-------    -1 |   -90 |   180 |   -45
- -0.5 |   -30 |   120 |
+ -0.5 |       |   120 |     0 |     0 |    90 |     0
-  0.5 |    30 |    60 |
+  0.5 |       |    60 |     1 |    90 |     0 |    45 (5 rows)

Any ideas?

On Tue, Apr 5, 2016 at 10:16 AM, Peter Eisentraut <peter_e@gmx.net> wrote:
> On 01/22/2016 03:46 PM, Tom Lane wrote:
>>
>> Add trigonometric functions that work in degrees.
>>
>> The implementations go to some lengths to deliver exact results for values
>> where an exact result can be expected, such as sind(30) = 0.5 exactly.
>
> I have a host here that is having regression test failures from this commit:
>
> --- src/test/regress/expected/float8.out
> +++ src/test/regress/results/float8.out
> @@ -490,9 +490,9 @@
>    x   | asind | acosd | atand
>  ------+-------+-------+-------
>     -1 |   -90 |   180 |   -45
> - -0.5 |   -30 |   120 |
> + -0.5 |       |   120 |
>      0 |     0 |    90 |     0
> -  0.5 |    30 |    60 |
> +  0.5 |       |    60 |
>      1 |    90 |     0 |    45
>  (5 rows)
>
> Any ideas?

Likely an oversight not tracked by the buildfarm. What are you using here?
-- 
Michael

Michael Paquier <michael.paquier@gmail.com> writes:
> On Tue, Apr 5, 2016 at 10:16 AM, Peter Eisentraut <peter_e@gmx.net> wrote:
>> On 01/22/2016 03:46 PM, Tom Lane wrote:
>>> Add trigonometric functions that work in degrees.

>> I have a host here that is having regression test failures from this commit:

> Likely an oversight not tracked by the buildfarm. What are you using here?

Indeed.  Also, I trust you're checking 00347575e or later, and not
e1bd684a3 itself?
        regards, tom lane

On 04/04/2016 09:20 PM, Michael Paquier wrote:
> Likely an oversight not tracked by the buildfarm. What are you using here?

It should be a rather unspectactular Debian system.  I have tried a 
bunch of different compilers and optimization options, without success.  I'll keep looking.

Peter Eisentraut <peter_e@gmx.net> writes:
> On 04/04/2016 09:20 PM, Michael Paquier wrote:
>> Likely an oversight not tracked by the buildfarm. What are you using here?

> It should be a rather unspectactular Debian system.  I have tried a 
> bunch of different compilers and optimization options, without success. 
>   I'll keep looking.

Perhaps it'd be worthwhile to check the behavior before and after each
of the portability patches that followed on to e1bd684a3, and see if any
of them change the behavior for you (and if so just how).  That might
give us a clue what's going on.
        regards, tom lane

Peter Eisentraut <peter_e@gmx.net> writes:
> On 01/22/2016 03:46 PM, Tom Lane wrote:
>> Add trigonometric functions that work in degrees.

> I have a host here that is having regression test failures from this commit:

> --- src/test/regress/expected/float8.out
> +++ src/test/regress/results/float8.out
> @@ -490,9 +490,9 @@
>     x   | asind | acosd | atand
>   ------+-------+-------+-------
>      -1 |   -90 |   180 |   -45
> - -0.5 |   -30 |   120 |
> + -0.5 |       |   120 |
>       0 |     0 |    90 |     0
> -  0.5 |    30 |    60 |
> +  0.5 |       |    60 |
>       1 |    90 |     0 |    45
>   (5 rows)

BTW ... looking closer at that, it appears to show asind(-0.5) and
asind(0.5) returning NULL.  Which makes no sense at all, because
there is no provision in dasind() for returning a null, regardless
of the input value.

So I'm pretty baffled.  Maybe you could step through this and figure
out where it's going off the rails?
        regards, tom lane

On Fri, Apr 8, 2016 at 5:48 PM, Tom Lane <tgl@sss.pgh.pa.us> wrote:
> Peter Eisentraut <peter_e@gmx.net> writes:
>> On 01/22/2016 03:46 PM, Tom Lane wrote:
>>> Add trigonometric functions that work in degrees.
>
>> I have a host here that is having regression test failures from this commit:
>
>> --- src/test/regress/expected/float8.out
>> +++ src/test/regress/results/float8.out
>> @@ -490,9 +490,9 @@
>>     x   | asind | acosd | atand
>>   ------+-------+-------+-------
>>      -1 |   -90 |   180 |   -45
>> - -0.5 |   -30 |   120 |
>> + -0.5 |       |   120 |
>>       0 |     0 |    90 |     0
>> -  0.5 |    30 |    60 |
>> +  0.5 |       |    60 |
>>       1 |    90 |     0 |    45
>>   (5 rows)
>
> BTW ... looking closer at that, it appears to show asind(-0.5) and
> asind(0.5) returning NULL.  Which makes no sense at all, because
> there is no provision in dasind() for returning a null, regardless
> of the input value.
>
> So I'm pretty baffled.  Maybe you could step through this and figure
> out where it's going off the rails?

Peter, are you going to look into this further?  This is on the open
items list, but there seems to be nothing that can be done about it by
anyone other than, maybe, you.

If you're not going to look into it, I think we should delete the open
item.  There's no point in tracking issues that aren't actionable.

-- 
Robert Haas
EnterpriseDB: http://www.enterprisedb.com
The Enterprise PostgreSQL Company

On 04/15/2016 02:04 PM, Robert Haas wrote:
> Peter, are you going to look into this further?  This is on the open
> items list, but there seems to be nothing that can be done about it by
> anyone other than, maybe, you.

I don't know if it's worth tracking this as an open item and thus
kind-of release blocker if no one else has the problem.  But I
definitely still have it.  My initial suspicion was that this had
something to do with a partial upgrade to gcc 6 (not yet released), in
other words a messed up system.  But I was able to reproduce it in a
freshly installed chroot.  It only happens with various versions of gcc,
but not with clang.  So I'm going to have to keep digging.

On Mon, Apr 18, 2016 at 12:31 PM, Peter Eisentraut <peter_e@gmx.net> wrote:
> I don't know if it's worth tracking this as an open item and thus
> kind-of release blocker if no one else has the problem.  But I
> definitely still have it.  My initial suspicion was that this had
> something to do with a partial upgrade to gcc 6 (not yet released), in
> other words a messed up system.  But I was able to reproduce it in a
> freshly installed chroot.  It only happens with various versions of gcc,
> but not with clang.  So I'm going to have to keep digging.

gcc is moving slowly but surely to have 6.0 in a released state btw:
https://gcc.gnu.org/ml/gcc/2016-04/msg00103.html
-- 
Michael

Michael Paquier <michael.paquier@gmail.com> writes:
> On Mon, Apr 18, 2016 at 12:31 PM, Peter Eisentraut <peter_e@gmx.net> wrote:
>> I don't know if it's worth tracking this as an open item and thus
>> kind-of release blocker if no one else has the problem.  But I
>> definitely still have it.  My initial suspicion was that this had
>> something to do with a partial upgrade to gcc 6 (not yet released), in
>> other words a messed up system.  But I was able to reproduce it in a
>> freshly installed chroot.  It only happens with various versions of gcc,
>> but not with clang.  So I'm going to have to keep digging.

> gcc is moving slowly but surely to have 6.0 in a released state btw:
> https://gcc.gnu.org/ml/gcc/2016-04/msg00103.html

Given that it's apparently showing the results of asind as NULL, the
theory that comes to mind is some sort of optimization issue affecting
the output tuple's null-flags.  I have no idea why only this test would
be affected, though.  Anyway, a good way to test that theory would be
to see if the -O level affects it.
        regards, tom lane

On Mon, Apr 18, 2016 at 09:17:46AM -0400, Tom Lane wrote:
> Michael Paquier <michael.paquier@gmail.com> writes:
> > On Mon, Apr 18, 2016 at 12:31 PM, Peter Eisentraut <peter_e@gmx.net> wrote:
> >> I don't know if it's worth tracking this as an open item and thus
> >> kind-of release blocker if no one else has the problem.  But I
> >> definitely still have it.  My initial suspicion was that this had
> >> something to do with a partial upgrade to gcc 6 (not yet released), in
> >> other words a messed up system.  But I was able to reproduce it in a
> >> freshly installed chroot.  It only happens with various versions of gcc,
> >> but not with clang.  So I'm going to have to keep digging.
> 
> > gcc is moving slowly but surely to have 6.0 in a released state btw:
> > https://gcc.gnu.org/ml/gcc/2016-04/msg00103.html
> 
> Given that it's apparently showing the results of asind as NULL, the
> theory that comes to mind is some sort of optimization issue affecting
> the output tuple's null-flags.  I have no idea why only this test would
> be affected, though.  Anyway, a good way to test that theory would be
> to see if the -O level affects it.

I doubt asind is returning NULL.  Here's the query, which uses a CASE to
report NULL if asind returns any value not on a whitelist:

SELECT x,      CASE WHEN asind(x) IN (-90,-30,0,30,90) THEN asind(x) END AS asind,      CASE WHEN acosd(x) IN
(0,60,90,120,180)THEN acosd(x) END AS acosd,      CASE WHEN atand(x) IN (-45,0,45) THEN atand(x) END AS atand
 
FROM (VALUES (-1), (-0.5), (0), (0.5), (1)) AS t(x);

I can see the benefit for atand(-0.5) and for atand(0.5), since those are
inexact.  Does the CASE gain us anything for asind or acosd?

Results under -O0 would be a helpful data point, nonetheless.

Noah Misch <noah@leadboat.com> writes:
> On Mon, Apr 18, 2016 at 09:17:46AM -0400, Tom Lane wrote:
>> Given that it's apparently showing the results of asind as NULL ...

> I doubt asind is returning NULL.  Here's the query, which uses a CASE to
> report NULL if asind returns any value not on a whitelist:

> SELECT x,
>        CASE WHEN asind(x) IN (-90,-30,0,30,90) THEN asind(x) END AS asind,
>        CASE WHEN acosd(x) IN (0,60,90,120,180) THEN acosd(x) END AS acosd,
>        CASE WHEN atand(x) IN (-45,0,45) THEN atand(x) END AS atand
> FROM (VALUES (-1), (-0.5), (0), (0.5), (1)) AS t(x);

Oh, duh --- should have checked the query.  Yes, the most probable theory
must be that it's returning something that's slightly off from the exact
value.

> I can see the benefit for atand(-0.5) and for atand(0.5), since those are
> inexact.  Does the CASE gain us anything for asind or acosd?

None of these are expected to be inexact.  The point of the CASE is to
make it obvious if what's returned isn't *exactly* what we expect.

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.
        regards, tom lane

On Mon, Apr 18, 2016 at 10:22:59PM -0400, Tom Lane wrote:
> Noah Misch <noah@leadboat.com> writes:
> > On Mon, Apr 18, 2016 at 09:17:46AM -0400, Tom Lane wrote:
> >> Given that it's apparently showing the results of asind as NULL ...
> 
> > I doubt asind is returning NULL.  Here's the query, which uses a CASE to
> > report NULL if asind returns any value not on a whitelist:
> 
> > SELECT x,
> >        CASE WHEN asind(x) IN (-90,-30,0,30,90) THEN asind(x) END AS asind,
> >        CASE WHEN acosd(x) IN (0,60,90,120,180) THEN acosd(x) END AS acosd,
> >        CASE WHEN atand(x) IN (-45,0,45) THEN atand(x) END AS atand
> > FROM (VALUES (-1), (-0.5), (0), (0.5), (1)) AS t(x);
> 
> Oh, duh --- should have checked the query.  Yes, the most probable theory
> must be that it's returning something that's slightly off from the exact
> value.
> 
> > I can see the benefit for atand(-0.5) and for atand(0.5), since those are
> > inexact.  Does the CASE gain us anything for asind or acosd?
> 
> None of these are expected to be inexact.  The point of the CASE is to
> make it obvious if what's returned isn't *exactly* what we expect.

Ah, got it.

> 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?

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.
        regards, tom lane

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?

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

Dean Rasheed <dean.a.rasheed@gmail.com> writes:
> 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:
>>> 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.

Yeah, what I was thinking of printing is something like
asind(x),asind(x) IN (-90,-30,0,30,90) AS asind_exact,...

with extra_float_digits=3.  The point of this is not necessarily to give
any extra information, though it might, but for failures to be more easily
interpretable.  If I'd forgotten how the test worked just a few months
after committing it, how likely is it that some random user faced with a
similar failure would understand what they were seeing?

Also, though I agree that it might not help much to know whether the
output is 45.0000000000000001 or 44.9999999999999999, our thoughts would
be trending in quite a different direction if it turns out that the
output is radically wrong, or even a NaN.  The existing test cannot
exclude that possibility.
        regards, tom lane

On 19 April 2016 at 14:38, Tom Lane <tgl@sss.pgh.pa.us> wrote:
> Yeah, what I was thinking of printing is something like
>
>         asind(x),
>         asind(x) IN (-90,-30,0,30,90) AS asind_exact,
>         ...
>
> with extra_float_digits=3.  The point of this is not necessarily to give
> any extra information, though it might, but for failures to be more easily
> interpretable.  If I'd forgotten how the test worked just a few months
> after committing it, how likely is it that some random user faced with a
> similar failure would understand what they were seeing?
>
> Also, though I agree that it might not help much to know whether the
> output is 45.0000000000000001 or 44.9999999999999999, our thoughts would
> be trending in quite a different direction if it turns out that the
> output is radically wrong, or even a NaN.  The existing test cannot
> exclude that possibility.
>

OK, that sounds like it would be a useful improvement to the tests.

Regards,
Dean

Dean Rasheed <dean.a.rasheed@gmail.com> writes:
> On 19 April 2016 at 14:38, Tom Lane <tgl@sss.pgh.pa.us> wrote:
>> Yeah, what I was thinking of printing is something like
>> 
>> asind(x),
>> asind(x) IN (-90,-30,0,30,90) AS asind_exact,
>> ...
>> 
>> with extra_float_digits=3.

> OK, that sounds like it would be a useful improvement to the tests.

Pushed.  Peter, what results do you get from these tests on your
problematic machine?
        regards, tom lane

On 04/19/2016 04:48 PM, Tom Lane wrote:
> Pushed.  Peter, what results do you get from these tests on your
> problematic machine?
       acosd(x),       acosd(x) IN (0,60,90,120,180) AS acosd_exactFROM (VALUES (-1), (-0.5), (0), (0.5), (1)) AS
t(x);
-  x   | asind | asind_exact | acosd | acosd_exact
-------+-------+-------------+-------+-------------
-   -1 |   -90 | t           |   180 | t
- -0.5 |   -30 | t           |   120 | t
-    0 |     0 | t           |    90 | t
-  0.5 |    30 | t           |    60 | t
-    1 |    90 | t           |     0 | t
+  x   |        asind         | asind_exact | acosd | acosd_exact
+------+----------------------+-------------+-------+-------------
+   -1 |                  -90 | t           |   180 | t
+ -0.5 | -29.9999999999999964 | f           |   120 | t
+    0 |                    0 | t           |    90 | t
+  0.5 |  29.9999999999999964 | f           |    60 | t
+    1 |                   90 | t           |     0 | t(5 rows)

This is the same under the default -O2 and under -O0.

Peter Eisentraut <peter_e@gmx.net> writes:
> On 04/19/2016 04:48 PM, Tom Lane wrote:
>> Pushed.  Peter, what results do you get from these tests on your
>> problematic machine?

>         acosd(x),
>         acosd(x) IN (0,60,90,120,180) AS acosd_exact
>  FROM (VALUES (-1), (-0.5), (0), (0.5), (1)) AS t(x);
> -  x   | asind | asind_exact | acosd | acosd_exact
> -------+-------+-------------+-------+-------------
> -   -1 |   -90 | t           |   180 | t
> - -0.5 |   -30 | t           |   120 | t
> -    0 |     0 | t           |    90 | t
> -  0.5 |    30 | t           |    60 | t
> -    1 |    90 | t           |     0 | t
> +  x   |        asind         | asind_exact | acosd | acosd_exact
> +------+----------------------+-------------+-------+-------------
> +   -1 |                  -90 | t           |   180 | t
> + -0.5 | -29.9999999999999964 | f           |   120 | t
> +    0 |                    0 | t           |    90 | t
> +  0.5 |  29.9999999999999964 | f           |    60 | t
> +    1 |                   90 | t           |     0 | t
>  (5 rows)

> This is the same under the default -O2 and under -O0.

Hm.  This seems to prove that we're not getting exactly 1.0 from
(asin(x) / asin_0_5) with x = 0.5, but I'm having a hard time guessing
why that might be so when all the other cases work.

Could you send along the assembler code generated by the compiler (-S
output) for float.c?  Maybe that would shed some light.  Probably the
-O0 version would be easier to read.
        regards, tom lane

Peter Eisentraut <peter_e@gmx.net> writes:
> On 04/21/2016 08:18 PM, Tom Lane wrote:
>> Could you send along the assembler code generated by the compiler (-S
>> output) for float.c?  Maybe that would shed some light.  Probably the
>> -O0 version would be easier to read.

> Attached is a smaller test program that prints 29.9999999999999964 (same
> as failing test result) as well as its assembler code.

This doesn't prove much.  It's clear from the assembly code that the
compiler is evaluating "asin(0.5)" at compile time (and, presumably,
getting a different answer than would be produced at runtime).  But
that's entirely unsurprising given this source text: you've got a
naked application of asin() to a constant.  What we need to know,
and what this doesn't quite prove, is whether the compiler is managing
to see through the tricks that float.c uses to prevent compile-time
calculation of that value.

The fact that this produces the same number as we see in the regression
test does point to the idea that that's what's happening.  But it's not
conclusive.  Part of the reason I'm not sold is that even if the compiler
can see through those tricks when optimizing, it hardly seems like it
should do so at -O0.

If that is the answer, then the next question is how we can put more
roadblocks in the way of compile-time evaluation of asin(0.5).  Dean
suggested that creative use of "volatile" might do it, and I agree
that that sounds like a promising thing to pursue.
        regards, tom lane

I wrote:
> If that is the answer, then the next question is how we can put more
> roadblocks in the way of compile-time evaluation of asin(0.5).  Dean
> suggested that creative use of "volatile" might do it, and I agree
> that that sounds like a promising thing to pursue.

It occurred to me that we don't actually need "volatile".  What we need
is a variable that the compiler is not allowed to assume is a compile-time
constant, and a plain global variable is sufficient for that.  In the
attached patch, we no longer need an assumption that init_degree_constants
doesn't get inlined; we only need to assume that the compiler can't prove
the variables degree_c_thirty etc to be immutable.  Which it cannot, even
if it does global optimization across the whole PG executable, because it
has to consider that loadable extensions might change them.

I'm going to go ahead and push this, because it seems clearly more robust
than what we have.  But I'd appreciate a report on whether it fixes your
issue.

            regards, tom lane

diff --git a/src/backend/utils/adt/float.c b/src/backend/utils/adt/float.c
index c7c0b58..a95ebe5 100644
*** a/src/backend/utils/adt/float.c
--- b/src/backend/utils/adt/float.c
*************** static float8 atan_1_0 = 0;
*** 77,91 ****
  static float8 tan_45 = 0;
  static float8 cot_45 = 0;

  /* Local function prototypes */
  static int    float4_cmp_internal(float4 a, float4 b);
  static int    float8_cmp_internal(float8 a, float8 b);
  static double sind_q1(double x);
  static double cosd_q1(double x);
!
! /* This is INTENTIONALLY NOT STATIC.  Don't "fix" it. */
! void init_degree_constants(float8 thirty, float8 forty_five, float8 sixty,
!                       float8 one_half, float8 one);

  #ifndef HAVE_CBRT
  /*
--- 77,100 ----
  static float8 tan_45 = 0;
  static float8 cot_45 = 0;

+ /*
+  * These are intentionally not static; don't "fix" them.  They will never
+  * be referenced by other files, much less changed; but we don't want the
+  * compiler to know that, else it might try to precompute expressions
+  * involving them.  See comments for init_degree_constants().
+  */
+ float8        degree_c_thirty = 30.0;
+ float8        degree_c_forty_five = 45.0;
+ float8        degree_c_sixty = 60.0;
+ float8        degree_c_one_half = 0.5;
+ float8        degree_c_one = 1.0;
+
  /* Local function prototypes */
  static int    float4_cmp_internal(float4 a, float4 b);
  static int    float8_cmp_internal(float8 a, float8 b);
  static double sind_q1(double x);
  static double cosd_q1(double x);
! static void init_degree_constants(void);

  #ifndef HAVE_CBRT
  /*
*************** dtan(PG_FUNCTION_ARGS)
*** 1814,1848 ****
   * compilers out there that will precompute expressions such as sin(constant)
   * using a sin() function different from what will be used at runtime.  If we
   * want exact results, we must ensure that none of the scaling constants used
!  * in the degree-based trig functions are computed that way.
!  *
!  * The whole approach fails if init_degree_constants() gets inlined into the
!  * call sites, since then constant-folding can happen anyway.  Currently it
!  * seems sufficient to declare it non-static to prevent that.  We have no
!  * expectation that other files will call this, but don't tell gcc that.
   *
   * Other hazards we are trying to forestall with this kluge include the
   * possibility that compilers will rearrange the expressions, or compute
   * some intermediate results in registers wider than a standard double.
   */
! void
! init_degree_constants(float8 thirty, float8 forty_five, float8 sixty,
!                       float8 one_half, float8 one)
  {
!     sin_30 = sin(thirty * RADIANS_PER_DEGREE);
!     one_minus_cos_60 = 1.0 - cos(sixty * RADIANS_PER_DEGREE);
!     asin_0_5 = asin(one_half);
!     acos_0_5 = acos(one_half);
!     atan_1_0 = atan(one);
!     tan_45 = sind_q1(forty_five) / cosd_q1(forty_five);
!     cot_45 = cosd_q1(forty_five) / sind_q1(forty_five);
      degree_consts_set = true;
  }

  #define INIT_DEGREE_CONSTANTS() \
  do { \
      if (!degree_consts_set) \
!         init_degree_constants(30.0, 45.0, 60.0, 0.5, 1.0); \
  } while(0)


--- 1823,1853 ----
   * compilers out there that will precompute expressions such as sin(constant)
   * using a sin() function different from what will be used at runtime.  If we
   * want exact results, we must ensure that none of the scaling constants used
!  * in the degree-based trig functions are computed that way.  To do so, we
!  * compute them from the variables degree_c_thirty etc, which are also really
!  * constants, but the compiler cannot assume that.
   *
   * Other hazards we are trying to forestall with this kluge include the
   * possibility that compilers will rearrange the expressions, or compute
   * some intermediate results in registers wider than a standard double.
   */
! static void
! init_degree_constants(void)
  {
!     sin_30 = sin(degree_c_thirty * RADIANS_PER_DEGREE);
!     one_minus_cos_60 = 1.0 - cos(degree_c_sixty * RADIANS_PER_DEGREE);
!     asin_0_5 = asin(degree_c_one_half);
!     acos_0_5 = acos(degree_c_one_half);
!     atan_1_0 = atan(degree_c_one);
!     tan_45 = sind_q1(degree_c_forty_five) / cosd_q1(degree_c_forty_five);
!     cot_45 = cosd_q1(degree_c_forty_five) / sind_q1(degree_c_forty_five);
      degree_consts_set = true;
  }

  #define INIT_DEGREE_CONSTANTS() \
  do { \
      if (!degree_consts_set) \
!         init_degree_constants(); \
  } while(0)

On 04/25/2016 03:09 PM, Tom Lane wrote:
> I'm going to go ahead and push this, because it seems clearly more robust
> than what we have.  But I'd appreciate a report on whether it fixes your
> issue.

6b1a213bbd6599228b2b67f7552ff7cc378797bf did not fix it.

Attached is the assembler output (-O0) from float.c as of that commit.

Peter Eisentraut <peter_e@gmx.net> writes:
> On 04/25/2016 03:09 PM, Tom Lane wrote:
>> I'm going to go ahead and push this, because it seems clearly more robust
>> than what we have.  But I'd appreciate a report on whether it fixes your
>> issue.

> 6b1a213bbd6599228b2b67f7552ff7cc378797bf did not fix it.

OK ... but it's still a good change, because it removes the assumption
that the compiler won't inline init_degree_constants().

> Attached is the assembler output (-O0) from float.c as of that commit.

Ah.  Here's the problem: line 1873 is
   return (asin(x) / asin_0_5) * 30.0;

and that compiles into
subl    $8, %esppushl    -12(%ebp)    ... push xpushl    -16(%ebp)call    asin        ... call asin()addl    $16,
%espfldl   asin_0_5    ... divide by asin_0_5fdivrp    %st, %st(1)fldt    .LC46        ... multiply by 30.0fmulp
%st,%st(1)fstpl    -24(%ebp)    ... round to 64 bitsfldl    -24(%ebp)
 

Evidently, asin() is returning an 80-bit result, and that's not
getting rounded to double width before we divide by asin_0_5.

The least ugly change that might fix this is to explicitly cast the
result of asin to double:
   return ((float8) asin(x) / asin_0_5) * 30.0;

However, I'm not certain that that would do anything; the compiler
might feel that the result of asin() is already double.  The next
messier answer is to explicitly store the function result in a local
variable:
   {       float8 asin_x = asin(x);       return (asin_x / asin_0_5) * 30.0;   }

A sufficiently cavalier compiler might choose to optimize that away,
too.  A look at gcc's manual suggests that we might need to use
the -ffloat-store option to guarantee it will work; which is ugly
and I'd prefer not to turn that on globally anyway.  If it comes to
that, probably the better answer is to turn asin_x into a global
variable, similarly to what we just did with the constants.

Can you try the above variants of line 1873 and see if either of them
fixes the issue for you?
        regards, tom lane

On 04/21/2016 08:18 PM, Tom Lane wrote:
> Hm.  This seems to prove that we're not getting exactly 1.0 from
> (asin(x) / asin_0_5) with x = 0.5, but I'm having a hard time guessing
> why that might be so when all the other cases work.
>
> Could you send along the assembler code generated by the compiler (-S
> output) for float.c?  Maybe that would shed some light.  Probably the
> -O0 version would be easier to read.

Attached is a smaller test program that prints 29.9999999999999964 (same
as failing test result) as well as its assembler code.

I wrote:
> Evidently, asin() is returning an 80-bit result, and that's not
> getting rounded to double width before we divide by asin_0_5.

I found that I could reproduce this type of assembly-code behavior
on dromedary (gcc version 4.2.1) by using -mfpmath=387.  That box
doesn't show the visible regression-test failure, but that must be
down to its version of asin() not producing the same low-order bits
that yours does.  It's clear from the assembly output that the code
*would* misbehave given an appropriate asin() library function.

With that version of gcc, just casting the function output to double
changes nothing.  The local-variable solution likewise.  I do get
a useful fix if I declare the asin_x local variable volatile:
       .loc 1 1873 0       fstpl   (%esp)       call    _asin
+       fstpl   -16(%ebp)
+ LVL107:
+       fldl    -16(%ebp)       fdivl   _asin_0_5-"L00000000019$pb"(%ebx)       fmuls   LC21-"L00000000019$pb"(%ebx)

Interestingly, declaring asin_x as global is not enough to fix it,
because it stores into the global but doesn't fetch back:
       .loc 1 1873 0       fstpl   (%esp)       call    _asin
+       movl    L_asin_x$non_lazy_ptr-"L00000000019$pb"(%ebx), %eax
+       fstl    (%eax)       fdivl   _asin_0_5-"L00000000019$pb"(%ebx)       fmuls   LC21-"L00000000019$pb"(%ebx)

Declaring asin_x as a volatile global does fix it:
       .loc 1 1873 0       fstpl   (%esp)       call    _asin
+       movl    L_asin_x$non_lazy_ptr-"L00000000019$pb"(%ebx), %eax
+       fstpl   (%eax)
+       fldl    (%eax)       fdivl   _asin_0_5-"L00000000019$pb"(%ebx)       fmuls   LC21-"L00000000019$pb"(%ebx)

but that seems like just being uglier without much redeeming value.

In short, these tests suggest that we need a coding pattern like
this:
volatile float8 asin_x = asin(x);return (asin_x / asin_0_5) * 30.0;

We could drop the "volatile" by adopting -ffloat-store, but that
doesn't seem like a better answer, because -ffloat-store would
pessimize a lot of code elsewhere.  (It causes a whole lot of
changes in float.c, for sure.)
        regards, tom lane

On 26 April 2016 at 04:25, Tom Lane <tgl@sss.pgh.pa.us> wrote:
> In short, these tests suggest that we need a coding pattern like
> this:
>
>         volatile float8 asin_x = asin(x);
>         return (asin_x / asin_0_5) * 30.0;
>
> We could drop the "volatile" by adopting -ffloat-store, but that
> doesn't seem like a better answer, because -ffloat-store would
> pessimize a lot of code elsewhere.  (It causes a whole lot of
> changes in float.c, for sure.)

Agreed. That looks like the least hacky way of solving the problem. I
think it's more readable when the logic is kept local, and it's
preferable to avoid any compiler-specific options or global flags that
would affect other code.

6b1a213bbd6599228b2b67f7552ff7cc378797bf by itself looks like a
worthwhile improvement that ought to be more robust, even though it
didn't fix this problem.

Regards,
Dean

Dean Rasheed <dean.a.rasheed@gmail.com> writes:
> On 26 April 2016 at 04:25, Tom Lane <tgl@sss.pgh.pa.us> wrote:
>> In short, these tests suggest that we need a coding pattern like
>> this:
>> volatile float8 asin_x = asin(x);
>> return (asin_x / asin_0_5) * 30.0;

> Agreed. That looks like the least hacky way of solving the problem. I
> think it's more readable when the logic is kept local, and it's
> preferable to avoid any compiler-specific options or global flags that
> would affect other code.

OK, I've pushed a change along these lines.  Peter, would you see whether
HEAD fixes it for you?

The next time somebody proposes that we can get exact results out of
floating-point arithmetic, I'm going to run away screaming.
        regards, tom lane

On Tue, Apr 26, 2016 at 11:26 AM, Tom Lane <tgl@sss.pgh.pa.us> wrote:
> The next time somebody proposes that we can get exact results out of
> floating-point arithmetic, I'm going to run away screaming.

And you wonder why I avoid floating point variables like the plague...!

-- 
Robert Haas
EnterpriseDB: http://www.enterprisedb.com
The Enterprise PostgreSQL Company

Robert Haas <robertmhaas@gmail.com> writes:
> On Tue, Apr 26, 2016 at 11:26 AM, Tom Lane <tgl@sss.pgh.pa.us> wrote:
>> The next time somebody proposes that we can get exact results out of
>> floating-point arithmetic, I'm going to run away screaming.

> And you wonder why I avoid floating point variables like the plague...!

Float variables are fine for their intended use-case.  The problem is
with trying to guarantee identical results across all platforms.
        regards, tom lane

On Tue, Apr 26, 2016 at 11:36 AM, Tom Lane <tgl@sss.pgh.pa.us> wrote:
> Robert Haas <robertmhaas@gmail.com> writes:
>> On Tue, Apr 26, 2016 at 11:26 AM, Tom Lane <tgl@sss.pgh.pa.us> wrote:
>>> The next time somebody proposes that we can get exact results out of
>>> floating-point arithmetic, I'm going to run away screaming.
>
>> And you wonder why I avoid floating point variables like the plague...!
>
> Float variables are fine for their intended use-case.  The problem is
> with trying to guarantee identical results across all platforms.

Yeah, yeah... :-)

-- 
Robert Haas
EnterpriseDB: http://www.enterprisedb.com
The Enterprise PostgreSQL Company

On 26 April 2016 at 16:26, Tom Lane <tgl@sss.pgh.pa.us> wrote:
> The next time somebody proposes that we can get exact results out of
> floating-point arithmetic, I'm going to run away screaming.
>

Yeah, I think I will have the same reaction.

Thanks for all your hard work getting this to work.

Regards,
Dean

On 04/26/2016 11:26 AM, Tom Lane wrote:
> OK, I've pushed a change along these lines.  Peter, would you see whether
> HEAD fixes it for you?

Yeah, it passes now.  Nobody touch anything! ;-)