Thread: Numeric x^y for negative x

Numeric x^y is supported for x < 0 if y is an integer, but this
currently fails if y is outside the range of an int32:

SELECT (-1.0) ^ 2147483647;
      ?column?
---------------------
 -1.0000000000000000
(1 row)

SELECT (-1.0) ^ 2147483648;
ERROR:  cannot take logarithm of a negative number

because only the power_var_int() code path in power_var() handles
negative bases correctly. Attached is a patch to fix that.

Regards,
Dean

On Tue, 29 Jun 2021 at 12:08, Dean Rasheed <dean.a.rasheed@gmail.com> wrote:
>
> Numeric x^y is supported for x < 0 if y is an integer, but this
> currently fails if y is outside the range of an int32
>

I've been doing some more testing of this, and I spotted another
problem with numeric_power().

This is what happens when raising 0.9999999999 to increasingly large
powers, which should decrease to zero:

   exp             0.9999999999^exp
10000000000     0.3678794411530483
100000000000    0.[4 zeros]4539992973978489
1000000000000   0.[43 zeros]3720075957420456
10000000000000  0.[434 zeros]5075958643751518
20000000000000  0.[868 zeros]2576535615307575
21000000000000  0.[912 zeros]9584908195943232
22000000000000  0.[955 zeros]3565658653381070
23000000000000  0.[998 zeros]13
23100000000000  0.[1000 zeros]
23200000000000  0.[1000 zeros]
23300000000000  1.[1000 zeros]  *** WRONG ***
30000000000000  1.[1000 zeros]  *** WRONG ***
40000000000000  1.[1000 zeros]  *** WRONG ***
50000000000000  1.[1000 zeros]  *** WRONG ***
60000000000000  1.[1000 zeros]  *** WRONG ***
70000000000000  ERROR:  value overflows numeric format

The cases where it returns 1 are a trivial logic bug in the
local_rscale calculation in power_var() -- when it computes
local_rscale from rscale and val, it needs to do so before clipping
rscale to NUMERIC_MAX_DISPLAY_SCALE, otherwise it ends up setting
local_rscale = 0, and loses all precision.

I also don't think it should be throwing an overflow error here. Some
code paths through numeric_power() catch cases that would underflow,
and return zero instead, but not all cases are caught. There's a
similar overflow error with numeric_exp() for large negative inputs
(-5999 returns 0, but -6000 overflows).

It's arguable though that numeric power() and exp() (and mul() for
that matter) should never return 0 for finite non-zero inputs, but
instead should throw underflow errors, which would make them
compatible with their floating-point counterparts. I don't think
that's useful though, and it's more likely to break people's code for
no real benefit. No other numeric code throws underflow errors.

So I think we should just attempt to avoid all such overflow errors,
that are actually underflows, and return zero instead.

Regards,
Dean

On Thu, 1 Jul 2021 at 14:17, Dean Rasheed <dean.a.rasheed@gmail.com> wrote:
>
> On Tue, 29 Jun 2021 at 12:08, Dean Rasheed <dean.a.rasheed@gmail.com> wrote:
> >
> > Numeric x^y is supported for x < 0 if y is an integer, but this
> > currently fails if y is outside the range of an int32
>
> I've been doing some more testing of this, and I spotted another
> problem with numeric_power().
>
> [loss of precision and overflow errors]
>
> I think we should attempt to avoid all such overflow errors,
> that are actually underflows, and return zero instead.
>

Finally getting back to this ... attached is an updated patch that now
includes a fix for the loss-of-precision bug and the overflow errors.
I don't think it's really worth trying to split these up, since
they're all somewhat interrelated.

Regards,
Dean

On Wed, 7 Jul 2021 at 18:57, Zhihong Yu <zyu@yugabyte.com> wrote:
>
> +                   (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
> +                    errmsg("value overflows numeric format")));
>
> Here is an example of existing error message which I think is more readable than 'overflows numeric format':
>
>                  errmsg("bigint out of range")));
>
> Maybe rephrase as: value is out of range
>

Hmm, I don't know. That's the error that has been thrown by lots of
numeric functions for a long time now, and it seems fine to me.

Regards,
Dean

On Wed, 7 Jul 2021 18:36:56 +0100
Dean Rasheed <dean.a.rasheed@gmail.com> wrote:

> On Thu, 1 Jul 2021 at 14:17, Dean Rasheed <dean.a.rasheed@gmail.com> wrote:
> >
> > On Tue, 29 Jun 2021 at 12:08, Dean Rasheed <dean.a.rasheed@gmail.com> wrote:
> > >
> > > Numeric x^y is supported for x < 0 if y is an integer, but this
> > > currently fails if y is outside the range of an int32
> >
> > I've been doing some more testing of this, and I spotted another
> > problem with numeric_power().
> >
> > [loss of precision and overflow errors]
> >
> > I think we should attempt to avoid all such overflow errors,
> > that are actually underflows, and return zero instead.
> >
> 
> Finally getting back to this ... attached is an updated patch that now
> includes a fix for the loss-of-precision bug and the overflow errors.
> I don't think it's really worth trying to split these up, since
> they're all somewhat interrelated.

The patch can be applied cleanly.
(Though, I need to remove lines "new file mode 100644" else I get an error
 "error: git apply: bad git-diff - expected /dev/null on line 4".)

Compilation succeeded, and all tests passed.

This patch fixes numeric_power() to handle negative bases correctly and not
to raise an error "cannot take logarithm of a negative number", as well as a
bug that a result whose values is almost zero is incorrectly returend as 1.
The previous behaviors are obvious strange, and these fixes seem to me reasonable.

Also, improper overflow errors are corrected in numeric_power() and
numeric_exp() to return 0 when it is underflow in fact.
I think it is no problem that these functions return zero instead of underflow
errors because power_var_int() already do the same.

The patch includes additional tests for checking negative bases cases and
underflow and rounding of the almost zero results. It seems good.

Let me just make one comment.

                (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
                 errmsg("zero raised to a negative power is undefined")));

-   if (sign1 < 0 && !numeric_is_integral(num2))
-       ereport(ERROR,
-               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
-                errmsg("a negative number raised to a non-integer power yields a complex result")));
-
    /*
     * Initialize things
     */

I don't think we need to move this check from numeric_power to power_var.
I noticed the following comment in a numeric_power(). 

    /*    
     * The SQL spec requires that we emit a particular SQLSTATE error code for
     * certain error conditions.  Specifically, we don't return a
     * divide-by-zero error code for 0 ^ -1.
     */

In the original code, two checks that could raise an error of
ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION are following the comment.
I think these check codes are placed together under this comment intentionally,
so I suggest not to move one of them.


Regards,
Yugo Nagata

-- 
Yugo NAGATA <nagata@sraoss.co.jp>

On Tue, 20 Jul 2021 at 10:17, Yugo NAGATA <nagata@sraoss.co.jp> wrote:
>
> This patch fixes numeric_power() to handle negative bases correctly and not
> to raise an error "cannot take logarithm of a negative number", as well as a
> bug that a result whose values is almost zero is incorrectly returend as 1.
> The previous behaviors are obvious strange, and these fixes seem to me reasonable.
>
> Also, improper overflow errors are corrected in numeric_power() and
> numeric_exp() to return 0 when it is underflow in fact.
> I think it is no problem that these functions return zero instead of underflow
> errors because power_var_int() already do the same.
>
> The patch includes additional tests for checking negative bases cases and
> underflow and rounding of the almost zero results. It seems good.

Thanks for the review!


> Let me just make one comment.
>
>                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
>                  errmsg("zero raised to a negative power is undefined")));
>
> -   if (sign1 < 0 && !numeric_is_integral(num2))
> -       ereport(ERROR,
> -               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
> -                errmsg("a negative number raised to a non-integer power yields a complex result")));
> -
>     /*
>      * Initialize things
>      */
>
> I don't think we need to move this check from numeric_power to power_var.

Moving it to power_var() means that it only needs to be checked in the
case of a negative base, together with an exponent that cannot be
handled by power_var_int(), which saves unnecessary checking. It isn't
necessary to do this test at all if the exponent is an integer small
enough to fit in a 32-bit int. And if it's not an integer, or it's a
larger integer than that, it seems more logical to do the test in
power_var() near to the other code handling that case.


> I noticed the following comment in a numeric_power().
>
>     /*
>      * The SQL spec requires that we emit a particular SQLSTATE error code for
>      * certain error conditions.  Specifically, we don't return a
>      * divide-by-zero error code for 0 ^ -1.
>      */
>
> In the original code, two checks that could raise an error of
> ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION are following the comment.
> I think these check codes are placed together under this comment intentionally,
> so I suggest not to move one of them.

Ah, that's a good point about the SQL spec. The comment only refers to
the case of 0 ^ -1, but the SQL spec does indeed say that a negative
number to a non-integer power should return the same error code.

Here is an updated patch with additional comments about the required
error code when raising a negative number to a non-integer power, and
where it is checked.

Regards,
Dean

On Wed, 21 Jul 2021 11:10:16 +0100
Dean Rasheed <dean.a.rasheed@gmail.com> wrote:

> On Tue, 20 Jul 2021 at 10:17, Yugo NAGATA <nagata@sraoss.co.jp> wrote:
> >
> > This patch fixes numeric_power() to handle negative bases correctly and not
> > to raise an error "cannot take logarithm of a negative number", as well as a
> > bug that a result whose values is almost zero is incorrectly returend as 1.
> > The previous behaviors are obvious strange, and these fixes seem to me reasonable.
> >
> > Also, improper overflow errors are corrected in numeric_power() and
> > numeric_exp() to return 0 when it is underflow in fact.
> > I think it is no problem that these functions return zero instead of underflow
> > errors because power_var_int() already do the same.
> >
> > The patch includes additional tests for checking negative bases cases and
> > underflow and rounding of the almost zero results. It seems good.
> 
> Thanks for the review!
> 
> 
> > Let me just make one comment.
> >
> >                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
> >                  errmsg("zero raised to a negative power is undefined")));
> >
> > -   if (sign1 < 0 && !numeric_is_integral(num2))
> > -       ereport(ERROR,
> > -               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
> > -                errmsg("a negative number raised to a non-integer power yields a complex result")));
> > -
> >     /*
> >      * Initialize things
> >      */
> >
> > I don't think we need to move this check from numeric_power to power_var.
> 
> Moving it to power_var() means that it only needs to be checked in the
> case of a negative base, together with an exponent that cannot be
> handled by power_var_int(), which saves unnecessary checking. It isn't
> necessary to do this test at all if the exponent is an integer small
> enough to fit in a 32-bit int. And if it's not an integer, or it's a
> larger integer than that, it seems more logical to do the test in
> power_var() near to the other code handling that case.

Indeed, I agree with that this change saves unnecessary checking.

> 
> > I noticed the following comment in a numeric_power().
> >
> >     /*
> >      * The SQL spec requires that we emit a particular SQLSTATE error code for
> >      * certain error conditions.  Specifically, we don't return a
> >      * divide-by-zero error code for 0 ^ -1.
> >      */
> >
> > In the original code, two checks that could raise an error of
> > ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION are following the comment.
> > I think these check codes are placed together under this comment intentionally,
> > so I suggest not to move one of them.
> 
> Ah, that's a good point about the SQL spec. The comment only refers to
> the case of 0 ^ -1, but the SQL spec does indeed say that a negative
> number to a non-integer power should return the same error code.
> 
> Here is an updated patch with additional comments about the required
> error code when raising a negative number to a non-integer power, and
> where it is checked.

Thank you for updating the patch. I am fine with the additional comments.
I don't think there is any other problem left, so I marked it Ready-for-Committers.

Regards,
Yugo Nagata

-- 
Yugo NAGATA <nagata@sraoss.co.jp>

On Thu, 22 Jul 2021 at 06:13, Yugo NAGATA <nagata@sraoss.co.jp> wrote:
>
> Thank you for updating the patch. I am fine with the additional comments.
> I don't think there is any other problem left, so I marked it Ready-for-Committers.
>

Thanks for looking. Barring any further comments, I'll push this in a few days.

Regards,
Dean

On Thu, 22 Jul 2021 at 16:19, Dean Rasheed <dean.a.rasheed@gmail.com> wrote:
>
> On Thu, 22 Jul 2021 at 06:13, Yugo NAGATA <nagata@sraoss.co.jp> wrote:
> >
> > Thank you for updating the patch. I am fine with the additional comments.
> > I don't think there is any other problem left, so I marked it Ready-for-Committers.
>
> Thanks for looking. Barring any further comments, I'll push this in a few days.
>

So I have been testing this a lot over the last couple of days, and I
have concluded that the patch works well as far as it goes, but I did
manage to construct another case where numeric_power() loses
precision. I think, though, that it would be better to tackle that as
a separate patch.

In writing the commit message for this patch, I realised that it was
possible to tidy up the local_rscale calculation part of it a little,
to make it more obvious what was going wrong, so attached is a
slightly tweaked version. I'll hold off on committing it for a few
more days in case anyone else wants to have a look. Tom?

The other issue I found is related to the first part of power_var(),
where it does a low-precision calculation to get an estimate for the
weight of the result. It occurred to me that for certain input bases,
that calculation could be made to be quite inaccurate, and therefore
lead to choosing the wrong rscale later. This is the test case I
constructed:

  (1-1.5123456789012345678e-1000) ^ 1.15e1003

Here, the base is a sliver under 1, and so ln(base) is approximately
-1.5e-1000, and ln_dweight is -1000 (the decimal weight of ln(base)).
The problem is that the local_rscale used for the first low-precision
calculation is limited to NUMERIC_MAX_DISPLAY_SCALE (which is 1000),
so we only compute ln_base to a scale of 1000 at that stage, and the
result is rounded to exactly 2e-1000, which is off by a factor of
around 1.333333.

That makes it think the result weight will be -998, when actually it's
-755, so it then chooses a local_rscale for the full calculation
that's far too small, and the result is very inaccurate.

To fix this, I think it's necessary to remove the line that limits the
initial local_rscale. I tried that in a debugger, and managed to get a
result that agreed with the result from "bc -l" with a scale of 2000.

The reason I think it will be OK to remove that line is that it only
ever comes into play when ln_dweight is a large negative number (and
the smallest it can be is -16383). But that can only happen in
instances like this, where the base is very very close to 1. In such
cases, the ln(base) calculation is very fast, because it basically
only has to do a couple of Taylor series terms, and it's done. This
will still only be a low-precision estimate of the result (about 8
significant digits, shifted down a long way).

It might also be necessary to re-think the choice of local_rscale for
the mul_var() that follows. If the weight of exp is much larger than
the weight of ln_base (or vice versa), it doesn't really make sense to
compute the product to the same local_rscale. That might be a source
of other inaccuracies. I'll try to investigate some more.

Anyway, I don't think any of that should get in the way of committing
the current patch.

Regards,
Dean

Dean Rasheed <dean.a.rasheed@gmail.com> writes:
> On Thu, 22 Jul 2021 at 16:19, Dean Rasheed <dean.a.rasheed@gmail.com> wrote:
>> Thanks for looking. Barring any further comments, I'll push this in a few days.

> So I have been testing this a lot over the last couple of days, and I
> have concluded that the patch works well as far as it goes, but I did
> manage to construct another case where numeric_power() loses
> precision. I think, though, that it would be better to tackle that as
> a separate patch.

It looks like castoroides is not happy with this patch:

https://buildfarm.postgresql.org/cgi-bin/show_log.pl?nm=castoroides&dt=2021-08-01%2008%3A52%3A43

Maybe there's some hidden C99 dependency in what you did?
Although pademelon, which is one of our other pre-C99
dinosaurs, doesn't seem to be unhappy.

            regards, tom lane

On Thu, 5 Aug 2021 at 17:04, Tom Lane <tgl@sss.pgh.pa.us> wrote:
>
> It looks like castoroides is not happy with this patch:
>
> https://buildfarm.postgresql.org/cgi-bin/show_log.pl?nm=castoroides&dt=2021-08-01%2008%3A52%3A43
>
> Maybe there's some hidden C99 dependency in what you did?
> Although pademelon, which is one of our other pre-C99
> dinosaurs, doesn't seem to be unhappy.
>

Hmm, there's something very weird going on there. The 0.9999999999 ^
70000000000000 test, for example, is one that would have thrown an
overflow error before, but it's not doing that.

So somehow, when it hits the overflow/underflow test, Abs(val) is not
greater than NUMERIC_MAX_RESULT_SCALE * 3.01, which is 6020. The thing
is, when I step through it, I get val = -7000, which should trigger
that comfortably. Even if I play with the return value from
estimate_ln_dweight(), which relies on some double precision
arithmetic, making it -11 or -9 instead of -10, I still get val =
-7000. And even if I force val to be -6000, or even val = 0, so that
it doesn't trigger the overflow/underflow test, it still returns zero
in the end. The end result in this case just isn't very sensitive to
changes in these values.

So I'm wondering if it's somehow not even getting that far. Maybe if
the earlier test to see if exp can be represented as an integer is
failing, it might be going through power_var_int() instead, which
would explain it returning a non-zero value. That hypothesis would be
easy to test, by changing the test to 0.9999999999 ^ 70000000000000.5.

In any case, it would be interesting to see what castoroides returns for

select 0.9999999999 ^ 70000000000000;

and

select 0.9999999999 ^ 70000000000000.5;

Regards,
Dean

Dean Rasheed <dean.a.rasheed@gmail.com> writes:
> On Thu, 5 Aug 2021 at 17:04, Tom Lane <tgl@sss.pgh.pa.us> wrote:
>> It looks like castoroides is not happy with this patch:
>> https://buildfarm.postgresql.org/cgi-bin/show_log.pl?nm=castoroides&dt=2021-08-01%2008%3A52%3A43

> Hmm, there's something very weird going on there.

Yeah.  I tried to reproduce this on the gcc compile farm's Solaris 10
machine, but the test passed fine for me.  The only obvious configuration
difference I can find is that that machine has

$ cc -V
cc: Sun C 5.10 SunOS_sparc Patch 141861-10 2012/11/07

whereas castorides' compiler seems to be a few years older.  So this
does seem like it could be a compiler bug.

            regards, tom lane

On Fri, 6 Aug 2021 at 03:58, Tom Lane <tgl@sss.pgh.pa.us> wrote:
>
> Dean Rasheed <dean.a.rasheed@gmail.com> writes:
> > On Thu, 5 Aug 2021 at 17:04, Tom Lane <tgl@sss.pgh.pa.us> wrote:
> >> It looks like castoroides is not happy with this patch:
> >> https://buildfarm.postgresql.org/cgi-bin/show_log.pl?nm=castoroides&dt=2021-08-01%2008%3A52%3A43
>
> > Hmm, there's something very weird going on there.
>
> Yeah.  I tried to reproduce this on the gcc compile farm's Solaris 10
> machine, but the test passed fine for me.  The only obvious configuration
> difference I can find is that that machine has
>
> $ cc -V
> cc: Sun C 5.10 SunOS_sparc Patch 141861-10 2012/11/07
>
> whereas castorides' compiler seems to be a few years older.  So this
> does seem like it could be a compiler bug.
>

Ah, so the latest test results from castoroides confirm my previous
hypothesis, that it isn't even reaching the new code in power_var():

  0.9999999999 ^ 23300000000000 returned 1.0199545627709647
  0.9999999999 ^ 70000000000000 returned 0.9396000441558118

which are actually the results you'd get if you just cast the exponent
to an int32, throwing away the top 32 bits and compute the results:

  0.9999999999 ^ -197580800 = 1.0199545627709647
  0.9999999999 ^ 623009792 = 0.9396000441558118

So the "test for overflow by reverse-conversion" obviously isn't
working as intended, and it's going through power_var_int() when it
shouldn't. I don't think there's anything wrong with that code, so I
think this is a compiler bug.

I guess the best thing to do is just test the value against
PG_INT32_MIN/MAX, which is what int84() does. There are 2 other places
in numeric.c that use similar code to check for int16/32 overflow, so
it's possible that they're broken in the same way on that platform,
but they aren't covered by the regression tests, so it's also possible
that they're OK. Anyway, something like the attached seems likely to
be safer.

Regards,
Dean

Dean Rasheed <dean.a.rasheed@gmail.com> writes:
> So the "test for overflow by reverse-conversion" obviously isn't
> working as intended, and it's going through power_var_int() when it
> shouldn't. I don't think there's anything wrong with that code, so I
> think this is a compiler bug.

Yeah, looks like one.

> I guess the best thing to do is just test the value against
> PG_INT32_MIN/MAX, which is what int84() does. There are 2 other places
> in numeric.c that use similar code to check for int16/32 overflow, so
> it's possible that they're broken in the same way on that platform,
> but they aren't covered by the regression tests, so it's also possible
> that they're OK. Anyway, something like the attached seems likely to
> be safer.

Looks plausible by eyeball (I've not tested).

            regards, tom lane

On Fri, 6 Aug 2021 at 17:15, Tom Lane <tgl@sss.pgh.pa.us> wrote:
>
> > I guess the best thing to do is just test the value against
> > PG_INT32_MIN/MAX, which is what int84() does. There are 2 other places
> > in numeric.c that use similar code to check for int16/32 overflow, so
> > it's possible that they're broken in the same way on that platform,
> > but they aren't covered by the regression tests, so it's also possible
> > that they're OK. Anyway, something like the attached seems likely to
> > be safer.
>
> Looks plausible by eyeball (I've not tested).
>

So, I have back-branch patches for this ready to go. The question is,
is it better to push now, or wait until after next week's releases?

Regards,
Dean

Dean Rasheed <dean.a.rasheed@gmail.com> writes:
> On Fri, 6 Aug 2021 at 17:15, Tom Lane <tgl@sss.pgh.pa.us> wrote:
>> Looks plausible by eyeball (I've not tested).

> So, I have back-branch patches for this ready to go. The question is,
> is it better to push now, or wait until after next week's releases?

I'd push now, given we have a failing buildfarm member.

Admittedly, there may be nobody else using that compiler out in
the real world, but we don't know that.

            regards, tom lane

On Fri, 6 Aug 2021 at 21:26, Tom Lane <tgl@sss.pgh.pa.us> wrote:
>
> Dean Rasheed <dean.a.rasheed@gmail.com> writes:
> > On Fri, 6 Aug 2021 at 17:15, Tom Lane <tgl@sss.pgh.pa.us> wrote:
> >> Looks plausible by eyeball (I've not tested).
>
> > So, I have back-branch patches for this ready to go. The question is,
> > is it better to push now, or wait until after next week's releases?
>
> I'd push now, given we have a failing buildfarm member.
>
> Admittedly, there may be nobody else using that compiler out in
> the real world, but we don't know that.
>

OK. Will do.

Regards,
Dean

On Fri, Aug 06, 2021 at 09:27:03PM +0100, Dean Rasheed wrote:
> On Fri, 6 Aug 2021 at 21:26, Tom Lane <tgl@sss.pgh.pa.us> wrote:
> >
> > Dean Rasheed <dean.a.rasheed@gmail.com> writes:
> > > On Fri, 6 Aug 2021 at 17:15, Tom Lane <tgl@sss.pgh.pa.us> wrote:
> > >> Looks plausible by eyeball (I've not tested).
> >
> > > So, I have back-branch patches for this ready to go. The question is,
> > > is it better to push now, or wait until after next week's releases?
> >
> > I'd push now, given we have a failing buildfarm member.
> >
> > Admittedly, there may be nobody else using that compiler out in
> > the real world, but we don't know that.
> >
> 
> OK. Will do.
> 

Hi Dean,

It seems you already committed this. But it's still as "Ready for
committer" in the commitfest app. 

Are we waiting for something else or we can mark it as committed?

-- 
Jaime Casanova
Director de Servicios Profesionales
SystemGuards - Consultores de PostgreSQL

On Thu, 2 Sept 2021 at 00:39, Jaime Casanova
<jcasanov@systemguards.com.ec> wrote:
>
> Hi Dean,
>
> It seems you already committed this. But it's still as "Ready for
> committer" in the commitfest app.
>
> Are we waiting for something else or we can mark it as committed?
>

It's mostly done, but there is one more corner case where it loses
precision. I'll post an update shortly.

Regards,
Dean

On Thu, Sep 02, 2021 at 07:27:09AM +0100, Dean Rasheed wrote:
> On Thu, 2 Sept 2021 at 00:39, Jaime Casanova
> <jcasanov@systemguards.com.ec> wrote:
> >
> > Hi Dean,
> >
> > It seems you already committed this. But it's still as "Ready for
> > committer" in the commitfest app.
> >
> > Are we waiting for something else or we can mark it as committed?
> >
> 
> It's mostly done, but there is one more corner case where it loses
> precision. I'll post an update shortly.
> 

Great! I'm marking this as "waiting on author".

-- 
Jaime Casanova
Director de Servicios Profesionales
SystemGuards - Consultores de PostgreSQL

On Thu, Sep 02, 2021 at 07:27:09AM +0100, Dean Rasheed wrote:
>
> It's mostly done, but there is one more corner case where it loses
> precision. I'll post an update shortly.
>

I spent some more time looking at this, testing a variety of edge
cases, and the only case I could find that produces inaccurate results
was the one I noted previously -- computing x^y when x is very close
to 1 (less than around 1e-1000 away from it, so that ln_dweight is
less than around -1000). In this case, it loses precision due to the
way local_rscale is set for the initial low-precision calculation:

    local_rscale = 8 - ln_dweight;
    local_rscale = Max(local_rscale, NUMERIC_MIN_DISPLAY_SCALE);
    local_rscale = Min(local_rscale, NUMERIC_MAX_DISPLAY_SCALE);

This needs to be allowed to be greater than NUMERIC_MAX_DISPLAY_SCALE
(1000), otherwise the approximate result will lose all precision,
leading to a poor choice of scale for the full-precision calculation.

So the fix is just to remove the upper bound on this local_rscale, as
we do for the full-precision calculation. This doesn't impact
performance, because it's only computing the logarithm to 8
significant digits at this stage, and when x is very close to 1 like
this, ln_var() has very little work to do -- there is no argument
reduction to do, and the Taylor series terminates on the second term,
since 1-x is so small.

Coming up with a test case that doesn't have thousands of digits is a
bit fiddly, so I chose one where most of the significant digits of the
result are a long way after the decimal point and shifted them up,
which makes the loss of precision in HEAD more obvious. The expected
result can be verified using bc with a scale of 2000.

Regards,
Dean

On 2021-Sep-12, Dean Rasheed wrote:

> So the fix is just to remove the upper bound on this local_rscale, as
> we do for the full-precision calculation. This doesn't impact
> performance, because it's only computing the logarithm to 8
> significant digits at this stage, and when x is very close to 1 like
> this, ln_var() has very little work to do -- there is no argument
> reduction to do, and the Taylor series terminates on the second term,
> since 1-x is so small.

I came here just to opine that there should be a comment about there not
being a clamp to the maximum scale.  For example, log_var says "Set the
scales .. so that they each have more digits ..." which seems clear
enough; I think the new comment is a bit on the short side.

> Coming up with a test case that doesn't have thousands of digits is a
> bit fiddly, so I chose one where most of the significant digits of the
> result are a long way after the decimal point and shifted them up,
> which makes the loss of precision in HEAD more obvious. The expected
> result can be verified using bc with a scale of 2000.

I couldn't get bc (version 1.07.1) to output the result; it says

Runtime warning (func=(main), adr=47): non-zero scale in exponent
Runtime error (func=(main), adr=47): exponent too large in raise

-- 
Álvaro Herrera         PostgreSQL Developer  —  https://www.EnterpriseDB.com/

On Mon, 13 Sept 2021 at 17:51, Alvaro Herrera <alvherre@alvh.no-ip.org> wrote:
>
> I came here just to opine that there should be a comment about there not
> being a clamp to the maximum scale.  For example, log_var says "Set the
> scales .. so that they each have more digits ..." which seems clear
> enough; I think the new comment is a bit on the short side.
>

OK, that's a fair point. Updated version attached.

> I couldn't get bc (version 1.07.1) to output the result; it says
>
> Runtime warning (func=(main), adr=47): non-zero scale in exponent
> Runtime error (func=(main), adr=47): exponent too large in raise
>

Ah yes, bc's "^" operator is a bit limited. It doesn't support
fractional powers for example, and evidently doesn't like powers that
large. I'm so used to not using it that I didn't notice - I always
just use exp() and ln() in bc to compute powers:

scale=2000
e(l(1 - 1.500012345678*10^-1000) * 1.45*10^1003) * 10^1000

Regards,
Dean

On Mon, Sep 13, 2021 at 07:29:13PM +0100, Dean Rasheed wrote:
> On Mon, 13 Sept 2021 at 17:51, Alvaro Herrera <alvherre@alvh.no-ip.org> wrote:
> >
> > I came here just to opine that there should be a comment about there not
> > being a clamp to the maximum scale.  For example, log_var says "Set the
> > scales .. so that they each have more digits ..." which seems clear
> > enough; I think the new comment is a bit on the short side.
> >
> 
> OK, that's a fair point. Updated version attached.
> 

Hi Dean,

Are you planning to commit this soon?

-- 
Jaime Casanova
Director de Servicios Profesionales
SystemGuards - Consultores de PostgreSQL

On Thu, 30 Sept 2021 at 18:25, Jaime Casanova
<jcasanov@systemguards.com.ec> wrote:
>
> Are you planning to commit this soon?
>

Yes, I'll take a look at it next week.

I think it's worth backpatching, despite the fact that it's a pretty
obscure corner case that probably isn't affecting anyone -- similar
fixes in this area have been backpatched, and keeping the code in the
back branches in sync will help with future maintenance and testing,
if any other bugs are found.

Regards,
Dean

On Fri, Oct 01, 2021 at 07:56:33AM +0100, Dean Rasheed wrote:
> On Thu, 30 Sept 2021 at 18:25, Jaime Casanova
> <jcasanov@systemguards.com.ec> wrote:
> >
> > Are you planning to commit this soon?
> >
> 
> Yes, I'll take a look at it next week.
> 

Hi Dean,

Great! I'll move the CF entry to the Next Commitfest so we can move to
closable state.

-- 
Jaime Casanova
Director de Servicios Profesionales
SystemGuards - Consultores de PostgreSQL