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.9999999999 ^ 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.000000000123 ^ (-2147483648);
select coalesce(nullif(0.9999999999 ^ 23300000000000, 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);
