Re: Numeric x^y for negative x

2021-10-01 Thread Jaime Casanova
On Fri, Oct 01, 2021 at 07:56:33AM +0100, Dean Rasheed wrote:
> On Thu, 30 Sept 2021 at 18:25, Jaime Casanova
>  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




Re: Numeric x^y for negative x

2021-09-30 Thread Dean Rasheed
On Thu, 30 Sept 2021 at 18:25, Jaime Casanova
 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




Re: Numeric x^y for negative x

2021-09-30 Thread Jaime Casanova
On Mon, Sep 13, 2021 at 07:29:13PM +0100, Dean Rasheed wrote:
> On Mon, 13 Sept 2021 at 17:51, Alvaro Herrera  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




Re: Numeric x^y for negative x

2021-09-13 Thread Dean Rasheed
On Mon, 13 Sept 2021 at 17:51, Alvaro Herrera  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
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
new file mode 100644
index 796f517..1de7448
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -10266,9 +10266,13 @@ power_var(const NumericVar *base, const
 	 */
 	ln_dweight = estimate_ln_dweight(base);
 
+	/*
+	 * Set the scale for the low-precision calculation, computing ln(base) to
+	 * around 8 significant digits.  Note that ln_dweight may be as small as
+	 * -SHRT_MAX, so the scale may exceed NUMERIC_MAX_DISPLAY_SCALE here.
+	 */
 	local_rscale = 8 - ln_dweight;
 	local_rscale = Max(local_rscale, NUMERIC_MIN_DISPLAY_SCALE);
-	local_rscale = Min(local_rscale, NUMERIC_MAX_DISPLAY_SCALE);
 
 	ln_var(base, &ln_base, local_rscale);
 
diff --git a/src/test/regress/expected/numeric.out b/src/test/regress/expected/numeric.out
new file mode 100644
index efbb22a..e224eff
--- a/src/test/regress/expected/numeric.out
+++ b/src/test/regress/expected/numeric.out
@@ -2483,6 +2483,12 @@ select coalesce(nullif(0.99 ^ 23
   0
 (1 row)
 
+select round(((1 - 1.500012345678e-1000) ^ 1.45e1003) * 1e1000);
+  round   
+--
+ 25218976308958387188077465658068501556514992509509282366
+(1 row)
+
 -- cases that used to error out
 select 0.12 ^ (-25);
  ?column?  
diff --git a/src/test/regress/sql/numeric.sql b/src/test/regress/sql/numeric.sql
new file mode 100644
index 0418ff0..eeca63d
--- a/src/test/regress/sql/numeric.sql
+++ b/src/test/regress/sql/numeric.sql
@@ -1148,6 +1148,7 @@ select 1.2 ^ 345;
 select 0.12 ^ (-20);
 select 1.0123 ^ (-2147483648);
 select coalesce(nullif(0.99 ^ 233000, 0), 0) as rounds_to_zero;
+select round(((1 - 1.500012345678e-1000) ^ 1.45e1003) * 1e1000);
 
 -- cases that used to error out
 select 0.12 ^ (-25);


Re: Numeric x^y for negative x

2021-09-13 Thread Alvaro Herrera
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/




Re: Numeric x^y for negative x

2021-09-12 Thread Dean Rasheed
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
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
new file mode 100644
index 796f517..65dc3bd
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -10266,9 +10266,9 @@ power_var(const NumericVar *base, const
 	 */
 	ln_dweight = estimate_ln_dweight(base);
 
+	/* scale for initial low-precision calculation (8 significant digits) */
 	local_rscale = 8 - ln_dweight;
 	local_rscale = Max(local_rscale, NUMERIC_MIN_DISPLAY_SCALE);
-	local_rscale = Min(local_rscale, NUMERIC_MAX_DISPLAY_SCALE);
 
 	ln_var(base, &ln_base, local_rscale);
 
diff --git a/src/test/regress/expected/numeric.out b/src/test/regress/expected/numeric.out
new file mode 100644
index efbb22a..e224eff
--- a/src/test/regress/expected/numeric.out
+++ b/src/test/regress/expected/numeric.out
@@ -2483,6 +2483,12 @@ select coalesce(nullif(0.99 ^ 23
   0
 (1 row)
 
+select round(((1 - 1.500012345678e-1000) ^ 1.45e1003) * 1e1000);
+  round   
+--
+ 25218976308958387188077465658068501556514992509509282366
+(1 row)
+
 -- cases that used to error out
 select 0.12 ^ (-25);
  ?column?  
diff --git a/src/test/regress/sql/numeric.sql b/src/test/regress/sql/numeric.sql
new file mode 100644
index 0418ff0..eeca63d
--- a/src/test/regress/sql/numeric.sql
+++ b/src/test/regress/sql/numeric.sql
@@ -1148,6 +1148,7 @@ select 1.2 ^ 345;
 select 0.12 ^ (-20);
 select 1.0123 ^ (-2147483648);
 select coalesce(nullif(0.99 ^ 233000, 0), 0) as rounds_to_zero;
+select round(((1 - 1.500012345678e-1000) ^ 1.45e1003) * 1e1000);
 
 -- cases that used to error out
 select 0.12 ^ (-25);


Re: Numeric x^y for negative x

2021-09-10 Thread Jaime Casanova
On Thu, Sep 02, 2021 at 07:27:09AM +0100, Dean Rasheed wrote:
> On Thu, 2 Sept 2021 at 00:39, Jaime Casanova
>  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




Re: Numeric x^y for negative x

2021-09-01 Thread Dean Rasheed
On Thu, 2 Sept 2021 at 00:39, Jaime Casanova
 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




Re: Numeric x^y for negative x

2021-09-01 Thread Jaime Casanova
On Fri, Aug 06, 2021 at 09:27:03PM +0100, Dean Rasheed wrote:
> On Fri, 6 Aug 2021 at 21:26, Tom Lane  wrote:
> >
> > Dean Rasheed  writes:
> > > On Fri, 6 Aug 2021 at 17:15, Tom Lane  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




Re: Numeric x^y for negative x

2021-08-06 Thread Dean Rasheed
On Fri, 6 Aug 2021 at 21:26, Tom Lane  wrote:
>
> Dean Rasheed  writes:
> > On Fri, 6 Aug 2021 at 17:15, Tom Lane  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




Re: Numeric x^y for negative x

2021-08-06 Thread Tom Lane
Dean Rasheed  writes:
> On Fri, 6 Aug 2021 at 17:15, Tom Lane  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




Re: Numeric x^y for negative x

2021-08-06 Thread Dean Rasheed
On Fri, 6 Aug 2021 at 17:15, Tom Lane  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




Re: Numeric x^y for negative x

2021-08-06 Thread Tom Lane
Dean Rasheed  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




Re: Numeric x^y for negative x

2021-08-06 Thread Dean Rasheed
On Fri, 6 Aug 2021 at 03:58, Tom Lane  wrote:
>
> Dean Rasheed  writes:
> > On Thu, 5 Aug 2021 at 17:04, Tom Lane  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.99 ^ 233000 returned 1.0199545627709647
  0.99 ^ 70 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.99 ^ -197580800 = 1.0199545627709647
  0.99 ^ 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
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
new file mode 100644
index 13bb968..0c210a3
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -4289,11 +4289,13 @@ numericvar_to_int32(const NumericVar *va
 	if (!numericvar_to_int64(var, &val))
 		return false;
 
+	if (unlikely(val < PG_INT32_MIN) || unlikely(val > PG_INT32_MAX))
+		return false;
+	
 	/* Down-convert to int4 */
 	*result = (int32) val;
 
-	/* Test for overflow by reverse-conversion. */
-	return ((int64) *result == val);
+	return true;
 }
 
 Datum
@@ -4373,15 +4375,14 @@ numeric_int2(PG_FUNCTION_ARGS)
 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
  errmsg("smallint out of range")));
 
-	/* Down-convert to int2 */
-	result = (int16) val;
-
-	/* Test for overflow by reverse-conversion. */
-	if ((int64) result != val)
+	if (unlikely(val < PG_INT16_MIN) || unlikely(val > PG_INT16_MAX))
 		ereport(ERROR,
 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
  errmsg("smallint out of range")));
 
+	/* Down-convert to int2 */
+	result = (int16) val;
+
 	PG_RETURN_INT16(result);
 }
 
@@ -10186,10 +10187,7 @@ power_var(const NumericVar *base, const
 
 		if (numericvar_to_int64(exp, &expval64))
 		{
-			int			expval = (int) expval64;
-
-			/* Test for overflow by reverse-conversion. */
-			if ((int64) expval == expval64)
+			if (expval64 >= PG_INT32_MIN && expval64 <= PG_INT32_MAX)
 			{
 /* Okay, select rscale */
 rscale = NUMERIC_MIN_SIG_DIGITS;
@@ -10197,7 +10195,7 @@ power_var(const NumericVar *base, const
 rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
 rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
 
-power_var_int(base, expval, result, rscale);
+power_var_int(base, (int) expval64, result, rscale);
 return;
 			}
 		}
diff --git a/src/test/regress/expected/numeric.out b/src/test/regress/expected/numeric.out
new file mode 100644
index 80b42f8..3707aff
--- a/src/test/regress/expected/numeric.out
+++ b/src/test/regress/expected/numeric.out
@@ -1154,6 +1154,55 @@ SELECT * FROM fract_only;
 (7 rows)
 
 DROP TABLE fract_only;
+-- Check conversion to integers
+SELECT (-9223372036854775808.5)::int8; -- should fail
+ERROR:  bigint out of range
+SELECT (-9223372036854775808.4)::int8; -- ok
+ int8 
+--
+ -9223372036854775808
+(1 row)
+
+SELECT 9223372036854775807.4::int8; -- ok
+int8 
+-
+ 9223372036854775807
+(1 row)
+
+SELECT 9223372036854775807.5::int8; -- should fail
+ERROR:  bigint out of range
+SELECT (-2147483648.5)::int4; -- should fail
+ERROR:  integer out of range
+SELECT (-2147483648.4)::int4; -- ok
+int4 
+-
+ -2147483648
+(1 row)
+
+SELECT 2147483647.4::int4; -- ok
+int4
+
+ 2147483647
+(1 row)
+
+SELECT 2147483647.5::int4; -- should fail
+ERROR:  integer out of range
+SELECT (-32768.5)::int2; -- should fail
+ERROR:  smallint out of range
+SELECT (-32768.4)::int2; -- ok
+  int2  
+
+ -32768
+(1 row)
+
+SELECT 32767.4::int2; -- ok
+ int2  
+---
+ 32767
+(1 row)
+
+SEL

Re: Numeric x^y for negative x

2021-08-05 Thread Tom Lane
Dean Rasheed  writes:
> On Thu, 5 Aug 2021 at 17:04, Tom Lane  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




Re: Numeric x^y for negative x

2021-08-05 Thread Dean Rasheed
On Thu, 5 Aug 2021 at 17:04, Tom Lane  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.99 ^
70 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.99 ^ 70.5.

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

select 0.99 ^ 70;

and

select 0.99 ^ 70.5;

Regards,
Dean




Re: Numeric x^y for negative x

2021-08-05 Thread Tom Lane
Dean Rasheed  writes:
> On Thu, 22 Jul 2021 at 16:19, Dean Rasheed  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




Re: Numeric x^y for negative x

2021-07-29 Thread Dean Rasheed
On Thu, 22 Jul 2021 at 16:19, Dean Rasheed  wrote:
>
> On Thu, 22 Jul 2021 at 06:13, Yugo NAGATA  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.33.

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
From 071d877af3ffcb8e2abd445cb5963d90a2e766bc Mon Sep 17 00:00:00 2001
From: Dean Rasheed 
Date: Thu, 29 Jul 2021 17:26:08 +0100
Subject: [PATCH] Fix corner-case errors and loss of precision in
 numeric_power().

This fixes a couple of related problems that arise when raising
numbers to very large powers.

Firstly, when raising a negative number to a very large integer power,
the result should be well-defined, but the previous code would only
cope if the exponent was small enough to go through power_var_int().
Otherwise it would throw an internal error, attempting to take the
logarithm of a negative number. Fix this by adding suitable handling
to the general case in power_var() to cope with negative bases,
checking for integer powers there.

Next, when raising a (positive or negative) number whose absolute
value is slightly less than 1 to a very large power, the result should
approach zero as the power is increased. However, in some cases, for
sufficiently large powers, this would lose all precision and return 1
instead. This was due to the way that the local_rscale was being
calculated for the final full-precision calculation:

  local_rscale = rscale + (int) val - ln_dweight + 8

The first two terms on the right hand side are meant to give the
number of significant digits required in the result ("val" being the
estimated result weight). However, this failed to account for the fact
that rscale is clipped to a maximum of NUMERIC_MAX_DISPLAY_SCALE
(1000), and the result weight might be less then -1000, causing their
sum to be negative, leading to a loss of precision. Fix this by
forcing the number of significant digits to be nonnegative. It's OK
for it to be zero (when the result weight is less than -1000), since
the local_rscale value then includes a few extra digits to ensure an

Re: Numeric x^y for negative x

2021-07-22 Thread Dean Rasheed
On Thu, 22 Jul 2021 at 06:13, Yugo NAGATA  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




Re: Numeric x^y for negative x

2021-07-21 Thread Yugo NAGATA
On Wed, 21 Jul 2021 11:10:16 +0100
Dean Rasheed  wrote:

> On Tue, 20 Jul 2021 at 10:17, Yugo NAGATA  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 




Re: Numeric x^y for negative x

2021-07-21 Thread Dean Rasheed
On Tue, 20 Jul 2021 at 10:17, Yugo NAGATA  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
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
new file mode 100644
index 2a0f68f..8a8822c
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -3935,7 +3935,9 @@ numeric_power(PG_FUNCTION_ARGS)
 	/*
 	 * 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.
+	 * divide-by-zero error code for 0 ^ -1.  Raising a negative number to a
+	 * non-integer power must produce the same error code, but that case is
+	 * handled in power_var().
 	 */
 	sign1 = numeric_sign_internal(num1);
 	sign2 = numeric_sign_internal(num2);
@@ -3945,11 +3947,6 @@ numeric_power(PG_FUNCTION_ARGS)
 (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
 	 */
@@ -9762,12 +9759,18 @@ exp_var(const NumericVar *arg, NumericVa
 	 */
 	val = numericvar_to_double_no_overflow(&x);
 
-	/* Guard against overflow */
+	/* Guard against overflow/underflow */
 	/* If you change this limit, see also power_var()'s limit */
 	if (Abs(val) >= NUMERIC_MAX_RESULT_SCALE * 3)
-		ereport(ERROR,
-(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
- errmsg("value overflows numeric format")));
+	{
+		if (val > 0)
+			ereport(ERROR,
+	(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+	 errmsg("value overflows numeric format")));
+		zero_var(result);
+		result->dscale = rscale;
+		return;
+	}
 
 	/* decimal weight = log10(e^x) = x * log10(e) */
 	dweight = (int) (val * 0.434294481903252);
@@ -10125,6 +10128,8 @@ log_var(const NumericVar *base, const Nu
 static void
 power_var(const NumericVar *base, const NumericVar *exp, NumericVar *result)
 {
+	int			res_sign;
+	NumericVar	abs_base;
 	NumericVar	ln_base;
 	NumericVar	ln_num;
 	int			ln_dweight;
@@ -10168,9 +10173,40 @@ power_var(const NumericVar *base, const
 		return;
 	}
 
+	init_var(&abs_base)

Re: Numeric x^y for negative x

2021-07-20 Thread Yugo NAGATA
On Wed, 7 Jul 2021 18:36:56 +0100
Dean Rasheed  wrote:

> On Thu, 1 Jul 2021 at 14:17, Dean Rasheed  wrote:
> >
> > On Tue, 29 Jun 2021 at 12:08, Dean Rasheed  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 




Re: Numeric x^y for negative x

2021-07-07 Thread Dean Rasheed
On Wed, 7 Jul 2021 at 18:57, Zhihong Yu  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




Re: Numeric x^y for negative x

2021-07-07 Thread Zhihong Yu
On Wed, Jul 7, 2021 at 10:37 AM Dean Rasheed 
wrote:

> On Thu, 1 Jul 2021 at 14:17, Dean Rasheed 
> wrote:
> >
> > On Tue, 29 Jun 2021 at 12:08, Dean Rasheed 
> 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
>
Hi,

+   (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

Cheers


Re: Numeric x^y for negative x

2021-07-07 Thread Dean Rasheed
On Thu, 1 Jul 2021 at 14:17, Dean Rasheed  wrote:
>
> On Tue, 29 Jun 2021 at 12:08, Dean Rasheed  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
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
new file mode 100644
index bc71326..1481538
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -3937,11 +3937,6 @@ numeric_power(PG_FUNCTION_ARGS)
 (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
 	 */
@@ -9754,12 +9749,18 @@ exp_var(const NumericVar *arg, NumericVa
 	 */
 	val = numericvar_to_double_no_overflow(&x);
 
-	/* Guard against overflow */
+	/* Guard against overflow/underflow */
 	/* If you change this limit, see also power_var()'s limit */
 	if (Abs(val) >= NUMERIC_MAX_RESULT_SCALE * 3)
-		ereport(ERROR,
-(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
- errmsg("value overflows numeric format")));
+	{
+		if (val > 0)
+			ereport(ERROR,
+	(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+	 errmsg("value overflows numeric format")));
+		zero_var(result);
+		result->dscale = rscale;
+		return;
+	}
 
 	/* decimal weight = log10(e^x) = x * log10(e) */
 	dweight = (int) (val * 0.434294481903252);
@@ -10117,6 +10118,8 @@ log_var(const NumericVar *base, const Nu
 static void
 power_var(const NumericVar *base, const NumericVar *exp, NumericVar *result)
 {
+	int			res_sign;
+	NumericVar	abs_base;
 	NumericVar	ln_base;
 	NumericVar	ln_num;
 	int			ln_dweight;
@@ -10160,9 +10163,37 @@ power_var(const NumericVar *base, const
 		return;
 	}
 
+	init_var(&abs_base);
 	init_var(&ln_base);
 	init_var(&ln_num);
 
+	/*
+	 * If base is negative, insist that exp be an integer.  The result is then
+	 * positive if exp is even and negative if exp is odd.
+	 */
+	if (base->sign == NUMERIC_NEG)
+	{
+		/* check that exp is an integer */
+		if (exp->ndigits > 0 && exp->ndigits > exp->weight + 1)
+			ereport(ERROR,
+	(errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
+	 errmsg("a negative number raised to a non-integer power yields a complex result")));
+
+		/* test if exp is odd or even */
+		if (exp->ndigits > 0 && exp->ndigits == exp->weight + 1 &&
+			(exp->digits[exp->ndigits - 1] & 1))
+			res_sign = NUMERIC_NEG;
+		else
+			res_sign = NUMERIC_POS;
+
+		/* then work with abs(base) below */
+		set_var_from_var(base, &abs_base);
+		abs_base.sign = NUMERIC_POS;
+		base = &abs_base;
+	}
+	else
+		res_sign = NUMERIC_POS;
+
 	/*--
 	 * Decide on the scale for the ln() calculation.  For this we need an
 	 * estimate of the weight of the result, which we obtain by doing an
@@ -10193,25 +10224,31 @@ power_var(const NumericVar *base, const
 
 	val = numericvar_to_double_no_overflow(&ln_num);
 
-	/* initial overflow test with fuzz factor */
+	/* initial overflow/underflow test with fuzz factor */
 	if (Abs(val) > NUMERIC_MAX_RESULT_SCALE * 3.01)
-		ereport(ERROR,
-(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
- errmsg("value overflows numeric format")));
+	{
+		if (val > 0)
+			ereport(ERROR,
+	(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+	 errmsg("value overflows numeric format")));
+		zero_var(result);
+		result->dscale = NUMERIC_MAX_DISPLAY_SCALE;
+		return;
+	}
 
 	val *= 0.434294481903252;	/* approximate decimal result weight */
 
-	/* choose the result scale */
+	/* choose the result scale and the scale for the real calculation */
 	rscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
 	rscale = Max(rscale, base->dscale);
 	rscale = Max(rscale, exp->dscale);
 	rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
-	rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
 
-	/* set the scale for the real exp * ln(base) calculation */
 	local_rscale = rscale + (int) val - ln_dweight + 8;
 	local_rscale = Max(local_rscale, NUMERIC_MIN_DISPLAY_SCALE);
 
+	rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
+
 	/* and do the real calculation */
 
 	ln_var(base, &ln_base, local_rscale);
@@ -10220,8 +10257,12 @@ power_var(const NumericVar *base, const
 
 	exp_var(&ln_num, result, rscale);
 
+	if (res_sign == NUMERIC_NEG && result->ndig

Re: Numeric x^y for negative x

2021-07-01 Thread Dean Rasheed
On Tue, 29 Jun 2021 at 12:08, Dean Rasheed  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.99 to increasingly large
powers, which should decrease to zero:

   exp 0.99^exp
100 0.3678794411530483
10000.[4 zeros]4539992973978489
1   0.[43 zeros]3720075957420456
10  0.[434 zeros]5075958643751518
20  0.[868 zeros]2576535615307575
21  0.[912 zeros]9584908195943232
22  0.[955 zeros]3565658653381070
23  0.[998 zeros]13
231000  0.[1000 zeros]
232000  0.[1000 zeros]
233000  1.[1000 zeros]  *** WRONG ***
30  1.[1000 zeros]  *** WRONG ***
40  1.[1000 zeros]  *** WRONG ***
50  1.[1000 zeros]  *** WRONG ***
60  1.[1000 zeros]  *** WRONG ***
70  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




Numeric x^y for negative x

2021-06-29 Thread Dean Rasheed
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.
(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
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
new file mode 100644
index eb78f0b..3eb0d90
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -3934,11 +3934,6 @@ numeric_power(PG_FUNCTION_ARGS)
 (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
 	 */
@@ -10138,6 +10133,8 @@ log_var(const NumericVar *base, const Nu
 static void
 power_var(const NumericVar *base, const NumericVar *exp, NumericVar *result)
 {
+	int			res_sign;
+	NumericVar	abs_base;
 	NumericVar	ln_base;
 	NumericVar	ln_num;
 	int			ln_dweight;
@@ -10181,9 +10178,37 @@ power_var(const NumericVar *base, const
 		return;
 	}
 
+	init_var(&abs_base);
 	init_var(&ln_base);
 	init_var(&ln_num);
 
+	/*
+	 * If base is negative, insist that exp be an integer.  The result is then
+	 * positive if exp is even and negative if exp is odd.
+	 */
+	if (base->sign == NUMERIC_NEG)
+	{
+		/* digits after the decimal point? */
+		if (exp->ndigits > 0 && exp->ndigits > exp->weight + 1)
+			ereport(ERROR,
+	(errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
+	 errmsg("a negative number raised to a non-integer power yields a complex result")));
+
+		/* odd exponent? */
+		if (exp->ndigits > 0 && exp->ndigits == exp->weight + 1 &&
+			(exp->digits[exp->ndigits - 1] & 1))
+			res_sign = NUMERIC_NEG;
+		else
+			res_sign = NUMERIC_POS;
+
+		/* work with abs(base) */
+		set_var_from_var(base, &abs_base);
+		abs_base.sign = NUMERIC_POS;
+		base = &abs_base;
+	}
+	else
+		res_sign = NUMERIC_POS;
+
 	/*--
 	 * Decide on the scale for the ln() calculation.  For this we need an
 	 * estimate of the weight of the result, which we obtain by doing an
@@ -10240,9 +10265,12 @@ power_var(const NumericVar *base, const
 	mul_var(&ln_base, exp, &ln_num, local_rscale);
 
 	exp_var(&ln_num, result, rscale);
+	if (res_sign == NUMERIC_NEG && result->ndigits > 0)
+		result->sign = NUMERIC_NEG;
 
 	free_var(&ln_num);
 	free_var(&ln_base);
+	free_var(&abs_base);
 }
 
 /*
diff --git a/src/test/regress/expected/numeric.out b/src/test/regress/expected/numeric.out
new file mode 100644
index 30a5642..94bc773
--- a/src/test/regress/expected/numeric.out
+++ b/src/test/regress/expected/numeric.out
@@ -2340,6 +2340,37 @@ select 0.5678 ^ (-85);
  782333637740774446257.7719390061997396
 (1 row)
 
+-- negative base to integer powers
+select (-1.0) ^ 2147483646;
+  ?column?  
+
+ 1.
+(1 row)
+
+select (-1.0) ^ 2147483647;
+  ?column?   
+-
+ -1.
+(1 row)
+
+select (-1.0) ^ 2147483648;
+  ?column?  
+
+ 1.
+(1 row)
+
+select (-1.0) ^ 1000;
+  ?column?  
+
+ 1.
+(1 row)
+
+select (-1.0) ^ 1001;
+  ?column?   
+-
+ -1.
+(1 row)
+
 --
 -- Tests for raising to non-integer powers
 --
diff --git a/src/test/regress/sql/numeric.sql b/src/test/regress/sql/numeric.sql
new file mode 100644
index db812c8..58c1c69
--- a/src/test/regress/sql/numeric.sql
+++ b/src/test/regress/sql/numeric.sql
@@ -1095,6 +1095,13 @@ select 1.0123 ^ (-2147483648);
 select 0.12 ^ (-25);
 select 0.5678 ^ (-85);
 
+-- negative base to integer powers
+select (-1.0) ^ 2147483646;
+select (-1.0) ^ 2147483647;
+select (-1.0) ^ 2147483648;
+select (-1.0) ^ 1000;
+select (-1.0) ^ 1001;
+
 --
 -- Tests for raising to non-integer powers
 --