Thread: More inaccurate results from numeric pow()
Doing some more testing of the numeric code patched in [1] I noticed
another case where the result is inaccurate -- computing 0.12 ^
-2345.6 gives a very large number containing 2162 digits, but only the
first 2006 correct, while the last 156 digits are wrong.
The reason is this code in power_var():
/* limit to something that won't cause integer overflow */
val = Max(val, -NUMERIC_MAX_RESULT_SCALE);
val = Min(val, NUMERIC_MAX_RESULT_SCALE);
where "val" is the approximate decimal result weight. Here
NUMERIC_MAX_RESULT_SCALE is 2000, so it's clamping the estimated
result weight to 2000, and therefore reducing the rscale in the
subsequent calculations, causing the loss of precision at around 2000
digits.
In fact it's possible to predict exactly how large we need to allow
"val" to become, since the final result is computed using exp_var(),
which accepts inputs up to 6000, so the result weight "val" can be up
to around log10(exp(6000)) ~= 2606 before the final result causes an
overflow.
The obvious fix would be to modify the clamping limits. I think a
better answer though is to replace the clamping code with an overflow
test, immediately throwing an error if "val" is outside the allowed
range, per the attached patch.
This has the advantage that it avoids some expensive computations in
the case where the result will end up overflowing, but more
importantly it means that power_var() isn't so critically dependent on
the limits of exp_var() -- if someone in the future increased the
limits of exp_var() without touching power_var(), and power_var()
clamped to the old range, the problem would resurface. But doing an
overflow test in power_var() instead of clamping "val", it would
either compute an accurate result, or throw an overflow error early
on. There should be no possibility of it returning an inaccurate
result.
Regards,
Dean
[1] http://www.postgresql.org/message-id/CAEZATCV7w+8iB=07dJ8Q0zihXQT1semcQuTeK+4_rogC_zq5Hw@mail.gmail.com
Attachment
Dean Rasheed <dean.a.rasheed@gmail.com> writes:
> In fact it's possible to predict exactly how large we need to allow
> "val" to become, since the final result is computed using exp_var(),
> which accepts inputs up to 6000, so the result weight "val" can be up
> to around log10(exp(6000)) ~= 2606 before the final result causes an
> overflow.
> The obvious fix would be to modify the clamping limits. I think a
> better answer though is to replace the clamping code with an overflow
> test, immediately throwing an error if "val" is outside the allowed
> range, per the attached patch.
I don't much care for the hardwired magic number here, especially since
exp_var() does not have its limit expressed as "6000" but as
"NUMERIC_MAX_RESULT_SCALE * 3". I think you should rephrase the limit
to use that expression, and also add something like this in exp_var():
val = numericvar_to_double_no_overflow(&x);
/* Guard against overflow */
+ /* 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")));
Seems like a reasonable idea otherwise.
regards, tom lane
On Mon, May 2, 2016 at 1:02 PM, Dean Rasheed <dean.a.rasheed@gmail.com> wrote: > Doing some more testing of the numeric code patched in [1] I noticed > another case where the result is inaccurate -- computing 0.12 ^ > -2345.6 gives a very large number containing 2162 digits, but only the > first 2006 correct, while the last 156 digits are wrong. Just out of curiosity, how can you tell? Where do you get alternate output to compare against? Also, I wonder what we think the contract with the user is in cases like this. Surely, if we were dealing with floating point numbers, nobody would expect a calculation like this to be accurate beyond the first n digits, where n is surely much less than 2006. I like the fact that numeric has a lot more precision than any built-in floating point type, but does it have to get every digit in front of the decimal point exactly right no matter how many there are? rhaas=# select tan(pi()::numeric/2), tan(pi()/2); tan | tan -----------------+----------------------618986325617924 | 1.63312393531954e+16 (1 row) -- Robert Haas EnterpriseDB: http://www.enterprisedb.com The Enterprise PostgreSQL Company
On 2 May 2016 at 19:40, Robert Haas <robertmhaas@gmail.com> wrote: > On Mon, May 2, 2016 at 1:02 PM, Dean Rasheed <dean.a.rasheed@gmail.com> wrote: >> Doing some more testing of the numeric code patched in [1] I noticed >> another case where the result is inaccurate -- computing 0.12 ^ >> -2345.6 gives a very large number containing 2162 digits, but only the >> first 2006 correct, while the last 156 digits are wrong. > > Just out of curiosity, how can you tell? Where do you get alternate > output to compare against? > The easiest way is to use bc, although it's pretty slow for this kind of thing. > Also, I wonder what we think the contract with the user is in cases > like this. My expectation is that the numeric computations should generally produce results that are correct in all, or nearly all, digits returned. This is commonly expressed in terms of ULP's -- https://en.wikipedia.org/wiki/Unit_in_the_last_place. An error of a few ULP's may be OK, but I don't think hundreds of ULP's is OK. Actually I think this particular issue is mostly of academic interest, but we went to some effort to get accurate results in the previous patch, and this is just closing a loophole to hopefully complete that work. Surely, if we were dealing with floating point numbers, > nobody would expect a calculation like this to be accurate beyond the > first n digits, where n is surely much less than 2006. I like the > fact that numeric has a lot more precision than any built-in floating > point type, but does it have to get every digit in front of the > decimal point exactly right no matter how many there are? > I would say it should come close. Otherwise we're just returning a lot of noise. Note that there is a limit to how many digits we will ever return, so this is manageable. > rhaas=# select tan(pi()::numeric/2), tan(pi()/2); > tan | tan > -----------------+---------------------- > 618986325617924 | 1.63312393531954e+16 > (1 row) > That doesn't prove anything, since we haven't implemented tan(numeric) (and I don't plan to), so the above is actually select tan((pi()::numeric/2)::float8), tan(pi()/2); and it's not surprising that the 2 results are wildly different (and infinitely far away from the correct answer). It's using tan(float8) in both cases, just with slightly different float8 inputs. There's not really any reason to expect particularly accurate results from float8 functions in cases like this. Regards, Dean
On 2 May 2016 at 18:38, Tom Lane <tgl@sss.pgh.pa.us> wrote: > I don't much care for the hardwired magic number here, especially since > exp_var() does not have its limit expressed as "6000" but as > "NUMERIC_MAX_RESULT_SCALE * 3". I think you should rephrase the limit > to use that expression, and also add something like this in exp_var(): OK, that makes sense. Done that way. Regards, Dean